source: git/factory/fieldGCD.cc @ c1b9927

spielwiese
Last change on this file since c1b9927 was 27bb97f, checked in by Hans Schönemann <hannes@…>, 14 years ago
*hannes: debug stuff: out_cf git-svn-id: file:///usr/local/Singular/svn/trunk@12273 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 7.0 KB
Line 
1#include <stdio.h>
2#include <config.h>
3#ifndef NOSTREAMIO
4#ifdef HAVE_IOSTREAM
5#include <iostream>
6#define OSTREAM std::ostream
7#elif defined(HAVE_IOSTREAM_H)
8#include <iostream.h>
9#define OSTREAM ostream
10#endif
11#endif /* NOSTREAMIO */
12
13#include "cf_defs.h"
14#include "canonicalform.h"
15#include "cf_iter.h"
16#include "cf_primes.h"
17#include "cf_algorithm.h"
18#include "algext.h"
19#include "fieldGCD.h"
20#include "cf_map.h"
21#include "cf_generator.h"
22
23void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
24
25
26CanonicalForm fieldGCD( const CanonicalForm & F, const CanonicalForm & G );
27void CRA(const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew);
28
29
30CanonicalForm fieldGCD( const CanonicalForm & F, const CanonicalForm & G )
31{// this is the modular method by Brown
32 // assume F,G are multivariate polys over Z/p for big prime p
33  if(F.isZero())
34  {
35    if(G.isZero())
36      return G; // G is zero
37    if(G.inCoeffDomain())
38      return CanonicalForm(1);
39    return G/lc(G); // return monic G
40  }
41  if(G.isZero()) // F is non-zero
42  {
43    if(F.inCoeffDomain())
44      return CanonicalForm(1);
45    return F/lc(F); // return monic F
46  }
47  if(F.inCoeffDomain() || G.inCoeffDomain())
48    return CanonicalForm(1);
49  //out_cf("F=",F,"\n");
50  //out_cf("G=",G,"\n");
51  CFMap MM,NN;
52  CFArray ps(1,2);
53  ps[1] = F;
54  ps[2] = G;
55  compress(ps,MM,NN); // maps MM, NN are created
56  CanonicalForm f=MM(F);
57  CanonicalForm g=MM(G);
58  // here: f, g are compressed
59  // compute largest variable in f or g (least one is Variable(1))
60  int mv = f.level();
61  if(g.level() > mv)
62    mv = g.level();
63  // here: mv is level of the largest variable in f, g
64  if(mv == 1) // f,g univariate
65    return NN(gcd( f, g )); // do not forget to map back
66  // here: mv > 1
67  CanonicalForm cf = vcontent(f, Variable(2)); // cf is univariate poly in var(1)
68  CanonicalForm cg = vcontent(g, Variable(2));
69  CanonicalForm c = gcd( cf, cg );
70  f/=cf;
71  g/=cg;
72  if(f.inCoeffDomain() || g.inCoeffDomain())
73  {
74    //printf("=============== inCoeffDomain\n");
75    return NN(c);
76  }
77  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
78  int *M = new int[mv+1];
79  for(int i=2; i<=mv; i++)
80    L[i] = M[i] = 0;
81  L = leadDeg(f, L);
82  M = leadDeg(g, M);
83  CanonicalForm gamma = gcd( firstLC(f), firstLC(g) );
84  for(int i=2; i<=mv; i++) // entry at i=1 is not visited
85    if(M[i] < L[i])
86      L[i] = M[i];
87  // L is now upper bound for leading degrees of gcd
88  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
89  for(int i=2; i<=mv; i++)
90    dg_im[i] = 0; // initialize
91  CanonicalForm gamma_image, m=1;
92  CanonicalForm gm=0;
93  CanonicalForm g_image, alpha, gnew, mnew;
94  FFGenerator gen = FFGenerator();
95  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
96  {
97    alpha = gen.item();
98    gamma_image = gamma(alpha, Variable(1)); // plug in alpha for var(1)
99    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
100      continue;
101    g_image = gcd( f(alpha, Variable(1)), g(alpha, Variable(1)) ); // recursive call with one var less
102    if(g_image.inCoeffDomain()) // early termination
103    {
104      //printf("================== inCoeffDomain\n");
105      return NN(c);
106    }
107    for(int i=2; i<=mv; i++)
108      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
109    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
110    if(isEqual(dg_im, L, 2, mv))
111    {
112      g_image /= lc(g_image); // make g_image monic
113      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
114      CRA( g_image, Variable(1)-alpha, gm, m, gnew, mnew );
115      // gnew = gm mod m
116      // gnew = g_image mod var(1)-alpha
117      // mnew = m * (var(1)-alpha)
118      m = mnew;
119      if(gnew == gm) // gnew did not change
120      {
121        g_image = gm / vcontent(gm, Variable(2));
122        //out_cf("=========== try ",g_image,"\n");
123        if(fdivides(g_image,f) && fdivides(g_image,g)) // trial division
124        {
125          //printf("=========== okay\n");
126          return NN(c*g_image);
127        }
128      }
129      gm = gnew;
130      continue;
131    }
132    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
133      continue;
134
135    // here: dg_im < L --> all previous points were unlucky
136    //printf("=========== reset\n");
137    m = CanonicalForm(1); // reset
138    gm = 0; // reset
139    for(int i=2; i<=mv; i++) // tighten bound
140      L[i] = dg_im[i];
141  }
142  // hopefully, we never reach this point
143  Off( SW_USE_fieldGCD );
144  g_image = gcd(f,g);
145  On( SW_USE_fieldGCD );
146  return g_image;
147}
148
149
150void CRA(const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew)
151{ // this is polynomial Chinese Remaindering. q1, q2 are assumed to be coprime.
152  // polynomials of level <= 1 are considered coefficients
153  // xnew = x1 mod q1 (coefficientwise in the above sense)
154  // xnew = x2 mod q2
155  // qnew = q1*q2
156  if(x1.level() <= 1 && x2.level() <= 1) // base case
157  {
158    (void) extgcd(q1,q2,xnew,qnew);
159    xnew = x1 + (x2-x1) * xnew * q1;
160    qnew = q1*q2;
161    xnew = mod(xnew,qnew);
162    return;
163  }
164  CanonicalForm tmp,tmp2;
165  xnew = 0;
166  qnew = q1 * q2;
167  // here: x1.level() > 1 || x2.level() > 1
168  if(x1.level() > x2.level())
169  {
170    for(CFIterator i=x1; i.hasTerms(); i++)
171    {
172      if(i.exp() == 0) // const. term
173      {
174        CRA(i.coeff(),q1,x2,q2,tmp,tmp2);
175        xnew += tmp;
176      }
177      else
178      {
179        CRA(i.coeff(),q1,0,q2,tmp,tmp2);
180        xnew += tmp * power(x1.mvar(),i.exp());
181      }
182    }
183    return;
184  }
185  // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
186  if(x2.level() > x1.level())
187  {
188    for(CFIterator j=x2; j.hasTerms(); j++)
189    {
190      if(j.exp() == 0) // const. term
191      {
192        CRA(x1,q1,j.coeff(),q2,tmp,tmp2);
193        xnew += tmp;
194      }
195      else
196      {
197        CRA(0,q1,j.coeff(),q2,tmp,tmp2);
198        xnew += tmp * power(x2.mvar(),j.exp());
199      }
200    }
201    return;
202  }
203  // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
204  CFIterator i = x1;
205  CFIterator j = x2;
206  while(i.hasTerms() || j.hasTerms())
207  {
208    if(i.hasTerms())
209    {
210      if(j.hasTerms())
211      {
212        if(i.exp() == j.exp())
213        {
214          CRA(i.coeff(),q1,j.coeff(),q2,tmp,tmp2);
215          xnew += tmp * power(x1.mvar(),i.exp());
216          i++; j++;
217        }
218        else
219        {
220          if(i.exp() < j.exp())
221          {
222            CRA(i.coeff(),q1,0,q2,tmp,tmp2);
223            xnew += tmp * power(x1.mvar(),i.exp());
224            i++;
225          }
226          else // i.exp() > j.exp()
227          {
228            CRA(0,q1,j.coeff(),q2,tmp,tmp2);
229            xnew += tmp * power(x1.mvar(),j.exp());
230            j++;
231          }
232        }
233      }
234      else // j is out of terms
235      {
236        CRA(i.coeff(),q1,0,q2,tmp,tmp2);
237        xnew += tmp * power(x1.mvar(),i.exp());
238        i++;
239      }
240    }
241    else // i is out of terms
242    {
243      CRA(0,q1,j.coeff(),q2,tmp,tmp2);
244      xnew += tmp * power(x1.mvar(),j.exp());
245      j++;
246    }
247  }
248}
Note: See TracBrowser for help on using the repository browser.