source: git/factory/cfModResultant.cc @ a37b34

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