My Project
Functions
cf_chinese.cc File Reference

algorithms for chinese remaindering and rational reconstruction More...

#include "config.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"
#include "cf_assert.h"
#include "debug.h"
#include "canonicalform.h"
#include "cf_iter.h"
#include "cf_util.h"

Go to the source code of this file.

Functions

void chineseRemainder (const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
 void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew ) More...
 
void chineseRemainder (const CFArray &x, const CFArray &q, CanonicalForm &xnew, CanonicalForm &qnew)
 void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew ) More...
 
CanonicalForm Farey (const CanonicalForm &f, const CanonicalForm &q)
 Farey rational reconstruction. More...
 
static CanonicalForm chin_mul_inv (const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
 
void out_cf (const char *s1, const CanonicalForm &f, const char *s2)
 cf_algorithm.cc - simple mathematical algorithms. More...
 
void chineseRemainderCached (const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
 
void chineseRemainderCached (const CanonicalForm &a, const CanonicalForm &q1, const CanonicalForm &b, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew, CFArray &inv)
 

Detailed Description

algorithms for chinese remaindering and rational reconstruction

Used by: cf_gcd.cc, cf_linsys.cc

Header file: cf_algorithm.h

Definition in file cf_chinese.cc.

Function Documentation

◆ chin_mul_inv()

static CanonicalForm chin_mul_inv ( const CanonicalForm  a,
const CanonicalForm  b,
int  ind,
CFArray inv 
)
inlinestatic

Definition at line 276 of file cf_chinese.cc.

277 {
278  if (inv[ind].isZero())
279  {
280  CanonicalForm s,dummy;
281  (void)bextgcd( a, b, s, dummy );
282  inv[ind]=s;
283  return s;
284  }
285  else
286  return inv[ind];
287 }
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
CanonicalForm b
Definition: cfModGcd.cc:4103
factory's main class
Definition: canonicalform.h:86
const CanonicalForm int s
Definition: facAbsFact.cc:51
bool isZero(const CFArray &A)
checks if entries of A are zero

◆ chineseRemainder() [1/2]

void chineseRemainder ( const CanonicalForm x1,
const CanonicalForm q1,
const CanonicalForm x2,
const CanonicalForm q2,
CanonicalForm xnew,
CanonicalForm qnew 
)

void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2) and qnew = q1*q2. q1 and q2 should be positive integers, pairwise prime, x1 and x2 should be polynomials with integer coefficients. If x1 and x2 are polynomials with positive coefficients, the result is guaranteed to have positive coefficients, too.

Note: This algorithm is optimized for the case q1>>q2.

This is a standard algorithm. See, for example, Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra', par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by Homomorphic Images' in B. Buchberger - 'Computer Algebra - Symbolic and Algebraic Computation'.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 57 of file cf_chinese.cc.

58 {
59  DEBINCLEVEL( cerr, "chineseRemainder" );
60 
61  DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
62  DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
63 
64  // We calculate xnew as follows:
65  // xnew = v1 + v2 * q1
66  // where
67  // v1 = x1 (mod q1)
68  // v2 = (x2-v1)/q1 (mod q2) (*)
69  //
70  // We do one extra test to check whether x2-v1 vanishes (mod
71  // q2) in (*) since it is not costly and may save us
72  // from calculating the inverse of q1 (mod q2).
73  //
74  // u: v1 (mod q2)
75  // d: x2-v1 (mod q2)
76  // s: 1/q1 (mod q2)
77  //
78  CanonicalForm v2, v1;
79  CanonicalForm u, d, s, dummy;
80 
81  v1 = mod( x1, q1 );
82  u = mod( v1, q2 );
83  d = mod( x2-u, q2 );
84  if ( d.isZero() )
85  {
86  xnew = v1;
87  qnew = q1 * q2;
88  DEBDECLEVEL( cerr, "chineseRemainder" );
89  return;
90  }
91  (void)bextgcd( q1, q2, s, dummy );
92  v2 = mod( d*s, q2 );
93  xnew = v1 + v2*q1;
94 
95  // After all, calculate new modulus. It is important that
96  // this is done at the very end of the algorithm, since q1
97  // and qnew may refer to the same object (same is true for x1
98  // and xnew).
99  qnew = q1 * q2;
100 
101  DEBDECLEVEL( cerr, "chineseRemainder" );
102 }
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CF_NO_INLINE bool isZero() const
int ilog2() const
int CanonicalForm::ilog2 () const
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45

◆ chineseRemainder() [2/2]

void chineseRemainder ( const CFArray x,
const CFArray q,
CanonicalForm xnew,
CanonicalForm qnew 
)

void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )

chineseRemainder - integer chinese remaindering.

Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the product of all q[i]. q[i] should be positive integers, pairwise prime. x[i] should be polynomials with integer coefficients. If all coefficients of all x[i] are positive integers, the result is guaranteed to have positive coefficients, too.

This is a standard algorithm, too, except for the fact that we use a divide-and-conquer method instead of a linear approach to calculate the remainder.

Note: Be sure you are calculating in Z, and not in Q!

Definition at line 124 of file cf_chinese.cc.

125 {
126  DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
127 
128  ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
129  CFArray X(x), Q(q);
130  int i, j, n = x.size(), start = x.min();
131 
132  DEBOUTLN( cerr, "array size = " << n );
133 
134  while ( n != 1 )
135  {
136  i = j = start;
137  while ( i < start + n - 1 )
138  {
139  // This is a little bit dangerous: X[i] and X[j] (and
140  // Q[i] and Q[j]) may refer to the same object. But
141  // xnew and qnew in the above function are modified
142  // at the very end of the function, so we do not
143  // modify x1 and q1, resp., by accident.
144  chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
145  i += 2;
146  j++;
147  }
148 
149  if ( n & 1 )
150  {
151  X[j] = X[i];
152  Q[j] = Q[i];
153  }
154  // Maybe we would get some memory back at this point if
155  // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
156  // at this point?
157 
158  n = ( n + 1) / 2;
159  }
160  xnew = X[start];
161  qnew = Q[q.min()];
162 
163  DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
164 }
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
#define ASSERT(expression, message)
Definition: cf_assert.h:99
void chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition: cf_chinese.cc:57
int size() const
Definition: ftmpl_array.cc:92
int min() const
Definition: ftmpl_array.cc:98
int j
Definition: facHensel.cc:110
STATIC_VAR jList * Q
Definition: janet.cc:30

◆ chineseRemainderCached() [1/2]

void chineseRemainderCached ( const CanonicalForm a,
const CanonicalForm q1,
const CanonicalForm b,
const CanonicalForm q2,
CanonicalForm xnew,
CanonicalForm qnew,
CFArray inv 
)

Definition at line 308 of file cf_chinese.cc.

309 {
310  CFArray A(2); A[0]=a; A[1]=b;
311  CFArray Q(2); Q[0]=q1; Q[1]=q2;
312  chineseRemainderCached(A,Q,xnew,qnew,inv);
313 }
void chineseRemainderCached(const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
Definition: cf_chinese.cc:290
#define A
Definition: sirandom.c:24

◆ chineseRemainderCached() [2/2]

void chineseRemainderCached ( const CFArray a,
const CFArray n,
CanonicalForm xnew,
CanonicalForm prod,
CFArray inv 
)

Definition at line 290 of file cf_chinese.cc.

291 {
292  CanonicalForm p, sum=0L; prod=1L;
293  int i;
294  int len=n.size();
295 
296  for (i = 0; i < len; i++) prod *= n[i];
297 
298  for (i = 0; i < len; i++)
299  {
300  p = prod / n[i];
301  sum += a[i] * chin_mul_inv(p, n[i], i, inv) * p;
302  }
303 
304  xnew = mod(sum , prod);
305 }
int p
Definition: cfModGcd.cc:4078
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:276
fq_nmod_poly_t prod
Definition: facHensel.cc:100

◆ Farey()

Farey rational reconstruction.

If NTL is available it uses the fast algorithm from NTL, i.e. Encarnacion, Collins.

Definition at line 202 of file cf_chinese.cc.

203 {
204  int is_rat=isOn(SW_RATIONAL);
205  Off(SW_RATIONAL);
206  Variable x = f.mvar();
207  CanonicalForm result = 0;
208  CanonicalForm c;
209  CFIterator i;
210 #ifdef HAVE_FLINT
211  fmpz_t FLINTq;
212  fmpz_init(FLINTq);
213  convertCF2initFmpz(FLINTq,q);
214  fmpz_t FLINTc;
215  fmpz_init(FLINTc);
216  fmpq_t FLINTres;
217  fmpq_init(FLINTres);
218 #elif defined(HAVE_NTL)
219  ZZ NTLq= convertFacCF2NTLZZ (q);
220  ZZ bound;
221  SqrRoot (bound, NTLq/2);
222 #else
223  factoryError("NTL/FLINT missing:Farey");
224 #endif
225  for ( i = f; i.hasTerms(); i++ )
226  {
227  c = i.coeff();
228  if ( c.inCoeffDomain())
229  {
230 #ifdef HAVE_FLINT
231  if (c.inZ())
232  {
233  convertCF2initFmpz(FLINTc,c);
234  fmpq_reconstruct_fmpz(FLINTres,FLINTc,FLINTq);
235  result += power (x, i.exp())*(convertFmpq2CF(FLINTres));
236  }
237 #elif defined(HAVE_NTL)
238  if (c.inZ())
239  {
240  ZZ NTLc= convertFacCF2NTLZZ (c);
241  bool lessZero= (sign (NTLc) == -1);
242  if (lessZero)
243  NTL::negate (NTLc, NTLc);
244  ZZ NTLnum, NTLden;
245  if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
246  {
247  if (lessZero)
248  NTL::negate (NTLnum, NTLnum);
249  CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
250  CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
251  On (SW_RATIONAL);
252  result += power (x, i.exp())*(num/den);
253  Off (SW_RATIONAL);
254  }
255  }
256 #else
257  if (c.inZ())
258  result += power (x, i.exp()) * Farey_n(c,q);
259 #endif
260  else
261  result += power( x, i.exp() ) * Farey(c,q);
262  }
263  else
264  result += power( x, i.exp() ) * Farey(c,q);
265  }
266  if (is_rat) On(SW_RATIONAL);
267 #ifdef HAVE_FLINT
268  fmpq_clear(FLINTres);
269  fmpz_clear(FLINTc);
270  fmpz_clear(FLINTq);
271 #endif
272  return result;
273 }
CanonicalForm convertFmpq2CF(const fmpq_t q)
conversion of a FLINT rational to CanonicalForm
void convertCF2initFmpz(fmpz_t result, const CanonicalForm &f)
conversion of a factory integer to fmpz_t(init.)
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:670
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:202
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
VAR void(* factoryError)(const char *s)
Definition: cf_util.cc:80
FILE * f
Definition: checklibs.c:9
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
bool inCoeffDomain() const
bool inZ() const
predicates
factory's class for variables
Definition: factory.h:127
return result
Definition: facAbsBiFact.cc:75
static int sign(int x)
Definition: ring.cc:3469

◆ out_cf()

void out_cf ( const char *  s1,
const CanonicalForm f,
const char *  s2 
)

cf_algorithm.cc - simple mathematical algorithms.

Hierarchy: mathematical algorithms on canonical forms

Developers note:

A "mathematical" algorithm is an algorithm which calculates some mathematical function in contrast to a "structural" algorithm which gives structural information on polynomials.

Compare these functions to the functions in ‘cf_ops.cc’, which are structural algorithms.

Definition at line 99 of file cf_factor.cc.

100 {
101  printf("%s",s1);
102  if (f.isZero()) printf("+0");
103  //else if (! f.inCoeffDomain() )
104  else if (! f.inBaseDomain() )
105  {
106  int l = f.level();
107  for ( CFIterator i = f; i.hasTerms(); i++ )
108  {
109  int e=i.exp();
110  if (i.coeff().isOne())
111  {
112  printf("+");
113  if (e==0) printf("1");
114  else
115  {
116  printf("%c",'a'+l-1);
117  if (e!=1) printf("^%d",e);
118  }
119  }
120  else
121  {
122  out_cf("+(",i.coeff(),")");
123  if (e!=0)
124  {
125  printf("*%c",'a'+l-1);
126  if (e!=1) printf("^%d",e);
127  }
128  }
129  }
130  }
131  else
132  {
133  if ( f.isImm() )
134  {
136  {
137  long a= imm2int (f.getval());
138  if ( a == gf_q )
139  printf ("+%ld", a);
140  else if ( a == 0L )
141  printf ("+1");
142  else if ( a == 1L )
143  printf ("+%c",gf_name);
144  else
145  {
146  printf ("+%c",gf_name);
147  printf ("^%ld",a);
148  }
149  }
150  else
151  {
152  long l=f.intval();
153  if (l<0) printf("%ld",l);
154  else printf("+%ld",l);
155  }
156  }
157  else
158  {
159  #ifdef NOSTREAMIO
160  if (f.inZ())
161  {
162  mpz_t m;
163  gmp_numerator(f,m);
164  char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
165  str = mpz_get_str( str, 10, m );
166  puts(str);
167  delete[] str;
168  mpz_clear(m);
169  }
170  else if (f.inQ())
171  {
172  mpz_t m;
173  gmp_numerator(f,m);
174  char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
175  str = mpz_get_str( str, 10, m );
176  while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
177  puts(str);putchar('/');
178  delete[] str;
179  mpz_clear(m);
181  str = new char[mpz_sizeinbase( m, 10 ) + 2];
182  str = mpz_get_str( str, 10, m );
183  while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
184  puts(str);
185  delete[] str;
186  mpz_clear(m);
187  }
188  #else
189  std::cout << f;
190  #endif
191  }
192  //if (f.inZ()) printf("(Z)");
193  //else if (f.inQ()) printf("(Q)");
194  //else if (f.inFF()) printf("(FF)");
195  //else if (f.inPP()) printf("(PP)");
196  //else if (f.inGF()) printf("(PP)");
197  //else
198  if (f.inExtension()) printf("E(%d)",f.level());
199  }
200  printf("%s",s2);
201 }
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
#define GaloisFieldDomain
Definition: cf_defs.h:18
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
cf_algorithm.cc - simple mathematical algorithms.
Definition: cf_factor.cc:99
static int gettype()
Definition: cf_factory.h:28
void FACTORY_PUBLIC gmp_numerator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:20
void FACTORY_PUBLIC gmp_denominator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:40
VAR int gf_q
Definition: gfops.cc:47
VAR char gf_name
Definition: gfops.cc:52
static long imm2int(const InternalCF *const imm)
Definition: imm.h:70
char * str(leftv arg)
Definition: shared.cc:704