Changeset 4e6908 in git for factory/fac_ezgcd.cc


Ignore:
Timestamp:
Jan 24, 2001, 7:09:41 PM (23 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
d83977461ed45e8d1baf489f2bd39e88bd603aa3
Parents:
4b24dc87e783c30424d898c8d71f88b84cb61baf
Message:
*hannes: even mor bug fixes to ezgcd


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

Legend:

Unmodified
Added
Removed

  • factory/fac_ezgcd.cc

    r4b24dc8 r4e6908  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: fac_ezgcd.cc,v 1.13 2001-01-23 17:41:52 Singular Exp $ */
     2/* $Id: fac_ezgcd.cc,v 1.14 2001-01-24 18:09:41 Singular Exp $ */
    33
    44#include <config.h>
    … …  
    5454    F = FF / f; G = GG / g;
    5555    if ( F.isUnivariate() && G.isUnivariate() )
    56         return d * gcd( F, G );
     56        return d * gcd( F, G );
    5757    else  if ( gcd_test_one( F, G, false ) )
    58         return d;
     58        return d;
    5959    lcF = LC( F, x ); lcG = LC( G, x );
    6060    lcD = gcd( lcF, lcG );
    … …  
    6464    bound = findBound( F, G, lcF, lcG, degF, degG );
    6565    if ( ! internal )
    66         b = REvaluation( 2, t, IntRandom( 100 ) );
     66        b = REvaluation( 2, t, IntRandom( 100 ) );
    6767    while ( ! gcdfound ) {
    68         /// ---> A2
    69         DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
    70         DEBOUTLN( cerr, "F = " << F );
    71         DEBOUTLN( cerr, "G = " << G );
    72         findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
    73         DEBOUTLN( cerr, "found evaluation b = " << b );
    74         DEBOUTLN( cerr, "F(b) = " << Fb );
    75         DEBOUTLN( cerr, "G(b) = " << Gb );
    76         DEBOUTLN( cerr, "D(b) = " << Db );
    77         delta = degree( Db );
    78         /// ---> A3
    79         if ( delta == 0 )
    80             return d;
    81         /// ---> A4
    82         deltaold = delta + 1;
    83         while ( deltaold != delta ) {
    84             bt = b;
    85             findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
    86             if ( degree( Dbt ) == 0 )
    87                 return d;
    88             if ( degree( Dbt ) == delta )
    89                 deltaold = delta;
    90             else  if ( degree( Dbt ) < delta ) {
    91                 deltaold = delta;
    92                 delta = degree( Dbt );
    93                 b = bt;
    94                 Db = Dbt; Fb = Fbt; Gb = Gbt;
     68        /// ---> A2
     69        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
     70        DEBOUTLN( cerr, "F = " << F );
     71        DEBOUTLN( cerr, "G = " << G );
     72        findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
     73        DEBOUTLN( cerr, "found evaluation b = " << b );
     74        DEBOUTLN( cerr, "F(b) = " << Fb );
     75        DEBOUTLN( cerr, "G(b) = " << Gb );
     76        DEBOUTLN( cerr, "D(b) = " << Db );
     77        delta = degree( Db );
     78        /// ---> A3
     79        if ( delta == 0 )
     80            return d;
     81        /// ---> A4
     82        deltaold = delta + 1;
     83        while ( deltaold != delta ) {
     84            bt = b;
     85            findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
     86            if ( degree( Dbt ) == 0 )
     87                return d;
     88            if ( degree( Dbt ) == delta )
     89                deltaold = delta;
     90            else  if ( degree( Dbt ) < delta ) {
     91                deltaold = delta;
     92                delta = degree( Dbt );
     93                b = bt;
     94                Db = Dbt; Fb = Fbt; Gb = Gbt;
     95            }
     96        }
     97        DEBOUTLN( cerr, "now after A4, delta = " << delta );
     98        /// ---> A5
     99        if ( degF <= degG && delta == degF && divides( F, G ) )
     100            return d*F;
     101        if ( degG < degF && delta == degG && divides( G, F ) )
     102            return d*G;
     103        if ( delta != degF && delta != degG ) {
     104            /// ---> A6
     105            CanonicalForm xxx;
     106            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
     107            xxx= gcd( (DD[1] = Fb / Db), Db );
     108            if (xxx.inCoeffDomain()) {
     109                B = F;
     110                DD[2] = Db;
     111                lcDD[1] = lcF;
     112                lcDD[2] = lcF;
     113                B *= lcF;
     114            }
     115            //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
     116            else {
     117            xxx=gcd( (DD[1] = Gb / Db), Db );
     118            if (xxx.inCoeffDomain()) {
     119                B = G;
     120                DD[2] = Db;
     121                lcDD[1] = lcG;
     122                lcDD[2] = lcG;
     123                B *= lcG;
     124            }
     125            else {
     126#ifdef DEBUGOUTPUT
     127                CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
     128                DEBDECLEVEL( cerr, "ezgcd" );
     129                return dummyres;
     130#else
     131                return d * ezgcd_specialcase( F, G, b, bound );
     132#endif
     133            }
    95134            }
    96         }
    97         DEBOUTLN( cerr, "now after A4, delta = " << delta );
    98         /// ---> A5
    99         if ( degF <= degG && delta == degF && divides( F, G ) )
    100             return d*F;
    101         if ( degG < degF && delta == degG && divides( G, F ) )
    102             return d*G;
    103         if ( delta != degF && delta != degG ) {
    104             /// ---> A6
    105             CanonicalForm xxx;
    106             if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
    107                 B = F;
    108                 DD[2] = Db;
    109                 lcDD[1] = lcF;
    110                 lcDD[2] = lcF;
    111                 B *= lcF;
    112             }
    113             else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
    114                 B = G;
    115                 DD[2] = Db;
    116                 lcDD[1] = lcG;
    117                 lcDD[2] = lcG;
    118                 B *= lcG;
    119             }
    120             else {
    121 #ifdef DEBUGOUTPUT
    122                 CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
    123                 DEBDECLEVEL( cerr, "ezgcd" );
    124                 return dummyres;
    125 #else
    126                 return d * ezgcd_specialcase( F, G, b, bound );
    127 #endif
    128             }
    129             /// ---> A7
    130             DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    131             DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    132             bound = enlargeBound( B, DD[2], DD[1], bound );
    133             DEBOUTLN( cerr, "(hensel) B    = " << B );
    134             DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
    135             DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
    136             DEBOUTLN( cerr, "(hensel) DD   = " << DD );
    137             DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
    138             gcdfound = Hensel( B, DD, lcDD, b, bound, x );
    139             DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
    140             /// ---> A8 (gcdfound)
    141         }
     135            /// ---> A7
     136            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     137            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     138            bound = enlargeBound( B, DD[2], DD[1], bound );
     139            DEBOUTLN( cerr, "(hensel) B    = " << B );
     140            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
     141            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
     142            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
     143            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
     144            gcdfound = Hensel( B, DD, lcDD, b, bound, x );
     145            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
     146            /// ---> A8 (gcdfound)
     147        }
    142148    }
    143149    /// ---> A9
    … …  
    152158    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt, d;
    153159    CFArray DD( 1, 2 ), lcDD( 1, 2 );
    154     Variable x = F.mvar();
     160    Variable x = Variable( 1 );
    155161    bool gcdfound;
    156162
    … …  
    160166    Ft = ezgcd( F, F.deriv( x ) );
    161167    if ( Ft.isOne() ) {
    162         // In this case F is squarefree and we came here by bad chance
    163         // (means: bad evaluation point).  Try again with another
    164         // evaluation point.  Bug fix (?) by JS.  The bad example was
     168        // In this case F is squarefree and we came here by bad chance
     169        // (means: bad evaluation point).  Try again with another
     170        // evaluation point.  Bug fix (?) by JS.  The bad example was
    165171        // gcd.debug -ocr /+USE_EZGCD/@12/CB \
    166         //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)' \
     172        //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)' \
    167173        //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
    168         b.nextpoint();
     174        b.nextpoint();
    169175        return ezgcd( F, G, b, true );
    170176    }
    171            
     177
    172178    DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
    173179    L = F / Ft;
    … …  
    182188    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
    183189    while ( true ) {
    184         /// ---> S3
    185         DEBOUTLN( cerr, "ezgcdspec: (S3)" );
    186         DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
    187         DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
    188         DD[1] = gcd( Lb, gcd( Fb, Gb ) );
    189         /// ---> S4
    190         DEBOUTLN( cerr, "ezgcdspec: (S4)" );
    191         if ( degree( DD[1] ) == 0 )
    192             return Ds;
    193         /// ---> S5
    194         DEBOUTLN( cerr, "ezgcdspec: (S5)" );
    195         DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
    196         DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
    197         Db = DD[1];
    198         if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
    199             DEBOUTLN( cerr, "ezgcdspec: (S55)" );
    200             lcDD[2] = LC( L, x );
    201             lcDD[1] = d;
    202             DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
    203             DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
    204             LL = L * d;
    205             modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
    206             DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
    207             DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
    208             DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
    209             DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
    210             DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
    211             DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
    212             DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
    213             DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
    214             gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
    215             ASSERT( gcdfound, "fatal error in algorithm" );
    216             Dt = pp( DD[1] );
    217         }
    218         DEBOUTLN( cerr, "ezgcdspec: (S7)" );
    219         Ds = Ds * Dt;
    220         Fb = Fb / Db;
    221         Gb = Gb / Db;
     190        /// ---> S3
     191        DEBOUTLN( cerr, "ezgcdspec: (S3)" );
     192        DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
     193        DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
     194        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
     195        /// ---> S4
     196        DEBOUTLN( cerr, "ezgcdspec: (S4)" );
     197        if ( degree( DD[1] ) == 0 )
     198            return Ds;
     199        /// ---> S5
     200        DEBOUTLN( cerr, "ezgcdspec: (S5)" );
     201        DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
     202        DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
     203        Db = DD[1];
     204        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
     205            DEBOUTLN( cerr, "ezgcdspec: (S55)" );
     206            lcDD[2] = LC( L, x );
     207            lcDD[1] = d;
     208            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
     209            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
     210            LL = L * d;
     211            modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
     212            DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
     213            DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
     214            DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
     215            DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
     216            DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
     217            DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
     218            DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
     219            DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
     220            gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
     221            ASSERT( gcdfound, "fatal error in algorithm" );
     222            Dt = pp( DD[1] );
     223        }
     224        DEBOUTLN( cerr, "ezgcdspec: (S7)" );
     225        Ds = Ds * Dt;
     226        Fb = Fb / Db;
     227        Gb = Gb / Db;
    222228    }
    223229    // this point is never reached, but to avoid compiler warnings let's return a value
    … …  
    230236    bool ok;
    231237    if ( delta != 0 )
    232         b.nextpoint();
     238        b.nextpoint();
    233239    do {
    234         DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
    235         Fb = b( F );
    236         ok = degree( Fb ) == degF;
    237         if ( ok ) {
    238             Gb = b( G );
    239             ok = degree( Gb ) == degG;
    240         }
    241         if ( ok ) {
    242             Db = gcd( Fb, Gb );
    243             if ( delta > 0 )
    244                 ok = degree( Db ) < delta;
    245         }
    246         if ( ! ok )
    247             b.nextpoint();
     240        DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
     241        Fb = b( F );
     242        ok = degree( Fb ) == degF;
     243        if ( ok ) {
     244            Gb = b( G );
     245            ok = degree( Gb ) == degG;
     246        }
     247        if ( ok ) {
     248            Db = gcd( Fb, Gb );
     249            if ( delta > 0 )
     250                ok = degree( Db ) < delta;
     251        }
     252        if ( ! ok )
     253            b.nextpoint();
    248254    } while ( ! ok );
    249255}
    … …  
    259265
    260266    CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
    261         tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
     267        tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
    262268    int p = pk.getp();
    263269    int k = pk.getk();
    264270    CanonicalForm bound = pk.getpk();
    265271    while ( bound < limit ) {
    266         k++;
    267         bound *= p;
     272        k++;
     273        bound *= p;
    268274    }
    269275    k *= 2;
    … …  
    276282{
    277283    CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
    278         gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
     284        gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
    279285    int p, i = 0;
    280286    do {
    281         p = cf_getBigPrime( i );
    282         i++;
     287        p = cf_getBigPrime( i );
     288        i++;
    283289    } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
    284290    CanonicalForm bound = p;
    285291    i = 1;
    286292    while ( bound < limit ) {
    287         i++;
    288         bound *= p;
     293        i++;
     294        bound *= p;
    289295    }
    290296    return modpk( p, i );
Note: See TracChangeset for help on using the changeset viewer.