source: git/factory/algext.cc @ 5079887

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