Changeset d57f19 in git for factory/cf_chinese.cc


Ignore:
Timestamp:
Feb 2, 1998, 9:57:29 AM (26 years ago)
Author:
Jens Schmidt <schmidt@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
f3b98ad1c11ea4bf6d3b9525f75480b2abf56cac
Parents:
9dac21e20830c8d04a26496a2d4e9328a8770887
Message:
	* cf_chinese.cc (chineseRemainder( const CanonicalForm & ...)):
	  completely rewritten

	* cf_chinese.cc (chineseRemainder): debug output added

	* cf_chinese.cc	(chineseRemainder( const CFArray & ...)): slightly
	  speeded up

	* cf_gcd.cc (iextgcd, igcd): functions removed.  All references
	  replaced by `bextgcd()' and `bgcd()', resp.


git-svn-id: file:///usr/local/Singular/svn/trunk@1112 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_chinese.cc

    r9dac21e rd57f19  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_chinese.cc,v 1.8 1997-11-19 17:42:15 schmidt Exp $ */
     2/* $Id: cf_chinese.cc,v 1.9 1998-02-02 08:57:29 schmidt Exp $ */
    33
    44//{{{ docu
    … …  
    1515
    1616#include "assert.h"
     17#include "debug.h"
    1718
    1819#include "canonicalform.h"
    … …  
    3031// coefficients, too.
    3132//
     33// Note: This algorithm is optimized for the case q1>>q2.
     34//
    3235// This is a standard algorithm.  See, for example,
    33 // Geddes/Czapor/Labahn - 'Alogorithms for Computer Algebra',
    34 // par. 5.6 and 5.8.
     36// Geddes/Czapor/Labahn - 'Algorithms for Computer Algebra',
     37// par. 5.6 and 5.8, or the article of M. Lauer - 'Computing by
     38// Homomorphic Images' in B. Buchberger - 'Computer Algebra -
     39// Symbolic and Algebraic Computation'.
    3540//
    36 // Note: be sure you are calculating in Z, and not in Q!
     41// Note: Be sure you are calculating in Z, and not in Q!
    3742//
    3843//}}}
    … …  
    4045chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew )
    4146{
    42     CanonicalForm a1, a2;
    43     (void)iextgcd( q1, q2, a1, a2 );
     47    DEBINCLEVEL( cerr, "chineseRemainder" );
     48
     49    DEBOUTLN( cerr, "log(q1) = " << q1.ilog2() );
     50    DEBOUTLN( cerr, "log(q2) = " << q2.ilog2() );
     51
     52    // We calculate xnew as follows:
     53    //     xnew = v1 + v2 * q1
     54    // where
     55    //     v1 = x1 (mod q1)
     56    //     v2 = (x2-v1)/q1 (mod q2)  (*)
     57    //
     58    // We do one extra test to check whether x2-v1 vanishes (mod
     59    // q2) in (*) since it is not costly and may save us
     60    // from calculating the inverse of q1 (mod q2).
     61    //
     62    // u: v1 (mod q2)
     63    // d: x2-v1 (mod q2)
     64    // s: 1/q1 (mod q2)
     65    //
     66    CanonicalForm v2, v1;
     67    CanonicalForm u, d, s, dummy;
     68
     69    v1 = mod( x1, q1 );
     70    u = mod( v1, q2 );
     71    d = mod( x2-u, q2 );
     72    if ( d.isZero() ) {
     73        xnew = v1;
     74        qnew = q1 * q2;
     75        return;
     76    }
     77    (void)bextgcd( q1, q2, s, dummy );
     78    v2 = mod( d*s, q2 );
     79    xnew = v1 + v2*q1;
     80
     81    // After all, calculate new modulus.  It is important that
     82    // this is done at the very end of the algorithm, since q1
     83    // and qnew may refer to the same object (same is true for x1
     84    // and xnew).
    4485    qnew = q1 * q2;
    45     xnew = mod( q1*a1*x2 + q2*a2*x1, qnew );
     86
     87    DEBDECLEVEL( cerr, "chineseRemainder" );
    4688}
    4789//}}}
    … …  
    63105// to calculate the remainder.
    64106//
    65 // Note: be sure you are calculating in Z, and not in Q!
     107// Note: Be sure you are calculating in Z, and not in Q!
    66108//
    67109//}}}
    … …  
    69111chineseRemainder ( const CFArray & x, const CFArray & q, CanonicalForm & xnew, CanonicalForm & qnew )
    70112{
     113    DEBINCLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
     114
    71115    ASSERT( x.min() == q.min() && x.size() == q.size(), "incompatible arrays" );
    72116    CFArray X(x), Q(q);
    73     int i, j, n = x.size();
     117    int i, j, n = x.size(), start = x.min();
     118
     119    DEBOUTLN( cerr, "array size = " << n );
    74120
    75121    while ( n != 1 ) {
    76         i = j = x.min();
    77         while ( i < x.min() + n - 1 ) {
     122        i = j = start;
     123        while ( i < start + n - 1 ) {
     124            // This is a little bit dangerous: X[i] and X[j] (and
     125            // Q[i] and Q[j]) may refer to the same object.  But
     126            // xnew and qnew in the above function are modified
     127            // at the very end of the function, so we do not
     128            // modify x1 and q1, resp., by accident.
    78129            chineseRemainder( X[i], Q[i], X[i+1], Q[i+1], X[j], Q[j] );
    79130            i += 2;
    80131            j++;
    81132        }
     133
    82134        if ( n & 1 ) {
    83135            X[j] = X[i];
    84136            Q[j] = Q[i];
    85137        }
     138        // Maybe we would get some memory back at this point if
     139        // we would set X[j+1, ..., n] and Q[j+1, ..., n] to zero
     140        // at this point?
     141
    86142        n = ( n + 1) / 2;
    87143    }
    88     xnew = X[x.min()];
     144    xnew = X[start];
    89145    qnew = Q[q.min()];
     146
     147    DEBDECLEVEL( cerr, "chineseRemainder( ... CFArray ... )" );
    90148}
    91149//}}}
Note: See TracChangeset for help on using the changeset viewer.