source: git/factory/cf_gcd.cc @ 1a2d66

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