Changeset 845303c in git


Ignore:
Timestamp:
Jul 8, 2020, 5:47:17 PM (4 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
afa77551770f2019ad99556a95d1327e8a25c616
Parents:
9bd7cb708540aed8dc6ea40c071e5b0ff24576e9
Message:
fix factory: make check w/o NTL, w.FINT2.5.2
File:
1 edited

Legend:

Unmodified
Added
Removed
  • factory/cf_gcd.cc

    r9bd7cb r845303c  
    2828#include "FLINTconvert.h"
    2929#include "facAlgFuncUtil.h"
     30#include "templates/ftmpl_functions.h"
    3031
    3132#ifdef HAVE_NTL
     
    7475{
    7576    return icontent( f, 0 );
     77}
     78
     79#ifdef HAVE_FLINT
     80static CanonicalForm
     81gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
     82{
     83  nmod_poly_t F1, G1;
     84  convertFacCF2nmod_poly_t (F1, F);
     85  convertFacCF2nmod_poly_t (G1, G);
     86  nmod_poly_gcd (F1, F1, G1);
     87  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
     88  nmod_poly_clear (F1);
     89  nmod_poly_clear (G1);
     90  return result;
     91}
     92
     93static CanonicalForm
     94gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
     95{
     96  fmpz_poly_t F1, G1;
     97  convertFacCF2Fmpz_poly_t(F1, F);
     98  convertFacCF2Fmpz_poly_t(G1, G);
     99  fmpz_poly_gcd (F1, F1, G1);
     100  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
     101  fmpz_poly_clear (F1);
     102  fmpz_poly_clear (G1);
     103  return result;
     104}
     105#endif
     106
     107#ifdef HAVE_NTL
     108static CanonicalForm
     109gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
     110{
     111    ZZX F1=convertFacCF2NTLZZX(F);
     112    ZZX G1=convertFacCF2NTLZZX(G);
     113    ZZX R=GCD(F1,G1);
     114    return convertNTLZZX2CF(R,F.mvar());
     115}
     116
     117static CanonicalForm
     118gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
     119{
     120  if (fac_NTL_char!=getCharacteristic())
     121  {
     122    fac_NTL_char=getCharacteristic();
     123    zz_p::init(getCharacteristic());
     124  }
     125  zz_pX F1=convertFacCF2NTLzzpX(F);
     126  zz_pX G1=convertFacCF2NTLzzpX(G);
     127  zz_pX R=GCD(F1,G1);
     128  return  convertNTLzzpX2CF(R,F.mvar());
     129}
     130#endif
     131
     132//{{{ static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q )
     133//{{{ docu
     134//
     135// balance_p() - map f from positive to symmetric representation
     136//   mod q.
     137//
     138// This makes sense for univariate polynomials over Z only.
     139// q should be an integer.
     140//
     141// Used by gcd_poly_univar0().
     142//
     143//}}}
     144
     145static CanonicalForm
     146balance_p ( const CanonicalForm & f, const CanonicalForm & q, const CanonicalForm & qh )
     147{
     148    Variable x = f.mvar();
     149    CanonicalForm result = 0;
     150    CanonicalForm c;
     151    CFIterator i;
     152    for ( i = f; i.hasTerms(); i++ )
     153    {
     154        c = i.coeff();
     155        if ( c.inCoeffDomain())
     156        {
     157          if ( c > qh )
     158            result += power( x, i.exp() ) * (c - q);
     159          else
     160            result += power( x, i.exp() ) * c;
     161        }
     162        else
     163          result += power( x, i.exp() ) * balance_p(c,q,qh);
     164    }
     165    return result;
     166}
     167
     168static CanonicalForm
     169balance_p ( const CanonicalForm & f, const CanonicalForm & q )
     170{
     171    CanonicalForm qh = q / 2;
     172    return balance_p (f, q, qh);
     173}
     174
     175static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
     176{
     177  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
     178  int p, i;
     179
     180  if ( primitive )
     181  {
     182    f = F;
     183    g = G;
     184    c = 1;
     185  }
     186  else
     187  {
     188    CanonicalForm cF = content( F ), cG = content( G );
     189    f = F / cF;
     190    g = G / cG;
     191    c = bgcd( cF, cG );
     192  }
     193  cg = gcd( f.lc(), g.lc() );
     194  cl = ( f.lc() / cg ) * g.lc();
     195//     B = 2 * cg * tmin(
     196//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
     197//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
     198//         )+1;
     199  M = tmin( maxNorm(f), maxNorm(g) );
     200  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
     201  q = 0;
     202  i = cf_getNumSmallPrimes() - 1;
     203  while ( true )
     204  {
     205    B = BB;
     206    while ( i >= 0 && q < B )
     207    {
     208      p = cf_getSmallPrime( i );
     209      i--;
     210      while ( i >= 0 && mod( cl, p ) == 0 )
     211      {
     212        p = cf_getSmallPrime( i );
     213        i--;
     214      }
     215      setCharacteristic( p );
     216      Dp = gcd( mapinto( f ), mapinto( g ) );
     217      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
     218      setCharacteristic( 0 );
     219      if ( Dp.degree() == 0 )
     220        return c;
     221      if ( q.isZero() )
     222      {
     223        D = mapinto( Dp );
     224        q = p;
     225        B = power(CanonicalForm(2),D.degree())*M+1;
     226      }
     227      else
     228      {
     229        if ( Dp.degree() == D.degree() )
     230        {
     231          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
     232          q = newq;
     233          D = newD;
     234        }
     235        else if ( Dp.degree() < D.degree() )
     236        {
     237          // all previous p's are bad primes
     238          q = p;
     239          D = mapinto( Dp );
     240          B = power(CanonicalForm(2),D.degree())*M+1;
     241        }
     242        // else p is a bad prime
     243      }
     244    }
     245    if ( i >= 0 )
     246    {
     247      // now balance D mod q
     248      D = pp( balance_p( D, q ) );
     249      if ( fdivides( D, f ) && fdivides( D, g ) )
     250        return D * c;
     251      else
     252        q = 0;
     253    }
     254    else
     255      return gcd_poly( F, G );
     256    DEBOUTLN( cerr, "another try ..." );
     257  }
     258}
     259
     260static CanonicalForm
     261gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
     262{
     263    if (f.inCoeffDomain() || g.inCoeffDomain()) //zero case should be caught by gcd
     264      return 1;
     265    CanonicalForm pi, pi1;
     266    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
     267    bool bpure, ezgcdon= isOn (SW_USE_EZGCD_P);
     268    int delta = degree( f ) - degree( g );
     269
     270    if ( delta >= 0 )
     271    {
     272        pi = f; pi1 = g;
     273    }
     274    else
     275    {
     276        pi = g; pi1 = f; delta = -delta;
     277    }
     278    if (pi.isUnivariate())
     279      Ci= 1;
     280    else
     281    {
     282      if (!ezgcdon)
     283        On (SW_USE_EZGCD_P);
     284      Ci = content( pi );
     285      if (!ezgcdon)
     286        Off (SW_USE_EZGCD_P);
     287      pi = pi / Ci;
     288    }
     289    if (pi1.isUnivariate())
     290      Ci1= 1;
     291    else
     292    {
     293      if (!ezgcdon)
     294        On (SW_USE_EZGCD_P);
     295      Ci1 = content( pi1 );
     296      if (!ezgcdon)
     297        Off (SW_USE_EZGCD_P);
     298      pi1 = pi1 / Ci1;
     299    }
     300    C = gcd( Ci, Ci1 );
     301    int d= 0;
     302    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
     303    {
     304        if ( gcd_test_one( pi1, pi, true, d ) )
     305        {
     306          C=abs(C);
     307          //out_cf("GCD:",C,"\n");
     308          return C;
     309        }
     310        bpure = false;
     311    }
     312    else
     313    {
     314        bpure = isPurePoly(pi) && isPurePoly(pi1);
     315#ifdef HAVE_FLINT
     316        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
     317          return gcd_univar_flintp(pi,pi1)*C;
     318#else
     319#ifdef HAVE_NTL
     320        if ( bpure && (CFFactory::gettype() != GaloisFieldDomain))
     321            return gcd_univar_ntlp(pi, pi1 ) * C;
     322#endif
     323#endif
     324    }
     325    Variable v = f.mvar();
     326    Hi = power( LC( pi1, v ), delta );
     327    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
     328
     329    if (!(pi.isUnivariate() && pi1.isUnivariate()))
     330    {
     331      if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
     332      {
     333        On (SW_USE_FF_MOD_GCD);
     334        C *= gcd (pi, pi1);
     335        Off (SW_USE_FF_MOD_GCD);
     336        return C;
     337      }
     338    }
     339
     340    if ( (delta+1) % 2 )
     341        bi = 1;
     342    else
     343        bi = -1;
     344    CanonicalForm oldPi= pi, oldPi1= pi1, powHi;
     345    while ( degree( pi1, v ) > 0 )
     346    {
     347        if (!(pi.isUnivariate() && pi1.isUnivariate()))
     348        {
     349          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
     350          {
     351            On (SW_USE_FF_MOD_GCD);
     352            C *= gcd (oldPi, oldPi1);
     353            Off (SW_USE_FF_MOD_GCD);
     354            return C;
     355          }
     356        }
     357        pi2 = psr( pi, pi1, v );
     358        pi2 = pi2 / bi;
     359        pi = pi1; pi1 = pi2;
     360        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
     361        if (!pi1.isUnivariate() && (size (pi1)/maxNumVars > 500))
     362        {
     363            On (SW_USE_FF_MOD_GCD);
     364            C *= gcd (oldPi, oldPi1);
     365            Off (SW_USE_FF_MOD_GCD);
     366            return C;
     367        }
     368        if ( degree( pi1, v ) > 0 )
     369        {
     370            delta = degree( pi, v ) - degree( pi1, v );
     371            powHi= power (Hi, delta-1);
     372            if ( (delta+1) % 2 )
     373                bi = LC( pi, v ) * powHi*Hi;
     374            else
     375                bi = -LC( pi, v ) * powHi*Hi;
     376            Hi = power( LC( pi1, v ), delta ) / powHi;
     377            if (!(pi.isUnivariate() && pi1.isUnivariate()))
     378            {
     379              if (size (Hi)*size (pi)/(maxNumVars*3) > 1500) //maybe this needs more tuning
     380              {
     381                On (SW_USE_FF_MOD_GCD);
     382                C *= gcd (oldPi, oldPi1);
     383                Off (SW_USE_FF_MOD_GCD);
     384                return C;
     385              }
     386            }
     387        }
     388    }
     389    if ( degree( pi1, v ) == 0 )
     390    {
     391      C=abs(C);
     392      //out_cf("GCD:",C,"\n");
     393      return C;
     394    }
     395    if (!pi.isUnivariate())
     396    {
     397      if (!ezgcdon)
     398        On (SW_USE_EZGCD_P);
     399      Ci= gcd (LC (oldPi,v), LC (oldPi1,v));
     400      pi /= LC (pi,v)/Ci;
     401      Ci= content (pi);
     402      pi /= Ci;
     403      if (!ezgcdon)
     404        Off (SW_USE_EZGCD_P);
     405    }
     406    if ( bpure )
     407        pi /= pi.lc();
     408    C=abs(C*pi);
     409    //out_cf("GCD:",C,"\n");
     410    return C;
     411}
     412
     413static CanonicalForm
     414gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
     415{
     416    CanonicalForm pi, pi1;
     417    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
     418    int delta = degree( f ) - degree( g );
     419
     420    if ( delta >= 0 )
     421    {
     422        pi = f; pi1 = g;
     423    }
     424    else
     425    {
     426        pi = g; pi1 = f; delta = -delta;
     427    }
     428    Ci = content( pi ); Ci1 = content( pi1 );
     429    pi1 = pi1 / Ci1; pi = pi / Ci;
     430    C = gcd( Ci, Ci1 );
     431    int d= 0;
     432    if ( pi.isUnivariate() && pi1.isUnivariate() )
     433    {
     434#ifdef HAVE_FLINT
     435        if (isPurePoly(pi) && isPurePoly(pi1) )
     436            return gcd_univar_flint0(pi, pi1 ) * C;
     437#else
     438#ifdef HAVE_NTL
     439        if ( isPurePoly(pi) && isPurePoly(pi1) )
     440            return gcd_univar_ntl0(pi, pi1 ) * C;
     441#endif
     442#endif
     443        return gcd_poly_univar0( pi, pi1, true ) * C;
     444    }
     445    else if ( gcd_test_one( pi1, pi, true, d ) )
     446      return C;
     447    Variable v = f.mvar();
     448    Hi = power( LC( pi1, v ), delta );
     449    if ( (delta+1) % 2 )
     450        bi = 1;
     451    else
     452        bi = -1;
     453    while ( degree( pi1, v ) > 0 )
     454    {
     455        pi2 = psr( pi, pi1, v );
     456        pi2 = pi2 / bi;
     457        pi = pi1; pi1 = pi2;
     458        if ( degree( pi1, v ) > 0 )
     459        {
     460            delta = degree( pi, v ) - degree( pi1, v );
     461            if ( (delta+1) % 2 )
     462                bi = LC( pi, v ) * power( Hi, delta );
     463            else
     464                bi = -LC( pi, v ) * power( Hi, delta );
     465            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
     466        }
     467    }
     468    if ( degree( pi1, v ) == 0 )
     469        return C;
     470    else
     471        return C * pp( pi );
    76472}
    77473
     
    130526    #endif
    131527    else
    132     fc = subResGCD_p( fc, gc );
     528    fc = gcd_poly_p( fc, gc );
    133529  }
    134530  else if (!fc_and_gc_Univariate) /* && char==0*/
     
    145541    if ( isOn( SW_USE_EZGCD ) )
    146542      fc= ezgcd (fc, gc);
    147     #endif 
     543    #endif
    148544    #ifdef HAVE_NTL
    149545    else if (isOn(SW_USE_CHINREM_GCD))
     
    152548    #endif
    153549    {
    154        fc = subResGCD_0( fc, gc );
     550       fc = gcd_poly_0( fc, gc );
    155551    }
    156552  }
    157553  else
    158554  {
    159     fc = subResGCD_0( fc, gc );
     555    fc = gcd_poly_0( fc, gc );
    160556  }
    161557  if ((getCharacteristic()>0)&&(!hasAlgVar(fc))) fc/=fc.lc();
Note: See TracChangeset for help on using the changeset viewer.