source: git/factory/cfModResultant.cc @ 80b8fe

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