source: git/factory/cf_chinese.cc @ 806c18

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