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 #include "cf_assert.h"
23 #include "debug.h"
24 
25 #include "canonicalform.h"
26 #include "cf_iter.h"
27 
28 
29 /** void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
30  *
31  * chineseRemainder - integer chinese remaindering.
32  *
33  * Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
34  * and qnew = q1*q2. q1 and q2 should be positive integers,
35  * pairwise prime, x1 and x2 should be polynomials with integer
36  * coefficients. If x1 and x2 are polynomials with positive
37  * coefficients, the result is guaranteed to have positive
38  * coefficients, too.
39  *
40  * Note: This algorithm is optimized for the case q1>>q2.
41  *
42  * This is a standard algorithm. See, for example,
43  * Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
44  * par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
45  * Homomorphic Images' in B. Buchberger - 'Computer Algebra -
46  * Symbolic and Algebraic Computation'.
47  *
48  * Note: Be sure you are calculating in Z, and not in Q!
49  *
50 **/
51 void
52 chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
53 {
54  DEBINCLEVEL( cerr, "chineseRemainder" );
55 
56  DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
57  DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
58 
59  // We calculate xnew as follows:
60  // xnew = v1 + v2 * q1
61  // where
62  // v1 = x1 (mod q1)
63  // v2 = (x2-v1)/q1 (mod q2) (*)
64  //
65  // We do one extra test to check whether x2-v1 vanishes (mod
66  // q2) in (*) since it is not costly and may save us
67  // from calculating the inverse of q1 (mod q2).
68  //
69  // u: v1 (mod q2)
70  // d: x2-v1 (mod q2)
71  // s: 1/q1 (mod q2)
72  //
73  CanonicalForm v2, v1;
74  CanonicalForm u, d, s, dummy;
75 
76  v1 = mod( x1, q1 );
77  u = mod( v1, q2 );
78  d = mod( x2-u, q2 );
79  if ( d.isZero() )
80  {
81  xnew = v1;
82  qnew = q1 * q2;
83  DEBDECLEVEL( cerr, "chineseRemainder" );
84  return;
85  }
86  (void)bextgcd( q1, q2, s, dummy );
87  v2 = mod( d*s, q2 );
88  xnew = v1 + v2*q1;
89 
90  // After all, calculate new modulus. It is important that
91  // this is done at the very end of the algorithm, since q1
92  // and qnew may refer to the same object (same is true for x1
93  // and xnew).
94  qnew = q1 * q2;
95 
96  DEBDECLEVEL( cerr, "chineseRemainder" );
97 }
98 //}}}
99 
100 /** void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
101  *
102  * chineseRemainder - integer chinese remaindering.
103  *
104  * Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
105  * product of all q[i]. q[i] should be positive integers,
106  * pairwise prime. x[i] should be polynomials with integer
107  * coefficients. If all coefficients of all x[i] are positive
108  * integers, the result is guaranteed to have positive
109  * coefficients, too.
110  *
111  * This is a standard algorithm, too, except for the fact that we
112  * use a divide-and-conquer method instead of a linear approach
113  * to calculate the remainder.
114  *
115  * Note: Be sure you are calculating in Z, and not in Q!
116  *
117 **/
118 void
119 chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
120 {
121  DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
122 
123  ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
124  CFArray X(x), Q(q);
125  int i, j, n = x.size(), start = x.min();
126 
127  DEBOUTLN( cerr, "array size = " << n );
128 
129  while ( n != 1 )
130  {
131  i = j = start;
132  while ( i < start + n - 1 )
133  {
134  // This is a little bit dangerous: X[i] and X[j] (and
135  // Q[i] and Q[j]) may refer to the same object. But
136  // xnew and qnew in the above function are modified
137  // at the very end of the function, so we do not
138  // modify x1 and q1, resp., by accident.
139  chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
140  i += 2;
141  j++;
142  }
143 
144  if ( n & 1 )
145  {
146  X[j] = X[i];
147  Q[j] = Q[i];
148  }
149  // Maybe we would get some memory back at this point if
150  // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
151  // at this point?
152 
153  n = ( n + 1) / 2;
154  }
155  xnew = X[start];
156  qnew = Q[q.min()];
157 
158  DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
159 }
160 
161 #ifndef HAVE_NTL
163 //"USAGE: Farey_n (N,P); P, N number;
164 //RETURN: a rational number a/b such that a/b=N mod P
165 // and |a|,|b|<(P/2)^{1/2}
166 {
167  //assume(P>0);
168  // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
169  if (N<0) N +=P;
170  CanonicalForm A,B,C,D,E;
171  E=P;
172  B=1;
173  while (!N.isZero())
174  {
175  if (2*N*N<P)
176  {
177  On(SW_RATIONAL);
178  N /=B;
179  Off(SW_RATIONAL);
180  return(N);
181  }
182  D=mod(E , N);
183  C=A-(E-mod(E , N))/N*B;
184  E=N;
185  N=D;
186  A=B;
187  B=C;
188  }
189  return(0);
190 }
191 #endif
192 
193 /**
194  * Farey rational reconstruction. If NTL is available it uses the fast algorithm
195  * from NTL, i.e. Encarnacion, Collins.
196 **/
198 {
199  int is_rat=isOn(SW_RATIONAL);
200  Off(SW_RATIONAL);
201  Variable x = f.mvar();
202  CanonicalForm result = 0;
203  CanonicalForm c;
204  CFIterator i;
205 #ifdef HAVE_NTL
206  ZZ NTLq= convertFacCF2NTLZZ (q);
207  ZZ bound;
208  SqrRoot (bound, NTLq/2);
209 #endif
210  for ( i = f; i.hasTerms(); i++ )
211  {
212  c = i.coeff();
213  if ( c.inCoeffDomain())
214  {
215 #ifdef HAVE_NTL
216  if (c.inZ())
217  {
218  ZZ NTLc= convertFacCF2NTLZZ (c);
219  bool lessZero= (sign (NTLc) == -1);
220  if (lessZero)
221  NTL::negate (NTLc, NTLc);
222  ZZ NTLnum, NTLden;
223  if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
224  {
225  if (lessZero)
226  NTL::negate (NTLnum, NTLnum);
227  CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
228  CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
229  On (SW_RATIONAL);
230  result += power (x, i.exp())*(num/den);
231  Off (SW_RATIONAL);
232  }
233  }
234  else
235  result += power( x, i.exp() ) * Farey(c,q);
236 #else
237  if (c.inZ())
238  result += power( x, i.exp() ) * Farey_n(c,q);
239  else
240  result += power( x, i.exp() ) * Farey(c,q);
241 #endif
242  }
243  else
244  result += power( x, i.exp() ) * Farey(c,q);
245  }
246  if (is_rat) On(SW_RATIONAL);
247  return result;
248 }
249 
250 // returns x where (a * x) % b == 1, inv is a cache
252 {
253  if (inv[ind].isZero())
254  {
255  CanonicalForm s,dummy;
256  (void)bextgcd( a, b, s, dummy );
257  inv[ind]=s;
258  return s;
259  }
260  else
261  return inv[ind];
262 }
263 
265 {
266  CanonicalForm p, sum = 0; prod=1;
267  int i;
268  int len=n.size();
269 
270  for (i = 0; i < len; i++) prod *= n[i];
271 
272  for (i = 0; i < len; i++)
273  {
274  p = prod / n[i];
275  sum += a[i] * chin_mul_inv(p, n[i], i, inv) * p;
276  }
277 
278  xnew = mod(sum , prod);
279 }
280 // http://rosettacode.org/wiki/Chinese_remainder_theorem#C
281 
int ilog2() const
int CanonicalForm::ilog2 () const
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
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:52
const CanonicalForm int s
Definition: facAbsFact.cc:55
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:667
#define D(A)
Definition: gentable.cc:123
Conversion to and from NTL.
const poly a
Definition: syzextra.cc:212
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
void Off(int sw)
switches
CanonicalForm num(const CanonicalForm &f)
functions to print debug output
return P p
Definition: myNF.cc:203
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
factory&#39;s class for variables
Definition: factory.h:115
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a...
factory&#39;s main class
Definition: canonicalform.h:75
assertions for Factory
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
#define Q
Definition: sirandom.c:25
void chineseRemainderCached(CFArray &a, CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
Definition: cf_chinese.cc:264
int min() const
Definition: ftmpl_array.cc:98
int size() const
Definition: ftmpl_array.cc:92
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:197
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:49
int j
Definition: myNF.cc:70
#define A
Definition: sirandom.c:23
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:28
bool isOn(int sw)
switches
void On(int sw)
switches
Iterators for CanonicalForm&#39;s.
FILE * f
Definition: checklibs.c:9
int i
Definition: cfEzgcd.cc:123
class to iterate through CanonicalForm&#39;s
Definition: cf_iter.h:44
REvaluation E(1, terms.length(), IntRandom(25))
CanonicalForm den(const CanonicalForm &f)
bool inZ() const
predicates
b *CanonicalForm B
Definition: facBivar.cc:51
fq_nmod_poly_t prod
Definition: facHensel.cc:95
static CanonicalForm chin_mul_inv(CanonicalForm a, CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:251
Variable x
Definition: cfModGcd.cc:4023
bool isZero(const CFArray &A)
checks if entries of A are zero
#define ASSERT(expression, message)
Definition: cf_assert.h:99
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
kBucketDestroy & P
Definition: myNF.cc:191
const poly b
Definition: syzextra.cc:213
static int sign(int x)
Definition: ring.cc:3342
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:281
return result
Definition: facAbsBiFact.cc:76
Header for factory&#39;s main class CanonicalForm.
bool inCoeffDomain() const