source: git/factory/cfGcdAlgExt.cc @ 517ffa

fieker-DuValspielwiese
Last change on this file since 517ffa was d0dbd7, checked in by Hans Schoenemann <hannes@…>, 6 years ago
chg: more NEW_ARRAY/DELETE_ARRY in factory
  • Property mode set to 100644
File size: 27.3 KB
Line 
1
2#include "config.h"
3
4
5#ifndef NOSTREAMIO
6#ifdef HAVE_CSTDIO
7#include <cstdio>
8#else
9#include <stdio.h>
10#endif
11#ifdef HAVE_IOSTREAM_H
12#include <iostream.h>
13#elif defined(HAVE_IOSTREAM)
14#include <iostream>
15#endif
16#endif
17
18#include "cf_assert.h"
19#include "timing.h"
20
21#include "templates/ftmpl_functions.h"
22#include "cf_defs.h"
23#include "canonicalform.h"
24#include "cf_iter.h"
25#include "cf_primes.h"
26#include "cf_algorithm.h"
27#include "cfGcdAlgExt.h"
28#include "cfUnivarGcd.h"
29#include "cf_map.h"
30#include "cf_generator.h"
31#include "facMul.h"
32#include "cfNTLzzpEXGCD.h"
33
34#ifdef HAVE_NTL
35#include "NTLconvert.h"
36#endif
37
38#ifdef HAVE_FLINT
39#include "FLINTconvert.h"
40#endif
41
42TIMING_DEFINE_PRINT(alg_content_p)
43TIMING_DEFINE_PRINT(alg_content)
44TIMING_DEFINE_PRINT(alg_compress)
45TIMING_DEFINE_PRINT(alg_termination)
46TIMING_DEFINE_PRINT(alg_termination_p)
47TIMING_DEFINE_PRINT(alg_reconstruction)
48TIMING_DEFINE_PRINT(alg_newton_p)
49TIMING_DEFINE_PRINT(alg_recursion_p)
50TIMING_DEFINE_PRINT(alg_gcd_p)
51TIMING_DEFINE_PRINT(alg_euclid_p)
52
53/// compressing two polynomials F and G, M is used for compressing,
54/// N to reverse the compression
55static int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
56                CFMap & N, bool topLevel)
57{
58  int n= tmax (F.level(), G.level());
59  int * degsf= NEW_ARRAY(int,n + 1);
60  int * degsg= NEW_ARRAY(int,n + 1);
61
62  for (int i = 0; i <= n; i++)
63    degsf[i]= degsg[i]= 0;
64
65  degsf= degrees (F, degsf);
66  degsg= degrees (G, degsg);
67
68  int both_non_zero= 0;
69  int f_zero= 0;
70  int g_zero= 0;
71  int both_zero= 0;
72  int Flevel=F.level();
73  int Glevel=G.level();
74
75  if (topLevel)
76  {
77    for (int i= 1; i <= n; i++)
78    {
79      if (degsf[i] != 0 && degsg[i] != 0)
80      {
81        both_non_zero++;
82        continue;
83      }
84      if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
85      {
86        f_zero++;
87        continue;
88      }
89      if (degsg[i] == 0 && degsf[i] && i <= Flevel)
90      {
91        g_zero++;
92        continue;
93      }
94    }
95
96    if (both_non_zero == 0)
97    {
98      DELETE_ARRAY(degsf);
99      DELETE_ARRAY(degsg);
100      return 0;
101    }
102
103    // map Variables which do not occur in both polynomials to higher levels
104    int k= 1;
105    int l= 1;
106    for (int i= 1; i <= n; i++)
107    {
108      if (degsf[i] != 0 && degsg[i] == 0 && i <= Flevel)
109      {
110        if (k + both_non_zero != i)
111        {
112          M.newpair (Variable (i), Variable (k + both_non_zero));
113          N.newpair (Variable (k + both_non_zero), Variable (i));
114        }
115        k++;
116      }
117      if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
118      {
119        if (l + g_zero + both_non_zero != i)
120        {
121          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
122          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
123        }
124        l++;
125      }
126    }
127
128    // sort Variables x_{i} in increasing order of
129    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
130    int m= tmax (Flevel, Glevel);
131    int min_max_deg;
132    k= both_non_zero;
133    l= 0;
134    int i= 1;
135    while (k > 0)
136    {
137      if (degsf [i] != 0 && degsg [i] != 0)
138        min_max_deg= tmax (degsf[i], degsg[i]);
139      else
140        min_max_deg= 0;
141      while (min_max_deg == 0)
142      {
143        i++;
144        min_max_deg= tmax (degsf[i], degsg[i]);
145        if (degsf [i] != 0 && degsg [i] != 0)
146          min_max_deg= tmax (degsf[i], degsg[i]);
147        else
148          min_max_deg= 0;
149      }
150      for (int j= i + 1; j <=  m; j++)
151      {
152        if (tmax (degsf[j],degsg[j]) <= min_max_deg && degsf[j] != 0 && degsg [j] != 0)
153        {
154          min_max_deg= tmax (degsf[j], degsg[j]);
155          l= j;
156        }
157      }
158      if (l != 0)
159      {
160        if (l != k)
161        {
162          M.newpair (Variable (l), Variable(k));
163          N.newpair (Variable (k), Variable(l));
164          degsf[l]= 0;
165          degsg[l]= 0;
166          l= 0;
167        }
168        else
169        {
170          degsf[l]= 0;
171          degsg[l]= 0;
172          l= 0;
173        }
174      }
175      else if (l == 0)
176      {
177        if (i != k)
178        {
179          M.newpair (Variable (i), Variable (k));
180          N.newpair (Variable (k), Variable (i));
181          degsf[i]= 0;
182          degsg[i]= 0;
183        }
184        else
185        {
186          degsf[i]= 0;
187          degsg[i]= 0;
188        }
189        i++;
190      }
191      k--;
192    }
193  }
194  else
195  {
196    //arrange Variables such that no gaps occur
197    for (int i= 1; i <= n; i++)
198    {
199      if (degsf[i] == 0 && degsg[i] == 0)
200      {
201        both_zero++;
202        continue;
203      }
204      else
205      {
206        if (both_zero != 0)
207        {
208          M.newpair (Variable (i), Variable (i - both_zero));
209          N.newpair (Variable (i - both_zero), Variable (i));
210        }
211      }
212    }
213  }
214
215  DELETE_ARRAY(degsf);
216  DELETE_ARRAY(degsg);
217
218  return 1;
219}
220
221void tryInvert( const CanonicalForm & F, const CanonicalForm & M, CanonicalForm & inv, bool & fail )
222{ // F, M are required to be "univariate" polynomials in an algebraic variable
223  // we try to invert F modulo M
224  if(F.inBaseDomain())
225  {
226    if(F.isZero())
227    {
228      fail = true;
229      return;
230    }
231    inv = 1/F;
232    return;
233  }
234  CanonicalForm b;
235  Variable a = M.mvar();
236  Variable x = Variable(1);
237  if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
238    fail = true;
239  else
240    inv = replacevar( inv, x, a ); // change back to alg var
241}
242
243#ifndef HAVE_NTL
244void tryDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
245                CanonicalForm& R, CanonicalForm& inv, const CanonicalForm& mipo,
246                bool& fail)
247{
248  if (F.inCoeffDomain())
249  {
250    Q= 0;
251    R= F;
252    return;
253  }
254
255  CanonicalForm A, B;
256  Variable x= F.mvar();
257  A= F;
258  B= G;
259  int degA= degree (A, x);
260  int degB= degree (B, x);
261
262  if (degA < degB)
263  {
264    R= A;
265    Q= 0;
266    return;
267  }
268
269  tryInvert (Lc (B), mipo, inv, fail);
270  if (fail)
271    return;
272
273  R= A;
274  Q= 0;
275  CanonicalForm Qi;
276  for (int i= degA -degB; i >= 0; i--)
277  {
278    if (degree (R, x) == i + degB)
279    {
280      Qi= Lc (R)*inv*power (x, i);
281      Qi= reduce (Qi, mipo);
282      R -= Qi*B;
283      R= reduce (R, mipo);
284      Q += Qi;
285    }
286  }
287}
288
289void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm & M, CanonicalForm & result, bool & fail )
290{
291  CanonicalForm P;
292  if(A.inCoeffDomain())
293  {
294    tryInvert( A, M, P, fail );
295    if(fail)
296      return;
297    result = 1;
298    return;
299  }
300  if(B.inCoeffDomain())
301  {
302    tryInvert( B, M, P, fail );
303    if(fail)
304      return;
305    result = 1;
306    return;
307  }
308  // here: both not inCoeffDomain
309  if( A.degree() > B.degree() )
310  {
311    P = A; result = B;
312  }
313  else
314  {
315    P = B; result = A;
316  }
317  CanonicalForm inv;
318  if( result.isZero() )
319  {
320    tryInvert( Lc(P), M, inv, fail );
321    if(fail)
322      return;
323    result = inv*P; // monify result (not reduced, yet)
324    result= reduce (result, M);
325    return;
326  }
327  Variable x = P.mvar();
328  CanonicalForm rem, Q;
329  // here: degree(P) >= degree(result)
330  while(true)
331  {
332    tryDivrem (P, result, Q, rem, inv, M, fail);
333    if (fail)
334      return;
335    if( rem.isZero() )
336    {
337      result *= inv;
338      result= reduce (result, M);
339      return;
340    }
341    if(result.degree(x) >= rem.degree(x))
342    {
343      P = result;
344      result = rem;
345    }
346    else
347      P = rem;
348  }
349}
350#endif
351
352static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
353static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
354static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
355
356static inline CanonicalForm
357tryNewtonInterp (const CanonicalForm & alpha, const CanonicalForm & u,
358              const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
359              const Variable & x, const CanonicalForm& M, bool& fail)
360{
361  CanonicalForm interPoly;
362
363  CanonicalForm inv;
364  tryInvert (newtonPoly (alpha, x), M, inv, fail);
365  if (fail)
366    return 0;
367
368  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
369  return interPoly;
370}
371
372void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
373{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
374  // M is assumed to be monic
375  if(F.isZero())
376  {
377    if(G.isZero())
378    {
379      result = G; // G is zero
380      return;
381    }
382    if(G.inCoeffDomain())
383    {
384      tryInvert(G,M,result,fail);
385      if(fail)
386        return;
387      result = 1;
388      return;
389    }
390    // try to make G monic modulo M
391    CanonicalForm inv;
392    tryInvert(Lc(G),M,inv,fail);
393    if(fail)
394      return;
395    result = inv*G;
396    result= reduce (result, M);
397    return;
398  }
399  if(G.isZero()) // F is non-zero
400  {
401    if(F.inCoeffDomain())
402    {
403      tryInvert(F,M,result,fail);
404      if(fail)
405        return;
406      result = 1;
407      return;
408    }
409    // try to make F monic modulo M
410    CanonicalForm inv;
411    tryInvert(Lc(F),M,inv,fail);
412    if(fail)
413      return;
414    result = inv*F;
415    result= reduce (result, M);
416    return;
417  }
418  // here: F,G both nonzero
419  if(F.inCoeffDomain())
420  {
421    tryInvert(F,M,result,fail);
422    if(fail)
423      return;
424    result = 1;
425    return;
426  }
427  if(G.inCoeffDomain())
428  {
429    tryInvert(G,M,result,fail);
430    if(fail)
431      return;
432    result = 1;
433    return;
434  }
435  TIMING_START (alg_compress)
436  CFMap MM,NN;
437  int lev= myCompress (F, G, MM, NN, topLevel);
438  if (lev == 0)
439  {
440    result= 1;
441    return;
442  }
443  CanonicalForm f=MM(F);
444  CanonicalForm g=MM(G);
445  TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
446  // here: f,g are compressed
447  // compute largest variable in f or g (least one is Variable(1))
448  int mv = f.level();
449  if(g.level() > mv)
450    mv = g.level();
451  // here: mv is level of the largest variable in f, g
452  Variable v1= Variable (1);
453#ifdef HAVE_NTL
454  Variable v= M.mvar();
455  if (fac_NTL_char != getCharacteristic())
456  {
457    fac_NTL_char= getCharacteristic();
458    zz_p::init (getCharacteristic());
459  }
460  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
461  zz_pE::init (NTLMipo);
462  zz_pEX NTLResult;
463  zz_pEX NTLF;
464  zz_pEX NTLG;
465#endif
466  if(mv == 1) // f,g univariate
467  {
468    TIMING_START (alg_euclid_p)
469#ifdef HAVE_NTL
470    NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
471    NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
472    tryNTLGCD (NTLResult, NTLF, NTLG, fail);
473    if (fail)
474      return;
475    result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
476#else
477    tryEuclid(f,g,M,result,fail);
478    if(fail)
479      return;
480#endif
481    TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
482    result= NN (reduce (result, M)); // do not forget to map back
483    return;
484  }
485  TIMING_START (alg_content_p)
486  // here: mv > 1
487  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
488  if(fail)
489    return;
490  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
491  if(fail)
492    return;
493  CanonicalForm c;
494#ifdef HAVE_NTL
495  NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
496  NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
497  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
498  if (fail)
499    return;
500  c= convertNTLzz_pEX2CF (NTLResult, v1, v);
501#else
502  tryEuclid(cf,cg,M,c,fail);
503  if(fail)
504    return;
505#endif
506  // f /= cf
507  f.tryDiv (cf, M, fail);
508  if(fail)
509    return;
510  // g /= cg
511  g.tryDiv (cg, M, fail);
512  if(fail)
513    return;
514  TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
515  if(f.inCoeffDomain())
516  {
517    tryInvert(f,M,result,fail);
518    if(fail)
519      return;
520    result = NN(c);
521    return;
522  }
523  if(g.inCoeffDomain())
524  {
525    tryInvert(g,M,result,fail);
526    if(fail)
527      return;
528    result = NN(c);
529    return;
530  }
531  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
532  int *N = new int[mv+1];
533  for(int i=2; i<=mv; i++)
534    L[i] = N[i] = 0;
535  L = leadDeg(f, L);
536  N = leadDeg(g, N);
537  CanonicalForm gamma;
538  TIMING_START (alg_euclid_p)
539#ifdef HAVE_NTL
540  NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
541  NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
542  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
543  if (fail)
544    return;
545  gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
546#else
547  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
548  if(fail)
549    return;
550#endif
551  TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
552  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
553    if(N[i] < L[i])
554      L[i] = N[i];
555  // L is now upper bound for degrees of gcd
556  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
557  for(int i=2; i<=mv; i++)
558    dg_im[i] = 0; // initialize
559  CanonicalForm gamma_image, m=1;
560  CanonicalForm gm=0;
561  CanonicalForm g_image, alpha, gnew;
562  FFGenerator gen = FFGenerator();
563  Variable x= Variable (1);
564  bool divides= true;
565  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
566  {
567    alpha = gen.item();
568    gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
569    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
570      continue;
571    TIMING_START (alg_recursion_p)
572    tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
573    TIMING_END_AND_PRINT (alg_recursion_p,
574                         "time for recursive calls in alg gcd mod p: ")
575    if(fail)
576      return;
577    g_image = reduce(g_image, M);
578    if(g_image.inCoeffDomain()) // early termination
579    {
580      tryInvert(g_image,M,result,fail);
581      if(fail)
582        return;
583      result = NN(c);
584      return;
585    }
586    for(int i=2; i<=mv; i++)
587      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
588    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
589    if(isEqual(dg_im, L, 2, mv))
590    {
591      CanonicalForm inv;
592      tryInvert (firstLC (g_image), M, inv, fail);
593      if (fail)
594        return;
595      g_image *= inv;
596      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
597      g_image= reduce (g_image, M);
598      TIMING_START (alg_newton_p)
599      gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
600      TIMING_END_AND_PRINT (alg_newton_p,
601                            "time for Newton interpolation in alg gcd mod p: ")
602      // gnew = gm mod m
603      // gnew = g_image mod var(1)-alpha
604      // mnew = m * (var(1)-alpha)
605      if(fail)
606        return;
607      m *= (x - alpha);
608      if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
609      {
610        TIMING_START (alg_termination_p)
611        cf = tryvcontent(gnew, Variable(2), M, fail);
612        if(fail)
613          return;
614        divides = true;
615        g_image= gnew;
616        g_image.tryDiv (cf, M, fail);
617        if(fail)
618          return;
619        divides= tryFdivides (g_image,f, M, fail); // trial division (f)
620        if(fail)
621          return;
622        if(divides)
623        {
624          bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
625          if(fail)
626            return;
627          if(divides2)
628          {
629            result = NN(reduce (c*g_image, M));
630            TIMING_END_AND_PRINT (alg_termination_p,
631                      "time for successful termination test in alg gcd mod p: ")
632            return;
633          }
634        }
635        TIMING_END_AND_PRINT (alg_termination_p,
636                    "time for unsuccessful termination test in alg gcd mod p: ")
637      }
638      gm = gnew;
639      continue;
640    }
641
642    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
643      continue;
644
645    // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
646    m = CanonicalForm(1); // reset
647    gm = 0; // reset
648    for(int i=2; i<=mv; i++) // tighten bound
649      L[i] = dg_im[i];
650  }
651  // we are out of evaluation points
652  fail = true;
653}
654
655static CanonicalForm
656myicontent ( const CanonicalForm & f, const CanonicalForm & c )
657{
658#ifdef HAVE_NTL
659    if (f.isOne() || c.isOne())
660      return 1;
661    if ( f.inBaseDomain() && c.inBaseDomain())
662    {
663      if (c.isZero()) return abs(f);
664      return bgcd( f, c );
665    }
666    else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
667              (f.inCoeffDomain() && c.inBaseDomain()) ||
668              (f.inBaseDomain() && c.inCoeffDomain()))
669    {
670      if (c.isZero()) return abs (f);
671#ifdef HAVE_FLINT
672      fmpz_poly_t FLINTf, FLINTc;
673      convertFacCF2Fmpz_poly_t (FLINTf, f);
674      convertFacCF2Fmpz_poly_t (FLINTc, c);
675      fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
676      CanonicalForm result;
677      if (f.inCoeffDomain())
678        result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
679      else
680        result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
681      fmpz_poly_clear (FLINTc);
682      fmpz_poly_clear (FLINTf);
683      return result;
684#else
685      ZZX NTLf= convertFacCF2NTLZZX (f);
686      ZZX NTLc= convertFacCF2NTLZZX (c);
687      NTLc= GCD (NTLc, NTLf);
688      if (f.inCoeffDomain())
689        return convertNTLZZX2CF(NTLc,f.mvar());
690      else
691        return convertNTLZZX2CF(NTLc,c.mvar());
692#endif
693    }
694    else
695    {
696        CanonicalForm g = c;
697        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
698            g = myicontent( i.coeff(), g );
699        return g;
700    }
701#else
702    return 1;
703#endif
704}
705
706static CanonicalForm myicontent ( const CanonicalForm & f )
707{
708#ifdef HAVE_NTL
709    return myicontent( f, 0 );
710#else
711    return 1;
712#endif
713}
714
715CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
716{ // f,g in Q(a)[x1,...,xn]
717  if(F.isZero())
718  {
719    if(G.isZero())
720      return G; // G is zero
721    if(G.inCoeffDomain())
722      return CanonicalForm(1);
723    CanonicalForm lcinv= 1/Lc (G);
724    return G*lcinv; // return monic G
725  }
726  if(G.isZero()) // F is non-zero
727  {
728    if(F.inCoeffDomain())
729      return CanonicalForm(1);
730    CanonicalForm lcinv= 1/Lc (F);
731    return F*lcinv; // return monic F
732  }
733  if(F.inCoeffDomain() || G.inCoeffDomain())
734    return CanonicalForm(1);
735  // here: both NOT inCoeffDomain
736  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
737  int p, i;
738  int *bound, *other; // degree vectors
739  bool fail;
740  bool off_rational=!isOn(SW_RATIONAL);
741  On( SW_RATIONAL ); // needed by bCommonDen
742  f = F * bCommonDen(F);
743  g = G * bCommonDen(G);
744  TIMING_START (alg_content)
745  CanonicalForm contf= myicontent (f);
746  CanonicalForm contg= myicontent (g);
747  f /= contf;
748  g /= contg;
749  CanonicalForm gcdcfcg= myicontent (contf, contg);
750  TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
751  Variable a, b;
752  if(hasFirstAlgVar(f,a))
753  {
754    if(hasFirstAlgVar(g,b))
755    {
756      if(b.level() > a.level())
757        a = b;
758    }
759  }
760  else
761  {
762    if(!hasFirstAlgVar(g,a))// both not in extension
763    {
764      Off( SW_RATIONAL );
765      Off( SW_USE_QGCD );
766      tmp = gcdcfcg*gcd( f, g );
767      On( SW_USE_QGCD );
768      if (off_rational) Off(SW_RATIONAL);
769      return tmp;
770    }
771  }
772  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
773  setReduce(a,false); // do not reduce expressions modulo mipo
774  tmp = getMipo(a);
775  M = tmp * bCommonDen(tmp);
776  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
777  Off( SW_RATIONAL ); // needed by mod
778  // calculate upper bound for degree vector of gcd
779  int mv = f.level(); i = g.level();
780  if(i > mv)
781    mv = i;
782  // here: mv is level of the largest variable in f, g
783  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
784  other = new int[mv+1];
785  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
786    bound[i] = other[i] = 0;
787  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
788  other = leadDeg(g,other);
789  for(int i=1; i<=mv; i++) // entry at i=0 not visited
790    if(other[i] < bound[i])
791      bound[i] = other[i];
792  // now 'bound' is the smaller vector
793  cl = lc(M) * lc(f) * lc(g);
794  q = 1;
795  D = 0;
796  CanonicalForm test= 0;
797  bool equal= false;
798  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
799  {
800    p = cf_getBigPrime(i);
801    if( mod( cl, p ).isZero() ) // skip lc-bad primes
802      continue;
803    fail = false;
804    setCharacteristic(p);
805    mipo = mapinto(M);
806    mipo /= mipo.lc();
807    // here: mipo is monic
808    TIMING_START (alg_gcd_p)
809    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
810    TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
811    if( fail ) // mipo splits in char p
812    {
813      setCharacteristic(0);
814      continue;
815    }
816    if( Dp.inCoeffDomain() ) // early termination
817    {
818      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
819      setCharacteristic(0);
820      if(fail)
821        continue;
822      setReduce(a,true);
823      if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
824      delete[] other;
825      delete[] bound;
826      return gcdcfcg;
827    }
828    setCharacteristic(0);
829    // here: Dp NOT inCoeffDomain
830    for(int i=1; i<=mv; i++)
831      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
832    other = leadDeg(Dp,other);
833
834    if(isEqual(bound, other, 1, mv)) // equal
835    {
836      chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
837      // tmp = Dp mod p
838      // tmp = D mod q
839      // newq = p*q
840      q = newq;
841      if( D != tmp )
842        D = tmp;
843      On( SW_RATIONAL );
844      TIMING_START (alg_reconstruction)
845      tmp = Farey( D, q ); // Farey
846      tmp *= bCommonDen (tmp);
847      TIMING_END_AND_PRINT (alg_reconstruction,
848                            "time for rational reconstruction in alg gcd: ")
849      setReduce(a,true); // reduce expressions modulo mipo
850      On( SW_RATIONAL ); // needed by fdivides
851      if (test != tmp)
852        test= tmp;
853      else
854        equal= true; // modular image did not add any new information
855      TIMING_START (alg_termination)
856#ifdef HAVE_NTL
857#ifdef HAVE_FLINT
858      if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
859          && f.level() == tmp.level() && tmp.level() == g.level())
860      {
861        CanonicalForm Q, R;
862        newtonDivrem (f, tmp, Q, R);
863        if (R.isZero())
864        {
865          newtonDivrem (g, tmp, Q, R);
866          if (R.isZero())
867          {
868            Off (SW_RATIONAL);
869            setReduce (a,true);
870            if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
871            TIMING_END_AND_PRINT (alg_termination,
872                                 "time for successful termination test in alg gcd: ")
873            delete[] other;
874            delete[] bound;
875            return tmp*gcdcfcg;
876          }
877        }
878      }
879      else
880#endif
881#endif
882      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
883      {
884        Off( SW_RATIONAL );
885        setReduce(a,true);
886        if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
887        TIMING_END_AND_PRINT (alg_termination,
888                            "time for successful termination test in alg gcd: ")
889        delete[] other;
890        delete[] bound;
891        return tmp*gcdcfcg;
892      }
893      TIMING_END_AND_PRINT (alg_termination,
894                          "time for unsuccessful termination test in alg gcd: ")
895      Off( SW_RATIONAL );
896      setReduce(a,false); // do not reduce expressions modulo mipo
897      continue;
898    }
899    if( isLess(bound, other, 1, mv) ) // current prime unlucky
900      continue;
901    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
902    q = p;
903    D = mapinto(Dp); // shortcut CRA // shortcut CRA
904    for(int i=1; i<=mv; i++) // tighten bound
905      bound[i] = other[i];
906  }
907  // hopefully, we never reach this point
908  setReduce(a,true);
909  delete[] other;
910  delete[] bound;
911  Off( SW_USE_QGCD );
912  D = gcdcfcg*gcd( f, g );
913  On( SW_USE_QGCD );
914  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
915  return D;
916}
917
918
919int * leadDeg(const CanonicalForm & f, int *degs)
920{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
921  // if f is in a coeff domain, the zero pointer is returned
922  // 'a' should point to an array of sufficient size level(f)+1
923  if(f.inCoeffDomain())
924    return 0;
925  CanonicalForm tmp = f;
926  do
927  {
928    degs[tmp.level()] = tmp.degree();
929    tmp = LC(tmp);
930  }
931  while(!tmp.inCoeffDomain());
932  return degs;
933}
934
935
936bool isLess(int *a, int *b, int lower, int upper)
937{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
938  for(int i=upper; i>=lower; i--)
939    if(a[i] == b[i])
940      continue;
941    else
942      return a[i] < b[i];
943  return true;
944}
945
946
947bool isEqual(int *a, int *b, int lower, int upper)
948{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
949  for(int i=lower; i<=upper; i++)
950    if(a[i] != b[i])
951      return false;
952  return true;
953}
954
955
956CanonicalForm firstLC(const CanonicalForm & f)
957{ // returns the leading coefficient (LC) of level <= 1
958  CanonicalForm ret = f;
959  while(ret.level() > 1)
960    ret = LC(ret);
961  return ret;
962}
963
964#ifndef HAVE_NTL
965void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
966{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
967  // F and G must have the same level AND level > 0
968  // we try to calculate gcd(F,G) = s*F + t*G
969  // if a zero divisor is encountered, 'fail' is set to one
970  // M is assumed to be monic
971  CanonicalForm P;
972  if(F.inCoeffDomain())
973  {
974    tryInvert( F, M, P, fail );
975    if(fail)
976      return;
977    result = 1;
978    s = P; t = 0;
979    return;
980  }
981  if(G.inCoeffDomain())
982  {
983    tryInvert( G, M, P, fail );
984    if(fail)
985      return;
986    result = 1;
987    s = 0; t = P;
988    return;
989  }
990  // here: both not inCoeffDomain
991  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
992  if( F.degree() > G.degree() )
993  {
994    P = F; result = G;  s=v=0; t=u=1;
995  }
996  else
997  {
998    P = G; result = F; s=v=1; t=u=0;
999  }
1000  Variable x = P.mvar();
1001  // here: degree(P) >= degree(result)
1002  while(true)
1003  {
1004    tryDivrem (P, result, q, rem, inv, M, fail);
1005    if(fail)
1006      return;
1007    if( rem.isZero() )
1008    {
1009      s*=inv;
1010      s= reduce (s, M);
1011      t*=inv;
1012      t= reduce (t, M);
1013      result *= inv; // monify result
1014      result= reduce (result, M);
1015      return;
1016    }
1017    sum += q;
1018    if(result.degree(x) >= rem.degree(x))
1019    {
1020      P=result;
1021      result=rem;
1022      tmp=u-sum*s;
1023      u=s;
1024      s=tmp;
1025      tmp=v-sum*t;
1026      v=t;
1027      t=tmp;
1028      sum = 0; // reset
1029    }
1030    else
1031      P = rem;
1032  }
1033}
1034#endif
1035
1036static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1037{ // as 'content', but takes care of zero divisors
1038  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1039  Variable y = f.mvar();
1040  if ( y == x )
1041    return trycf_content( f, 0, M, fail );
1042  if ( y < x )
1043     return f;
1044  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1045}
1046
1047
1048static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1049{ // as vcontent, but takes care of zero divisors
1050  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1051  if ( f.mvar() <= x )
1052    return trycontent( f, x, M, fail );
1053  CFIterator i;
1054  CanonicalForm d = 0, e, ret;
1055  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1056  {
1057    e = tryvcontent( i.coeff(), x, M, fail );
1058    if(fail)
1059      break;
1060    tryBrownGCD( d, e, M, ret, fail );
1061    d = ret;
1062  }
1063  return d;
1064}
1065
1066
1067static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
1068{ // as cf_content, but takes care of zero divisors
1069  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1070  {
1071    CFIterator i = f;
1072    CanonicalForm tmp = g, result;
1073    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1074    {
1075      tryBrownGCD( i.coeff(), tmp, M, result, fail );
1076      tmp = result;
1077      i++;
1078    }
1079    return result;
1080  }
1081  return abs( f );
1082}
1083
Note: See TracBrowser for help on using the repository browser.