source: git/factory/cf_chinese.cc @ 16f511

spielwiese
Last change on this file since 16f511 was 16f511, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fixed the usage of "config.h" (if defined HAVE_CONFIG_H)
  • Property mode set to 100644
File size: 6.8 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3//{{{ docu
4//
5// cf_chinese.cc - algorithms for chinese remaindering.
6//
7// Used by: cf_gcd.cc, cf_linsys.cc
8//
9// Header file: cf_algorithm.h
10//
11//}}}
12
13#ifdef HAVE_CONFIG_H
14#include "config.h"
15#endif /* HAVE_CONFIG_H */
16
17#ifdef HAVE_NTL
18#include "NTLconvert.h"
19#endif
20
21#include "cf_assert.h"
22#include "debug.h"
23
24#include "canonicalform.h"
25#include "cf_iter.h"
26
27
28//{{{ void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
29//{{{ docu
30//
31// chineseRemainder - integer chinese remaindering.
32//
33// Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
34// and qnew = q1*q2.  q1 and q2 should be positive integers,
35// pairwise prime, x1 and x2 should be polynomials with integer
36// coefficients.  If x1 and x2 are polynomials with positive
37// coefficients, the result is guaranteed to have positive
38// coefficients, too.
39//
40// Note: This algorithm is optimized for the case q1>>q2.
41//
42// This is a standard algorithm.  See, for example,
43// Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
44// par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
45// Homomorphic Images' in B. Buchberger - 'Computer Algebra -
46// Symbolic and Algebraic Computation'.
47//
48// Note: Be sure you are calculating in Z, and not in Q!
49//
50//}}}
51void
52chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
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}
98//}}}
99
100//{{{ void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
101//{{{ docu
102//
103// chineseRemainder - integer chinese remaindering.
104//
105// Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
106// product of all q[i].  q[i] should be positive integers,
107// pairwise prime.  x[i] should be polynomials with integer
108// coefficients.  If all coefficients of all x[i] are positive
109// integers, the result is guaranteed to have positive
110// coefficients, too.
111//
112// This is a standard algorithm, too, except for the fact that we
113// use a divide-and-conquer method instead of a linear approach
114// to calculate the remainder.
115//
116// Note: Be sure you are calculating in Z, and not in Q!
117//
118//}}}
119void
120chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
121{
122    DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
123
124    ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
125    CFArray X(x), Q(q);
126    int i, j, n = x.size(), start = x.min();
127
128    DEBOUTLN( cerr, "array size = " << n );
129
130    while ( n != 1 )
131    {
132        i = j = start;
133        while ( i < start + n - 1 )
134        {
135            // This is a little bit dangerous: X[i] and X[j] (and
136            // Q[i] and Q[j]) may refer to the same object.  But
137            // xnew and qnew in the above function are modified
138            // at the very end of the function, so we do not
139            // modify x1 and q1, resp., by accident.
140            chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
141            i += 2;
142            j++;
143        }
144
145        if ( n & 1 )
146        {
147            X[j] = X[i];
148            Q[j] = Q[i];
149        }
150        // Maybe we would get some memory back at this point if
151        // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
152        // at this point?
153
154        n = ( n + 1) / 2;
155    }
156    xnew = X[start];
157    qnew = Q[q.min()];
158
159    DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
160}
161//}}}
162
163CanonicalForm Farey_n (CanonicalForm N, const CanonicalForm P)
164//"USAGE:  Farey_n (N,P); P, N number;
165//RETURN:  a rational number a/b such that a/b=N mod P
166//         and |a|,|b|<(P/2)^{1/2}
167{
168   //assume(P>0);
169   // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
170   if (N<0) N +=P;
171   CanonicalForm A,B,C,D,E;
172   E=P;
173   B=1;
174   while (!N.isZero())
175   {
176        if (2*N*N<P)
177        {
178           On(SW_RATIONAL);
179           N /=B;
180           Off(SW_RATIONAL);
181           return(N);
182        }
183        D=mod(E , N);
184        C=A-(E-mod(E , N))/N*B;
185        E=N;
186        N=D;
187        A=B;
188        B=C;
189   }
190   return(0);
191}
192
193CanonicalForm Farey ( const CanonicalForm & f, const CanonicalForm & q )
194{
195    int is_rat=isOn(SW_RATIONAL);
196    Off(SW_RATIONAL);
197    Variable x = f.mvar();
198    CanonicalForm result = 0;
199    CanonicalForm c;
200    CFIterator i;
201#ifdef HAVE_NTL
202    ZZ NTLq= convertFacCF2NTLZZ (q);
203    ZZ bound;
204    SqrRoot (bound, NTLq/2);
205#endif
206    for ( i = f; i.hasTerms(); i++ )
207    {
208        c = i.coeff();
209        if ( c.inCoeffDomain())
210        {
211#ifdef HAVE_NTL
212          if (c.inZ() && isOn (SW_USE_NTL))
213          {
214            ZZ NTLc= convertFacCF2NTLZZ (c);
215            bool lessZero= (sign (NTLc) == -1);
216            if (lessZero)
217              NTL::negate (NTLc, NTLc);
218            ZZ NTLnum, NTLden;
219            if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
220            {
221              if (lessZero)
222                NTL::negate (NTLnum, NTLnum);
223              CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
224              CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
225              On (SW_RATIONAL);
226              result += power (x, i.exp())*(num/den);
227              Off (SW_RATIONAL);
228            }
229          }
230          else
231            result += power( x, i.exp() ) * Farey(c,q);
232#else
233          if (c.inZ())
234            result += power( x, i.exp() ) * Farey_n(c,q);
235          else
236            result += power( x, i.exp() ) * Farey(c,q);
237#endif
238        }
239        else
240          result += power( x, i.exp() ) * Farey(c,q);
241    }
242    if (is_rat) On(SW_RATIONAL);
243    return result;
244}
245
Note: See TracBrowser for help on using the repository browser.