Changeset c244d04 in git for factory/fac_ezgcd.cc


Ignore:
Timestamp:
Apr 24, 2012, 1:01:15 PM (11 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
Children:
9f84ad2f3376850c4d78697edb203d32c07d3094
Parents:
f7a4e9271997208f4ee40a14fd1e1b5c6c0c1ffc
git-author:
Martin Lee <martinlee84@web.de>2012-04-24 13:01:15+02:00
git-committer:
Martin Lee <martinlee84@web.de>2012-05-07 14:15:53+02:00
Message:
chg: deleted unused old functions + editing chg: added description
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/fac_ezgcd.cc

    rf7a4e9 rc244d04  
    1 /* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id$ */
     1/*****************************************************************************\
     2 * Computer Algebra System SINGULAR
     3\*****************************************************************************/
     4/** @file fac_ezgcd.cc
     5 *
     6 * This file implements the GCD of two multivariate polynomials over Q using
     7 * EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor,
     8 * Labahnn
     9 *
     10 * @author Martin Lee
     11 *
     12 **/
     13/*****************************************************************************/
    314
    415#include "config.h"
    … …  
    920#include "cf_defs.h"
    1021#include "canonicalform.h"
    11 #include "fac_util.h"
     22#include "fac_util.h" // HEADER for this file
    1223#include "cf_algorithm.h"
    1324#include "cf_reval.h"
    … …  
    331342}
    332343
    333 static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG );
    334 
    335 static CanonicalForm ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal );
    336 
    337 static CanonicalForm ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound );
    338 
    339 static modpk findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG );
    340 
    341 static modpk enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk );
     344static CanonicalForm
     345ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
     346        bool internal )
     347{
     348  bool isRat= isOn (SW_RATIONAL);
     349  if (!isRat)
     350    On (SW_RATIONAL);
     351  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
     352                lcD, cand, result;
     353  CFArray DD( 1, 2 ), lcDD( 1, 2 );
     354  int degF, degG, delta, t;
     355  REvaluation bt;
     356  bool gcdfound = false;
     357  Variable x = Variable(1);
     358  modpk bound;
     359
     360  F= FF;
     361  G= GG;
     362
     363  CFMap M,N;
     364  int smallestDegLev;
     365  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
     366
     367  if (best_level == 0)
     368  {
     369    if (!isRat)
     370      Off (SW_RATIONAL);
     371    return G.genOne();
     372  }
     373
     374  F= M (F);
     375  G= M (G);
     376
     377
     378  DEBINCLEVEL( cerr, "ezgcd" );
     379  DEBOUTLN( cerr, "FF = " << FF );
     380  DEBOUTLN( cerr, "GG = " << GG );
     381  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
     382  DEBOUTLN( cerr, "f = " << f );
     383  DEBOUTLN( cerr, "g = " << g );
     384  F /= f; G /= g;
     385  if ( F.isUnivariate() && G.isUnivariate() )
     386  {
     387    DEBDECLEVEL( cerr, "ezgcd" );
     388    if (!isRat)
     389      Off (SW_RATIONAL);
     390    if(F.mvar()==G.mvar())
     391      d*=gcd(F,G);
     392    return N (d);
     393  }
     394  else  if ( gcd_test_one( F, G, false ) )
     395  {
     396    DEBDECLEVEL( cerr, "ezgcd" );
     397    if (!isRat)
     398      Off (SW_RATIONAL);
     399    return N (d);
     400  }
     401  lcF = LC( F, x ); lcG = LC( G, x );
     402  lcD = gcd( lcF, lcG );
     403  delta = 0;
     404  degF = degree( F, x ); degG = degree( G, x );
     405  t = tmax( F.level(), G.level() );
     406  //bound = findBound( F, G, lcF, lcG, degF, degG );
     407  if ( ! internal )
     408    b = REvaluation( 2, t, IntRandom( 25 ) );
     409  while ( ! gcdfound )
     410  {
     411    /// ---> A2
     412    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
     413    DEBOUTLN( cerr, "F = " << F );
     414    DEBOUTLN( cerr, "G = " << G );
     415    findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
     416    DEBOUTLN( cerr, "found evaluation b = " << b );
     417    DEBOUTLN( cerr, "F(b) = " << Fb );
     418    DEBOUTLN( cerr, "G(b) = " << Gb );
     419    DEBOUTLN( cerr, "D(b) = " << Db );
     420    delta = degree( Db );
     421    /// ---> A3
     422    if ( delta == 0 )
     423    {
     424      DEBDECLEVEL( cerr, "ezgcd" );
     425      if (!isRat)
     426        Off (SW_RATIONAL);
     427      return N (d);
     428    }
     429    /// ---> A4
     430    //deltaold = delta;
     431    while ( 1 )
     432    {
     433      bt = b;
     434      findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
     435      int dd=degree( Dbt );
     436      if ( dd /*degree( Dbt )*/ == 0 )
     437      {
     438        DEBDECLEVEL( cerr, "ezgcd" );
     439        if (!isRat)
     440          Off (SW_RATIONAL);
     441        return N (d);
     442      }
     443      if ( dd /*degree( Dbt )*/ == delta )
     444        break;
     445      else  if ( dd /*degree( Dbt )*/ < delta )
     446      {
     447        delta = dd /*degree( Dbt )*/;
     448        b = bt;
     449        Db = Dbt; Fb = Fbt; Gb = Gbt;
     450      }
     451    }
     452    DEBOUTLN( cerr, "now after A4, delta = " << delta );
     453    /// ---> A5
     454    if (delta == degF)
     455    {
     456      if (degF <= degG  && fdivides (F, G))
     457      {
     458        DEBDECLEVEL( cerr, "ezgcd" );
     459        if (!isRat)
     460          Off (SW_RATIONAL);
     461        return N (d*F);
     462      }
     463      else
     464        delta--;
     465    }
     466    else if (delta == degG)
     467    {
     468      if (degG <= degF && fdivides( G, F ))
     469      {
     470        DEBDECLEVEL( cerr, "ezgcd" );
     471        if (!isRat)
     472          Off (SW_RATIONAL);
     473        return N (d*G);
     474      }
     475      else
     476        delta--;
     477    }
     478    if ( delta != degF && delta != degG )
     479    {
     480      /// ---> A6
     481      bool B_is_F;
     482      CanonicalForm xxx1, xxx2;
     483      CanonicalForm buf;
     484      DD[1] = Fb / Db;
     485      buf= Gb/Db;
     486      xxx1 = gcd( DD[1], Db );
     487      xxx2 = gcd( buf, Db );
     488      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     489          (size (F) <= size (G)))
     490            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
     491      {
     492        B = F;
     493        DD[2] = Db;
     494        lcDD[1] = lcF;
     495        lcDD[2] = lcD;
     496        B_is_F = true;
     497      }
     498      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     499                (size (G) < size (F)))
     500                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
     501      {
     502        DD[1] = buf;
     503        B = G;
     504        DD[2] = Db;
     505        lcDD[1] = lcG;
     506        lcDD[2] = lcD;
     507        B_is_F = false;
     508      }
     509      else
     510      {
     511        //special case
     512        Off(SW_USE_EZGCD);
     513        result=gcd( F, G );
     514        On(SW_USE_EZGCD);
     515        if (!isRat)
     516          Off (SW_RATIONAL);
     517        return N (d*result);
     518      }
     519      /// ---> A7
     520      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     521      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     522      DEBOUTLN( cerr, "(hensel) B    = " << B );
     523      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
     524      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
     525      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
     526      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
     527      gcdfound= Hensel (B*lcD, DD, b, lcDD);
     528      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
     529
     530      if (gcdfound == -1)
     531      {
     532        Off (SW_USE_EZGCD);
     533        result= gcd (F,G);
     534        On (SW_USE_EZGCD);
     535        if (!isRat)
     536          Off (SW_RATIONAL);
     537        return N (d*result);
     538      }
     539
     540      if (gcdfound)
     541      {
     542        cand = DD[2] / content( DD[2], Variable(1) );
     543        gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
     544      }
     545      /// ---> A8 (gcdfound)
     546    }
     547    delta--;
     548  }
     549  /// ---> A9
     550  DEBDECLEVEL( cerr, "ezgcd" );
     551  if (!isRat)
     552    Off (SW_RATIONAL);
     553  return N (d*cand);
     554}
    342555
    343556CanonicalForm
    344557ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
    345558{
    346     REvaluation b;
    347     return ezgcd( FF, GG, b, false );
    348 }
    349 
    350 static CanonicalForm
    351 ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal )
    352 {
    353     CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD, cand, result;
    354     CFArray DD( 1, 2 ), lcDD( 1, 2 );
    355     int degF, degG, delta, t;
    356     REvaluation bt;
    357     bool gcdfound = false;
    358     Variable x = Variable(1);
    359     modpk bound;
    360 
    361     F= FF;
    362     G= GG;
    363 
    364     CFMap M,N;
    365     int smallestDegLev;
    366     int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
    367 
    368     if (best_level == 0) return G.genOne();
    369 
    370     F= M (F);
    371     G= M (G);
    372 
    373 
    374     DEBINCLEVEL( cerr, "ezgcd" );
    375     DEBOUTLN( cerr, "FF = " << FF );
    376     DEBOUTLN( cerr, "GG = " << GG );
    377     f = content( F, x ); g = content( G, x ); d = gcd( f, g );
    378     DEBOUTLN( cerr, "f = " << f );
    379     DEBOUTLN( cerr, "g = " << g );
    380     F /= f; G /= g;
    381     if ( F.isUnivariate() && G.isUnivariate() )
    382     {
    383         DEBDECLEVEL( cerr, "ezgcd" );
    384         if(F.mvar()==G.mvar())
    385           d*=gcd(F,G);
    386         return N (d);
    387     }
    388     else  if ( gcd_test_one( F, G, false ) )
    389     {
    390         DEBDECLEVEL( cerr, "ezgcd" );
    391         return N (d);
    392     }
    393     lcF = LC( F, x ); lcG = LC( G, x );
    394     lcD = gcd( lcF, lcG );
    395     delta = 0;
    396     degF = degree( F, x ); degG = degree( G, x );
    397     t = tmax( F.level(), G.level() );
    398     //bound = findBound( F, G, lcF, lcG, degF, degG );
    399     if ( ! internal )
    400         b = REvaluation( 2, t, IntRandom( 25 ) );
    401     while ( ! gcdfound )
    402     {
    403         /// ---> A2
    404         DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
    405         DEBOUTLN( cerr, "F = " << F );
    406         DEBOUTLN( cerr, "G = " << G );
    407         findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
    408         DEBOUTLN( cerr, "found evaluation b = " << b );
    409         DEBOUTLN( cerr, "F(b) = " << Fb );
    410         DEBOUTLN( cerr, "G(b) = " << Gb );
    411         DEBOUTLN( cerr, "D(b) = " << Db );
    412         delta = degree( Db );
    413         /// ---> A3
    414         if ( delta == 0 )
    415         {
    416           DEBDECLEVEL( cerr, "ezgcd" );
    417           return N (d);
    418         }
    419         /// ---> A4
    420         //deltaold = delta;
    421         while ( 1 ) {
    422             bt = b;
    423             findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
    424             int dd=degree( Dbt );
    425             if ( dd /*degree( Dbt )*/ == 0 )
    426             {
    427                 DEBDECLEVEL( cerr, "ezgcd" );
    428                 return N (d);
    429             }
    430             if ( dd /*degree( Dbt )*/ == delta )
    431                 break;
    432             else  if ( dd /*degree( Dbt )*/ < delta ) {
    433                 delta = dd /*degree( Dbt )*/;
    434                 b = bt;
    435                 Db = Dbt; Fb = Fbt; Gb = Gbt;
    436             }
    437         }
    438         DEBOUTLN( cerr, "now after A4, delta = " << delta );
    439         /// ---> A5
    440         if (delta == degF)
    441         {
    442           if (degF <= degG  && fdivides (F, G))
    443           {
    444             DEBDECLEVEL( cerr, "ezgcd" );
    445             return N (d*F);
    446           }
    447           else
    448             delta--;
    449         }
    450         else if (delta == degG)
    451         {
    452           if (degG <= degF && fdivides( G, F ))
    453           {
    454             DEBDECLEVEL( cerr, "ezgcd" );
    455             return N (d*G);
    456           }
    457           else
    458             delta--;
    459         }
    460         if ( delta != degF && delta != degG )
    461         {
    462             /// ---> A6
    463             bool B_is_F;
    464             CanonicalForm xxx1, xxx2;
    465             CanonicalForm buf;
    466             DD[1] = Fb / Db;
    467             buf= Gb/Db;
    468             xxx1 = gcd( DD[1], Db );
    469             xxx2 = gcd( buf, Db );
    470             if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
    471                 (size (F) <= size (G)))
    472                  || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
    473             {
    474               B = F;
    475               DD[2] = Db;
    476               lcDD[1] = lcF;
    477               lcDD[2] = lcD;
    478               B_is_F = true;
    479             }
    480             else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
    481                      (size (G) < size (F)))
    482                      || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
    483             {
    484               DD[1] = buf;
    485               B = G;
    486               DD[2] = Db;
    487               lcDD[1] = lcG;
    488               lcDD[2] = lcD;
    489               B_is_F = false;
    490             }
    491             else
    492             {
    493               Off(SW_USE_EZGCD);
    494               result=gcd( F, G );
    495               On(SW_USE_EZGCD);
    496               return N (d*result);
    497             }
    498             /*CanonicalForm xxx;
    499             //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
    500             DD[1] = Fb / Db;
    501             xxx= gcd( DD[1], Db );
    502             DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
    503             DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
    504             DEBOUTLN( cerr, "Db = " << Db );
    505             if (xxx.inCoeffDomain()) {
    506                 B = F;
    507                 DD[2] = Db;
    508                 lcDD[1] = lcF;
    509                 lcDD[2] = lcF;
    510                 B *= lcF;
    511                 B_is_F=true;
    512             }
    513             //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
    514             else
    515             {
    516               DD[1] = Gb / Db;
    517               xxx=gcd( DD[1], Db );
    518               DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
    519               DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
    520               DEBOUTLN( cerr, "Db = " << Db );
    521               if (xxx.inCoeffDomain())
    522               {
    523                 B = G;
    524                 DD[2] = Db;
    525                 lcDD[1] = lcG;
    526                 lcDD[2] = lcG;
    527                 B *= lcG;
    528                 B_is_F=false;
    529               }
    530               else
    531               {
    532 #ifdef DEBUGOUTPUT
    533                 CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
    534                 DEBDECLEVEL( cerr, "ezgcd" );
    535                 return dummyres;
    536 #else
    537                 return d * ezgcd_specialcase( F, G, b, bound );
    538 #endif
    539               }
    540             }*/
    541             /// ---> A7
    542             DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    543             DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    544             //bound = enlargeBound( B, DD[2], DD[1], bound );
    545             DEBOUTLN( cerr, "(hensel) B    = " << B );
    546             DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
    547             DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
    548             DEBOUTLN( cerr, "(hensel) DD   = " << DD );
    549             DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
    550             gcdfound= Hensel (B*lcD, DD, b, lcDD);
    551             DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
    552 
    553             if (gcdfound == -1)
    554             {
    555               Off (SW_USE_EZGCD);
    556               result= gcd (F,G);
    557               On (SW_USE_EZGCD_P);
    558               return N (d*result);
    559             }
    560 
    561             if (gcdfound)
    562             {
    563               cand = DD[2] / content( DD[2], Variable(1) );
    564               gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
    565             }
    566             /// ---> A8 (gcdfound)
    567         }
    568         delta--;
    569     }
    570     /// ---> A9
    571     DEBDECLEVEL( cerr, "ezgcd" );
    572     return N (d*cand);
    573 }
    574 
    575 static CanonicalForm
    576 ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & /*b*/, const modpk & /*bound*/ )
    577 {
    578     CanonicalForm d;
    579 #if 1
    580     Off(SW_USE_EZGCD);
    581     //bool ntl0=isOn(SW_USE_NTL_GCD_0);
    582     //Off(SW_USE_NTL_GCD_0);
    583     //bool ntlp=isOn(SW_USE_NTL_GCD_P);
    584     //Off(SW_USE_NTL_GCD_P);
    585     d=gcd( F, G );
    586     //if (ntl0) On(SW_USE_NTL_GCD_0);
    587     //if (ntlp) On(SW_USE_NTL_GCD_P);
    588     On(SW_USE_EZGCD);
    589     return d;
    590 #else
    591     /*
    592      * I am not sure, if these lines are much of use.
    593      * Remove?
    594      */
    595     DEBOUTLN( cerr, "ezgcd: special case" );
    596     CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt;
    597     CFArray DD( 1, 2 ), lcDD( 1, 2 );
    598     Variable x = Variable( 1 );
    599     bool gcdfound;
    600 
    601     Dt = 1;
    602     /// ---> S1
    603     DEBOUTLN( cerr, "ezgcdspec: (S1)" );
    604     Ft = ezgcd( F, F.deriv( x ) );
    605     if ( Ft.isOne() )
    606     {
    607         // In this case F is squarefree and we came here by bad chance
    608         // (means: bad evaluation point).  Try again with another
    609         // evaluation point.  Bug fix (?) by JS.  The bad example was
    610         // gcd.debug -ocr /+USE_EZGCD/@12/CB
    611         //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)'
    612         //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
    613         b.nextpoint();
    614         return ezgcd( F, G, b, true );
    615     }
    616 
    617     DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
    618     L = F / Ft;
    619     /// ---> S2
    620     DEBOUTLN( cerr, "ezgcdspec: (S2)" );
    621 
    622     L = ezgcd( L, G, b, true );
    623     DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds );
    624     Gt = G / L;
    625     Ds = L; Dt = L;
    626     Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
    627     d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
    628     while ( true ) {
    629         /// ---> S3
    630         DEBOUTLN( cerr, "ezgcdspec: (S3)" );
    631         DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
    632         DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
    633         DD[1] = gcd( Lb, gcd( Fb, Gb ) );
    634         /// ---> S4
    635         DEBOUTLN( cerr, "ezgcdspec: (S4)" );
    636         if ( degree( DD[1] ) == 0 )
    637             return Ds;
    638         /// ---> S5
    639         DEBOUTLN( cerr, "ezgcdspec: (S5)" );
    640         DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
    641         DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
    642         Db = DD[1];
    643         if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
    644             DEBOUTLN( cerr, "ezgcdspec: (S55)" );
    645             lcDD[2] = LC( L, x );
    646             lcDD[1] = d;
    647             DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
    648             DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
    649             LL = L * d;
    650             modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
    651             DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
    652             DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
    653             DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
    654             DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
    655             DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
    656             DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
    657             DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
    658             DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
    659             gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
    660             ASSERT( gcdfound, "fatal error in algorithm" );
    661             Dt = pp( DD[1] );
    662         }
    663         DEBOUTLN( cerr, "ezgcdspec: (S7)" );
    664         Ds = Ds * Dt;
    665         Fb = Fb / Db;
    666         Gb = Gb / Db;
    667     }
    668     // this point is never reached, but to avoid compiler warnings let's return a value
    669     return 0;
    670 #endif
    671 }
    672 
    673 static void
    674 findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG )
    675 {
    676     bool ok;
    677     if ( delta != 0 )
    678         b.nextpoint();
    679     DEBOUTLN( cerr, "ezgcd: (findeval) F = " << F  <<", G="<< G);
    680     DEBOUTLN( cerr, "ezgcd: (findeval) degF = " << degF << ", degG="<<degG );
    681     do {
    682         DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
    683         Fb = b( F );
    684         ok = degree( Fb ) == degF;
    685         if ( ok ) {
    686             Gb = b( G );
    687             ok = degree( Gb ) == degG;
    688         }
    689 
    690         if ( ok )
    691         {
    692             Db = gcd( Fb, Gb );
    693             if ( delta > 0 )
    694               ok = degree( Db ) < delta;
    695         }
    696         if ( ! ok )
    697         {
    698             b.nextpoint();
    699         }
    700     } while ( ! ok );
    701 }
    702 
    703 static modpk
    704 enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk )
    705 {
    706     DEBOUTLN( cerr, "ezgcd: (enlarge bound) F      = " << F );
    707     DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb     = " << Lb );
    708     DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db     = " << Db );
    709     DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) );
    710     DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) );
    711 
    712     CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
    713         tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
    714     int p = pk.getp();
    715     int k = pk.getk();
    716     CanonicalForm bound = pk.getpk();
    717     while ( bound < limit ) {
    718         k++;
    719         bound *= p;
    720     }
    721     k *= 2;
    722     DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) );
    723     return modpk( p, k );
    724 }
    725 
    726 static modpk
    727 findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG )
    728 {
    729     CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
    730         gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
    731     int p, i = 0;
    732     do {
    733         p = cf_getBigPrime( i );
    734         i++;
    735     } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
    736     CanonicalForm bound = p;
    737     i = 1;
    738     while ( bound < limit ) {
    739         i++;
    740         bound *= p;
    741     }
    742     return modpk( p, i );
    743 }
     559  REvaluation b;
     560  return ezgcd( FF, GG, b, false );
     561}
     562
Note: See TracChangeset for help on using the changeset viewer.