My Project
Loading...
Searching...
No Matches
Functions
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 ) 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 (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. More...
 
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:4103
factory's main class
Definition: canonicalform.h:86
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}
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
#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
Definition: ftmpl_array.cc:92
int min() const
Definition: ftmpl_array.cc:98
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)
Definition: cf_chinese.cc:290
#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:4078
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:276
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)
Definition: NTLconvert.cc:285
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:670
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.
Definition: cf_chinese.cc:202
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: variable.h:33
return result
Definition: facAbsBiFact.cc:75
static int sign(int x)
Definition: ring.cc:3427

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

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