source: git/factory/algext.cc @ 4a05ed

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