source: git/factory/cfModResultant.cc @ fd80670

spielwiese
Last change on this file since fd80670 was bffe62d, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: skip zz_p::init() if it is already correctly initialized
  • 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  if (fac_NTL_char != getCharacteristic())
261  {
262    fac_NTL_char= getCharacteristic();
263    zz_p::init (getCharacteristic());
264  }
265  zz_pX NTLF= convertFacCF2NTLzzpX (F);
266  zz_pX NTLG= convertFacCF2NTLzzpX (G);
267
268  zz_p NTLResult= resultant (NTLF, NTLG);
269
270  return CanonicalForm (to_long (rep (NTLResult)));
271#endif
272#else
273  return resultant (F, G, F.mvar());
274#endif
275}
276
277static inline
278void evalPoint (const CanonicalForm& F, const CanonicalForm& G,
279                CanonicalForm& FEval, CanonicalForm& GEval, int& evalPoint)
280{
281  int degF, degG;
282  Variable x= Variable (1);
283  degF= degree (F, x);
284  degG= degree (G, x);
285  do
286  {
287    evalPoint++;
288    if (evalPoint >= getCharacteristic())
289      break;
290    FEval= F (evalPoint, 2);
291    GEval= G (evalPoint, 2);
292    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
293      continue;
294    else
295      return;
296  }
297  while (evalPoint < getCharacteristic());
298}
299
300static inline CanonicalForm
301newtonInterp (const CanonicalForm alpha, const CanonicalForm u,
302              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
303              const Variable & x)
304{
305  CanonicalForm interPoly;
306
307  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
308                            *newtonPoly;
309  return interPoly;
310}
311
312CanonicalForm
313resultantFp (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
314             bool prob)
315{
316  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
317
318  if (A.isZero() || B.isZero())
319    return 0;
320
321  int degAx= degree (A, x);
322  int degBx= degree (B, x);
323  if (A.level() < x.level())
324    return power (A, degBx);
325  if (B.level() < x.level())
326    return power (B, degAx);
327
328  if (degAx == 0)
329    return power (A, degBx);
330  else if (degBx == 0)
331    return power (B, degAx);
332
333  if (A.isUnivariate() && B.isUnivariate() && A.level() == B.level())
334    return uniResultant (A, B);
335
336  CanonicalForm F= A;
337  CanonicalForm G= B;
338
339  CFMap M, N;
340  myCompress (F, G, M, N, x);
341
342  F= M (F);
343  G= M (G);
344
345  Variable y= Variable (2);
346
347  int i= -1;
348  CanonicalForm GEval, FEval, recResult, H;
349  CanonicalForm newtonPoly= 1;
350  CanonicalForm modResult= 0;
351
352  Variable z= Variable (1);
353  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
354
355  int count= 0;
356  int equalCount= 0;
357  do
358  {
359    evalPoint (F, G, FEval, GEval, i);
360
361    ASSERT (i < getCharacteristic(), "ran out of points");
362
363    recResult= resultantFp (FEval, GEval, z);
364
365    H= newtonInterp (i, recResult, newtonPoly, modResult, y);
366
367    if (H == modResult)
368      equalCount++;
369    else
370      equalCount= 0;
371
372    count++;
373    if (count > bound || (prob && equalCount == 2))
374      break;
375
376    modResult= H;
377    newtonPoly *= (y - i);
378  } while (1);
379
380  return N (H);
381}
382
383static inline
384CanonicalForm
385balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
386{
387    Variable x = f.mvar();
388    CanonicalForm result = 0, qh = q / 2;
389    CanonicalForm c;
390    CFIterator i;
391    for ( i = f; i.hasTerms(); i++ ) {
392        c = mod( i.coeff(), q );
393        if ( c > qh )
394            result += power( x, i.exp() ) * (c - q);
395        else
396            result += power( x, i.exp() ) * c;
397    }
398    return result;
399}
400
401static inline
402CanonicalForm
403symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
404{
405  CanonicalForm result= 0;
406  if (f.isUnivariate() || f.inCoeffDomain())
407    return balanceUni (f, q);
408  else
409  {
410    Variable x= f.mvar();
411    for (CFIterator i= f; i.hasTerms(); i++)
412      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
413  }
414  return result;
415}
416
417CanonicalForm
418resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
419            bool prob)
420{
421  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
422
423  int degAx= degree (A, x);
424  int degBx= degree (B, x);
425  if (A.level() < x.level())
426    return power (A, degBx);
427  if (B.level() < x.level())
428    return power (B, degAx);
429
430  if (degAx == 0)
431    return power (A, degBx);
432  else if (degBx == 0)
433    return power (B, degAx);
434
435  CanonicalForm F= A;
436  CanonicalForm G= B;
437
438  Variable X= x;
439  if (F.level() != x.level() || G.level() != x.level())
440  {
441    if (F.level() > G.level())
442      X= F.mvar();
443    else
444      X= G.mvar();
445    F= swapvar (F, X, x);
446    G= swapvar (G, X, x);
447  }
448
449  // now X is the main variable
450
451  CanonicalForm d= 0;
452  CanonicalForm dd= 0;
453  CanonicalForm buf;
454  for (CFIterator i= F; i.hasTerms(); i++)
455  {
456    buf= oneNorm (i.coeff());
457    d= (buf > d) ? buf : d;
458  }
459  CanonicalForm e= 0, ee= 0;
460  for (CFIterator i= G; i.hasTerms(); i++)
461  {
462    buf= oneNorm (i.coeff());
463    e= (buf > e) ? buf : e;
464  }
465  d= power (d, degBx);
466  e= power (e, degAx);
467  CanonicalForm bound= 1;
468  for (int i= degBx + degAx; i > 1; i--)
469    bound *= i;
470  bound *= d*e;
471  bound *= 2;
472
473  bool onRational= isOn (SW_RATIONAL);
474  if (onRational)
475    Off (SW_RATIONAL);
476  int i = cf_getNumBigPrimes() - 1;
477  int p;
478  CanonicalForm l= lc (F)*lc(G);
479  CanonicalForm resultModP, q (0), newResult, newQ;
480  CanonicalForm result;
481  int equalCount= 0;
482  CanonicalForm test, newTest;
483  int count= 0;
484  do
485  {
486    p = cf_getBigPrime( i );
487    i--;
488    while ( i >= 0 && mod( l, p ) == 0)
489    {
490      p = cf_getBigPrime( i );
491      i--;
492    }
493
494    ASSERT (i >= 0, "ran out of primes"); //sic
495
496    setCharacteristic (p);
497
498    TIMING_START (fac_resultant_p);
499    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
500    TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
501
502    setCharacteristic (0);
503
504    count++;
505    if ( q.isZero() )
506    {
507      result= mapinto(resultModP);
508      q= p;
509    }
510    else
511    {
512      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
513      q= newQ;
514      result= newResult;
515      test= symmetricRemainder (result,q);
516      if (test != newTest)
517      {
518        newTest= test;
519        equalCount= 0;
520      }
521      else
522        equalCount++;
523      if (newQ > bound || (prob && equalCount == 2))
524      {
525        result= test;
526        break;
527      }
528    }
529  } while (1);
530
531  if (onRational)
532    On (SW_RATIONAL);
533  return swapvar (result, X, x);
534}
535
Note: See TracBrowser for help on using the repository browser.