My Project
Loading...
Searching...
No Matches
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 )
 
void chineseRemainder (const CFArray &x, const CFArray &q, CanonicalForm &xnew, CanonicalForm &qnew)
 void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
 
CanonicalForm Farey (const CanonicalForm &f, const CanonicalForm &q)
 Farey rational reconstruction.
 
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.
 
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:4111
factory's main class
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}
Array< CanonicalForm > CFArray
int i
Definition cfEzgcd.cc:132
Variable x
Definition cfModGcd.cc:4090
#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
int min() const
int j
Definition facHensel.cc:110
#define Q
Definition sirandom.c:26

◆ 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)
#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:4086
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
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);
206 Variable x = f.mvar();
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);
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)
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
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.
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
static int sign(int x)
Definition ring.cc:3503

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

104{
105 printf("%s",s1);
106 if (f.isZero()) printf("+0");
107 //else if (! f.inCoeffDomain() )
108 else if (! f.inBaseDomain() )
109 {
110 int l = f.level();
111 for ( CFIterator i = f; i.hasTerms(); i++ )
112 {
113 int e=i.exp();
114 if (i.coeff().isOne())
115 {
116 printf("+");
117 if (e==0) printf("1");
118 else
119 {
120 printf("%c",'a'+l-1);
121 if (e!=1) printf("^%d",e);
122 }
123 }
124 else
125 {
126 out_cf("+(",i.coeff(),")");
127 if (e!=0)
128 {
129 printf("*%c",'a'+l-1);
130 if (e!=1) printf("^%d",e);
131 }
132 }
133 }
134 }
135 else
136 {
137 if ( f.isImm() )
138 {
140 {
141 long a= imm2int (f.getval());
142 if ( a == gf_q )
143 printf ("+%ld", a);
144 else if ( a == 0L )
145 printf ("+1");
146 else if ( a == 1L )
147 printf ("+%c",gf_name);
148 else
149 {
150 printf ("+%c",gf_name);
151 printf ("^%ld",a);
152 }
153 }
154 else
155 {
156 long l=f.intval();
157 if (l<0) printf("%ld",l);
158 else printf("+%ld",l);
159 }
160 }
161 else
162 {
163 #ifdef NOSTREAMIO
164 if (f.inZ())
165 {
166 mpz_t m;
168 char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
169 str = mpz_get_str( str, 10, m );
170 puts(str);
171 delete[] str;
172 mpz_clear(m);
173 }
174 else if (f.inQ())
175 {
176 mpz_t m;
178 char * str = new char[mpz_sizeinbase( m, 10 ) + 2];
179 str = mpz_get_str( str, 10, m );
180 while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
181 puts(str);putchar('/');
182 delete[] str;
183 mpz_clear(m);
185 str = new char[mpz_sizeinbase( m, 10 ) + 2];
186 str = mpz_get_str( str, 10, m );
187 while(str[strlen(str)]<' ') { str[strlen(str)]='\0'; }
188 puts(str);
189 delete[] str;
190 mpz_clear(m);
191 }
192 #else
193 std::cout << f;
194 #endif
195 }
196 //if (f.inZ()) printf("(Z)");
197 //else if (f.inQ()) printf("(Q)");
198 //else if (f.inFF()) printf("(FF)");
199 //else if (f.inPP()) printf("(PP)");
200 //else if (f.inGF()) printf("(PP)");
201 //else
202 if (f.inExtension()) printf("E(%d)",f.level());
203 }
204 printf("%s",s2);
205}
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:103
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:699