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

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