source: git/factory/cfModResultant.cc @ 5ab7d6

spielwiese
Last change on this file since 5ab7d6 was 72bfc8, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: deleted @internal
  • Property mode set to 100644
File size: 11.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file cfModResultant.cc
5 *
6 * modular resultant algorithm
7 *
8 * @author Martin Lee
9 *
10 **/
11/*****************************************************************************/
12
13#include "config.h"
14
15#include "cf_assert.h"
16
17#include "cfModResultant.h"
18#include "cf_primes.h"
19#include "templates/ftmpl_functions.h"
20#include "cf_map.h"
21#include "cf_algorithm.h"
22#include "cf_iter.h"
23
24#ifdef HAVE_NTL
25#include "NTLconvert.h"
26#endif
27
28//TODO arrange by bound= deg (F,xlevel)*deg (G,i)+deg (G,xlevel)*deg (F, i)
29static inline
30void myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
31                  CFMap & N, const Variable& x)
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  int degsfx, degsgx;
48
49  if (x.level() != 1)
50  {
51    int xlevel= x.level();
52
53    for (int i= 1; i <= n; i++)
54    {
55      if (degsf[i] != 0 && degsg[i] != 0)
56      {
57        both_non_zero++;
58        continue;
59      }
60      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
61      {
62        f_zero++;
63        continue;
64      }
65      if (degsg[i] == 0 && degsf[i] && i <= F.level())
66      {
67        g_zero++;
68        continue;
69      }
70    }
71
72    M.newpair (Variable (xlevel), Variable (1));
73    N.newpair (Variable (1), Variable (xlevel));
74    degsfx= degsf [xlevel];
75    degsgx= degsg [xlevel];
76    degsf [xlevel]= 0;
77    degsg [xlevel]= 0;
78    if (getNumVars (F) == 2 || getNumVars (G) == 2)
79    {
80      int pos= 2;
81      for (int i= 1; i <= n; i++)
82      {
83        if (i != xlevel)
84        {
85          if (i != pos && (degsf[i] != 0 || degsg [i] != 0))
86          {
87            M.newpair (Variable (i), Variable (pos));
88            N.newpair (Variable (pos), Variable (i));
89            pos++;
90          }
91        }
92      }
93      delete [] degsf;
94      delete [] degsg;
95      return;
96    }
97
98    if (both_non_zero == 0)
99    {
100      delete [] degsf;
101      delete [] degsg;
102      return;
103    }
104
105    // map Variables which do not occur in both polynomials to higher levels
106    int k= 1;
107    int l= 1;
108    for (int i= 1; i <= n; i++)
109    {
110      if (i == xlevel)
111        continue;
112      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
113      {
114        if (k + both_non_zero != i)
115        {
116          M.newpair (Variable (i), Variable (k + both_non_zero));
117          N.newpair (Variable (k + both_non_zero), Variable (i));
118        }
119        k++;
120      }
121      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
122      {
123        if (l + g_zero + both_non_zero != i)
124        {
125          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
126          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
127        }
128        l++;
129      }
130    }
131
132    int m= tmax (F.level(), G.level());
133    int min_max_deg;
134    k= both_non_zero;
135    l= 0;
136    int i= 1;
137    while (k > 0)
138    {
139      if (degsf [i] != 0 && degsg [i] != 0)
140        min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
141      else
142        min_max_deg= 0;
143      while (min_max_deg == 0)
144      {
145        i++;
146        if (degsf [i] != 0 && degsg [i] != 0)
147          min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
148        else
149          min_max_deg= 0;
150      }
151      for (int j= i + 1; j <=  m; j++)
152      {
153        if (degsgx*degsf[j] + degsfx*degsg[j] <= min_max_deg &&
154            degsf[j] != 0 && degsg [j] != 0)
155        {
156          min_max_deg= degsgx*degsf[j] + degsfx*degsg[j];
157          l= j;
158        }
159      }
160      if (l != 0)
161      {
162        if (l != k && l != xlevel && k != 1)
163        {
164          M.newpair (Variable (l), Variable(k));
165          N.newpair (Variable (k), Variable(l));
166          degsf[l]= 0;
167          degsg[l]= 0;
168          l= 0;
169        }
170        else
171        {
172          degsf[l]= 0;
173          degsg[l]= 0;
174          l= 0;
175        }
176      }
177      else if (l == 0)
178      {
179        if (i != k && i != xlevel && k != 1)
180        {
181          M.newpair (Variable (i), Variable (k));
182          N.newpair (Variable (k), Variable (i));
183          degsf[i]= 0;
184          degsg[i]= 0;
185        }
186        else
187        {
188          degsf[i]= 0;
189          degsg[i]= 0;
190        }
191        i++;
192      }
193      k--;
194    }
195  }
196  else
197  {
198    //arrange Variables such that no gaps occur
199    for (int i= 1; i <= n; i++)
200    {
201      if (degsf[i] == 0 && degsg[i] == 0)
202      {
203        both_zero++;
204        continue;
205      }
206      else
207      {
208        if (both_zero != 0 && i != 1)
209        {
210          M.newpair (Variable (i), Variable (i - both_zero));
211          N.newpair (Variable (i - both_zero), Variable (i));
212        }
213      }
214    }
215  }
216
217  delete [] degsf;
218  delete [] degsg;
219}
220
221static inline
222CanonicalForm oneNorm (const CanonicalForm& F)
223{
224  if (F.inZ())
225    return abs (F);
226
227  CanonicalForm result= 0;
228  for (CFIterator i= F; i.hasTerms(); i++)
229    result += oneNorm (i.coeff());
230  return result;
231}
232
233// if F and G are both non constant, make sure their level is equal
234static inline
235CanonicalForm uniResultant (const CanonicalForm& F, const CanonicalForm& G)
236{
237#ifdef HAVE_NTL
238  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
239  if (F.inCoeffDomain() && G.inCoeffDomain())
240    return 1;
241  if (F.isZero() || G.isZero())
242    return 0;
243
244  zz_pBak bak;
245  bak.save();
246  zz_p::init (getCharacteristic());
247  zz_pX NTLF= convertFacCF2NTLzzpX (F);
248  zz_pX NTLG= convertFacCF2NTLzzpX (G);
249
250  zz_p NTLResult= resultant (NTLF, NTLG);
251
252  bak.restore();
253  return CanonicalForm (to_long (rep (NTLResult)));
254#else
255  return resultant (F, G, F.mvar());
256#endif
257}
258
259static inline
260void evalPoint (const CanonicalForm& F, const CanonicalForm& G,
261                CanonicalForm& FEval, CanonicalForm& GEval, int& evalPoint)
262{
263  int degF, degG;
264  Variable x= Variable (1);
265  degF= degree (F, x);
266  degG= degree (G, x);
267  do
268  {
269    evalPoint++;
270    if (evalPoint >= getCharacteristic())
271      break;
272    FEval= F (evalPoint, 2);
273    GEval= G (evalPoint, 2);
274    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
275      continue;
276    else
277      return;
278  }
279  while (evalPoint < getCharacteristic());
280}
281
282static inline CanonicalForm
283newtonInterp (const CanonicalForm alpha, const CanonicalForm u,
284              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
285              const Variable & x)
286{
287  CanonicalForm interPoly;
288
289  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
290                            *newtonPoly;
291  return interPoly;
292}
293
294CanonicalForm
295resultantFp (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
296             bool prob)
297{
298  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
299
300  if (A.isZero() || B.isZero())
301    return 0;
302
303  int degAx= degree (A, x);
304  int degBx= degree (B, x);
305  if (A.level() < x.level())
306    return power (A, degBx);
307  if (B.level() < x.level())
308    return power (B, degAx);
309
310  if (degAx == 0)
311    return power (A, degBx);
312  else if (degBx == 0)
313    return power (B, degAx);
314
315  CanonicalForm F= A;
316  CanonicalForm G= B;
317
318  CFMap M, N;
319  myCompress (F, G, M, N, x);
320
321  F= M (F);
322  G= M (G);
323
324  Variable y= Variable (2);
325
326  if (F.isUnivariate() && G.isUnivariate() && F.level() == G.level())
327    return N(uniResultant (F, G));
328
329  int i= -1;
330  CanonicalForm GEval, FEval, recResult, H;
331  CanonicalForm newtonPoly= 1;
332  CanonicalForm modResult= 0;
333
334  Variable z= Variable (1);
335  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
336
337  int count= 0;
338  int equalCount= 0;
339  do
340  {
341    evalPoint (F, G, FEval, GEval, i);
342
343    ASSERT (i < getCharacteristic(), "ran out of points");
344
345    recResult= resultantFp (FEval, GEval, z);
346
347    H= newtonInterp (i, recResult, newtonPoly, modResult, y);
348
349    if (H == modResult)
350      equalCount++;
351
352    count++;
353    if (count > bound || (prob && equalCount == 2))
354      break;
355
356    modResult= H;
357    newtonPoly *= (y - i);
358  } while (1);
359
360  return N (H);
361}
362
363static inline
364CanonicalForm
365balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
366{
367    Variable x = f.mvar();
368    CanonicalForm result = 0, qh = q / 2;
369    CanonicalForm c;
370    CFIterator i;
371    for ( i = f; i.hasTerms(); i++ ) {
372        c = mod( i.coeff(), q );
373        if ( c > qh )
374            result += power( x, i.exp() ) * (c - q);
375        else
376            result += power( x, i.exp() ) * c;
377    }
378    return result;
379}
380
381static inline
382CanonicalForm
383symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
384{
385  CanonicalForm result= 0;
386  if (f.isUnivariate() || f.inCoeffDomain())
387    return balanceUni (f, q);
388  else
389  {
390    Variable x= f.mvar();
391    for (CFIterator i= f; i.hasTerms(); i++)
392      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
393  }
394  return result;
395}
396
397CanonicalForm
398resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
399            bool prob)
400{
401  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
402
403  int degAx= degree (A, x);
404  int degBx= degree (B, x);
405  if (A.level() < x.level())
406    return power (A, degBx);
407  if (B.level() < x.level())
408    return power (B, degAx);
409
410  if (degAx == 0)
411    return power (A, degBx);
412  else if (degBx == 0)
413    return power (B, degAx);
414
415  CanonicalForm F= A;
416  CanonicalForm G= B;
417
418  Variable X= x;
419  if (F.level() != x.level() || G.level() != x.level())
420  {
421    if (F.level() > G.level())
422      X= F.mvar();
423    else
424      X= G.mvar();
425    F= swapvar (F, X, x);
426    G= swapvar (G, X, x);
427  }
428
429  // now X is the main variable
430
431  CanonicalForm d= 0;
432  CanonicalForm dd= 0;
433  CanonicalForm buf;
434  for (CFIterator i= F; i.hasTerms(); i++)
435  {
436    buf= oneNorm (i.coeff());
437    d= (buf > d) ? buf : d;
438  }
439  CanonicalForm e= 0, ee= 0;
440  for (CFIterator i= G; i.hasTerms(); i++)
441  {
442    buf= oneNorm (i.coeff());
443    e= (buf > e) ? buf : e;
444  }
445  d= power (d, degBx);
446  e= power (e, degAx);
447  CanonicalForm bound= 1;
448  for (int i= degBx + degAx; i > 1; i--)
449    bound *= i;
450  bound *= d*e;
451  bound *= 2;
452
453  bool onRational= isOn (SW_RATIONAL);
454  if (onRational)
455    Off (SW_RATIONAL);
456  int i = cf_getNumBigPrimes() - 1;
457  int p;
458  CanonicalForm l= lc (F)*lc(G);
459  CanonicalForm resultModP, q (0), newResult, newQ;
460  CanonicalForm result;
461  do
462  {
463    p = cf_getBigPrime( i );
464    i--;
465    while ( i >= 0 && mod( l, p ) == 0)
466    {
467      p = cf_getBigPrime( i );
468      i--;
469    }
470
471    ASSERT (i >= 0, "ran out of primes"); //sic
472
473    setCharacteristic (p);
474
475    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
476
477    setCharacteristic (0);
478
479    if ( q.isZero() )
480    {
481      result= mapinto(resultModP);
482      q= p;
483    }
484    else
485    {
486      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
487      q= newQ;
488      result= newResult;
489      if (newQ > bound)
490      {
491        result= symmetricRemainder (result, q);
492        break;
493      }
494    }
495  } while (1);
496
497  if (onRational)
498    On (SW_RATIONAL);
499  return swapvar (result, X, x);
500}
501
Note: See TracBrowser for help on using the repository browser.