source: git/factory/algext.cc @ c8f95d

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