source: git/factory/algext.cc @ 5df7d0

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