# source:git/factory/cf_chinese.cc@498648

spielwiese
Last change on this file since 498648 was 498648, checked in by Martin Lee <martinlee84@…>, 12 years ago
rational reconstruction via NTL git-svn-id: file:///usr/local/Singular/svn/trunk@14308 2c84dea3-7e68-4137-9b89-c4e89433aadc
• 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, sm_util.cc
9//
10// Header file: cf_algorithm.h
11//
12//}}}
13
14#include <config.h>
15
16#include "assert.h"
17#include "debug.h"
18
19#include "canonicalform.h"
20#include "cf_iter.h"
21
22#ifdef HAVE_NTL
23#include "NTLconvert.h"
24#endif
25
26//{{{ void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
27//{{{ docu
28//
29// chineseRemainder - integer chinese remaindering.
30//
31// Calculate xnew such that xnew=x1 (mod q1) and xnew=x2 (mod q2)
32// and qnew = q1*q2.  q1 and q2 should be positive integers,
33// pairwise prime, x1 and x2 should be polynomials with integer
34// coefficients.  If x1 and x2 are polynomials with positive
35// coefficients, the result is guaranteed to have positive
36// coefficients, too.
37//
38// Note: This algorithm is optimized for the case q1>>q2.
39//
40// This is a standard algorithm.  See, for example,
41// Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
42// par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
43// Homomorphic Images' in B. Buchberger - 'Computer Algebra -
44// Symbolic and Algebraic Computation'.
45//
46// Note: Be sure you are calculating in Z, and not in Q!
47//
48//}}}
49void
50chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
51{
52    DEBINCLEVEL( cerr, "chineseRemainder" );
53
54    DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
55    DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
56
57    // We calculate xnew as follows:
58    //     xnew = v1 + v2 * q1
59    // where
60    //     v1 = x1 (mod q1)
61    //     v2 = (x2-v1)/q1 (mod q2)  (*)
62    //
63    // We do one extra test to check whether x2-v1 vanishes (mod
64    // q2) in (*) since it is not costly and may save us
65    // from calculating the inverse of q1 (mod q2).
66    //
67    // u: v1 (mod q2)
68    // d: x2-v1 (mod q2)
69    // s: 1/q1 (mod q2)
70    //
71    CanonicalForm v2, v1;
72    CanonicalForm u, d, s, dummy;
73
74    v1 = mod( x1, q1 );
75    u = mod( v1, q2 );
76    d = mod( x2-u, q2 );
77    if ( d.isZero() )
78    {
79        xnew = v1;
80        qnew = q1 * q2;
81        DEBDECLEVEL( cerr, "chineseRemainder" );
82        return;
83    }
84    (void)bextgcd( q1, q2, s, dummy );
85    v2 = mod( d*s, q2 );
86    xnew = v1 + v2*q1;
87
88    // After all, calculate new modulus.  It is important that
89    // this is done at the very end of the algorithm, since q1
90    // and qnew may refer to the same object (same is true for x1
91    // and xnew).
92    qnew = q1 * q2;
93
94    DEBDECLEVEL( cerr, "chineseRemainder" );
95}
96//}}}
97
98//{{{ void chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
99//{{{ docu
100//
101// chineseRemainder - integer chinese remaindering.
102//
103// Calculate xnew such that xnew=x[i] (mod q[i]) and qnew is the
104// product of all q[i].  q[i] should be positive integers,
105// pairwise prime.  x[i] should be polynomials with integer
106// coefficients.  If all coefficients of all x[i] are positive
107// integers, the result is guaranteed to have positive
108// coefficients, too.
109//
110// This is a standard algorithm, too, except for the fact that we
111// use a divide-and-conquer method instead of a linear approach
112// to calculate the remainder.
113//
114// Note: Be sure you are calculating in Z, and not in Q!
115//
116//}}}
117void
118chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
119{
120    DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
121
122    ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
123    CFArray X(x), Q(q);
124    int i, j, n = x.size(), start = x.min();
125
126    DEBOUTLN( cerr, "array size = " << n );
127
128    while ( n != 1 )
129    {
130        i = j = start;
131        while ( i < start + n - 1 )
132        {
133            // This is a little bit dangerous: X[i] and X[j] (and
134            // Q[i] and Q[j]) may refer to the same object.  But
135            // xnew and qnew in the above function are modified
136            // at the very end of the function, so we do not
137            // modify x1 and q1, resp., by accident.
138            chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
139            i += 2;
140            j++;
141        }
142
143        if ( n & 1 )
144        {
145            X[j] = X[i];
146            Q[j] = Q[i];
147        }
148        // Maybe we would get some memory back at this point if
149        // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
150        // at this point?
151
152        n = ( n + 1) / 2;
153    }
154    xnew = X[start];
155    qnew = Q[q.min()];
156
157    DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
158}
159//}}}
160
161CanonicalForm Farey_n (CanonicalForm N, const CanonicalForm P)
162//"USAGE:  Farey_n (N,P); P, N number;
163//RETURN:  a rational number a/b such that a/b=N mod P
164//         and |a|,|b|<(P/2)^{1/2}
165{
166   //assume(P>0);
167   // assume !isOn(SW_RATIONAL): mod is a no-op otherwise
168   if (N<0) N +=P;
169   CanonicalForm A,B,C,D,E;
170   E=P;
171   B=1;
172   while (!N.isZero())
173   {
174        if (2*N*N<P)
175        {
176           On(SW_RATIONAL);
177           N /=B;
178           Off(SW_RATIONAL);
179           return(N);
180        }
181        D=mod(E , N);
182        C=A-(E-mod(E , N))/N*B;
183        E=N;
184        N=D;
185        A=B;
186        B=C;
187   }
188   return(0);
189}
190
191CanonicalForm Farey ( const CanonicalForm & f, const CanonicalForm & q )
192{
193    int is_rat=isOn(SW_RATIONAL);
194    Off(SW_RATIONAL);
195    Variable x = f.mvar();
196    CanonicalForm result = 0;
197    CanonicalForm c;
198    CFIterator i;
199    for ( i = f; i.hasTerms(); i++ )
200    {
201        c = i.coeff();
202        if ( c.inCoeffDomain())
203        {
204#ifdef HAVE_NTL
205          if (c.inZ() && isOn (SW_USE_NTL))
206          {
207            ZZ NTLc= convertFacCF2NTLZZ (c);
208            bool lessZero= (sign (NTLc) == -1);
209            if (lessZero)
210              negate (NTLc, NTLc);
211            ZZ NTLq= convertFacCF2NTLZZ (q);
212            ZZ bound;
213            SqrRoot (bound, NTLq/2);
214            ZZ NTLnum, NTLden;
215            if (ReconstructRational (NTLnum, NTLden, NTLc, NTLq, bound, bound))
216            {
217              if (lessZero)
218                negate (NTLnum, NTLnum);
219              CanonicalForm num= convertNTLZZX2CF (to_ZZX (NTLnum), Variable (1));
220              CanonicalForm den= convertNTLZZX2CF (to_ZZX (NTLden), Variable (1));
221              On (SW_RATIONAL);
222              result += power (x, i.exp())*(num/den);
223              Off (SW_RATIONAL);
224            }
225          }
226          else
227#endif
228            result += power( x, i.exp() ) * Farey_n(c,q);
229        }
230        else
231          result += power( x, i.exp() ) * Farey(c,q);
232    }
233    if (is_rat) On(SW_RATIONAL);
234    return result;
235}
236
Note: See TracBrowser for help on using the repository browser.