source: git/factory/cf_chinese.cc @ e88604

jengelh-datetimespielwiese
Last change on this file since e88604 was e88604, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: deleted sparse modular gcd from sm_sparsemod.*
  • Property mode set to 100644
File size: 6.7 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4//{{{ docu
5//
6// cf_chinese.cc - algorithms for chinese remaindering.
7//
8// Used by: cf_gcd.cc, cf_linsys.cc
9//
10// Header file: cf_algorithm.h
11//
12//}}}
13
14#include "config.h"
15
16#ifdef HAVE_NTL
17#include "NTLconvert.h"
18#endif
19
20#include "cf_assert.h"
21#include "debug.h"
22
23#include "canonicalform.h"
24#include "cf_iter.h"
25
26
27//{{{ void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
28//{{{ docu
29//
30// chineseRemainder - integer chinese remaindering.
31//
32// Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
33// and qnew = q1*q2.  q1 and q2 should be positive integers,
34// pairwise prime, x1 and x2 should be polynomials with integer
35// coefficients.  If x1 and x2 are polynomials with positive
36// coefficients, the result is guaranteed to have positive
37// coefficients, too.
38//
39// Note: This algorithm is optimized for the case q1>>q2.
40//
41// This is a standard algorithm.  See, for example,
42// Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
43// par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
44// Homomorphic Images' in B. Buchberger - 'Computer Algebra -
45// Symbolic and Algebraic Computation'.
46//
47// Note: Be sure you are calculating in Z, and not in Q!
48//
49//}}}
50void
51chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
52{
53    DEBINCLEVEL( cerr, "chineseRemainder" );
54
55    DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
56    DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
57
58    // We calculate xnew as follows:
59    //     xnew = v1 + v2 * q1
60    // where
61    //     v1 = x1 (mod q1)
62    //     v2 = (x2-v1)/q1 (mod q2)  (*)
63    //
64    // We do one extra test to check whether x2-v1 vanishes (mod
65    // q2) in (*) since it is not costly and may save us
66    // from calculating the inverse of q1 (mod q2).
67    //
68    // u: v1 (mod q2)
69    // d: x2-v1 (mod q2)
70    // s: 1/q1 (mod q2)
71    //
72    CanonicalForm v2, v1;
73    CanonicalForm u, d, s, dummy;
74
75    v1 = mod( x1, q1 );
76    u = mod( v1, q2 );
77    d = mod( x2-u, q2 );
78    if ( d.isZero() )
79    {
80        xnew = v1;
81        qnew = q1 * q2;
82        DEBDECLEVEL( cerr, "chineseRemainder" );
83        return;
84    }
85    (void)bextgcd( q1, q2, s, dummy );
86    v2 = mod( d*s, q2 );
87    xnew = v1 + v2*q1;
88
89    // After all, calculate new modulus.  It is important that
90    // this is done at the very end of the algorithm, since q1
91    // and qnew may refer to the same object (same is true for x1
92    // and xnew).
93    qnew = q1 * q2;
94
95    DEBDECLEVEL( cerr, "chineseRemainder" );
96}
97//}}}
98
99//{{{ void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
100//{{{ docu
101//
102// chineseRemainder - integer chinese remaindering.
103//
104// Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
105// product of all q[i].  q[i] should be positive integers,
106// pairwise prime.  x[i] should be polynomials with integer
107// coefficients.  If all coefficients of all x[i] are positive
108// integers, the result is guaranteed to have positive
109// coefficients, too.
110//
111// This is a standard algorithm, too, except for the fact that we
112// use a divide-and-conquer method instead of a linear approach
113// to calculate the remainder.
114//
115// Note: Be sure you are calculating in Z, and not in Q!
116//
117//}}}
118void
119chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
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}
160//}}}
161
162CanonicalForm Farey_n (CanonicalForm N, const CanonicalForm P)
163//"USAGE:  Farey_n (N,P); P, N number;
164//RETURN:  a rational number a/b such that a/b=N mod P
165//         and |a|,|b|<(P/2)^{1/2}
166{
167   //assume(P>0);
168   // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
169   if (N<0) N +=P;
170   CanonicalForm A,B,C,D,E;
171   E=P;
172   B=1;
173   while (!N.isZero())
174   {
175        if (2*N*N<P)
176        {
177           On(SW_RATIONAL);
178           N /=B;
179           Off(SW_RATIONAL);
180           return(N);
181        }
182        D=mod(E , N);
183        C=A-(E-mod(E , N))/N*B;
184        E=N;
185        N=D;
186        A=B;
187        B=C;
188   }
189   return(0);
190}
191
192CanonicalForm Farey ( const CanonicalForm & f, const CanonicalForm & q )
193{
194    int is_rat=isOn(SW_RATIONAL);
195    Off(SW_RATIONAL);
196    Variable x = f.mvar();
197    CanonicalForm result = 0;
198    CanonicalForm c;
199    CFIterator i;
200#ifdef HAVE_NTL
201    ZZ NTLq= convertFacCF2NTLZZ (q);
202    ZZ bound;
203    SqrRoot (bound, NTLq/2);
204#endif
205    for ( i = f; i.hasTerms(); i++ )
206    {
207        c = i.coeff();
208        if ( c.inCoeffDomain())
209        {
210#ifdef HAVE_NTL
211          if (c.inZ() && isOn (SW_USE_NTL))
212          {
213            ZZ NTLc= convertFacCF2NTLZZ (c);
214            bool lessZero= (sign (NTLc) == -1);
215            if (lessZero)
216              NTL::negate (NTLc, NTLc);
217            ZZ NTLnum, NTLden;
218            if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
219            {
220              if (lessZero)
221                NTL::negate (NTLnum, NTLnum);
222              CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
223              CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
224              On (SW_RATIONAL);
225              result += power (x, i.exp())*(num/den);
226              Off (SW_RATIONAL);
227            }
228          }
229          else
230#endif
231            result += power( x, i.exp() ) * Farey_n(c,q);
232        }
233        else
234          result += power( x, i.exp() ) * Farey(c,q);
235    }
236    if (is_rat) On(SW_RATIONAL);
237    return result;
238}
239
Note: See TracBrowser for help on using the repository browser.