source: git/factory/algext.cc @ 1682691

fieker-DuValspielwiese
Last change on this file since 1682691 was 1682691, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: more replacement of divisions by multiplications
  • Property mode set to 100644
File size: 29.0 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    CanonicalForm lcinv= 1/Lc (G);
647    return G*lcinv; // return monic G
648  }
649  if(G.isZero()) // F is non-zero
650  {
651    if(F.inCoeffDomain())
652      return CanonicalForm(1);
653    CanonicalForm lcinv= 1/Lc (F);
654    return F*lcinv; // return monic F
655  }
656  if(F.inCoeffDomain() || G.inCoeffDomain())
657    return CanonicalForm(1);
658  // here: both NOT inCoeffDomain
659  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
660  int p, i;
661  int *bound, *other; // degree vectors
662  bool fail;
663  bool off_rational=!isOn(SW_RATIONAL);
664  On( SW_RATIONAL ); // needed by bCommonDen
665  f = F * bCommonDen(F);
666  g = G * bCommonDen(G);
667  CanonicalForm contf= myicontent (f);
668  CanonicalForm contg= myicontent (g);
669  f /= contf;
670  g /= contg;
671  CanonicalForm gcdcfcg= myicontent (contf, contg);
672  Variable a, b;
673  if(hasFirstAlgVar(f,a))
674  {
675    if(hasFirstAlgVar(g,b))
676    {
677      if(b.level() > a.level())
678        a = b;
679    }
680  }
681  else
682  {
683    if(!hasFirstAlgVar(g,a))// both not in extension
684    {
685      Off( SW_RATIONAL );
686      Off( SW_USE_QGCD );
687      tmp = gcdcfcg*gcd( f, g );
688      On( SW_USE_QGCD );
689      if (off_rational) Off(SW_RATIONAL);
690      return tmp;
691    }
692  }
693  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
694  // (in the sequel b is used to swap alg/poly vars)
695  setReduce(a,false); // do not reduce expressions modulo mipo
696  tmp = getMipo(a);
697  M = tmp * bCommonDen(tmp);
698  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
699  Off( SW_RATIONAL ); // needed by mod
700  // calculate upper bound for degree vector of gcd
701  int mv = f.level(); i = g.level();
702  if(i > mv)
703    mv = i;
704  // here: mv is level of the largest variable in f, g
705  b = Variable(mv+1);
706  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
707  other = new int[mv+1];
708  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
709    bound[i] = other[i] = 0;
710  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
711  other = leadDeg(g,other);
712  for(int i=1; i<=mv; i++) // entry at i=0 not visited
713    if(other[i] < bound[i])
714      bound[i] = other[i];
715  // now 'bound' is the smaller vector
716  cl = lc(M) * lc(f) * lc(g);
717  q = 1;
718  D = 0;
719  CanonicalForm test= 0;
720  bool equal= false;
721  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
722  {
723    p = cf_getBigPrime(i);
724    if( mod( cl, p ).isZero() ) // skip lc-bad primes
725      continue;
726    fail = false;
727    setCharacteristic(p);
728    mipo = mapinto(M);
729    mipo /= mipo.lc();
730    // here: mipo is monic
731    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
732    if( fail ) // mipo splits in char p
733      continue;
734    if( Dp.inCoeffDomain() ) // early termination
735    {
736      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
737      if(fail)
738        continue;
739      setReduce(a,true);
740      if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
741      setCharacteristic(0);
742      return gcdcfcg;
743    }
744    setCharacteristic(0);
745    // here: Dp NOT inCoeffDomain
746    for(int i=1; i<=mv; i++)
747      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
748    other = leadDeg(Dp,other);
749
750    if(isEqual(bound, other, 1, mv)) // equal
751    {
752      chineseRemainder( D, q, replacevar( mapinto(Dp), a, b ), p, tmp, newq );
753      // tmp = Dp mod p
754      // tmp = D mod q
755      // newq = p*q
756      q = newq;
757      if( D != tmp )
758        D = tmp;
759      On( SW_RATIONAL );
760      tmp = replacevar( Farey( D, q ), b, a ); // Farey and switch back to alg var
761      setReduce(a,true); // reduce expressions modulo mipo
762      On( SW_RATIONAL ); // needed by fdivides
763      if (test != tmp)
764        test= tmp;
765      else
766        equal= true; // modular image did not add any new information
767      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
768      {
769        Off( SW_RATIONAL );
770        setReduce(a,true);
771        if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
772        return tmp*gcdcfcg;
773      }
774      Off( SW_RATIONAL );
775      setReduce(a,false); // do not reduce expressions modulo mipo
776      continue;
777    }
778    if( isLess(bound, other, 1, mv) ) // current prime unlucky
779      continue;
780    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
781    q = p;
782    D = replacevar( mapinto(Dp), a, b ); // shortcut CRA // shortcut CRA
783    for(int i=1; i<=mv; i++) // tighten bound
784      bound[i] = other[i];
785  }
786  // hopefully, we never reach this point
787  setReduce(a,true);
788  Off( SW_USE_QGCD );
789  D = gcdcfcg*gcd( f, g );
790  On( SW_USE_QGCD );
791  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
792  return D;
793}
794
795
796int * leadDeg(const CanonicalForm & f, int *degs)
797{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
798  // if f is in a coeff domain, the zero pointer is returned
799  // 'a' should point to an array of sufficient size level(f)+1
800  if(f.inCoeffDomain())
801    return 0;
802  CanonicalForm tmp = f;
803  do
804  {
805    degs[tmp.level()] = tmp.degree();
806    tmp = LC(tmp);
807  }
808  while(!tmp.inCoeffDomain());
809  return degs;
810}
811
812
813bool isLess(int *a, int *b, int lower, int upper)
814{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
815  for(int i=upper; i>=lower; i--)
816    if(a[i] == b[i])
817      continue;
818    else
819      return a[i] < b[i];
820  return true;
821}
822
823
824bool isEqual(int *a, int *b, int lower, int upper)
825{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
826  for(int i=lower; i<=upper; i++)
827    if(a[i] != b[i])
828      return false;
829  return true;
830}
831
832
833CanonicalForm firstLC(const CanonicalForm & f)
834{ // returns the leading coefficient (LC) of level <= 1
835  CanonicalForm ret = f;
836  while(ret.level() > 1)
837    ret = LC(ret);
838  return ret;
839}
840
841void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
842{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
843  // F and G must have the same level AND level > 0
844  // we try to calculate gcd(F,G) = s*F + t*G
845  // if a zero divisor is encontered, 'fail' is set to one
846  // M is assumed to be monic
847  CanonicalForm P;
848  if(F.inCoeffDomain())
849  {
850    tryInvert( F, M, P, fail );
851    if(fail)
852      return;
853    result = 1;
854    s = P; t = 0;
855    return;
856  }
857  if(G.inCoeffDomain())
858  {
859    tryInvert( G, M, P, fail );
860    if(fail)
861      return;
862    result = 1;
863    s = 0; t = P;
864    return;
865  }
866  // here: both not inCoeffDomain
867  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
868  if( F.degree() > G.degree() )
869  {
870    P = F; result = G;  s=v=0; t=u=1;
871  }
872  else
873  {
874    P = G; result = F; s=v=1; t=u=0;
875  }
876  Variable x = P.mvar();
877  // here: degree(P) >= degree(result)
878  while(true)
879  {
880    tryDivrem (P, result, q, rem, inv, M, fail);
881    if(fail)
882      return;
883    if( rem.isZero() )
884    {
885      s*=inv;
886      s= reduce (s, M);
887      t*=inv;
888      t= reduce (t, M);
889      result *= inv; // monify result
890      result= reduce (result, M);
891      return;
892    }
893    sum += q;
894    if(result.degree(x) >= rem.degree(x))
895    {
896      P=result;
897      result=rem;
898      tmp=u-sum*s;
899      u=s;
900      s=tmp;
901      tmp=v-sum*t;
902      v=t;
903      t=tmp;
904      sum = 0; // reset
905    }
906    else
907      P = rem;
908  }
909}
910
911
912static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
913{ // as 'content', but takes care of zero divisors
914  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
915  Variable y = f.mvar();
916  if ( y == x )
917    return trycf_content( f, 0, M, fail );
918  if ( y < x )
919     return f;
920  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
921}
922
923
924static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
925{ // as vcontent, but takes care of zero divisors
926  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
927  if ( f.mvar() <= x )
928    return trycontent( f, x, M, fail );
929  CFIterator i;
930  CanonicalForm d = 0, e, ret;
931  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
932  {
933    e = tryvcontent( i.coeff(), x, M, fail );
934    if(fail)
935      break;
936    tryBrownGCD( d, e, M, ret, fail );
937    d = ret;
938  }
939  return d;
940}
941
942
943static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
944{ // as cf_content, but takes care of zero divisors
945  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
946  {
947    CFIterator i = f;
948    CanonicalForm tmp = g, result;
949    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
950    {
951      tryBrownGCD( i.coeff(), tmp, M, result, fail );
952      tmp = result;
953      i++;
954    }
955    return result;
956  }
957  return abs( f );
958}
959
960
961static void tryDivide( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, CanonicalForm & result, bool & divides, bool & fail )
962{ // M "univariate" monic polynomial
963  // f, g polynomials with coeffs modulo M.
964  // if f is divisible by g, 'divides' is set to 1 and 'result' == f/g mod M coefficientwise.
965  // 'fail' is set to 1, iff a zero divisor is encountered.
966  // divides==1 implies fail==0
967  // required: getReduce(M.mvar())==0
968  if(g.inBaseDomain())
969  {
970    result = f/g;
971    divides = true;
972    return;
973  }
974  if(g.inCoeffDomain())
975  {
976    tryInvert(g,M,result,fail);
977    if(fail)
978      return;
979    result = reduce(f*result, M);
980    divides = true;
981    return;
982  }
983  // here: g NOT inCoeffDomain
984  Variable x = g.mvar();
985  if(f.degree(x) < g.degree(x))
986  {
987    divides = false;
988    return;
989  }
990  // here: f.degree(x) > 0 and f.degree(x) >= g.degree(x)
991  CanonicalForm F = f;
992  CanonicalForm q, leadG = LC(g);
993  result = 0;
994  while(!F.isZero())
995  {
996    tryDivide(F.LC(x),leadG,M,q,divides,fail);
997    if(fail || !divides)
998      return;
999    if(F.degree(x)<g.degree(x))
1000    {
1001      divides = false;
1002      return;
1003    }
1004    q *= power(x,F.degree(x)-g.degree(x));
1005    result += q;
1006    F = reduce(F-q*g, M);
1007  }
1008  result = reduce(result, M);
1009  divides = true;
1010}
1011
1012void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
1013{
1014  // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
1015  // F and G must have the same level AND level > 0
1016  // we try to calculate gcd(f,g) = s*F + t*G
1017  // if a zero divisor is encontered, 'fail' is set to one
1018  Variable a, b;
1019  if( !hasFirstAlgVar(F,a) && !hasFirstAlgVar(G,b) ) // note lazy evaluation
1020  {
1021    result = extgcd( F, G, s, t ); // no zero divisors possible
1022    return;
1023  }
1024  if( b.level() > a.level() )
1025    a = b;
1026  // here: a is the biggest alg. var in F and G
1027  CanonicalForm M = getMipo(a);
1028  CanonicalForm P;
1029  if( degree(F) > degree(G) )
1030  {
1031    P=F; result=G; s=0; t=1;
1032  }
1033  else
1034  {
1035    P=G; result=F; s=1; t=0;
1036  }
1037  CanonicalForm inv, rem, q, u, v;
1038  // here: degree(P) >= degree(result)
1039  while(true)
1040  {
1041    tryInvert( Lc(result), M, inv, fail );
1042    if(fail)
1043      return;
1044    // here: Lc(result) is invertible modulo M
1045    q = Lc(P)*inv * power( P.mvar(), degree(P)-degree(result) );
1046    rem = P - q*result;
1047    // here: s*F + t*G = result
1048    if( rem.isZero() )
1049    {
1050      s*=inv;
1051      t*=inv;
1052      result *= inv; // monify result
1053      return;
1054    }
1055    P=result;
1056    result=rem;
1057    rem=u-q*s;
1058    u=s;
1059    s=rem;
1060    rem=v-q*t;
1061    v=t;
1062    t=rem;
1063  }
1064}
1065
1066void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
1067{ // polys of level <= 1 are considered coefficients. q1,q2 are assumed to be coprime
1068  // xnew = x1 mod q1 (coefficientwise in the above sense)
1069  // xnew = x2 mod q2
1070  // qnew = q1*q2
1071  CanonicalForm tmp;
1072  if(x1.level() <= 1 && x2.level() <= 1) // base case
1073  {
1074    tryExtgcd(q1,q2,tmp,xnew,qnew,fail);
1075    if(fail)
1076      return;
1077    xnew = x1 + (x2-x1) * xnew * q1;
1078    qnew = q1*q2;
1079    xnew = mod(xnew,qnew);
1080    return;
1081  }
1082  CanonicalForm tmp2;
1083  xnew = 0;
1084  qnew = q1 * q2;
1085  // here: x1.level() > 1 || x2.level() > 1
1086  if(x1.level() > x2.level())
1087  {
1088    for(CFIterator i=x1; i.hasTerms(); i++)
1089    {
1090      if(i.exp() == 0) // const. term
1091      {
1092        tryCRA(i.coeff(),q1,x2,q2,tmp,tmp2,fail);
1093        if(fail)
1094          return;
1095        xnew += tmp;
1096      }
1097      else
1098      {
1099        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1100        if(fail)
1101          return;
1102        xnew += tmp * power(x1.mvar(),i.exp());
1103      }
1104    }
1105    return;
1106  }
1107  // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
1108  if(x2.level() > x1.level())
1109  {
1110    for(CFIterator j=x2; j.hasTerms(); j++)
1111    {
1112      if(j.exp() == 0) // const. term
1113      {
1114        tryCRA(x1,q1,j.coeff(),q2,tmp,tmp2,fail);
1115        if(fail)
1116          return;
1117        xnew += tmp;
1118      }
1119      else
1120      {
1121        tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1122        if(fail)
1123          return;
1124        xnew += tmp * power(x2.mvar(),j.exp());
1125      }
1126    }
1127    return;
1128  }
1129  // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
1130  CFIterator i = x1;
1131  CFIterator j = x2;
1132  while(i.hasTerms() || j.hasTerms())
1133  {
1134    if(i.hasTerms())
1135    {
1136      if(j.hasTerms())
1137      {
1138        if(i.exp() == j.exp())
1139        {
1140          tryCRA(i.coeff(),q1,j.coeff(),q2,tmp,tmp2,fail);
1141          if(fail)
1142            return;
1143          xnew += tmp * power(x1.mvar(),i.exp());
1144          i++; j++;
1145        }
1146        else
1147        {
1148          if(i.exp() < j.exp())
1149          {
1150            tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1151            if(fail)
1152              return;
1153            xnew += tmp * power(x1.mvar(),i.exp());
1154            i++;
1155          }
1156          else // i.exp() > j.exp()
1157          {
1158            tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1159            if(fail)
1160              return;
1161            xnew += tmp * power(x1.mvar(),j.exp());
1162            j++;
1163          }
1164        }
1165      }
1166      else // j is out of terms
1167      {
1168        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1169        if(fail)
1170          return;
1171        xnew += tmp * power(x1.mvar(),i.exp());
1172        i++;
1173      }
1174    }
1175    else // i is out of terms
1176    {
1177      tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1178      if(fail)
1179        return;
1180      xnew += tmp * power(x1.mvar(),j.exp());
1181      j++;
1182    }
1183  }
1184}
1185
Note: See TracBrowser for help on using the repository browser.