source: git/factory/cfModResultant.cc @ f6237dd

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