source: git/factory/cf_gcd.cc @ a209e1d

spielwiese
Last change on this file since a209e1d was 5b2d2b, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: gcd_test_one now checks for coprimality over extension chg: gcd_test_one returns degree of gcd used to check for coprimality
  • Property mode set to 100644
File size: 34.5 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2
3#include "config.h"
4
5#include "timing.h"
6#include "cf_assert.h"
7#include "debug.h"
8
9#include "cf_defs.h"
10#include "canonicalform.h"
11#include "cf_iter.h"
12#include "cf_reval.h"
13#include "cf_primes.h"
14#include "cf_algorithm.h"
15#include "cf_factory.h"
16#include "fac_util.h"
17#include "templates/ftmpl_functions.h"
18#include "algext.h"
19#include "cf_gcd_smallp.h"
20#include "cf_map_ext.h"
21#include "cf_util.h"
22#include "gfops.h"
23
24#ifdef HAVE_NTL
25#include <NTL/ZZX.h>
26#include "NTLconvert.h"
27bool isPurePoly(const CanonicalForm & );
28#ifndef HAVE_FLINT
29static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
30static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
31#endif
32#endif
33
34#ifdef HAVE_FLINT
35#include "FLINTconvert.h"
36static CanonicalForm gcd_univar_flint0 (const CanonicalForm &, const CanonicalForm &);
37static CanonicalForm gcd_univar_flintp (const CanonicalForm &, const CanonicalForm &);
38#endif
39
40static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
41
42void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
43
44CanonicalForm chinrem_gcd(const CanonicalForm & FF,const CanonicalForm & GG);
45
46bool
47gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap, int & d )
48{
49    d= 0;
50    int count = 0;
51    // assume polys have same level;
52
53    Variable v= Variable (1);
54    bool algExtension= (hasFirstAlgVar (f, v) || hasFirstAlgVar (g, v));
55    CanonicalForm lcf, lcg;
56    if ( swap )
57    {
58        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
59        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
60    }
61    else
62    {
63        lcf = LC( f, Variable(1) );
64        lcg = LC( g, Variable(1) );
65    }
66
67    CanonicalForm F, G;
68    if ( swap )
69    {
70        F=swapvar( f, Variable(1), f.mvar() );
71        G=swapvar( g, Variable(1), g.mvar() );
72    }
73    else
74    {
75        F = f;
76        G = g;
77    }
78
79    #define TEST_ONE_MAX 50
80    int p= getCharacteristic();
81    bool passToGF= false;
82    int k= 1;
83    if (p > 0 && p < TEST_ONE_MAX && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
84    {
85      if (p == 2)
86        setCharacteristic (2, 6, 'Z');
87      else if (p == 3)
88        setCharacteristic (3, 4, 'Z');
89      else if (p == 5 || p == 7)
90        setCharacteristic (p, 3, 'Z');
91      else
92        setCharacteristic (p, 2, 'Z');
93      passToGF= true;
94    }
95    else if (p > 0 && CFFactory::gettype() == GaloisFieldDomain && ipower (p , getGFDegree()) < TEST_ONE_MAX)
96    {
97      k= getGFDegree();
98      if (ipower (p, 2*k) > TEST_ONE_MAX)
99        setCharacteristic (p, 2*k, gf_name);
100      else
101        setCharacteristic (p, 3*k, gf_name);
102      F= GFMapUp (F, k);
103      G= GFMapUp (G, k);
104      lcf= GFMapUp (lcf, k);
105      lcg= GFMapUp (lcg, k);
106    }
107    else if (p > 0 && p < TEST_ONE_MAX && algExtension)
108    {
109      bool extOfExt= false;
110#ifdef HAVE_NTL
111      int d= degree (getMipo (v));
112      CFList source, dest;
113      Variable v2;
114      CanonicalForm primElem, imPrimElem;
115      if (p == 2 && d < 6)
116      {
117        zz_p::init (p);
118        bool primFail= false;
119        Variable vBuf;
120        primElem= primitiveElement (v, vBuf, primFail);
121        ASSERT (!primFail, "failure in integer factorizer");
122        if (d < 3)
123        {
124          zz_pX NTLIrredpoly;
125          BuildIrred (NTLIrredpoly, d*3);
126          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
127          v2= rootOf (newMipo);
128        }
129        else
130        {
131          zz_pX NTLIrredpoly;
132          BuildIrred (NTLIrredpoly, d*2);
133          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
134          v2= rootOf (newMipo);
135        }
136        imPrimElem= mapPrimElem (primElem, v, v2);
137        extOfExt= true;
138      }
139      else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
140      {
141        zz_p::init (p);
142        bool primFail= false;
143        Variable vBuf;
144        primElem= primitiveElement (v, vBuf, primFail);
145        ASSERT (!primFail, "failure in integer factorizer");
146        zz_pX NTLIrredpoly;
147        BuildIrred (NTLIrredpoly, d*2);
148        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
149        v2= rootOf (newMipo);
150        imPrimElem= mapPrimElem (primElem, v, v2);
151        extOfExt= true;
152      }
153      if (extOfExt)
154      {
155        F= mapUp (F, v, v2, primElem, imPrimElem, source, dest);
156        G= mapUp (G, v, v2, primElem, imPrimElem, source, dest);
157        lcf= mapUp (lcf, v, v2, primElem, imPrimElem, source, dest);
158        lcg= mapUp (lcg, v, v2, primElem, imPrimElem, source, dest);
159        v= v2;
160      }
161#endif
162    }
163
164    CFRandom * sample;
165    if ((!algExtension && p > 0) || p == 0)
166      sample = CFRandomFactory::generate();
167    else
168      sample = AlgExtRandomF (v).clone();
169
170    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
171    delete sample;
172
173    if (passToGF)
174    {
175      lcf= lcf.mapinto();
176      lcg= lcg.mapinto();
177    }
178
179    CanonicalForm eval1, eval2;
180    if (passToGF)
181    {
182      eval1= e (lcf);
183      eval2= e (lcg);
184    }
185    else
186    {
187      eval1= e (lcf);
188      eval2= e (lcg);
189    }
190
191    while ( ( eval1.isZero() || eval2.isZero() ) && count < TEST_ONE_MAX )
192    {
193        e.nextpoint();
194        count++;
195        eval1= e (lcf);
196        eval2= e (lcg);
197    }
198    if ( count >= TEST_ONE_MAX )
199    {
200        if (passToGF)
201          setCharacteristic (p);
202        if (k > 1)
203          setCharacteristic (p, k, gf_name);
204        return false;
205    }
206
207
208    if (passToGF)
209    {
210      F= F.mapinto();
211      G= G.mapinto();
212      eval1= e (F);
213      eval2= e (G);
214    }
215    else
216    {
217      eval1= e (F);
218      eval2= e (G);
219    }
220
221    CanonicalForm c= gcd (eval1, eval2);
222    d= c.degree();
223    bool result= d < 1;
224    if (d < 0)
225      d= 0;
226
227    if (passToGF)
228      setCharacteristic (p);
229    if (k > 1)
230      setCharacteristic (p, k, gf_name);
231    return result;
232}
233
234//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
235//{{{ docu
236//
237// icontent() - return gcd of c and all coefficients of f which
238//   are in a coefficient domain.
239//
240// Used by icontent().
241//
242//}}}
243static CanonicalForm
244icontent ( const CanonicalForm & f, const CanonicalForm & c )
245{
246    if ( f.inBaseDomain() )
247    {
248      if (c.isZero()) return abs(f);
249      return bgcd( f, c );
250    }
251    //else if ( f.inCoeffDomain() )
252    //   return gcd(f,c);
253    else
254    {
255        CanonicalForm g = c;
256        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
257            g = icontent( i.coeff(), g );
258        return g;
259    }
260}
261//}}}
262
263//{{{ CanonicalForm icontent ( const CanonicalForm & f )
264//{{{ docu
265//
266// icontent() - return gcd over all coefficients of f which are
267//   in a coefficient domain.
268//
269//}}}
270CanonicalForm
271icontent ( const CanonicalForm & f )
272{
273    return icontent( f, 0 );
274}
275//}}}
276
277//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
278//{{{ docu
279//
280// extgcd() - returns polynomial extended gcd of f and g.
281//
282// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
283// The gcd is calculated using an extended euclidean polynomial
284// remainder sequence, so f and g should be polynomials over an
285// euclidean domain.  Normalizes result.
286//
287// Note: be sure that f and g have the same level!
288//
289//}}}
290CanonicalForm
291extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
292{
293  if (f.isZero())
294  {
295    a= 0;
296    b= 1;
297    return g;
298  }
299  else if (g.isZero())
300  {
301    a= 1;
302    b= 0;
303    return f;
304  }
305#ifdef HAVE_NTL
306#ifdef HAVE_FLINT
307  if (( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain)
308  &&  (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
309  {
310    nmod_poly_t F1, G1, A, B, R;
311    convertFacCF2nmod_poly_t (F1, f);
312    convertFacCF2nmod_poly_t (G1, g);
313    nmod_poly_init (R, getCharacteristic());
314    nmod_poly_init (A, getCharacteristic());
315    nmod_poly_init (B, getCharacteristic());
316    nmod_poly_xgcd (R, A, B, F1, G1);
317    a= convertnmod_poly_t2FacCF (A, f.mvar());
318    b= convertnmod_poly_t2FacCF (B, f.mvar());
319    CanonicalForm r= convertnmod_poly_t2FacCF (R, f.mvar());
320    nmod_poly_clear (F1);
321    nmod_poly_clear (G1);
322    nmod_poly_clear (A);
323    nmod_poly_clear (B);
324    nmod_poly_clear (R);
325    return r;
326  }
327#else
328  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain)
329  &&  (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
330  {
331    if (fac_NTL_char!=getCharacteristic())
332    {
333      fac_NTL_char=getCharacteristic();
334      #ifdef NTL_ZZ
335      ZZ r;
336      r=getCharacteristic();
337      ZZ_pContext ccc(r);
338      #else
339      zz_pContext ccc(getCharacteristic());
340      #endif
341      ccc.restore();
342      #ifdef NTL_ZZ
343      ZZ_p::init(r);
344      #else
345      zz_p::init(getCharacteristic());
346      #endif
347    }
348    #ifdef NTL_ZZ
349    ZZ_pX F1=convertFacCF2NTLZZpX(f);
350    ZZ_pX G1=convertFacCF2NTLZZpX(g);
351    ZZ_pX R;
352    ZZ_pX A,B;
353    XGCD(R,A,B,F1,G1);
354    a=convertNTLZZpX2CF(A,f.mvar());
355    b=convertNTLZZpX2CF(B,f.mvar());
356    return convertNTLZZpX2CF(R,f.mvar());
357    #else
358    zz_pX F1=convertFacCF2NTLzzpX(f);
359    zz_pX G1=convertFacCF2NTLzzpX(g);
360    zz_pX R;
361    zz_pX A,B;
362    XGCD(R,A,B,F1,G1);
363    a=convertNTLzzpX2CF(A,f.mvar());
364    b=convertNTLzzpX2CF(B,f.mvar());
365    return convertNTLzzpX2CF(R,f.mvar());
366    #endif
367  }
368#endif
369#ifdef HAVE_FLINT
370  if (( getCharacteristic() ==0) && (f.level()==g.level())
371       && isPurePoly(f) && isPurePoly(g))
372  {
373    fmpq_poly_t F1, G1;
374    convertFacCF2Fmpq_poly_t (F1, f);
375    convertFacCF2Fmpq_poly_t (G1, g);
376    fmpq_poly_t R, A, B;
377    fmpq_poly_init (R);
378    fmpq_poly_init (A);
379    fmpq_poly_init (B);
380    fmpq_poly_xgcd (R, A, B, F1, G1);
381    a= convertFmpq_poly_t2FacCF (A, f.mvar());
382    b= convertFmpq_poly_t2FacCF (B, f.mvar());
383    CanonicalForm r= convertFmpq_poly_t2FacCF (R, f.mvar());
384    fmpq_poly_clear (F1);
385    fmpq_poly_clear (G1);
386    fmpq_poly_clear (A);
387    fmpq_poly_clear (B);
388    fmpq_poly_clear (R);
389    return r;
390  }
391#else
392  if (isOn(SW_USE_NTL_GCD_0) && ( getCharacteristic() ==0)
393  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
394  {
395    CanonicalForm fc=bCommonDen(f);
396    CanonicalForm gc=bCommonDen(g);
397    ZZX F1=convertFacCF2NTLZZX(f*fc);
398    ZZX G1=convertFacCF2NTLZZX(g*gc);
399    ZZX R=GCD(F1,G1);
400    CanonicalForm r=convertNTLZZX2CF(R,f.mvar());
401    ZZ RR;
402    ZZX A,B;
403    if (r.inCoeffDomain())
404    {
405      XGCD(RR,A,B,F1,G1,1);
406      CanonicalForm rr=convertZZ2CF(RR);
407      ASSERT (!rr.isZero(), "NTL:XGCD failed");
408      a=convertNTLZZX2CF(A,f.mvar())*fc/rr;
409      b=convertNTLZZX2CF(B,f.mvar())*gc/rr;
410      return CanonicalForm(1);
411    }
412    else
413    {
414      fc=bCommonDen(f);
415      gc=bCommonDen(g);
416      F1=convertFacCF2NTLZZX(f*fc/r);
417      G1=convertFacCF2NTLZZX(g*gc/r);
418      XGCD(RR,A,B,F1,G1,1);
419      a=convertNTLZZX2CF(A,f.mvar())*fc;
420      b=convertNTLZZX2CF(B,f.mvar())*gc;
421      CanonicalForm rr=convertZZ2CF(RR);
422      ASSERT (!rr.isZero(), "NTL:XGCD failed");
423      a /= rr;
424      b /= rr;
425      return r;
426    }
427  }
428#endif
429#endif
430  // may contain bug in the co-factors, see track 107
431  CanonicalForm contf = content( f );
432  CanonicalForm contg = content( g );
433
434  CanonicalForm p0 = f / contf, p1 = g / contg;
435  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
436
437  while ( ! p1.isZero() )
438  {
439      divrem( p0, p1, q, r );
440      p0 = p1; p1 = r;
441      r = g0 - g1 * q;
442      g0 = g1; g1 = r;
443      r = f0 - f1 * q;
444      f0 = f1; f1 = r;
445  }
446  CanonicalForm contp0 = content( p0 );
447  a = f0 / ( contf * contp0 );
448  b = g0 / ( contg * contp0 );
449  p0 /= contp0;
450  if ( p0.sign() < 0 )
451  {
452      p0 = -p0;
453      a = -a;
454      b = -b;
455  }
456  return p0;
457}
458//}}}
459
460//{{{ static CanonicalForm balance_p ( const CanonicalForm & f, const CanonicalForm & q )
461//{{{ docu
462//
463// balance_p() - map f from positive to symmetric representation
464//   mod q.
465//
466// This makes sense for univariate polynomials over Z only.
467// q should be an integer.
468//
469// Used by gcd_poly_univar0().
470//
471//}}}
472static CanonicalForm
473balance_p ( const CanonicalForm & f, const CanonicalForm & q )
474{
475    Variable x = f.mvar();
476    CanonicalForm result = 0, qh = q / 2;
477    CanonicalForm c;
478    CFIterator i;
479    for ( i = f; i.hasTerms(); i++ )
480    {
481        c = i.coeff();
482        if ( c.inCoeffDomain())
483        {
484          if ( c > qh )
485            result += power( x, i.exp() ) * (c - q);
486          else
487            result += power( x, i.exp() ) * c;
488        }
489        else
490          result += power( x, i.exp() ) * balance_p(c,q);
491    }
492    return result;
493}
494
495/*static CanonicalForm
496balance ( const CanonicalForm & f, const CanonicalForm & q )
497{
498    Variable x = f.mvar();
499    CanonicalForm result = 0, qh = q / 2;
500    CanonicalForm c;
501    CFIterator i;
502    for ( i = f; i.hasTerms(); i++ ) {
503        c = mod( i.coeff(), q );
504        if ( c > qh )
505            result += power( x, i.exp() ) * (c - q);
506        else
507            result += power( x, i.exp() ) * c;
508    }
509    return result;
510}*/
511//}}}
512
513static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
514{
515  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
516  int p, i;
517
518  if ( primitive )
519  {
520    f = F;
521    g = G;
522    c = 1;
523  }
524  else
525  {
526    CanonicalForm cF = content( F ), cG = content( G );
527    f = F / cF;
528    g = G / cG;
529    c = bgcd( cF, cG );
530  }
531  cg = gcd( f.lc(), g.lc() );
532  cl = ( f.lc() / cg ) * g.lc();
533//     B = 2 * cg * tmin(
534//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
535//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
536//         )+1;
537  M = tmin( maxNorm(f), maxNorm(g) );
538  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
539  q = 0;
540  i = cf_getNumSmallPrimes() - 1;
541  while ( true )
542  {
543    B = BB;
544    while ( i >= 0 && q < B )
545    {
546      p = cf_getSmallPrime( i );
547      i--;
548      while ( i >= 0 && mod( cl, p ) == 0 )
549      {
550        p = cf_getSmallPrime( i );
551        i--;
552      }
553      setCharacteristic( p );
554      Dp = gcd( mapinto( f ), mapinto( g ) );
555      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
556      setCharacteristic( 0 );
557      if ( Dp.degree() == 0 )
558        return c;
559      if ( q.isZero() )
560      {
561        D = mapinto( Dp );
562        q = p;
563        B = power(CanonicalForm(2),D.degree())*M+1;
564      }
565      else
566      {
567        if ( Dp.degree() == D.degree() )
568        {
569          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
570          q = newq;
571          D = newD;
572        }
573        else if ( Dp.degree() < D.degree() )
574        {
575          // all previous p's are bad primes
576          q = p;
577          D = mapinto( Dp );
578          B = power(CanonicalForm(2),D.degree())*M+1;
579        }
580        // else p is a bad prime
581      }
582    }
583    if ( i >= 0 )
584    {
585      // now balance D mod q
586      D = pp( balance_p( D, q ) );
587      if ( fdivides( D, f ) && fdivides( D, g ) )
588        return D * c;
589      else
590        q = 0;
591    }
592    else
593      return gcd_poly( F, G );
594    DEBOUTLN( cerr, "another try ..." );
595  }
596}
597
598static CanonicalForm
599gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
600{
601    CanonicalForm pi, pi1;
602    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
603    bool bpure;
604    int delta = degree( f ) - degree( g );
605
606    if ( delta >= 0 )
607    {
608        pi = f; pi1 = g;
609    }
610    else
611    {
612        pi = g; pi1 = f; delta = -delta;
613    }
614    Ci = content( pi ); Ci1 = content( pi1 );
615    pi1 = pi1 / Ci1; pi = pi / Ci;
616    C = gcd( Ci, Ci1 );
617    int d= 0;
618    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
619    {
620        if ( gcd_test_one( pi1, pi, true, d ) )
621        {
622          C=abs(C);
623          //out_cf("GCD:",C,"\n");
624          return C;
625        }
626        bpure = false;
627    }
628    else
629    {
630        bpure = isPurePoly(pi) && isPurePoly(pi1);
631#ifdef HAVE_FLINT
632        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
633          return gcd_univar_flintp(pi,pi1)*C;
634#else
635#ifdef HAVE_NTL
636        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
637            return gcd_univar_ntlp(pi, pi1 ) * C;
638#endif
639#endif
640    }
641    Variable v = f.mvar();
642    Hi = power( LC( pi1, v ), delta );
643    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
644
645    if (!(pi.isUnivariate() && pi1.isUnivariate()))
646    {
647      if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
648      {
649        On (SW_USE_FF_MOD_GCD);
650        C *= gcd (pi, pi1);
651        Off (SW_USE_FF_MOD_GCD);
652        return C;
653      }
654    }
655
656    if ( (delta+1) % 2 )
657        bi = 1;
658    else
659        bi = -1;
660    CanonicalForm oldPi= pi, oldPi1= pi1, powHi;
661    while ( degree( pi1, v ) > 0 )
662    {
663        if (!(pi.isUnivariate() && pi1.isUnivariate()))
664        {
665          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
666          {
667            On (SW_USE_FF_MOD_GCD);
668            C *= gcd (oldPi, oldPi1);
669            Off (SW_USE_FF_MOD_GCD);
670            return C;
671          }
672        }
673        pi2 = psr( pi, pi1, v );
674        pi2 = pi2 / bi;
675        pi = pi1; pi1 = pi2;
676        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
677        if (!(pi1.isUnivariate()) && (size (pi1)/maxNumVars > 500))
678        {
679            On (SW_USE_FF_MOD_GCD);
680            C *= gcd (oldPi, oldPi1);
681            Off (SW_USE_FF_MOD_GCD);
682            return C;
683        }
684        if ( degree( pi1, v ) > 0 )
685        {
686            delta = degree( pi, v ) - degree( pi1, v );
687            powHi= power (Hi, delta-1);
688            if ( (delta+1) % 2 )
689                bi = LC( pi, v ) * powHi*Hi;
690            else
691                bi = -LC( pi, v ) * powHi*Hi;
692            Hi = power( LC( pi1, v ), delta ) / powHi;
693            if (!(pi.isUnivariate() && pi1.isUnivariate()))
694            {
695              if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
696              {
697                On (SW_USE_FF_MOD_GCD);
698                C *= gcd (oldPi, oldPi1);
699                Off (SW_USE_FF_MOD_GCD);
700                return C;
701              }
702            }
703        }
704    }
705    if ( degree( pi1, v ) == 0 )
706    {
707      C=abs(C);
708      //out_cf("GCD:",C,"\n");
709      return C;
710    }
711    pi /= content( pi );
712    if ( bpure )
713        pi /= pi.lc();
714    C=abs(C*pi);
715    //out_cf("GCD:",C,"\n");
716    return C;
717}
718
719static CanonicalForm
720gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
721{
722    CanonicalForm pi, pi1;
723    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
724    int delta = degree( f ) - degree( g );
725
726    if ( delta >= 0 )
727    {
728        pi = f; pi1 = g;
729    }
730    else
731    {
732        pi = g; pi1 = f; delta = -delta;
733    }
734    Ci = content( pi ); Ci1 = content( pi1 );
735    pi1 = pi1 / Ci1; pi = pi / Ci;
736    C = gcd( Ci, Ci1 );
737    int d= 0;
738    if ( pi.isUnivariate() && pi1.isUnivariate() )
739    {
740#ifdef HAVE_FLINT
741        if (isPurePoly(pi) && isPurePoly(pi1) )
742            return gcd_univar_flint0(pi, pi1 ) * C;
743#else
744#ifdef HAVE_NTL
745        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
746            return gcd_univar_ntl0(pi, pi1 ) * C;
747#endif
748#endif
749        return gcd_poly_univar0( pi, pi1, true ) * C;
750    }
751    else if ( gcd_test_one( pi1, pi, true, d ) )
752      return C;
753    Variable v = f.mvar();
754    Hi = power( LC( pi1, v ), delta );
755    if ( (delta+1) % 2 )
756        bi = 1;
757    else
758        bi = -1;
759    while ( degree( pi1, v ) > 0 )
760    {
761        pi2 = psr( pi, pi1, v );
762        pi2 = pi2 / bi;
763        pi = pi1; pi1 = pi2;
764        if ( degree( pi1, v ) > 0 )
765        {
766            delta = degree( pi, v ) - degree( pi1, v );
767            if ( (delta+1) % 2 )
768                bi = LC( pi, v ) * power( Hi, delta );
769            else
770                bi = -LC( pi, v ) * power( Hi, delta );
771            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
772        }
773    }
774    if ( degree( pi1, v ) == 0 )
775        return C;
776    else
777        return C * pp( pi );
778}
779
780//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
781//{{{ docu
782//
783// gcd_poly() - calculate polynomial gcd.
784//
785// This is the dispatcher for polynomial gcd calculation.  We call either
786// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
787// characteristic and settings of SW_USE_EZGCD.
788//
789// Used by gcd() and gcd_poly_univar0().
790//
791//}}}
792#if 0
793int si_factor_reminder=1;
794#endif
795CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
796{
797  CanonicalForm fc, gc, d1;
798  bool fc_isUnivariate=f.isUnivariate();
799  bool gc_isUnivariate=g.isUnivariate();
800  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
801  fc = f;
802  gc = g;
803  if ( getCharacteristic() != 0 )
804  {
805    #ifdef HAVE_NTL
806    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
807    {
808      fc= EZGCD_P (fc, gc);
809    }
810    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
811    {
812      Variable a;
813      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
814        fc=GCD_Fp_extension (fc, gc, a);
815      else if (CFFactory::gettype() == GaloisFieldDomain)
816        fc=GCD_GF (fc, gc);
817      else
818        fc=GCD_small_p (fc, gc);
819    }
820    else
821    #endif
822    fc = gcd_poly_p( fc, gc );
823  }
824  else if (!fc_and_gc_Univariate)
825  {
826    if ( isOn( SW_USE_EZGCD ) )
827      fc= ezgcd (fc, gc);
828    else if (isOn(SW_USE_CHINREM_GCD))
829      fc = chinrem_gcd( fc, gc);
830    else
831    {
832       fc = gcd_poly_0( fc, gc );
833    }
834  }
835  else
836  {
837    fc = gcd_poly_0( fc, gc );
838  }
839  if ( d1.degree() > 0 )
840    fc *= d1;
841  return fc;
842}
843//}}}
844
845//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
846//{{{ docu
847//
848// cf_content() - return gcd(g, content(f)).
849//
850// content(f) is calculated with respect to f's main variable.
851//
852// Used by gcd(), content(), content( CF, Variable ).
853//
854//}}}
855static CanonicalForm
856cf_content ( const CanonicalForm & f, const CanonicalForm & g )
857{
858    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
859    {
860        CFIterator i = f;
861        CanonicalForm result = g;
862        while ( i.hasTerms() && ! result.isOne() )
863        {
864            result = gcd( i.coeff(), result );
865            i++;
866        }
867        return result;
868    }
869    else
870        return abs( f );
871}
872//}}}
873
874//{{{ CanonicalForm content ( const CanonicalForm & f )
875//{{{ docu
876//
877// content() - return content(f) with respect to main variable.
878//
879// Normalizes result.
880//
881//}}}
882CanonicalForm
883content ( const CanonicalForm & f )
884{
885    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
886    {
887        CFIterator i = f;
888        CanonicalForm result = abs( i.coeff() );
889        i++;
890        while ( i.hasTerms() && ! result.isOne() )
891        {
892            result = gcd( i.coeff(), result );
893            i++;
894        }
895        return result;
896    }
897    else
898        return abs( f );
899}
900//}}}
901
902//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
903//{{{ docu
904//
905// content() - return content(f) with respect to x.
906//
907// x should be a polynomial variable.
908//
909//}}}
910CanonicalForm
911content ( const CanonicalForm & f, const Variable & x )
912{
913    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
914    Variable y = f.mvar();
915
916    if ( y == x )
917        return cf_content( f, 0 );
918    else  if ( y < x )
919        return f;
920    else
921        return swapvar( content( swapvar( f, y, x ), y ), y, x );
922}
923//}}}
924
925//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
926//{{{ docu
927//
928// vcontent() - return content of f with repect to variables >= x.
929//
930// The content is recursively calculated over all coefficients in
931// f having level less than x.  x should be a polynomial
932// variable.
933//
934//}}}
935CanonicalForm
936vcontent ( const CanonicalForm & f, const Variable & x )
937{
938    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
939
940    if ( f.mvar() <= x )
941        return content( f, x );
942    else {
943        CFIterator i;
944        CanonicalForm d = 0;
945        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
946            d = gcd( d, vcontent( i.coeff(), x ) );
947        return d;
948    }
949}
950//}}}
951
952//{{{ CanonicalForm pp ( const CanonicalForm & f )
953//{{{ docu
954//
955// pp() - return primitive part of f.
956//
957// Returns zero if f equals zero, otherwise f / content(f).
958//
959//}}}
960CanonicalForm
961pp ( const CanonicalForm & f )
962{
963    if ( f.isZero() )
964        return f;
965    else
966        return f / content( f );
967}
968//}}}
969
970CanonicalForm
971gcd ( const CanonicalForm & f, const CanonicalForm & g )
972{
973    bool b = f.isZero();
974    if ( b || g.isZero() )
975    {
976        if ( b )
977            return abs( g );
978        else
979            return abs( f );
980    }
981    if ( f.inPolyDomain() || g.inPolyDomain() )
982    {
983        if ( f.mvar() != g.mvar() )
984        {
985            if ( f.mvar() > g.mvar() )
986                return cf_content( f, g );
987            else
988                return cf_content( g, f );
989        }
990        if (isOn(SW_USE_QGCD))
991        {
992          Variable m;
993          if (
994          (getCharacteristic() == 0) &&
995          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
996          )
997          {
998            bool on_rational = isOn(SW_RATIONAL);
999            CanonicalForm r=QGCD(f,g);
1000            On(SW_RATIONAL);
1001            CanonicalForm cdF = bCommonDen( r );
1002            if (!on_rational) Off(SW_RATIONAL);
1003            return cdF*r;
1004          }
1005        }
1006
1007        if ( f.inExtension() && getReduce( f.mvar() ) )
1008            return CanonicalForm(1);
1009        else
1010        {
1011            if ( fdivides( f, g ) )
1012                return abs( f );
1013            else  if ( fdivides( g, f ) )
1014                return abs( g );
1015            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
1016            {
1017                CanonicalForm d;
1018                d = gcd_poly( f, g );
1019                return abs( d );
1020            }
1021            else
1022            {
1023                //printf ("here\n");
1024                CanonicalForm cdF = bCommonDen( f );
1025                CanonicalForm cdG = bCommonDen( g );
1026                Off( SW_RATIONAL );
1027                CanonicalForm l = lcm( cdF, cdG );
1028                On( SW_RATIONAL );
1029                CanonicalForm F = f * l, G = g * l;
1030                Off( SW_RATIONAL );
1031                l = gcd_poly( F, G );
1032                On( SW_RATIONAL );
1033                return abs( l );
1034            }
1035        }
1036    }
1037    if ( f.inBaseDomain() && g.inBaseDomain() )
1038        return bgcd( f, g );
1039    else
1040        return 1;
1041}
1042
1043//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1044//{{{ docu
1045//
1046// lcm() - return least common multiple of f and g.
1047//
1048// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1049//
1050// Returns zero if one of f or g equals zero.
1051//
1052//}}}
1053CanonicalForm
1054lcm ( const CanonicalForm & f, const CanonicalForm & g )
1055{
1056    if ( f.isZero() || g.isZero() )
1057        return 0;
1058    else
1059        return ( f / gcd( f, g ) ) * g;
1060}
1061//}}}
1062
1063#ifdef HAVE_NTL
1064#ifndef HAVE_FLINT
1065static CanonicalForm
1066gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1067{
1068    ZZX F1=convertFacCF2NTLZZX(F);
1069    ZZX G1=convertFacCF2NTLZZX(G);
1070    ZZX R=GCD(F1,G1);
1071    return convertNTLZZX2CF(R,F.mvar());
1072}
1073
1074static CanonicalForm
1075gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1076{
1077  if (fac_NTL_char!=getCharacteristic())
1078  {
1079    fac_NTL_char=getCharacteristic();
1080    #ifdef NTL_ZZ
1081    ZZ r;
1082    r=getCharacteristic();
1083    ZZ_pContext ccc(r);
1084    #else
1085    zz_pContext ccc(getCharacteristic());
1086    #endif
1087    ccc.restore();
1088    #ifdef NTL_ZZ
1089    ZZ_p::init(r);
1090    #else
1091    zz_p::init(getCharacteristic());
1092    #endif
1093  }
1094  #ifdef NTL_ZZ
1095  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1096  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1097  ZZ_pX R=GCD(F1,G1);
1098  return  convertNTLZZpX2CF(R,F.mvar());
1099  #else
1100  zz_pX F1=convertFacCF2NTLzzpX(F);
1101  zz_pX G1=convertFacCF2NTLzzpX(G);
1102  zz_pX R=GCD(F1,G1);
1103  return  convertNTLzzpX2CF(R,F.mvar());
1104  #endif
1105}
1106#endif
1107#endif
1108
1109#ifdef HAVE_FLINT
1110static CanonicalForm
1111gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
1112{
1113  nmod_poly_t F1, G1;
1114  convertFacCF2nmod_poly_t (F1, F);
1115  convertFacCF2nmod_poly_t (G1, G);
1116  nmod_poly_gcd (F1, F1, G1);
1117  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
1118  nmod_poly_clear (F1);
1119  nmod_poly_clear (G1);
1120  return result;
1121}
1122
1123static CanonicalForm
1124gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
1125{
1126  fmpz_poly_t F1, G1;
1127  convertFacCF2Fmpz_poly_t(F1, F);
1128  convertFacCF2Fmpz_poly_t(G1, G);
1129  fmpz_poly_gcd (F1, F1, G1);
1130  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
1131  fmpz_poly_clear (F1);
1132  fmpz_poly_clear (G1);
1133  return result;
1134}
1135#endif
1136
1137
1138/*
1139*  compute positions p1 and pe of optimal variables:
1140*    pe is used in "ezgcd" and
1141*    p1 in "gcd_poly1"
1142*/
1143/*static
1144void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1145{
1146    int i, o1, oe;
1147    if ( df[n] > dg[n] )
1148    {
1149        o1 = df[n]; oe = dg[n];
1150    }
1151    else
1152    {
1153        o1 = dg[n]; oe = df[n];
1154    }
1155    i = n-1;
1156    while ( i > 0 )
1157    {
1158        if ( df[i] != 0 )
1159        {
1160            if ( df[i] > dg[i] )
1161            {
1162                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1163                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1164            }
1165            else
1166            {
1167                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1168                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1169            }
1170        }
1171        i--;
1172    }
1173}*/
1174
1175/*
1176*  make some changes of variables, see optvalues
1177*/
1178/*static void
1179cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1180{
1181    int i, k, n;
1182    n = f.level();
1183    cc = 0;
1184    p1 = pe = n;
1185    if ( n == 1 )
1186        return;
1187    int * degsf = new int[n+1];
1188    int * degsg = new int[n+1];
1189    for ( i = n; i > 0; i-- )
1190    {
1191        degsf[i] = degsg[i] = 0;
1192    }
1193    degsf = degrees( f, degsf );
1194    degsg = degrees( g, degsg );
1195
1196    k = 0;
1197    for ( i = n-1; i > 0; i-- )
1198    {
1199        if ( degsf[i] == 0 )
1200        {
1201            if ( degsg[i] != 0 )
1202            {
1203                cc = -i;
1204                break;
1205            }
1206        }
1207        else
1208        {
1209            if ( degsg[i] == 0 )
1210            {
1211                cc = i;
1212                break;
1213            }
1214            else k++;
1215        }
1216    }
1217
1218    if ( ( cc == 0 ) && ( k != 0 ) )
1219        optvalues( degsf, degsg, n, p1, pe );
1220    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1221        pe = -pe;
1222
1223    delete [] degsf;
1224    delete [] degsg;
1225}*/
1226
1227TIMING_DEFINE_PRINT(chinrem_termination)
1228TIMING_DEFINE_PRINT(chinrem_recursion)
1229TIMING_DEFINE_PRINT(chinrem_reconstruction)
1230
1231CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1232{
1233  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq;
1234  int p, i, dp_deg, d_deg=-1;
1235
1236  CanonicalForm cd ( bCommonDen( FF ));
1237  f=cd*FF;
1238  Variable x= Variable (1);
1239  CanonicalForm cf, cg;
1240  cf= icontent (f);
1241  f /= cf;
1242  //cd = bCommonDen( f ); f *=cd;
1243  //f /=vcontent(f,Variable(1));
1244
1245  cd = bCommonDen( GG );
1246  g=cd*GG;
1247  cg= icontent (g);
1248  g /= cg;
1249  //cd = bCommonDen( g ); g *=cd;
1250  //g /=vcontent(g,Variable(1));
1251
1252  CanonicalForm Dn, test= 0;
1253  cl =  gcd (f.lc(),g.lc());
1254  CanonicalForm gcdcfcg= gcd (cf, cg);
1255  CanonicalForm fp, gp;
1256  CanonicalForm b= 1;
1257  int minCommonDeg= 0;
1258  for (i= tmax (f.level(), g.level()); i > 0; i--)
1259  {
1260    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1261      continue;
1262    else
1263    {
1264      minCommonDeg= tmin (degree (g, i), degree (f, i));
1265      break;
1266    }
1267  }
1268  if (i == 0)
1269    return gcdcfcg;
1270  for (; i > 0; i--)
1271  {
1272    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1273      continue;
1274    else
1275      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1276  }
1277  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
1278     power (CanonicalForm (2), minCommonDeg);
1279  bool equal= false;
1280  i = cf_getNumBigPrimes() - 1;
1281
1282  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
1283  int maxNumVars= tmax (getNumVars (f), getNumVars (g));
1284  //Off (SW_RATIONAL);
1285  while ( true )
1286  {
1287    p = cf_getBigPrime( i );
1288    i--;
1289    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
1290    {
1291      p = cf_getBigPrime( i );
1292      i--;
1293    }
1294    //printf("try p=%d\n",p);
1295    setCharacteristic( p );
1296    fp= mapinto (f);
1297    gp= mapinto (g);
1298    TIMING_START (chinrem_recursion)
1299#ifdef HAVE_NTL
1300    if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
1301      Dp = GCD_small_p (fp, gp, cofp, cogp);
1302    else
1303    {
1304      Dp= gcd_poly (fp, gp);
1305      cofp= fp/Dp;
1306      cogp= gp/Dp;
1307    }
1308#else
1309    Dp= gcd_poly (fp, gp);
1310    cofp= fp/Dp;
1311    cogp= gp/Dp;
1312#endif
1313    TIMING_END_AND_PRINT (chinrem_recursion,
1314                          "time for gcd mod p in modular gcd: ");
1315    Dp /=Dp.lc();
1316    cofp /= lc (cofp);
1317    cogp /= lc (cogp);
1318    setCharacteristic( 0 );
1319    dp_deg=totaldegree(Dp);
1320    if ( dp_deg == 0 )
1321    {
1322      //printf(" -> 1\n");
1323      return CanonicalForm(gcdcfcg);
1324    }
1325    if ( q.isZero() )
1326    {
1327      D = mapinto( Dp );
1328      cof= mapinto (cofp);
1329      cog= mapinto (cogp);
1330      d_deg=dp_deg;
1331      q = p;
1332    }
1333    else
1334    {
1335      if ( dp_deg == d_deg )
1336      {
1337        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1338        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
1339        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
1340        cof= newCof;
1341        cog= newCog;
1342        q = newq;
1343        D = newD;
1344      }
1345      else if ( dp_deg < d_deg )
1346      {
1347        //printf(" were all bad, try more\n");
1348        // all previous p's are bad primes
1349        q = p;
1350        D = mapinto( Dp );
1351        cof= mapinto (cof);
1352        cog= mapinto (cog);
1353        d_deg=dp_deg;
1354        test= 0;
1355        equal= false;
1356      }
1357      else
1358      {
1359        //printf(" was bad, try more\n");
1360      }
1361      //else dp_deg > d_deg: bad prime
1362    }
1363    if ( i >= 0 )
1364    {
1365      TIMING_START (chinrem_reconstruction);
1366      Dn= Farey(D,q);
1367      cofn= Farey(cof,q);
1368      cogn= Farey(cog,q);
1369      TIMING_END_AND_PRINT (chinrem_reconstruction,
1370                           "time for rational reconstruction in modular gcd: ");
1371      int is_rat= isOn (SW_RATIONAL);
1372      On (SW_RATIONAL);
1373      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1374      cofn *= bCommonDen (cofn);
1375      cogn *= bCommonDen (cogn);
1376      if (!is_rat)
1377        Off (SW_RATIONAL);
1378      Dn *=cd;
1379      if (test != Dn)
1380        test= Dn;
1381      else
1382        equal= true;
1383      //Dn /=vcontent(Dn,Variable(1));
1384      TIMING_START (chinrem_termination);
1385      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
1386          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
1387      {
1388        TIMING_END_AND_PRINT (chinrem_termination,
1389                            "time for successful termination in modular gcd: ");
1390        //printf(" -> success\n");
1391        return Dn*gcdcfcg;
1392      }
1393      TIMING_END_AND_PRINT (chinrem_termination,
1394                          "time for unsuccessful termination in modular gcd: ");
1395      equal= false;
1396      //else: try more primes
1397    }
1398    else
1399    { // try other method
1400      //printf("try other gcd\n");
1401      Off(SW_USE_CHINREM_GCD);
1402      D=gcd_poly( f, g );
1403      On(SW_USE_CHINREM_GCD);
1404      return D*gcdcfcg;
1405    }
1406  }
1407}
1408
1409
Note: See TracBrowser for help on using the repository browser.