Functions
cf_chinese.cc File Reference

algorithms for chinese remaindering and rational reconstruction More...

#include "config.h"
#include "NTLconvert.h"
#include "cf_assert.h"
#include "debug.h"
#include "canonicalform.h"
#include "cf_iter.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 (CanonicalForm a, CanonicalForm b, int ind, CFArray &inv)
 
void chineseRemainderCached (CFArray &a, CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, 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

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

Definition at line 251 of file cf_chinese.cc.

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 }
const CanonicalForm int s
Definition: facAbsFact.cc:55
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a...
factory's main class
Definition: canonicalform.h:75
bool isZero(const CFArray &A)
checks if entries of A are zero
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 52 of file cf_chinese.cc.

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 }
const CanonicalForm int s
Definition: facAbsFact.cc:55
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
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
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
int ilog2() const
int CanonicalForm::ilog2 () const
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
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 119 of file cf_chinese.cc.

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 }
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
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
int size() const
Definition: ftmpl_array.cc:92
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
#define Q
Definition: sirandom.c:25
int j
Definition: myNF.cc:70
int i
Definition: cfEzgcd.cc:123
int min() const
Definition: ftmpl_array.cc:98
#define ASSERT(expression, message)
Definition: cf_assert.h:99
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
void chineseRemainderCached ( CFArray a,
CFArray n,
CanonicalForm xnew,
CanonicalForm prod,
CFArray inv 
)

Definition at line 264 of file cf_chinese.cc.

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 }
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
return P p
Definition: myNF.cc:203
int size() const
Definition: ftmpl_array.cc:92
factory&#39;s main class
Definition: canonicalform.h:75
int i
Definition: cfEzgcd.cc:123
static CanonicalForm chin_mul_inv(CanonicalForm a, CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:251

Farey rational reconstruction.

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

Definition at line 197 of file cf_chinese.cc.

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 }
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:666
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
void Off(int sw)
switches
bool inCoeffDomain() const
CanonicalForm num(const CanonicalForm &f)
factory&#39;s class for variables
Definition: variable.h:32
factory&#39;s main class
Definition: canonicalform.h:75
bool inZ() const
predicates
CF_NO_INLINE int hasTerms() const
check if iterator has reached < the end of CanonicalForm
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:197
CF_NO_INLINE CanonicalForm coeff() const
get the current coefficient
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
int i
Definition: cfEzgcd.cc:123
class to iterate through CanonicalForm&#39;s
Definition: cf_iter.h:44
CF_NO_INLINE int exp() const
get the current exponent
CanonicalForm den(const CanonicalForm &f)
Variable x
Definition: cfModGcd.cc:4023
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:281
return result
Definition: facAbsBiFact.cc:76
int sign(const CanonicalForm &a)