source: git/factory/algext.cc @ 4782bc

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