source: git/factory/cf_gcd.cc @ 2488dc3

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