Changeset edb4893 in git


Ignore:
Timestamp:
Jun 3, 1996, 10:32:56 AM (28 years ago)
Author:
Rüdiger Stobbe <stobbe@…>
Branches:
(u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
Children:
f59b88fb5075f45eedd00f61e2cb4e7c27bacdce
Parents:
3a37ecd689343f0dd008afb140c8ef3fe9474428
Message:
"gcd_poly: now uses new function gcd_poly_univar0 to compute univariate
          polynomial gcd's over Z.
gcd_poly_univar0: computes univariate polynomial gcd's in characteristic 0
                  via chinese remaindering.
maxnorm: computes the maximum norm of all coefficients of a polynomial.
"


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

Legend:

Unmodified
Added
Removed
  • factory/cf_gcd.cc

    r3a37ec redb4893  
    11// emacs edit mode for this file is -*- C++ -*-
    2 // $Id: cf_gcd.cc,v 1.0 1996-05-17 11:56:37 stobbe Exp $
     2// $Id: cf_gcd.cc,v 1.1 1996-06-03 08:32:56 stobbe Exp $
    33
    44/*
    55$Log: not supported by cvs2svn $
     6Revision 1.0  1996/05/17 11:56:37  stobbe
     7Initial revision
     8
    69*/
    710
     
    1114#include "cf_iter.h"
    1215#include "cf_reval.h"
    13 
    14 static void indent ( int i )
    15 {
    16     for ( int j = 0; j < i; j++ )
    17         cerr << "   ";
     16#include "cf_primes.h"
     17#include "cf_chinese.h"
     18#include "templates/functions.h"
     19
     20static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
     21
     22
     23static int
     24isqrt ( int a )
     25{
     26    int h, x0, x1 = a;
     27    do {
     28        x0 = x1;
     29        h = x0 * x0 + a - 1;
     30        if ( h % (2 * x0) == 0 )
     31            x1 = h / (2 * x0);
     32        else
     33            x1 = (h - 1)  / (2 * x0);
     34    } while ( x1 < x0 );
     35    return x1;
    1836}
    1937
     
    4159}
    4260
     61static CanonicalForm
     62maxnorm ( const CanonicalForm & f )
     63{
     64    CanonicalForm m = 0, h;
     65    CFIterator i;
     66    for ( i = f; i.hasTerms(); i++ )
     67        m = tmax( m, abs( i.coeff() ) );
     68    return m;
     69}
     70
     71static void
     72chinesePoly ( const CanonicalForm & f1, const CanonicalForm & q1, const CanonicalForm & f2, const CanonicalForm & q2, CanonicalForm & f, CanonicalForm & q )
     73{
     74    CFIterator i1 = f1, i2 = f2;
     75    CanonicalForm c;
     76    Variable x = f1.mvar();
     77    f = 0;
     78    while ( i1.hasTerms() && i2.hasTerms() ) {
     79        if ( i1.exp() == i2.exp() ) {
     80            chineseRemainder( i1.coeff(), q1, i2.coeff(), q2, c, q );
     81            f += power( x, i1.exp() ) * c;
     82            i1++; i2++;
     83        }
     84        else if ( i1.exp() > i2.exp() ) {
     85            chineseRemainder( 0, q1, i2.coeff(), q2, c, q );
     86            f += power( x, i2.exp() ) * c;
     87            i2++;
     88        }
     89        else {
     90            chineseRemainder( i1.coeff(), q1, 0, q2, c, q );
     91            f += power( x, i1.exp() ) * c;
     92            i1++;
     93        }
     94    }
     95    while ( i1.hasTerms() ) {
     96        chineseRemainder( i1.coeff(), q1, 0, q2, c, q );
     97        f += power( x, i1.exp() ) * c;
     98        i1++;
     99    }
     100    while ( i2.hasTerms() ) {
     101        chineseRemainder( 0, q1, i2.coeff(), q2, c, q );
     102        f += power( x, i2.exp() ) * c;
     103        i2++;
     104    }
     105}
     106
     107static CanonicalForm
     108balance ( const CanonicalForm & f, const CanonicalForm & q )
     109{
     110    CFIterator i;
     111    CanonicalForm result = 0, qh = q / 2;
     112    Variable x = f.mvar();
     113    for ( i = f; i.hasTerms(); i++ ) {
     114        if ( i.coeff() > qh )
     115            result += power( x, i.exp() ) * (i.coeff() - q);
     116        else
     117            result += power( x, i.exp() ) * i.coeff();
     118    }
     119    return result;
     120}
     121
    43122CanonicalForm
    44123igcd ( const CanonicalForm & f, const CanonicalForm & g )
     
    132211}
    133212
    134 CanonicalForm
    135 gcd_poly( const CanonicalForm & f, const CanonicalForm& g )
     213static CanonicalForm
     214gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
     215{
     216    CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
     217    int p, i, n;
     218
     219    if ( primitive ) {
     220        f = F;
     221        g = G;
     222        c = 1;
     223    }
     224    else {
     225        CanonicalForm cF = content( F ), cG = content( G );
     226        f = F / cF;
     227        g = G / cG;
     228        c = igcd( cF, cG );
     229    }
     230    cg = gcd( f.lc(), g.lc() );
     231    cl = ( f.lc() / cg ) * g.lc();
     232//     B = 2 * cg * tmin(
     233//      maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
     234//      maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
     235//      )+1;
     236    M = tmin( maxnorm(f), maxnorm(g) );
     237    BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
     238    q = 0;
     239    i = 1;
     240    n = cf_getNumBigPrimes();
     241    while ( true ) {
     242        B = BB;
     243        while ( i < n && q < B ) {
     244            p = cf_getBigPrime( i );
     245            i++;
     246            while ( i < n && mod( cl, p ) == 0 ) {
     247                p = cf_getBigPrime( i );
     248                i++;
     249            }
     250            setCharacteristic( p );
     251            Dp = gcd( mapinto( f ), mapinto( g ) );
     252            Dp = ( Dp / Dp.lc() ) * mapinto( cg );
     253            setCharacteristic( 0 );
     254            if ( Dp.degree() == 0 ) return c;
     255            if ( q.isZero() ) {
     256                D = mapinto( Dp );
     257                q = p;
     258                B = power(CanonicalForm(2),D.degree())*M+1;
     259            }
     260            else {
     261                if ( Dp.degree() == D.degree() ) {
     262                    chinesePoly( D, q, mapinto( Dp ), p, newD, newq );
     263                    q = newq;
     264                    D = newD;
     265                }
     266                else if ( Dp.degree() < D.degree() ) {
     267                    // all previous p's are bad primes
     268                    q = p;
     269                    D = mapinto( Dp );
     270                    B = power(CanonicalForm(2),D.degree())*M+1;
     271                }
     272                // else p is a bad prime
     273            }
     274        }
     275        if ( i < n ) {
     276            // now balance D mod q
     277            D = pp( balance( cg * D, q ) );
     278            if ( divides( D, f ) && divides( D, g ) )
     279                return D * c;
     280            else
     281                q = 0;
     282        }
     283        else {
     284            return gcd_poly( F, G, false );
     285        }
     286    }
     287}
     288
     289
     290static CanonicalForm
     291gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag )
    136292{
    137293    CanonicalForm C, Ci, Ci1, Hi, bi, pi, pi1, pi2;
     
    148304    C = gcd( Ci, Ci1 );
    149305    pi1 = pi1 / Ci1; pi = pi / Ci;
    150     if ( ! pi.isUnivariate() )
     306    if ( ! pi.isUnivariate() ) {
    151307        if ( gcd_test_one( pi1, pi ) )
    152308            return C;
     309    }
     310    else {
     311        // pi is univariate
     312        if ( modularflag )
     313            return gcd_poly_univar0( pi, pi1, true ) * C;
     314    }
    153315    delta = degree( pi, v ) - degree( pi1, v );
    154316    Hi = power( LC( pi1, v ), delta );
     
    261423                CanonicalForm F = f * l, G = g * l;
    262424                Off( SW_RATIONAL );
    263                 l = gcd_poly( F, G );
     425                l = gcd_poly( F, G, true );
    264426                On( SW_RATIONAL );
    265427                return l;
    266428            }
    267429            else
    268                 return gcd_poly( f, g );
     430                return gcd_poly( f, g, getCharacteristic()==0 );
    269431        }
    270432    else  if ( f.mvar() > g.mvar() )
Note: See TracChangeset for help on using the changeset viewer.