Changeset 3fd1ff in git


Ignore:
Timestamp:
Jul 1, 2008, 3:45:44 PM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
be03c4a08b21f0e77bc9586601f803a77c5bf3e6
Parents:
7ff04ad0f18519b82845de025fb9ed6890b0b366
Message:
*hannes: ezgcd_p


git-svn-id: file:///usr/local/Singular/svn/trunk@10821 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_factor.cc

    r7ff04a r3fd1ff  
    11/* emacs edit mode for this file is -*- C++ -*- */
    2 /* $Id: cf_factor.cc,v 1.40 2008-03-17 17:44:04 Singular Exp $ */
     2/* $Id: cf_factor.cc,v 1.41 2008-07-01 13:45:44 Singular Exp $ */
    33
    44//{{{ docu
    … …  
    120120      printf("+%d",f.intval());
    121121    }
    122     else
     122    else
     123    {
    123124    #ifdef NOSTREAMIO
     125      #ifdef SINGULAR
     126      if (f.inZ())
     127      {
     128        MP_INT m=gmp_numerator(f);
     129        char * str = new char[mpz_sizeinbase( &m, 10 ) + 2];
     130        str = mpz_get_str( str, 10, &m );
     131        printf("%s",str);
     132        delete[] str;
     133      }
     134      else if (f.inQ())
     135      {
     136        MP_INT m=gmp_numerator(f);
     137        char * str = new char[mpz_sizeinbase( &m, 10 ) + 2];
     138        str = mpz_get_str( str, 10, &m );
     139        printf("%s/",str);
     140        delete[] str;
     141        m=gmp_denominator(f);
     142        str = new char[mpz_sizeinbase( &m, 10 ) + 2];
     143        str = mpz_get_str( str, 10, &m );
     144        printf("%s",str);
     145        delete[] str;
     146      }
     147      #else
    124148      printf("+...");
     149      #endif
    125150    #else
    126151       std::cout << f;
    127152    #endif
     153    }
    128154    //if (f.inZ()) printf("(Z)");
    129155    //else if (f.inQ()) printf("(Q)");
    … …  
    254280  CFIterator i;
    255281  CanonicalForm elem, result(0);
    256  
     282
    257283  for (i=f; i.hasTerms(); i++)
    258284  {
    … …  
    261287    if ( deg < maxdeg )
    262288      result += elem * power(x,maxdeg-deg);
    263     else 
     289    else
    264290      result+=elem;
    265291  }
    … …  
    294320  CFIterator i;
    295321  CanonicalForm elem, result(0);
    296  
     322
    297323  for (i=f; i.hasTerms(); i++)
    298324  {
    … …  
    301327    if ( deg < maxdeg )
    302328      result += elem * power(x,maxdeg-deg);
    303     else 
     329    else
    304330      result+=elem;
    305331  }
    … …  
    373399      }
    374400      if ( d_xn != 0 ) // have to append xn^(d_xn)
    375         Unhomoglist.append(CFFactor(CanonicalForm(xn),d_xn)); 
    376       if(isOn(SW_USE_NTL_SORT)) Unhomoglist.sort(cmpCF); 
     401        Unhomoglist.append(CFFactor(CanonicalForm(xn),d_xn));
     402      if(isOn(SW_USE_NTL_SORT)) Unhomoglist.sort(cmpCF);
    377403      return Unhomoglist;
    378404    }
    … …  
    404430      if (getCharacteristic()!=2)
    405431      {
    406         // set remainder
    407432        if (fac_NTL_char!=getCharacteristic())
    408433        {
    409434          fac_NTL_char=getCharacteristic();
     435          #ifndef NTL_ZZ
     436          if (fac_NTL_char >NTL_SP_BOUND)
     437          {
     438            ZZ r;
     439            r=getCharacteristic();
     440            ZZ_pContext ccc(r);
     441            ccc.restore();
     442            ZZ_p::init(r);
     443          }
     444          else
     445          #endif
     446          {
     447            #ifdef NTL_ZZ
     448            ZZ r;
     449            r=getCharacteristic();
     450            ZZ_pContext ccc(r);
     451            #else
     452            zz_pContext ccc(getCharacteristic());
     453            #endif
     454            ccc.restore();
     455            #ifdef NTL_ZZ
     456            ZZ_p::init(r);
     457            #else
     458            zz_p::init(getCharacteristic());
     459            #endif
     460          }
     461        }
     462       #ifndef NTL_ZZ
     463       if (fac_NTL_char >NTL_SP_BOUND)
     464       {
     465          // convert to NTL
     466          ZZ_pX f1=convertFacCF2NTLZZpX(f);
     467          ZZ_p leadcoeff = LeadCoeff(f1);
     468          //make monic
     469          f1=f1 / LeadCoeff(f1);
     470          // factorize
     471          vec_pair_ZZ_pX_long factors;
     472          CanZass(factors,f1);
     473          // convert back to factory
     474          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
     475        }
     476        else
     477        #endif
     478        {
     479          // convert to NTL
    410480          #ifdef NTL_ZZ
    411           ZZ r;
    412           r=getCharacteristic();
    413           ZZ_pContext ccc(r);
     481          ZZ_pX f1=convertFacCF2NTLZZpX(f);
     482          ZZ_p leadcoeff = LeadCoeff(f1);
    414483          #else
    415           zz_pContext ccc(getCharacteristic());
     484          zz_pX f1=convertFacCF2NTLzzpX(f);
     485          zz_p leadcoeff = LeadCoeff(f1);
    416486          #endif
    417           ccc.restore();
     487          //make monic
     488          f1=f1 / LeadCoeff(f1);
     489          // factorize
    418490          #ifdef NTL_ZZ
    419           ZZ_p::init(r);
     491          vec_pair_ZZ_pX_long factors;
    420492          #else
    421           zz_p::init(getCharacteristic());
     493          vec_pair_zz_pX_long factors;
     494          #endif
     495          CanZass(factors,f1);
     496          // convert back to factory
     497          #ifdef NTL_ZZ
     498          F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
     499          #else
     500          F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
    422501          #endif
    423502        }
    424         // convert to NTL
    425         #ifdef NTL_ZZ
    426         ZZ_pX f1=convertFacCF2NTLZZpX(f);
    427         ZZ_p leadcoeff = LeadCoeff(f1);
    428         #else
    429         zz_pX f1=convertFacCF2NTLzzpX(f);
    430         zz_p leadcoeff = LeadCoeff(f1);
    431         #endif
    432         //make monic
    433         f1=f1 / LeadCoeff(f1);
    434 
    435         // factorize
    436         #ifdef NTL_ZZ
    437         vec_pair_ZZ_pX_long factors;
    438         #else
    439         vec_pair_zz_pX_long factors;
    440         #endif
    441         CanZass(factors,f1);
    442 
    443         // convert back to factory
    444         #ifdef NTL_ZZ
    445         F=convertNTLvec_pair_ZZpX_long2FacCFFList(factors,leadcoeff,f.mvar());
    446         #else
    447         F=convertNTLvec_pair_zzpX_long2FacCFFList(factors,leadcoeff,f.mvar());
    448         #endif
    449503        //test_cff(F,f);
    450504      }
    … …  
    565619  }
    566620  //out_cff(F);
    567   if(isOn(SW_USE_NTL_SORT)) F.sort(cmpCF); 
     621  if(isOn(SW_USE_NTL_SORT)) F.sort(cmpCF);
    568622  return F;
    569623}

  • factory/ffreval.cc

    r7ff04a r3fd1ff  
    1414#include "ftmpl_functions.h"
    1515#include "ffreval.h"
    16 
    17 void FFREvaluation::init()
    18 {
    19   int n = values.max();
    20 
    21   for( int i=values.min(); i<=n; i++ )
    22   {
    23     CanonicalForm tmp=gen->generate();
    24     values[i] = tmp; // generate random point
    25     offset[i] = tmp; // generate random point
    26   }
    27 }
    28 
    29 bool FFREvaluation::step()
     16#include "cf_binom.h"
     17#include "fac_iterfor.h"
     18#include "cf_iter.h"
     19
     20void FFREvaluation::nextpoint()
    3021{
    3122  // enumerates the points stored in values
    3223  int n = values.max();
    3324  int p = getCharacteristic();
    34 
    35   int i;
    36   for(i=values.min(); i<=n; i++ )
     25  for( int i=values.min(); i<=n; i++ )
    3726  {
    3827    if( values[i] != p-1 )
    … …  
    4332    values[i] += 1; // becomes 0
    4433  }
    45   for(i=values.min();i<=n;i++)
    46   {
    47     if(values[i]!=offset[i]) return true;
     34}
     35
     36bool FFREvaluation::hasNext()
     37{
     38  int n = values.max();
     39
     40  for( int i=values.min(); i<=n; i++ )
     41  {
     42    if( values[i]!=start[i] )
     43      return true;
    4844  }
    4945  return false;
    … …  
    5753      delete gen;
    5854    values = e.values;
    59     offset = e.offset;
     55    start = e.start;
    6056    if( e.gen == NULL )
    6157      gen = NULL;
    … …  
    6864/* ------------------------------------------------------------------------ */
    6965/* forward declarations: fin_ezgcd stuff*/
    70 static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal );
    71 static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG );
    72 static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound );
    73 
    74 CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
    75 {
    76     FFREvaluation b;
    77     return fin_ezgcd( FF, GG, b, false );
    78 }
    79 
    80 static CanonicalForm fin_ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, FFREvaluation & b, bool internal )
    81 {
    82     CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
    83     CFArray DD( 1, 2 ), lcDD( 1, 2 );
    84     int degF, degG, delta, t;
    85     FFREvaluation bt;
    86     bool gcdfound = false;
    87     Variable x = Variable(1);
    88     modpk p = modpk(getCharacteristic(), 1); // k = 1
    89     DEBINCLEVEL( cerr, "fin_ezgcd" );
    90     DEBOUTLN( cerr, "FF = " << FF );
    91     DEBOUTLN( cerr, "GG = " << GG );
    92     f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
    93     DEBOUTLN( cerr, "f = " << f );
    94     DEBOUTLN( cerr, "g = " << g );
    95     F = FF / f; G = GG / g;
    96     if ( F.isUnivariate() && G.isUnivariate() )
    97     {
    98         DEBDECLEVEL( cerr, "fin_ezgcd" );
    99         if(F.mvar()==G.mvar())
    100           d*=gcd_poly(F,G);
     66static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int bound, int & counter );
     67static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a );
     68static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a );
     69static CanonicalForm checkCombination_P ( const CanonicalForm & W, const Evaluation & a, const IteratedFor & e, int k );
     70static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n );
     71static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h );
     72static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x );
     73
     74CanonicalForm fin_ezgcd( const CanonicalForm & FF, const CanonicalForm & GG )
     75{
     76  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
     77  CFArray DD( 1, 2 ), lcDD( 1, 2 );
     78  int degF, degG, delta, count, maxeval;
     79  count = 0; // number of eval. used
     80  maxeval = getCharacteristic(); // bound on the number of eval. to use
     81  REvaluation b, bt;
     82  bool gcdfound = false;
     83  Variable x = Variable(1);
     84  f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
     85  F = FF / f; G = GG / g;
     86  if( F.isUnivariate() && G.isUnivariate() )
     87  {
     88    if( F.mvar() == G.mvar() )
     89      d *= gcd( F, G );
     90    return d;
     91  }
     92  if( gcd_test_one( F, G, false ) )
     93    return d;
     94
     95  lcF = LC( F, x ); lcG = LC( G, x );
     96  lcD = gcd( lcF, lcG );
     97  delta = 0;
     98  degF = degree( F, x ); degG = degree( G, x );
     99  b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
     100  while( !gcdfound )
     101  {
     102    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
     103    { // too many eval. used --> try another method
     104      Off( SW_USE_EZGCD_P );
     105      d *= gcd( F, G );
     106      On( SW_USE_EZGCD_P );
     107      return d;
     108    }
     109    delta = degree( Db );
     110    if( delta == 0 )
     111      return d;
     112    while( true )
     113    {
     114      bt = b;
     115      if( !findeval_P( F, G, Fbt, Gbt, Dbt, bt, delta+1, degF, degG, maxeval, count ))
     116      { // too many eval. used --> try another method
     117        Off( SW_USE_EZGCD_P );
     118        d *= gcd( F, G );
     119        On( SW_USE_EZGCD_P );
    101120        return d;
    102     }
    103     else  if ( gcd_test_one( F, G, false ) )
    104     {
    105         DEBDECLEVEL( cerr, "fin_ezgcd" );
     121      }
     122      int dd = degree( Dbt );
     123      if( dd == 0 )
    106124        return d;
    107     }
    108     lcF = LC( F, x ); lcG = LC( G, x );
    109     lcD = gcd( lcF, lcG );
    110     delta = 0;
    111     degF = degree( F, x ); degG = degree( G, x );
    112     t = tmax( F.level(), G.level() );
    113     if ( ! internal )
    114     {
    115         b = FFREvaluation( 2, t, FFRandom() );
    116         b.init(); // choose random point
    117     }
    118     while ( ! gcdfound )
    119     {
    120         /// ---> A2
    121         DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
    122         DEBOUTLN( cerr, "F = " << F );
    123         DEBOUTLN( cerr, "G = " << G );
    124         if( ! fin_findeval( F, G, Fb, Gb, Db, b, delta, degF, degG ) )
     125      if( dd == delta )
     126        break;
     127      if( dd < delta )
     128      {
     129        delta = dd;
     130        b = bt;
     131        Db = Dbt; Fb = Fbt; Gb = Gbt;
     132      }
     133    }
     134    if( degF <= degG && delta == degF && fdivides( F, G ) )
     135      return d*F;
     136    if( degG < degF && delta == degG && fdivides( G, F ) )
     137      return d*G;
     138    if( delta != degF && delta != degG )
     139    {
     140      bool B_is_F;
     141      CanonicalForm xxx;
     142      DD[1] = Fb / Db;
     143      xxx = gcd( DD[1], Db );
     144      if( xxx.inCoeffDomain() )
     145      {
     146        B = F;
     147        DD[2] = Db;
     148        lcDD[1] = lcF;
     149        lcDD[2] = lcF;
     150        B *= lcF;
     151        B_is_F = true;
     152      }
     153      else
     154      {
     155        DD[1] = Gb / Db;
     156        xxx = gcd( DD[1], Db );
     157        if( xxx.inCoeffDomain() )
    125158        {
    126            Off(SW_USE_EZGCD_P);
    127            d *= gcd_poly(F,G);
    128            On(SW_USE_EZGCD_P);
    129            return d;
     159          B = G;
     160          DD[2] = Db;
     161          lcDD[1] = lcG;
     162          lcDD[2] = lcG;
     163          B *= lcG;
     164          B_is_F = false;
    130165        }
    131 
    132         DEBOUTLN( cerr, "found evaluation b = " << b );
    133         DEBOUTLN( cerr, "F(b) = " << Fb );
    134         DEBOUTLN( cerr, "G(b) = " << Gb );
    135         DEBOUTLN( cerr, "D(b) = " << Db );
    136         delta = degree( Db );
    137         /// ---> A3
    138         if ( delta == 0 )
     166        else // special case handling
    139167        {
    140           DEBDECLEVEL( cerr, "fin_ezgcd" );
     168          Off( SW_USE_EZGCD_P );
     169          d *= gcd( F, G ); // try another method
     170          On( SW_USE_EZGCD_P );
    141171          return d;
    142172        }
    143         /// ---> A4
    144         //deltaold = delta;
    145         while ( 1 )
     173      }
     174      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     175      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     176
     177      gcdfound = Hensel_P( B, DD, lcDD, b, x );
     178
     179      if( gcdfound )
     180      {
     181        CanonicalForm cand = DD[2] / content( DD[2], Variable(1) );
     182        if( B_is_F )
     183          gcdfound = fdivides( cand, G );
     184        else
     185          gcdfound = fdivides( cand, F );
     186      }
     187    }
     188    delta++;
     189  }
     190  return d * DD[2] / content( DD[2],Variable(1) );
     191}
     192
     193
     194static bool findeval_P( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG, int maxeval, int & count )
     195{
     196  if( delta != 0 )
     197  {
     198    b.nextpoint();
     199    if( count++ > maxeval )
     200      return false; // too many eval. used
     201  }
     202  while( true )
     203  {
     204    Fb = b( F );
     205    if( degree( Fb ) == degF )
     206    {
     207      Gb = b( G );
     208      if( degree( Gb ) == degG )
     209      {
     210        Db = gcd( Fb, Gb );
     211        if( delta > 0 )
    146212        {
    147             bt = b;
    148             if( ! fin_findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG ) )
     213          if( degree( Db ) < delta )
     214            return true;
     215        }
     216        else
     217          return true;
     218      }
     219    }
     220    b.nextpoint();
     221    if( count++ > maxeval ) // too many eval. used
     222      return false;
     223  }
     224}
     225
     226
     227static void solveF_P ( const CFArray & P, const CFArray & Q, const CFArray & S, const CFArray & T, const CanonicalForm & C, int r, CFArray & a )
     228{
     229  CanonicalForm bb, b = C;
     230  for ( int j = 1; j < r; j++ )
     231  {
     232    divrem( S[j] * b, Q[j], a[j], bb );
     233    a[j] = T[j] * b + a[j] * P[j];
     234    b = bb;
     235  }
     236  a[r] = b;
     237}
     238
     239
     240static CanonicalForm derivAndEval_P ( const CanonicalForm & f, int n, const Variable & x, const CanonicalForm & a )
     241{
     242  if ( n == 0 )
     243    return f( a, x );
     244  if ( f.degree( x ) < n )
     245    return 0;
     246  CFIterator i;
     247  CanonicalForm sum = 0, fact;
     248  int min, j;
     249  Variable v = Variable( f.level() + 1 );
     250  for ( i = swapvar( f, x, v); i.hasTerms() && i.exp() >= n; i++ )
     251  {
     252    fact = 1;
     253    min = i.exp() - n;
     254    for ( j = i.exp(); j > min; j-- )
     255      fact *= j;
     256    sum += fact * i.coeff() * power( v, min );
     257  }
     258  return sum( a, v );
     259}
     260
     261
     262static CFArray findCorrCoeffs_P ( const CFArray & P, const CFArray & Q, const CFArray & P0, const CFArray & Q0, const CFArray & S, const CFArray & T, const CanonicalForm & C, const Evaluation & I, int r, int k, int h, int * n )
     263{
     264  bool what;
     265  int i, j, m;
     266  CFArray A(1,r), a(1,r);
     267  CanonicalForm C0, Dm, g, prodcomb;
     268  C0 = I( C, 2, k-1 );
     269  solveF_P( P0, Q0, S, T, 1, r, a );
     270  for ( i = 1; i <= r; i++ )
     271    A[i] = (a[i] * C0) % P0[i];
     272  for ( m = 0; m <= h && ( m == 0 || Dm != 0 ); m++ )
     273  {
     274    Dm = -C;
     275    prodcomb = 1;
     276    for ( i = 1; i <= r; i++ )
     277    {
     278      Dm += prodcomb * A[i] * Q[i];
     279      prodcomb *= P[i];
     280    }
     281    if ( Dm != 0 )
     282    {
     283      if ( k == 2 )
     284      {
     285        solveF_P( P0, Q0, S, T, Dm, r, a );
     286        for ( i = 1; i <= r; i++ )
     287          A[i] -= a[i];
     288      }
     289      else
     290      {
     291        IteratedFor e( 2, k-1, m+1 );
     292        while ( e.iterations_left() )
     293        {
     294          j = 0;
     295          what = true;
     296          for ( i = 2; i < k; i++ )
     297          {
     298            j += e[i];
     299            if ( e[i] > n[i] )
    149300            {
    150               Off(SW_USE_EZGCD_P);
    151               d *= gcd_poly(F,G);
    152               On(SW_USE_EZGCD_P);
    153               return d;
     301              what = false;
     302              break;
    154303            }
    155 
    156             int dd=degree( Dbt );
    157             if ( dd /*degree( Dbt )*/ == 0 )
     304          }
     305          if ( what && j == m+1 )
     306          {
     307            g = Dm;
     308            for ( i = k-1; i >= 2 && ! g.isZero(); i-- )
     309            g = derivAndEval_P( g, e[i], Variable( i ), I[i] );
     310            if ( ! g.isZero() )
    158311            {
    159                 DEBDECLEVEL( cerr, "fin_ezgcd" );
    160                 return d;
    161             }
    162             if ( dd /*degree( Dbt )*/ == delta )
    163                 break;
    164             else  if ( dd /*degree( Dbt )*/ < delta )
    165             {
    166                 delta = dd /*degree( Dbt )*/;
    167                 b = bt;
    168                 Db = Dbt; Fb = Fbt; Gb = Gbt;
    169             }
    170         }
    171         DEBOUTLN( cerr, "now after A4, delta = " << delta );
    172         /// ---> A5
    173         if ( degF <= degG && delta == degF && fdivides( F, G ) )
    174         {
    175             DEBDECLEVEL( cerr, "fin_ezgcd" );
    176             return d*F;
    177         }
    178         if ( degG < degF && delta == degG && fdivides( G, F ) )
    179         {
    180             DEBDECLEVEL( cerr, "fin_ezgcd" );
    181             return d*G;
    182         }
    183         if ( delta != degF && delta != degG )
    184         {
    185             bool B_is_F;
    186             /// ---> A6
    187             CanonicalForm xxx;
    188             //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
    189             DD[1] = Fb / Db;
    190             xxx= gcd( DD[1], Db );
    191             DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
    192             DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
    193             DEBOUTLN( cerr, "Db = " << Db );
    194             if (xxx.inCoeffDomain())
    195             {
    196                 B = F;
    197                 DD[2] = Db;
    198                 lcDD[1] = lcF;
    199                 lcDD[2] = lcF;
    200                 B *= lcF;
    201                 B_is_F=true;
    202             }
    203             //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
    204             else
    205             {
    206               DD[1] = Gb / Db;
    207               xxx=gcd( DD[1], Db );
    208               DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
    209               DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
    210               DEBOUTLN( cerr, "Db = " << Db );
    211               if (xxx.inCoeffDomain())
     312              prodcomb = 1;
     313              for ( i = 2; i < k; i++ )
     314                for ( j = 2; j <= e[i]; j++ )
     315                  prodcomb *= j;
     316              g /= prodcomb;
     317              if( ! (g.mvar() > Variable(1)) )
    212318              {
    213                 B = G;
    214                 DD[2] = Db;
    215                 lcDD[1] = lcG;
    216                 lcDD[2] = lcG;
    217                 B *= lcG;
    218                 B_is_F=false;
    219               }
    220               else
    221               {
    222 #ifdef DEBUGOUTPUT
    223                 CanonicalForm dummyres = d * fin_ezgcd_specialcase( F, G, b, p );
    224                 DEBDECLEVEL( cerr, "fin_ezgcd" );
    225                 return dummyres;
    226 #else
    227                 return d * fin_ezgcd_specialcase( F, G, b, p );
    228 #endif
     319                prodcomb = 1;
     320                for ( i = k-1; i > 1; i-- )
     321                  prodcomb *= binomialpower( Variable(i), -I[i], e[i] );
     322                solveF_P( P0, Q0, S, T, g, r, a );
     323                for ( i = 1; i <= r; i++ )
     324                  A[i] -= a[i] * prodcomb;
    229325              }
    230326            }
    231             /// ---> A7
    232             DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    233             DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    234             DEBOUTLN( cerr, "(hensel) B    = " << B );
    235             DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
    236             DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
    237             DEBOUTLN( cerr, "(hensel) DD   = " << DD );
    238             DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
    239             gcdfound = Hensel( B, DD, lcDD, b, p, x );
    240             DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
    241             if (gcdfound)
    242             {
    243               CanonicalForm xxx=content(DD[2],Variable(1));
    244               CanonicalForm cand=DD[2] / xxx; //content(DD[2],Variable(1));
    245 #if 0
    246               gcdfound= fdivides(cand,G) &&  fdivides(cand,F);
    247 #else
    248               if (B_is_F /*B==F*lcF*/)
    249               {
    250                 DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
    251                 gcdfound= fdivides(cand,G);
    252               }
    253               else
    254               {
    255                 DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
    256                 gcdfound= fdivides(cand,F);
    257               }
    258 #endif
    259             }
    260             /// ---> A8 (gcdfound)
     327          }
     328          e++;
    261329        }
    262         delta++;
    263     }
    264     /// ---> A9
    265     DEBDECLEVEL( cerr, "fin_ezgcd" );
    266     return d * DD[2] / content(DD[2],Variable(1));
    267 }
    268 
    269 static CanonicalForm fin_ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, FFREvaluation & b, const modpk & bound )
    270 {
    271     CanonicalForm d;
    272     Off(SW_USE_EZGCD_P);
    273     d=gcd_poly( F, G );
    274     On(SW_USE_EZGCD_P);
    275     return d;
    276 }
    277 
    278 static bool fin_findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, FFREvaluation & b, int delta, int degF, int degG )
    279 {
    280     int i;
    281     bool ok;
    282     DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) F = " << F  <<", G="<< G);
    283     DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) degF = " << degF << ", degG="<<degG );
    284     while(1)
    285     {
    286         DEBOUTLN( cerr, "fin_ezgcd: (fin_findeval) b = " << b );
    287         Fb = b( F );
    288         ok = degree( Fb ) == degF;
    289         if ( ok )
    290         {
    291             Gb = b( G );
    292             ok = degree( Gb ) == degG;
    293         }
    294         if ( ok )
    295         {
    296             Db = gcd( Fb, Gb );
    297             if ( delta > 0 )
    298               ok = degree( Db ) < delta;
    299         }
    300         if ( ok )
    301             break;
    302 
    303         if( b.step() ) // can choose valid point
    304            continue;
    305 
    306         return false;
    307     }
    308     return true;
    309 }
     330      }
     331    }
     332  }
     333  return A;
     334}
     335
     336
     337static bool liftStep_P ( CFArray & P, int k, int r, int t, const Evaluation & A, const CFArray & lcG, const CanonicalForm & Uk, int * n, int h )
     338{
     339  CFArray K( 1, r ), Q( 1, r ), Q0( 1, r ), P0( 1, r ), S( 1, r ), T( 1, r ), alpha( 1, r );
     340  CanonicalForm Rm, C, D, xa = Variable(k) - A[k];
     341  CanonicalForm xa1 = xa, xa2 = xa*xa;
     342  int i, m;
     343  for ( i = 1; i <= r; i++ )
     344  {
     345    Variable vm = Variable( t + 1 );
     346    Variable v1 = Variable(1);
     347    K[i] = swapvar( replaceLc( swapvar( P[i], v1, vm ), swapvar( A( lcG[i], k+1, t ), v1, vm ) ), v1, vm );
     348    P[i] = A( K[i], k, t );
     349  }
     350  Q[r] = 1;
     351  for ( i = r; i > 1; i-- )
     352  {
     353    Q[i-1] = Q[i] * P[i];
     354    P0[i] = A( P[i], 2, k-1 );
     355    Q0[i] = A( Q[i], 2, k-1 );
     356    (void) extgcd( P0[i], Q0[i], S[i], T[i] );
     357  }
     358  P0[1] = A( P[1], 2, k-1 );
     359  Q0[1] = A( Q[1], 2, k-1 );
     360  (void) extgcd( P0[1], Q0[1], S[1], T[1] );
     361  for ( m = 1; m <= n[k]+1; m++ )
     362  {
     363    Rm = prod( K ) - Uk;
     364    for ( i = 2; i < k; i++ )
     365      Rm.mod( binomialpower( Variable(i), -A[i], n[i]+1 ) );
     366    if ( mod( Rm, xa2 ) != 0 )
     367    {
     368      C = derivAndEval_P( Rm, m, Variable( k ), A[k] );
     369      D = 1;
     370      for ( i = 2; i <= m; i++ ) D *= i;
     371        C /= D;
     372      alpha = findCorrCoeffs_P( P, Q, P0, Q0, S, T, C, A, r, k, h, n );
     373      for ( i = 1; i <= r; i++ )
     374        K[i] -= alpha[i] * xa1;
     375    }
     376    xa1 = xa2;
     377    xa2 *= xa;
     378  }
     379  for ( i = 1; i <= r; i++ )
     380    P[i] = K[i];
     381  return prod( K ) - Uk == 0;
     382}
     383
     384
     385static bool Hensel_P ( const CanonicalForm & U, CFArray & G, const CFArray & lcG, const Evaluation & A, const Variable & x )
     386{
     387  int k, i, h, t = A.max();
     388  bool goodeval = true;
     389  CFArray Uk( A.min(), A.max() );
     390  int * n = new int[t+1];
     391  Uk[t] = U;
     392  for ( k = t-1; k > 1; k-- )
     393  {
     394    Uk[k] = Uk[k+1]( A[k+1], Variable( k+1 ) );
     395    n[k] = degree( Uk[k], Variable( k ) );
     396  }
     397  for ( k = A.min(); goodeval && (k <= t); k++ )
     398  {
     399    h = totaldegree( Uk[k], Variable(A.min()), Variable(k-1) );
     400    for ( i = A.min(); i <= k; i++ )
     401      n[i] = degree( Uk[k], Variable(i) );
     402    goodeval = liftStep_P( G, k, G.max(), t, A, lcG, Uk[k], n, h );
     403  }
     404  delete[] n;
     405  return goodeval;
     406}
     407

  • factory/ffreval.h

    r7ff04a r3fd1ff  
    1111{
    1212  private:
    13     CFArray offset; // random point
     13    CFArray start; // random point
    1414  public:
    15     FFREvaluation( ) : REvaluation() {}
    16     FFREvaluation( int min, int max ) : offset( min, max ) {}
    17     FFREvaluation( int min, int max, const FFRandom & sample ) : REvaluation( min, max, sample ) {}
    18     void init();
    19     bool step();
     15    FFREvaluation( ) : REvaluation(), start() {}
     16    FFREvaluation( int min, int max, const FFRandom & sample ) : REvaluation( min, max, sample ), start( min, max )
     17    {
     18      for( int i=min; i<=max; i++ )
     19        values[i] = start[i] = gen->generate();  //generate random point
     20
     21      nextpoint();
     22    }
    2023    FFREvaluation& operator= ( const FFREvaluation & e );
     24    void nextpoint();
     25    bool hasNext();
    2126};
    2227
Note: See TracChangeset for help on using the changeset viewer.