My Project
Loading...
Searching...
No Matches
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#ifdef HAVE_FLINT
23#include "FLINTconvert.h"
24#endif
25
26#include "cf_assert.h"
27#include "debug.h"
28
29#include "canonicalform.h"
30#include "cf_iter.h"
31#include "cf_util.h"
32
33
34/** void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
35 *
36 * chineseRemainder - integer chinese remaindering.
37 *
38 * Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
39 * and qnew = q1*q2. q1 and q2 should be positive integers,
40 * pairwise prime, x1 and x2 should be polynomials with integer
41 * coefficients. If x1 and x2 are polynomials with positive
42 * coefficients, the result is guaranteed to have positive
43 * coefficients, too.
44 *
45 * Note: This algorithm is optimized for the case q1>>q2.
46 *
47 * This is a standard algorithm. See, for example,
48 * Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
49 * par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
50 * Homomorphic Images' in B. Buchberger - 'Computer Algebra -
51 * Symbolic and Algebraic Computation'.
52 *
53 * Note: Be sure you are calculating in Z, and not in Q!
54 *
55**/
56void
57chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
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}
103//}}}
104
105/** void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
106 *
107 * chineseRemainder - integer chinese remaindering.
108 *
109 * Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
110 * product of all q[i]. q[i] should be positive integers,
111 * pairwise prime. x[i] should be polynomials with integer
112 * coefficients. If all coefficients of all x[i] are positive
113 * integers, the result is guaranteed to have positive
114 * coefficients, too.
115 *
116 * This is a standard algorithm, too, except for the fact that we
117 * use a divide-and-conquer method instead of a linear approach
118 * to calculate the remainder.
119 *
120 * Note: Be sure you are calculating in Z, and not in Q!
121 *
122**/
123void
124chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
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}
165
166#if !defined(HAVE_NTL) && !defined(HAVE_FLINT)
168//"USAGE: Farey_n (N,P); P, N number;
169//RETURN: a rational number a/b such that a/b=N mod P
170// and |a|,|b|<(P/2)^{1/2}
171{
172 //assume(P>0);
173 // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
174 if (N<0) N +=P;
175 CanonicalForm A,B,C,D,E;
176 E=P;
177 B=1;
178 while (!N.isZero())
179 {
180 if (2*N*N<P)
181 {
183 N /=B;
185 return(N);
186 }
187 D=mod(E , N);
188 C=A-(E-mod(E , N))/N*B;
189 E=N;
190 N=D;
191 A=B;
192 B=C;
193 }
194 return(0);
195}
196#endif
197
198/**
199 * Farey rational reconstruction. If NTL is available it uses the fast algorithm
200 * from NTL, i.e. Encarnacion, Collins.
201**/
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}
274
275// returns x where (a * x) % b == 1, inv is a cache
276static inline CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
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}
288
289void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
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}
306// http://rosettacode.org/wiki/Chinese_remainder_theorem#C
307
308void chineseRemainderCached ( const CanonicalForm & a, const CanonicalForm &q1, const CanonicalForm & b, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew,CFArray &inv )
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}
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.)
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:670
Conversion to and from NTL.
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm bextgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm bextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CanonicalForm num(const CanonicalForm &f)
CanonicalForm den(const CanonicalForm &f)
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int i
Definition: cfEzgcd.cc:132
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
CanonicalForm b
Definition: cfModGcd.cc:4103
assertions for Factory
#define ASSERT(expression, message)
Definition: cf_assert.h:99
CanonicalForm Farey(const CanonicalForm &f, const CanonicalForm &q)
Farey rational reconstruction.
Definition: cf_chinese.cc:202
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
static CanonicalForm chin_mul_inv(const CanonicalForm a, const CanonicalForm b, int ind, CFArray &inv)
Definition: cf_chinese.cc:276
void out_cf(const char *s1, const CanonicalForm &f, const char *s2)
cf_algorithm.cc - simple mathematical algorithms.
Definition: cf_factor.cc:99
void chineseRemainderCached(const CFArray &a, const CFArray &n, CanonicalForm &xnew, CanonicalForm &prod, CFArray &inv)
Definition: cf_chinese.cc:290
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
Iterators for CanonicalForm's.
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
int size() const
Definition: ftmpl_array.cc:92
int min() const
Definition: ftmpl_array.cc:98
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
bool inCoeffDomain() const
bool inZ() const
predicates
int ilog2() const
int CanonicalForm::ilog2 () const
factory's class for variables
Definition: variable.h:33
functions to print debug output
#define DEBINCLEVEL(stream, msg)
Definition: debug.h:44
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
#define DEBDECLEVEL(stream, msg)
Definition: debug.h:45
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
REvaluation E(1, terms.length(), IntRandom(25))
b *CanonicalForm B
Definition: facBivar.cc:52
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
bool isZero(const CFArray &A)
checks if entries of A are zero
#define D(A)
Definition: gentable.cc:131
static int sign(int x)
Definition: ring.cc:3427
#define A
Definition: sirandom.c:24
#define Q
Definition: sirandom.c:26