My Project
cf_chinese.cc
Go to the documentation of this file.
1 /* emacs edit mode for this file is -*- C++ -*- */
2 
3 /**
4  * @file cf_chinese.cc
5  *
6  * algorithms for chinese remaindering and rational reconstruction
7  *
8  * Used by: cf_gcd.cc, cf_linsys.cc
9  *
10  * Header file: cf_algorithm.h
11  *
12 **/
13 
14 
15 #include "config.h"
16 
17 
18 #ifdef HAVE_NTL
19 #include "NTLconvert.h"
20 #endif
21 
22 #ifdef HAVE_FLINT
23 #include "FLINTconvert.h"
24 #endif
25 
26 #include "cf_assert.h"
27 #include "debug.h"
28 
29 #include "canonicalform.h"
30 #include "cf_iter.h"
31 #include "cf_util.h"
32 
33 
34 /** void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
35  *
36  * chineseRemainder - integer chinese remaindering.
37  *
38  * Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
39  * and qnew = q1*q2. q1 and q2 should be positive integers,
40  * pairwise prime, x1 and x2 should be polynomials with integer
41  * coefficients. If x1 and x2 are polynomials with positive
42  * coefficients, the result is guaranteed to have positive
43  * coefficients, too.
44  *
45  * Note: This algorithm is optimized for the case q1>>q2.
46  *
47  * This is a standard algorithm. See, for example,
48  * Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
49  * par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
50  * Homomorphic Images' in B. Buchberger - 'Computer Algebra -
51  * Symbolic and Algebraic Computation'.
52  *
53  * Note: Be sure you are calculating in Z, and not in Q!
54  *
55 **/
56 void
57 chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
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 }
103 //}}}
104 
105 /** void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
106  *
107  * chineseRemainder - integer chinese remaindering.
108  *
109  * Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
110  * product of all q[i]. q[i] should be positive integers,
111  * pairwise prime. x[i] should be polynomials with integer
112  * coefficients. If all coefficients of all x[i] are positive
113  * integers, the result is guaranteed to have positive
114  * coefficients, too.
115  *
116  * This is a standard algorithm, too, except for the fact that we
117  * use a divide-and-conquer method instead of a linear approach
118  * to calculate the remainder.
119  *
120  * Note: Be sure you are calculating in Z, and not in Q!
121  *
122 **/
123 void
124 chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
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 }
165 
166 #if !defined(HAVE_NTL) && !defined(HAVE_FLINT)
167 CanonicalForm Farey_n (CanonicalForm N, const CanonicalForm P)
168 //"USAGE: Farey_n (N,P); P, N number;
169 //RETURN: a rational number a/b such that a/b=N mod P
170 // and |a|,|b|<(P/2)^{1/2}
171 {
172  //assume(P>0);
173  // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
174  if (N<0) N +=P;
175  CanonicalForm A,B,C,D,E;
176  E=P;
177  B=1;
178  while (!N.isZero())
179  {
180  if (2*N*N<P)
181  {
182  On(SW_RATIONAL);
183  N /=B;
184  Off(SW_RATIONAL);
185  return(N);
186  }
187  D=mod(E , N);
188  C=A-(E-mod(E , N))/N*B;
189  E=N;
190  N=D;
191  A=B;
192  B=C;
193  }
194  return(0);
195 }
196 #endif
197 
198 /**
199  * Farey rational reconstruction. If NTL is available it uses the fast algorithm
200  * from NTL, i.e. Encarnacion, Collins.
201 **/
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 }
274 
275 // returns x where (a * x) % b == 1, inv is a cache
276 static inline CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
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 }
288 
289 void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
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 }
306 // http://rosettacode.org/wiki/Chinese_remainder_theorem#C
307 
308 void chineseRemainderCached ( const CanonicalForm & a, const CanonicalForm &q1, const CanonicalForm & b, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew,CFArray &inv )
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 }
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.)
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:670
Conversion to and from NTL.
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4084
int p
Definition: cfModGcd.cc:4080
CanonicalForm b
Definition: cfModGcd.cc:4105
assertions for Factory
#define ASSERT(expression, message)
Definition: cf_assert.h:99
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:202
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
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:276
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
cf_algorithm.cc - simple mathematical algorithms.
Definition: cf_factor.cc:99
void chineseRemainderCached(const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
Definition: cf_chinese.cc:290
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:30
Iterators for CanonicalForm's.
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
int size() const
Definition: ftmpl_array.cc:92
int min() const
Definition: ftmpl_array.cc:98
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
bool inCoeffDomain() const
bool inZ() const
predicates
int ilog2() const
int CanonicalForm::ilog2 () const
factory's class for variables
Definition: factory.h:134
functions to print debug output
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
REvaluation E(1, terms.length(), IntRandom(25))
b *CanonicalForm B
Definition: facBivar.cc:52
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
bool isZero(const CFArray &A)
checks if entries of A are zero
#define D(A)
Definition: gentable.cc:131
STATIC_VAR jList * Q
Definition: janet.cc:30
static int sign(int x)
Definition: ring.cc:3380
#define A
Definition: sirandom.c:24