My Project  UNKNOWN_GIT_VERSION
Functions
cf_chinese.cc File Reference
#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 out_cf (const char *s1, const CanonicalForm &f, const char *s2)
 cf_algorithm.cc - simple mathematical algorithms. More...
 
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

◆ chin_mul_inv()

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

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 }

◆ 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 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 }

◆ 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 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 }

◆ chineseRemainderCached()

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

Definition at line 265 of file cf_chinese.cc.

266 {
267  CanonicalForm p, sum=0L; prod=1L;
268  int i;
269  int len=n.size();
270 
271  for (i = 0; i < len; i++) prod *= n[i];
272 
273  for (i = 0; i < len; i++)
274  {
275  p = prod / n[i];
276  sum += a[i] * chin_mul_inv(p, n[i], i, inv) * p;
277  }
278 
279  xnew = mod(sum , prod);
280 }

◆ Farey()

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 }

◆ 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 90 of file cf_factor.cc.

91 {
92  printf("%s",s1);
93  if (f.isZero()) printf("+0");
94  //else if (! f.inCoeffDomain() )
95  else if (! f.inBaseDomain() )
96  {
97  int l = f.level();
98  for ( CFIterator i = f; i.hasTerms(); i++ )
99  {
100  int e=i.exp();
101  if (i.coeff().isOne())
102  {
103  printf("+");
104  if (e==0) printf("1");
105  else
106  {
107  printf("v(%d)",l);
108  if (e!=1) printf("^%d",e);
109  }
110  }
111  else
112  {
113  out_cf("+(",i.coeff(),")");
114  if (e!=0)
115  {
116  printf("*v(%d)",l);
117  if (e!=1) printf("^%d",e);
118  }
119  }
120  }
121  }
122  else
123  {
124  if ( f.isImm() )
125  {
127  {
128  long a= imm2int (f.getval());
129  if ( a == gf_q )
130  printf ("+%ld", a);
131  else if ( a == 0L )
132  printf ("+1");
133  else if ( a == 1L )
134  printf ("+%c",gf_name);
135  else
136  {
137  printf ("+%c",gf_name);
138  printf ("^%ld",a);
139  }
140  }
141  else
142  printf("+%ld",f.intval());
143  }
144  else
145  {
146  #ifdef NOSTREAMIO
147  if (f.inZ())
148  {
149  mpz_t m;
150  gmp_numerator(f,m);
151  char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
152  str = mpz_get_str( str, 10, m );
153  puts(str);
154  delete[] str;
155  mpz_clear(m);
156  }
157  else if (f.inQ())
158  {
159  mpz_t m;
160  gmp_numerator(f,m);
161  char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
162  str = mpz_get_str( str, 10, m );
163  puts(str);putchar('/');
164  delete[] str;
165  mpz_clear(m);
167  str = new char[mpz_sizeinbase( m, 10 ) + 2];
168  str = mpz_get_str( str, 10, m );
169  puts(str);
170  delete[] str;
171  mpz_clear(m);
172  }
173  #else
174  std::cout << f;
175  #endif
176  }
177  //if (f.inZ()) printf("(Z)");
178  //else if (f.inQ()) printf("(Q)");
179  //else if (f.inFF()) printf("(FF)");
180  //else if (f.inPP()) printf("(PP)");
181  //else if (f.inGF()) printf("(PP)");
182  //else
183  if (f.inExtension()) printf("E(%d)",f.level());
184  }
185  printf("%s",s2);
186 }
SW_RATIONAL
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:28
isOn
bool isOn(int sw)
switches
Definition: canonicalform.cc:1912
isZero
bool isZero(const CFArray &A)
checks if entries of A are zero
Definition: facSparseHensel.h:468
j
int j
Definition: facHensel.cc:105
f
FILE * f
Definition: checklibs.c:9
CFIterator
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
x
Variable x
Definition: cfModGcd.cc:4023
DEBOUTLN
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
result
return result
Definition: facAbsBiFact.cc:76
num
CanonicalForm num(const CanonicalForm &f)
Definition: canonicalform.h:330
prod
fq_nmod_poly_t prod
Definition: facHensel.cc:95
CFFactory::gettype
static int gettype()
Definition: cf_factory.h:28
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition: canonicalform.cc:1837
sign
static int sign(int x)
Definition: ring.cc:3346
chin_mul_inv
static CanonicalForm chin_mul_inv(CanonicalForm a, CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:251
mod
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
gf_name
char gf_name
Definition: gfops.cc:52
CanonicalForm::ilog2
int ilog2() const
int CanonicalForm::ilog2 () const
Definition: canonicalform.cc:1352
b
CanonicalForm b
Definition: cfModGcd.cc:4044
CanonicalForm
factory's main class
Definition: canonicalform.h:83
GaloisFieldDomain
#define GaloisFieldDomain
Definition: cf_defs.h:22
Array::min
int min() const
Definition: ftmpl_array.cc:98
i
int i
Definition: cfEzgcd.cc:125
Array
Definition: ftmpl_array.h:17
out_cf
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
cf_algorithm.cc - simple mathematical algorithms.
Definition: cf_factor.cc:90
ASSERT
#define ASSERT(expression, message)
Definition: cf_assert.h:99
Array::size
int size() const
Definition: ftmpl_array.cc:92
chineseRemainder
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
DEBINCLEVEL
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
DEBDECLEVEL
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
gf_q
int gf_q
Definition: gfops.cc:47
den
CanonicalForm den(const CanonicalForm &f)
Definition: canonicalform.h:333
Off
void Off(int sw)
switches
Definition: canonicalform.cc:1905
CanonicalForm::inZ
bool inZ() const
predicates
Definition: canonicalform.cc:66
CanonicalForm::inCoeffDomain
bool inCoeffDomain() const
Definition: canonicalform.cc:119
imm2int
static long imm2int(const InternalCF *const imm)
Definition: imm.h:70
Variable
factory's class for variables
Definition: factory.h:118
bound
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
Farey
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:197
m
int m
Definition: cfEzgcd.cc:121
l
int l
Definition: cfEzgcd.cc:93
gmp_denominator
void gmp_denominator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:40
gmp_numerator
void gmp_numerator(const CanonicalForm &f, mpz_ptr result)
Definition: singext.cc:20
convertNTLZZX2CF
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:278
p
int p
Definition: cfModGcd.cc:4019
s
const CanonicalForm int s
Definition: facAbsFact.cc:55
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
Q
#define Q
Definition: sirandom.c:25
convertFacCF2NTLZZ
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:664
bextgcd
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition: canonicalform.cc:1663
On
void On(int sw)
switches
Definition: canonicalform.cc:1898