source: git/factory/cfModResultant.cc @ a9c298

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