/* emacs edit mode for this file is -*- C++ -*- */ #include "config.h" #include "timing.h" #include "cf_assert.h" #include "debug.h" #include "cf_defs.h" #include "canonicalform.h" #include "cf_iter.h" #include "cf_reval.h" #include "cf_primes.h" #include "cf_algorithm.h" #include "cf_factory.h" #include "fac_util.h" #include "templates/ftmpl_functions.h" #include "algext.h" #include "cf_gcd_smallp.h" #include "cf_map_ext.h" #include "cf_util.h" #include "gfops.h" #ifdef HAVE_NTL #include #include "NTLconvert.h" bool isPurePoly(const CanonicalForm & ); #ifndef HAVE_FLINT static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & ); static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & ); #endif #endif #ifdef HAVE_FLINT #include "FLINTconvert.h" static CanonicalForm gcd_univar_flint0 (const CanonicalForm &, const CanonicalForm &); static CanonicalForm gcd_univar_flintp (const CanonicalForm &, const CanonicalForm &); #endif static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & ); void out_cf(const char *s1,const CanonicalForm &f,const char *s2); CanonicalForm chinrem_gcd(const CanonicalForm & FF,const CanonicalForm & GG); bool gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap, int & d ) { d= 0; int count = 0; // assume polys have same level; Variable v= Variable (1); bool algExtension= (hasFirstAlgVar (f, v) || hasFirstAlgVar (g, v)); CanonicalForm lcf, lcg; if ( swap ) { lcf = swapvar( LC( f ), Variable(1), f.mvar() ); lcg = swapvar( LC( g ), Variable(1), f.mvar() ); } else { lcf = LC( f, Variable(1) ); lcg = LC( g, Variable(1) ); } CanonicalForm F, G; if ( swap ) { F=swapvar( f, Variable(1), f.mvar() ); G=swapvar( g, Variable(1), g.mvar() ); } else { F = f; G = g; } #define TEST_ONE_MAX 50 int p= getCharacteristic(); bool passToGF= false; int k= 1; if (p > 0 && p < TEST_ONE_MAX && CFFactory::gettype() != GaloisFieldDomain && !algExtension) { if (p == 2) setCharacteristic (2, 6, 'Z'); else if (p == 3) setCharacteristic (3, 4, 'Z'); else if (p == 5 || p == 7) setCharacteristic (p, 3, 'Z'); else setCharacteristic (p, 2, 'Z'); passToGF= true; } else if (p > 0 && CFFactory::gettype() == GaloisFieldDomain && ipower (p , getGFDegree()) < TEST_ONE_MAX) { k= getGFDegree(); if (ipower (p, 2*k) > TEST_ONE_MAX) setCharacteristic (p, 2*k, gf_name); else setCharacteristic (p, 3*k, gf_name); F= GFMapUp (F, k); G= GFMapUp (G, k); lcf= GFMapUp (lcf, k); lcg= GFMapUp (lcg, k); } else if (p > 0 && p < TEST_ONE_MAX && algExtension) { bool extOfExt= false; #ifdef HAVE_NTL int d= degree (getMipo (v)); CFList source, dest; Variable v2; CanonicalForm primElem, imPrimElem; if (p == 2 && d < 6) { if (fac_NTL_char != 2) { fac_NTL_char= 2; zz_p::init (p); } bool primFail= false; Variable vBuf; primElem= primitiveElement (v, vBuf, primFail); ASSERT (!primFail, "failure in integer factorizer"); if (d < 3) { zz_pX NTLIrredpoly; BuildIrred (NTLIrredpoly, d*3); CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); v2= rootOf (newMipo); } else { zz_pX NTLIrredpoly; BuildIrred (NTLIrredpoly, d*2); CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); v2= rootOf (newMipo); } imPrimElem= mapPrimElem (primElem, v, v2); extOfExt= true; } else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3)) { if (fac_NTL_char != p) { fac_NTL_char= p; zz_p::init (p); } bool primFail= false; Variable vBuf; primElem= primitiveElement (v, vBuf, primFail); ASSERT (!primFail, "failure in integer factorizer"); zz_pX NTLIrredpoly; BuildIrred (NTLIrredpoly, d*2); CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1)); v2= rootOf (newMipo); imPrimElem= mapPrimElem (primElem, v, v2); extOfExt= true; } if (extOfExt) { F= mapUp (F, v, v2, primElem, imPrimElem, source, dest); G= mapUp (G, v, v2, primElem, imPrimElem, source, dest); lcf= mapUp (lcf, v, v2, primElem, imPrimElem, source, dest); lcg= mapUp (lcg, v, v2, primElem, imPrimElem, source, dest); v= v2; } #endif } CFRandom * sample; if ((!algExtension && p > 0) || p == 0) sample = CFRandomFactory::generate(); else sample = AlgExtRandomF (v).clone(); REvaluation e( 2, tmax( f.level(), g.level() ), *sample ); delete sample; if (passToGF) { lcf= lcf.mapinto(); lcg= lcg.mapinto(); } CanonicalForm eval1, eval2; if (passToGF) { eval1= e (lcf); eval2= e (lcg); } else { eval1= e (lcf); eval2= e (lcg); } while ( ( eval1.isZero() || eval2.isZero() ) && count < TEST_ONE_MAX ) { e.nextpoint(); count++; eval1= e (lcf); eval2= e (lcg); } if ( count >= TEST_ONE_MAX ) { if (passToGF) setCharacteristic (p); if (k > 1) setCharacteristic (p, k, gf_name); return false; } if (passToGF) { F= F.mapinto(); G= G.mapinto(); eval1= e (F); eval2= e (G); } else { eval1= e (F); eval2= e (G); } CanonicalForm c= gcd (eval1, eval2); d= c.degree(); bool result= d < 1; if (d < 0) d= 0; if (passToGF) setCharacteristic (p); if (k > 1) setCharacteristic (p, k, gf_name); return result; } //{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c ) //{{{ docu // // icontent() - return gcd of c and all coefficients of f which // are in a coefficient domain. // // Used by icontent(). // //}}} static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c ) { if ( f.inBaseDomain() ) { if (c.isZero()) return abs(f); return bgcd( f, c ); } //else if ( f.inCoeffDomain() ) // return gcd(f,c); else { CanonicalForm g = c; for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ ) g = icontent( i.coeff(), g ); return g; } } //}}} //{{{ CanonicalForm icontent ( const CanonicalForm & f ) //{{{ docu // // icontent() - return gcd over all coefficients of f which are // in a coefficient domain. // //}}} CanonicalForm icontent ( const CanonicalForm & f ) { return icontent( f, 0 ); } //}}} //{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b ) //{{{ docu // // extgcd() - returns polynomial extended gcd of f and g. // // Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g). // The gcd is calculated using an extended euclidean polynomial // remainder sequence, so f and g should be polynomials over an // euclidean domain. Normalizes result. // // Note: be sure that f and g have the same level! // //}}} CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b ) { if (f.isZero()) { a= 0; b= 1; return g; } else if (g.isZero()) { a= 1; b= 0; return f; } #ifdef HAVE_NTL #ifdef HAVE_FLINT if (( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain) && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g)) { nmod_poly_t F1, G1, A, B, R; convertFacCF2nmod_poly_t (F1, f); convertFacCF2nmod_poly_t (G1, g); nmod_poly_init (R, getCharacteristic()); nmod_poly_init (A, getCharacteristic()); nmod_poly_init (B, getCharacteristic()); nmod_poly_xgcd (R, A, B, F1, G1); a= convertnmod_poly_t2FacCF (A, f.mvar()); b= convertnmod_poly_t2FacCF (B, f.mvar()); CanonicalForm r= convertnmod_poly_t2FacCF (R, f.mvar()); nmod_poly_clear (F1); nmod_poly_clear (G1); nmod_poly_clear (A); nmod_poly_clear (B); nmod_poly_clear (R); return r; } #else if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain) && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g)) { if (fac_NTL_char!=getCharacteristic()) { fac_NTL_char=getCharacteristic(); #ifdef NTL_ZZ ZZ r; r=getCharacteristic(); ZZ_pContext ccc(r); #else zz_pContext ccc(getCharacteristic()); #endif ccc.restore(); #ifdef NTL_ZZ ZZ_p::init(r); #else zz_p::init(getCharacteristic()); #endif } #ifdef NTL_ZZ ZZ_pX F1=convertFacCF2NTLZZpX(f); ZZ_pX G1=convertFacCF2NTLZZpX(g); ZZ_pX R; ZZ_pX A,B; XGCD(R,A,B,F1,G1); a=convertNTLZZpX2CF(A,f.mvar()); b=convertNTLZZpX2CF(B,f.mvar()); return convertNTLZZpX2CF(R,f.mvar()); #else zz_pX F1=convertFacCF2NTLzzpX(f); zz_pX G1=convertFacCF2NTLzzpX(g); zz_pX R; zz_pX A,B; XGCD(R,A,B,F1,G1); a=convertNTLzzpX2CF(A,f.mvar()); b=convertNTLzzpX2CF(B,f.mvar()); return convertNTLzzpX2CF(R,f.mvar()); #endif } #endif #ifdef HAVE_FLINT if (( getCharacteristic() ==0) && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g)) { fmpq_poly_t F1, G1; convertFacCF2Fmpq_poly_t (F1, f); convertFacCF2Fmpq_poly_t (G1, g); fmpq_poly_t R, A, B; fmpq_poly_init (R); fmpq_poly_init (A); fmpq_poly_init (B); fmpq_poly_xgcd (R, A, B, F1, G1); a= convertFmpq_poly_t2FacCF (A, f.mvar()); b= convertFmpq_poly_t2FacCF (B, f.mvar()); CanonicalForm r= convertFmpq_poly_t2FacCF (R, f.mvar()); fmpq_poly_clear (F1); fmpq_poly_clear (G1); fmpq_poly_clear (A); fmpq_poly_clear (B); fmpq_poly_clear (R); return r; } #else if (isOn(SW_USE_NTL_GCD_0) && ( getCharacteristic() ==0) && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g)) { CanonicalForm fc=bCommonDen(f); CanonicalForm gc=bCommonDen(g); ZZX F1=convertFacCF2NTLZZX(f*fc); ZZX G1=convertFacCF2NTLZZX(g*gc); ZZX R=GCD(F1,G1); CanonicalForm r=convertNTLZZX2CF(R,f.mvar()); ZZ RR; ZZX A,B; if (r.inCoeffDomain()) { XGCD(RR,A,B,F1,G1,1); CanonicalForm rr=convertZZ2CF(RR); ASSERT (!rr.isZero(), "NTL:XGCD failed"); a=convertNTLZZX2CF(A,f.mvar())*fc/rr; b=convertNTLZZX2CF(B,f.mvar())*gc/rr; return CanonicalForm(1); } else { fc=bCommonDen(f); gc=bCommonDen(g); F1=convertFacCF2NTLZZX(f*fc/r); G1=convertFacCF2NTLZZX(g*gc/r); XGCD(RR,A,B,F1,G1,1); a=convertNTLZZX2CF(A,f.mvar())*fc; b=convertNTLZZX2CF(B,f.mvar())*gc; CanonicalForm rr=convertZZ2CF(RR); ASSERT (!rr.isZero(), "NTL:XGCD failed"); a /= rr; b /= rr; return r; } } #endif #endif // may contain bug in the co-factors, see track 107 CanonicalForm contf = content( f ); CanonicalForm contg = content( g ); CanonicalForm p0 = f / contf, p1 = g / contg; CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r; while ( ! p1.isZero() ) { divrem( p0, p1, q, r ); p0 = p1; p1 = r; r = g0 - g1 * q; g0 = g1; g1 = r; r = f0 - f1 * q; f0 = f1; f1 = r; } CanonicalForm contp0 = content( p0 ); a = f0 / ( contf * contp0 ); b = g0 / ( contg * contp0 ); p0 /= contp0; if ( p0.sign() < 0 ) { p0 = -p0; a = -a; b = -b; } return p0; } //}}} //{{{ static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q ) //{{{ docu // // balance_p() - map f from positive to symmetric representation // mod q. // // This makes sense for univariate polynomials over Z only. // q should be an integer. // // Used by gcd_poly_univar0(). // //}}} static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q ) { Variable x = f.mvar(); CanonicalForm result = 0, qh = q / 2; CanonicalForm c; CFIterator i; for ( i = f; i.hasTerms(); i++ ) { c = i.coeff(); if ( c.inCoeffDomain()) { if ( c > qh ) result += power( x, i.exp() ) * (c - q); else result += power( x, i.exp() ) * c; } else result += power( x, i.exp() ) * balance_p(c,q); } return result; } /*static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q ) { Variable x = f.mvar(); CanonicalForm result = 0, qh = q / 2; CanonicalForm c; CFIterator i; for ( i = f; i.hasTerms(); i++ ) { c = mod( i.coeff(), q ); if ( c > qh ) result += power( x, i.exp() ) * (c - q); else result += power( x, i.exp() ) * c; } return result; }*/ //}}} static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive ) { CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq; int p, i; if ( primitive ) { f = F; g = G; c = 1; } else { CanonicalForm cF = content( F ), cG = content( G ); f = F / cF; g = G / cG; c = bgcd( cF, cG ); } cg = gcd( f.lc(), g.lc() ); cl = ( f.lc() / cg ) * g.lc(); // B = 2 * cg * tmin( // maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1), // maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1) // )+1; M = tmin( maxNorm(f), maxNorm(g) ); BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M; q = 0; i = cf_getNumSmallPrimes() - 1; while ( true ) { B = BB; while ( i >= 0 && q < B ) { p = cf_getSmallPrime( i ); i--; while ( i >= 0 && mod( cl, p ) == 0 ) { p = cf_getSmallPrime( i ); i--; } setCharacteristic( p ); Dp = gcd( mapinto( f ), mapinto( g ) ); Dp = ( Dp / Dp.lc() ) * mapinto( cg ); setCharacteristic( 0 ); if ( Dp.degree() == 0 ) return c; if ( q.isZero() ) { D = mapinto( Dp ); q = p; B = power(CanonicalForm(2),D.degree())*M+1; } else { if ( Dp.degree() == D.degree() ) { chineseRemainder( D, q, mapinto( Dp ), p, newD, newq ); q = newq; D = newD; } else if ( Dp.degree() < D.degree() ) { // all previous p's are bad primes q = p; D = mapinto( Dp ); B = power(CanonicalForm(2),D.degree())*M+1; } // else p is a bad prime } } if ( i >= 0 ) { // now balance D mod q D = pp( balance_p( D, q ) ); if ( fdivides( D, f ) && fdivides( D, g ) ) return D * c; else q = 0; } else return gcd_poly( F, G ); DEBOUTLN( cerr, "another try ..." ); } } static CanonicalForm gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g ) { CanonicalForm pi, pi1; CanonicalForm C, Ci, Ci1, Hi, bi, pi2; bool bpure; int delta = degree( f ) - degree( g ); if ( delta >= 0 ) { pi = f; pi1 = g; } else { pi = g; pi1 = f; delta = -delta; } On (SW_USE_EZGCD_P); Ci = content( pi ); Ci1 = content( pi1 ); Off (SW_USE_EZGCD_P); pi1 = pi1 / Ci1; pi = pi / Ci; C = gcd( Ci, Ci1 ); int d= 0; if ( !( pi.isUnivariate() && pi1.isUnivariate() ) ) { if ( gcd_test_one( pi1, pi, true, d ) ) { C=abs(C); //out_cf("GCD:",C,"\n"); return C; } bpure = false; } else { bpure = isPurePoly(pi) && isPurePoly(pi1); #ifdef HAVE_FLINT if (bpure && (CFFactory::gettype() != GaloisFieldDomain)) return gcd_univar_flintp(pi,pi1)*C; #else #ifdef HAVE_NTL if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain)) return gcd_univar_ntlp(pi, pi1 ) * C; #endif #endif } Variable v = f.mvar(); Hi = power( LC( pi1, v ), delta ); int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1)); if (!(pi.isUnivariate() && pi1.isUnivariate())) { if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning { On (SW_USE_FF_MOD_GCD); C *= gcd (pi, pi1); Off (SW_USE_FF_MOD_GCD); return C; } } if ( (delta+1) % 2 ) bi = 1; else bi = -1; CanonicalForm oldPi= pi, oldPi1= pi1, powHi; while ( degree( pi1, v ) > 0 ) { if (!(pi.isUnivariate() && pi1.isUnivariate())) { if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500) { On (SW_USE_FF_MOD_GCD); C *= gcd (oldPi, oldPi1); Off (SW_USE_FF_MOD_GCD); return C; } } pi2 = psr( pi, pi1, v ); pi2 = pi2 / bi; pi = pi1; pi1 = pi2; maxNumVars= tmax (getNumVars (pi), getNumVars (pi1)); if (!(pi1.isUnivariate()) && (size (pi1)/maxNumVars > 500)) { On (SW_USE_FF_MOD_GCD); C *= gcd (oldPi, oldPi1); Off (SW_USE_FF_MOD_GCD); return C; } if ( degree( pi1, v ) > 0 ) { delta = degree( pi, v ) - degree( pi1, v ); powHi= power (Hi, delta-1); if ( (delta+1) % 2 ) bi = LC( pi, v ) * powHi*Hi; else bi = -LC( pi, v ) * powHi*Hi; Hi = power( LC( pi1, v ), delta ) / powHi; if (!(pi.isUnivariate() && pi1.isUnivariate())) { if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning { On (SW_USE_FF_MOD_GCD); C *= gcd (oldPi, oldPi1); Off (SW_USE_FF_MOD_GCD); return C; } } } } if ( degree( pi1, v ) == 0 ) { C=abs(C); //out_cf("GCD:",C,"\n"); return C; } pi /= content( pi ); if ( bpure ) pi /= pi.lc(); C=abs(C*pi); //out_cf("GCD:",C,"\n"); return C; } static CanonicalForm gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g ) { CanonicalForm pi, pi1; CanonicalForm C, Ci, Ci1, Hi, bi, pi2; int delta = degree( f ) - degree( g ); if ( delta >= 0 ) { pi = f; pi1 = g; } else { pi = g; pi1 = f; delta = -delta; } Ci = content( pi ); Ci1 = content( pi1 ); pi1 = pi1 / Ci1; pi = pi / Ci; C = gcd( Ci, Ci1 ); int d= 0; if ( pi.isUnivariate() && pi1.isUnivariate() ) { #ifdef HAVE_FLINT if (isPurePoly(pi) && isPurePoly(pi1) ) return gcd_univar_flint0(pi, pi1 ) * C; #else #ifdef HAVE_NTL if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) ) return gcd_univar_ntl0(pi, pi1 ) * C; #endif #endif return gcd_poly_univar0( pi, pi1, true ) * C; } else if ( gcd_test_one( pi1, pi, true, d ) ) return C; Variable v = f.mvar(); Hi = power( LC( pi1, v ), delta ); if ( (delta+1) % 2 ) bi = 1; else bi = -1; while ( degree( pi1, v ) > 0 ) { pi2 = psr( pi, pi1, v ); pi2 = pi2 / bi; pi = pi1; pi1 = pi2; if ( degree( pi1, v ) > 0 ) { delta = degree( pi, v ) - degree( pi1, v ); if ( (delta+1) % 2 ) bi = LC( pi, v ) * power( Hi, delta ); else bi = -LC( pi, v ) * power( Hi, delta ); Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 ); } } if ( degree( pi1, v ) == 0 ) return C; else return C * pp( pi ); } //{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g ) //{{{ docu // // gcd_poly() - calculate polynomial gcd. // // This is the dispatcher for polynomial gcd calculation. We call either // ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current // characteristic and settings of SW_USE_EZGCD. // // Used by gcd() and gcd_poly_univar0(). // //}}} #if 0 int si_factor_reminder=1; #endif CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g ) { CanonicalForm fc, gc, d1; bool fc_isUnivariate=f.isUnivariate(); bool gc_isUnivariate=g.isUnivariate(); bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate; fc = f; gc = g; if ( getCharacteristic() != 0 ) { #ifdef HAVE_NTL if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P ))) { fc= EZGCD_P (fc, gc); } else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate) { Variable a; if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a)) fc=GCD_Fp_extension (fc, gc, a); else if (CFFactory::gettype() == GaloisFieldDomain) fc=GCD_GF (fc, gc); else fc=GCD_small_p (fc, gc); } else #endif fc = gcd_poly_p( fc, gc ); } else if (!fc_and_gc_Univariate) { if ( isOn( SW_USE_EZGCD ) ) fc= ezgcd (fc, gc); else if (isOn(SW_USE_CHINREM_GCD)) fc = chinrem_gcd( fc, gc); else { fc = gcd_poly_0( fc, gc ); } } else { fc = gcd_poly_0( fc, gc ); } if ( d1.degree() > 0 ) fc *= d1; return fc; } //}}} //{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g ) //{{{ docu // // cf_content() - return gcd(g, content(f)). // // content(f) is calculated with respect to f's main variable. // // Used by gcd(), content(), content( CF, Variable ). // //}}} static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g ) { if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) { CFIterator i = f; CanonicalForm result = g; while ( i.hasTerms() && ! result.isOne() ) { result = gcd( i.coeff(), result ); i++; } return result; } else return abs( f ); } //}}} //{{{ CanonicalForm content ( const CanonicalForm & f ) //{{{ docu // // content() - return content(f) with respect to main variable. // // Normalizes result. // //}}} CanonicalForm content ( const CanonicalForm & f ) { if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) { CFIterator i = f; CanonicalForm result = abs( i.coeff() ); i++; while ( i.hasTerms() && ! result.isOne() ) { result = gcd( i.coeff(), result ); i++; } return result; } else return abs( f ); } //}}} //{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x ) //{{{ docu // // content() - return content(f) with respect to x. // // x should be a polynomial variable. // //}}} CanonicalForm content ( const CanonicalForm & f, const Variable & x ) { ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" ); Variable y = f.mvar(); if ( y == x ) return cf_content( f, 0 ); else if ( y < x ) return f; else return swapvar( content( swapvar( f, y, x ), y ), y, x ); } //}}} //{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x ) //{{{ docu // // vcontent() - return content of f with repect to variables >= x. // // The content is recursively calculated over all coefficients in // f having level less than x. x should be a polynomial // variable. // //}}} CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x ) { ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" ); if ( f.mvar() <= x ) return content( f, x ); else { CFIterator i; CanonicalForm d = 0; for ( i = f; i.hasTerms() && ! d.isOne(); i++ ) d = gcd( d, vcontent( i.coeff(), x ) ); return d; } } //}}} //{{{ CanonicalForm pp ( const CanonicalForm & f ) //{{{ docu // // pp() - return primitive part of f. // // Returns zero if f equals zero, otherwise f / content(f). // //}}} CanonicalForm pp ( const CanonicalForm & f ) { if ( f.isZero() ) return f; else return f / content( f ); } //}}} CanonicalForm gcd ( const CanonicalForm & f, const CanonicalForm & g ) { bool b = f.isZero(); if ( b || g.isZero() ) { if ( b ) return abs( g ); else return abs( f ); } if ( f.inPolyDomain() || g.inPolyDomain() ) { if ( f.mvar() != g.mvar() ) { if ( f.mvar() > g.mvar() ) return cf_content( f, g ); else return cf_content( g, f ); } if (isOn(SW_USE_QGCD)) { Variable m; if ( (getCharacteristic() == 0) && (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m)) ) { bool on_rational = isOn(SW_RATIONAL); CanonicalForm r=QGCD(f,g); On(SW_RATIONAL); CanonicalForm cdF = bCommonDen( r ); if (!on_rational) Off(SW_RATIONAL); return cdF*r; } } if ( f.inExtension() && getReduce( f.mvar() ) ) return CanonicalForm(1); else { if ( fdivides( f, g ) ) return abs( f ); else if ( fdivides( g, f ) ) return abs( g ); if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) ) { CanonicalForm d; d = gcd_poly( f, g ); return abs( d ); } else { //printf ("here\n"); CanonicalForm cdF = bCommonDen( f ); CanonicalForm cdG = bCommonDen( g ); Off( SW_RATIONAL ); CanonicalForm l = lcm( cdF, cdG ); On( SW_RATIONAL ); CanonicalForm F = f * l, G = g * l; Off( SW_RATIONAL ); l = gcd_poly( F, G ); On( SW_RATIONAL ); return abs( l ); } } } if ( f.inBaseDomain() && g.inBaseDomain() ) return bgcd( f, g ); else return 1; } //{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g ) //{{{ docu // // lcm() - return least common multiple of f and g. // // The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g). // // Returns zero if one of f or g equals zero. // //}}} CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g ) { if ( f.isZero() || g.isZero() ) return 0; else return ( f / gcd( f, g ) ) * g; } //}}} #ifdef HAVE_NTL #ifndef HAVE_FLINT static CanonicalForm gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G ) { ZZX F1=convertFacCF2NTLZZX(F); ZZX G1=convertFacCF2NTLZZX(G); ZZX R=GCD(F1,G1); return convertNTLZZX2CF(R,F.mvar()); } static CanonicalForm gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G ) { if (fac_NTL_char!=getCharacteristic()) { fac_NTL_char=getCharacteristic(); #ifdef NTL_ZZ ZZ r; r=getCharacteristic(); ZZ_pContext ccc(r); #else zz_pContext ccc(getCharacteristic()); #endif ccc.restore(); #ifdef NTL_ZZ ZZ_p::init(r); #else zz_p::init(getCharacteristic()); #endif } #ifdef NTL_ZZ ZZ_pX F1=convertFacCF2NTLZZpX(F); ZZ_pX G1=convertFacCF2NTLZZpX(G); ZZ_pX R=GCD(F1,G1); return convertNTLZZpX2CF(R,F.mvar()); #else zz_pX F1=convertFacCF2NTLzzpX(F); zz_pX G1=convertFacCF2NTLzzpX(G); zz_pX R=GCD(F1,G1); return convertNTLzzpX2CF(R,F.mvar()); #endif } #endif #endif #ifdef HAVE_FLINT static CanonicalForm gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G) { nmod_poly_t F1, G1; convertFacCF2nmod_poly_t (F1, F); convertFacCF2nmod_poly_t (G1, G); nmod_poly_gcd (F1, F1, G1); CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar()); nmod_poly_clear (F1); nmod_poly_clear (G1); return result; } static CanonicalForm gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G ) { fmpz_poly_t F1, G1; convertFacCF2Fmpz_poly_t(F1, F); convertFacCF2Fmpz_poly_t(G1, G); fmpz_poly_gcd (F1, F1, G1); CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar()); fmpz_poly_clear (F1); fmpz_poly_clear (G1); return result; } #endif /* * compute positions p1 and pe of optimal variables: * pe is used in "ezgcd" and * p1 in "gcd_poly1" */ /*static void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe ) { int i, o1, oe; if ( df[n] > dg[n] ) { o1 = df[n]; oe = dg[n]; } else { o1 = dg[n]; oe = df[n]; } i = n-1; while ( i > 0 ) { if ( df[i] != 0 ) { if ( df[i] > dg[i] ) { if ( o1 > df[i]) { o1 = df[i]; p1 = i; } if ( oe <= dg[i]) { oe = dg[i]; pe = i; } } else { if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; } if ( oe <= df[i]) { oe = df[i]; pe = i; } } } i--; } }*/ /* * make some changes of variables, see optvalues */ /*static void cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe ) { int i, k, n; n = f.level(); cc = 0; p1 = pe = n; if ( n == 1 ) return; int * degsf = new int[n+1]; int * degsg = new int[n+1]; for ( i = n; i > 0; i-- ) { degsf[i] = degsg[i] = 0; } degsf = degrees( f, degsf ); degsg = degrees( g, degsg ); k = 0; for ( i = n-1; i > 0; i-- ) { if ( degsf[i] == 0 ) { if ( degsg[i] != 0 ) { cc = -i; break; } } else { if ( degsg[i] == 0 ) { cc = i; break; } else k++; } } if ( ( cc == 0 ) && ( k != 0 ) ) optvalues( degsf, degsg, n, p1, pe ); if ( ( pe != 1 ) && ( degsf[1] != 0 ) ) pe = -pe; delete [] degsf; delete [] degsg; }*/ TIMING_DEFINE_PRINT(chinrem_termination) TIMING_DEFINE_PRINT(chinrem_recursion) TIMING_DEFINE_PRINT(chinrem_reconstruction) CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG ) { CanonicalForm f, g, cl, q(0), Dp, newD, D, newq; int p, i, dp_deg, d_deg=-1; CanonicalForm cd ( bCommonDen( FF )); f=cd*FF; Variable x= Variable (1); CanonicalForm cf, cg; cf= icontent (f); f /= cf; //cd = bCommonDen( f ); f *=cd; //f /=vcontent(f,Variable(1)); cd = bCommonDen( GG ); g=cd*GG; cg= icontent (g); g /= cg; //cd = bCommonDen( g ); g *=cd; //g /=vcontent(g,Variable(1)); CanonicalForm Dn, test= 0; cl = gcd (f.lc(),g.lc()); CanonicalForm gcdcfcg= gcd (cf, cg); CanonicalForm fp, gp; CanonicalForm b= 1; int minCommonDeg= 0; for (i= tmax (f.level(), g.level()); i > 0; i--) { if (degree (f, i) <= 0 || degree (g, i) <= 0) continue; else { minCommonDeg= tmin (degree (g, i), degree (f, i)); break; } } if (i == 0) return gcdcfcg; for (; i > 0; i--) { if (degree (f, i) <= 0 || degree (g, i) <= 0) continue; else minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i))); } b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)* power (CanonicalForm (2), minCommonDeg); bool equal= false; i = cf_getNumBigPrimes() - 1; CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn; int maxNumVars= tmax (getNumVars (f), getNumVars (g)); //Off (SW_RATIONAL); while ( true ) { p = cf_getBigPrime( i ); i--; while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 ) { p = cf_getBigPrime( i ); i--; } //printf("try p=%d\n",p); setCharacteristic( p ); fp= mapinto (f); gp= mapinto (g); TIMING_START (chinrem_recursion) #ifdef HAVE_NTL if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500) Dp = GCD_small_p (fp, gp, cofp, cogp); else { Dp= gcd_poly (fp, gp); cofp= fp/Dp; cogp= gp/Dp; } #else Dp= gcd_poly (fp, gp); cofp= fp/Dp; cogp= gp/Dp; #endif TIMING_END_AND_PRINT (chinrem_recursion, "time for gcd mod p in modular gcd: "); Dp /=Dp.lc(); cofp /= lc (cofp); cogp /= lc (cogp); setCharacteristic( 0 ); dp_deg=totaldegree(Dp); if ( dp_deg == 0 ) { //printf(" -> 1\n"); return CanonicalForm(gcdcfcg); } if ( q.isZero() ) { D = mapinto( Dp ); cof= mapinto (cofp); cog= mapinto (cogp); d_deg=dp_deg; q = p; } else { if ( dp_deg == d_deg ) { chineseRemainder( D, q, mapinto( Dp ), p, newD, newq ); chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq); chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq); cof= newCof; cog= newCog; q = newq; D = newD; } else if ( dp_deg < d_deg ) { //printf(" were all bad, try more\n"); // all previous p's are bad primes q = p; D = mapinto( Dp ); cof= mapinto (cof); cog= mapinto (cog); d_deg=dp_deg; test= 0; equal= false; } else { //printf(" was bad, try more\n"); } //else dp_deg > d_deg: bad prime } if ( i >= 0 ) { TIMING_START (chinrem_reconstruction); Dn= Farey(D,q); cofn= Farey(cof,q); cogn= Farey(cog,q); TIMING_END_AND_PRINT (chinrem_reconstruction, "time for rational reconstruction in modular gcd: "); int is_rat= isOn (SW_RATIONAL); On (SW_RATIONAL); cd = bCommonDen( Dn ); // we need On(SW_RATIONAL) cofn *= bCommonDen (cofn); cogn *= bCommonDen (cogn); if (!is_rat) Off (SW_RATIONAL); Dn *=cd; if (test != Dn) test= Dn; else equal= true; //Dn /=vcontent(Dn,Variable(1)); TIMING_START (chinrem_termination); if ((terminationTest (f,g, cofn, cogn, Dn)) || ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g))) { TIMING_END_AND_PRINT (chinrem_termination, "time for successful termination in modular gcd: "); //printf(" -> success\n"); return Dn*gcdcfcg; } TIMING_END_AND_PRINT (chinrem_termination, "time for unsuccessful termination in modular gcd: "); equal= false; //else: try more primes } else { // try other method //printf("try other gcd\n"); Off(SW_USE_CHINREM_GCD); D=gcd_poly( f, g ); On(SW_USE_CHINREM_GCD); return D*gcdcfcg; } } }