source: git/factory/cf_gcd.cc @ 9dbf5c

spielwiese
Last change on this file since 9dbf5c was e16f7d, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: deleted fieldGCD
  • Property mode set to 100644
File size: 33.0 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    #endif
741    else 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 b= 1;
1235  int minCommonDeg= 0;
1236  for (i= tmax (f.level(), g.level()); i > 0; i--)
1237  {
1238    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1239      continue;
1240    else
1241    {
1242      minCommonDeg= tmin (degree (g, i), degree (f, i));
1243      break;
1244    }
1245  }
1246  if (i == 0)
1247    return gcdcfcg;
1248  for (; i > 0; i--)
1249  {
1250    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1251      continue;
1252    else
1253      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1254  }
1255  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
1256     power (CanonicalForm (2), minCommonDeg);
1257  bool equal= false;
1258  i = cf_getNumBigPrimes() - 1;
1259
1260  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
1261  //Off (SW_RATIONAL);
1262  while ( true )
1263  {
1264    p = cf_getBigPrime( i );
1265    i--;
1266    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
1267    {
1268      p = cf_getBigPrime( i );
1269      i--;
1270    }
1271    //printf("try p=%d\n",p);
1272    setCharacteristic( p );
1273#ifdef HAVE_NTL
1274    Dp = GCD_small_p (mapinto (f), mapinto (g), cofp, cogp);
1275#else
1276    Dp= gcd_poly (mapinto (f), mapinto (g));
1277    cofp= mapinto (f)/Dp;
1278    cogp= mapinto (g)/Dp;
1279#endif
1280    Dp /=Dp.lc();
1281    cofp /= lc (cofp);
1282    cogp /= lc (cogp);
1283    setCharacteristic( 0 );
1284    dp_deg=totaldegree(Dp);
1285    if ( dp_deg == 0 )
1286    {
1287      //printf(" -> 1\n");
1288      return CanonicalForm(gcdcfcg);
1289    }
1290    if ( q.isZero() )
1291    {
1292      D = mapinto( Dp );
1293      cof= mapinto (cofp);
1294      cog= mapinto (cogp);
1295      d_deg=dp_deg;
1296      q = p;
1297    }
1298    else
1299    {
1300      if ( dp_deg == d_deg )
1301      {
1302        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1303        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
1304        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
1305        cof= newCof;
1306        cog= newCog;
1307        q = newq;
1308        D = newD;
1309      }
1310      else if ( dp_deg < d_deg )
1311      {
1312        //printf(" were all bad, try more\n");
1313        // all previous p's are bad primes
1314        q = p;
1315        D = mapinto( Dp );
1316        cof= mapinto (cof);
1317        cog= mapinto (cog);
1318        d_deg=dp_deg;
1319        test= 0;
1320        equal= false;
1321      }
1322      else
1323      {
1324        //printf(" was bad, try more\n");
1325      }
1326      //else dp_deg > d_deg: bad prime
1327    }
1328    if ( i >= 0 )
1329    {
1330      Dn= Farey(D,q);
1331      cofn= Farey(cof,q);
1332      cogn= Farey(cog,q);
1333      int is_rat= isOn (SW_RATIONAL);
1334      On (SW_RATIONAL);
1335      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1336      cofn *= bCommonDen (cofn);
1337      cogn *= bCommonDen (cogn);
1338      if (!is_rat)
1339        Off (SW_RATIONAL);
1340      Dn *=cd;
1341      if (test != Dn)
1342        test= Dn;
1343      else
1344        equal= true;
1345      //Dn /=vcontent(Dn,Variable(1));
1346      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
1347          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
1348      {
1349        //printf(" -> success\n");
1350        return Dn*gcdcfcg;
1351      }
1352      equal= false;
1353      //else: try more primes
1354    }
1355    else
1356    { // try other method
1357      //printf("try other gcd\n");
1358      Off(SW_USE_CHINREM_GCD);
1359      D=gcd_poly( f, g );
1360      On(SW_USE_CHINREM_GCD);
1361      return D*gcdcfcg;
1362    }
1363  }
1364}
1365
1366
Note: See TracBrowser for help on using the repository browser.