source: git/factory/cfModResultant.cc @ 36dfca

fieker-DuValspielwiese
Last change on this file since 36dfca was 36dfca, checked in by Martin Lee <martinlee84@…>, 10 years ago
fix for #628
  • Property mode set to 100644
File size: 16.5 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#include "cf_irred.h"
27#include "cf_generator.h"
28#include "cf_random.h"
29#include "cf_map_ext.h"
30#include "facFqBivarUtil.h"
31
32#ifdef HAVE_NTL
33#include "NTLconvert.h"
34#endif
35
36#ifdef HAVE_FLINT
37#include "FLINTconvert.h"
38#endif
39
40#ifdef HAVE_NTL
41TIMING_DEFINE_PRINT(fac_resultant_p)
42
43//TODO arrange by bound= deg (F,xlevel)*deg (G,i)+deg (G,xlevel)*deg (F, i)
44static inline
45void myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
46                  CFMap & N, const Variable& x)
47{
48  int n= tmax (F.level(), G.level());
49  int * degsf= new int [n + 1];
50  int * degsg= new int [n + 1];
51
52  for (int i = 0; i <= n; i++)
53    degsf[i]= degsg[i]= 0;
54
55  degsf= degrees (F, degsf);
56  degsg= degrees (G, degsg);
57
58  int both_non_zero= 0;
59  int f_zero= 0;
60  int g_zero= 0;
61  int both_zero= 0;
62  int degsfx, degsgx;
63
64  if (x.level() != 1)
65  {
66    int xlevel= x.level();
67
68    for (int i= 1; i <= n; i++)
69    {
70      if (degsf[i] != 0 && degsg[i] != 0)
71      {
72        both_non_zero++;
73        continue;
74      }
75      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
76      {
77        f_zero++;
78        continue;
79      }
80      if (degsg[i] == 0 && degsf[i] && i <= F.level())
81      {
82        g_zero++;
83        continue;
84      }
85    }
86
87    M.newpair (Variable (xlevel), Variable (1));
88    N.newpair (Variable (1), Variable (xlevel));
89    degsfx= degsf [xlevel];
90    degsgx= degsg [xlevel];
91    degsf [xlevel]= 0;
92    degsg [xlevel]= 0;
93    if ((getNumVars (F) == 2 && getNumVars (G) == 1) ||
94        (getNumVars (G) == 2 && getNumVars (F) == 1) ||
95        (getNumVars (F) == 2 && getNumVars (F) == getNumVars (G)
96         && getVars (F) == getVars (G)))
97    {
98      int pos= 2;
99      for (int i= 1; i <= n; i++)
100      {
101        if (i != xlevel)
102        {
103          if (i != pos && (degsf[i] != 0 || degsg [i] != 0))
104          {
105            M.newpair (Variable (i), Variable (pos));
106            N.newpair (Variable (pos), Variable (i));
107            pos++;
108          }
109        }
110      }
111      delete [] degsf;
112      delete [] degsg;
113      return;
114    }
115
116    if (both_non_zero == 0)
117    {
118      delete [] degsf;
119      delete [] degsg;
120      return;
121    }
122
123    // map Variables which do not occur in both polynomials to higher levels
124    int k= 1;
125    int l= 1;
126    for (int i= 1; i <= n; i++)
127    {
128      if (i == xlevel)
129        continue;
130      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
131      {
132        if (k + both_non_zero != i)
133        {
134          M.newpair (Variable (i), Variable (k + both_non_zero));
135          N.newpair (Variable (k + both_non_zero), Variable (i));
136        }
137        k++;
138      }
139      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
140      {
141        if (l + g_zero + both_non_zero != i)
142        {
143          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
144          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
145        }
146        l++;
147      }
148    }
149
150    int m= n;
151    int min_max_deg;
152    k= both_non_zero;
153    l= 0;
154    int i= 1;
155    while (k > 0)
156    {
157      if (degsf [i] != 0 && degsg [i] != 0)
158        min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
159      else
160        min_max_deg= 0;
161      while (min_max_deg == 0 && i < m + 1)
162      {
163        i++;
164        if (degsf [i] != 0 && degsg [i] != 0)
165          min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
166        else
167          min_max_deg= 0;
168      }
169      for (int j= i + 1; j <= m; j++)
170      {
171        if (degsgx*degsf[j] + degsfx*degsg[j] <= min_max_deg &&
172            degsf[j] != 0 && degsg [j] != 0)
173        {
174          min_max_deg= degsgx*degsf[j] + degsfx*degsg[j];
175          l= j;
176        }
177      }
178      if (l != 0)
179      {
180        if (l != k && l != xlevel && k != 1)
181        {
182          M.newpair (Variable (l), Variable(k));
183          N.newpair (Variable (k), Variable(l));
184          degsf[l]= 0;
185          degsg[l]= 0;
186          l= 0;
187        }
188        else if (l < m + 1)
189        {
190          degsf[l]= 0;
191          degsg[l]= 0;
192          l= 0;
193        }
194      }
195      else if (l == 0)
196      {
197        if (i != k && i != xlevel && k != 1)
198        {
199          M.newpair (Variable (i), Variable (k));
200          N.newpair (Variable (k), Variable (i));
201          degsf[i]= 0;
202          degsg[i]= 0;
203        }
204        else if (i < m + 1)
205        {
206          degsf[i]= 0;
207          degsg[i]= 0;
208        }
209        i++;
210      }
211      k--;
212    }
213  }
214  else
215  {
216    //arrange Variables such that no gaps occur
217    for (int i= 1; i <= n; i++)
218    {
219      if (degsf[i] == 0 && degsg[i] == 0)
220      {
221        both_zero++;
222        continue;
223      }
224      else
225      {
226        if (both_zero != 0 && i != 1)
227        {
228          M.newpair (Variable (i), Variable (i - both_zero));
229          N.newpair (Variable (i - both_zero), Variable (i));
230        }
231      }
232    }
233  }
234
235  delete [] degsf;
236  delete [] degsg;
237}
238
239static inline
240CanonicalForm oneNorm (const CanonicalForm& F)
241{
242  if (F.inZ())
243    return abs (F);
244
245  CanonicalForm result= 0;
246  for (CFIterator i= F; i.hasTerms(); i++)
247    result += oneNorm (i.coeff());
248  return result;
249}
250
251// if F and G are both non constant, make sure their level is equal
252static inline
253CanonicalForm uniResultant (const CanonicalForm& F, const CanonicalForm& G)
254{
255#ifdef HAVE_NTL
256  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
257  if (F.inCoeffDomain() && G.inCoeffDomain())
258    return 1;
259  if (F.isZero() || G.isZero())
260    return 0;
261  Variable alpha;
262
263#ifdef HAVE_FLINT
264  if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha))
265  {
266    nmod_poly_t FLINTF, FLINTG;
267    convertFacCF2nmod_poly_t (FLINTF, F);
268    convertFacCF2nmod_poly_t (FLINTG, G);
269    mp_limb_t FLINTresult= nmod_poly_resultant (FLINTF, FLINTG);
270    nmod_poly_clear (FLINTF);
271    nmod_poly_clear (FLINTG);
272    return CanonicalForm ((long) FLINTresult);
273  }
274#else
275  if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha))
276  {
277    if (fac_NTL_char != getCharacteristic())
278    {
279      fac_NTL_char= getCharacteristic();
280      zz_p::init (getCharacteristic());
281    }
282    zz_pX NTLF= convertFacCF2NTLzzpX (F);
283    zz_pX NTLG= convertFacCF2NTLzzpX (G);
284
285    zz_p NTLResult= resultant (NTLF, NTLG);
286
287    return CanonicalForm (to_long (rep (NTLResult)));
288  }
289#endif
290  //at this point F or G has an algebraic var.
291  if (fac_NTL_char != getCharacteristic())
292  {
293    fac_NTL_char= getCharacteristic();
294    zz_p::init (getCharacteristic());
295  }
296  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
297  zz_pE::init (NTLMipo);
298  zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
299  zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
300  zz_pE NTLResult= resultant (NTLF, NTLG);
301
302  return convertNTLzzpE2CF (NTLResult, alpha);
303#else
304  return resultant (F, G, F.mvar());
305#endif
306}
307
308static inline
309void evalPoint (const CanonicalForm& F, const CanonicalForm& G,
310                CanonicalForm& FEval, CanonicalForm& GEval,
311                CFGenerator& evalPoint)
312{
313  int degF, degG;
314  Variable x= Variable (1);
315  degF= degree (F, x);
316  degG= degree (G, x);
317  do
318  {
319    if (!evalPoint.hasItems())
320      break;
321    FEval= F (evalPoint.item(), 2);
322    GEval= G (evalPoint.item(), 2);
323    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
324    {
325      evalPoint.next();
326      continue;
327    }
328    else
329      return;
330  }
331  while (evalPoint.hasItems());
332}
333
334static inline CanonicalForm
335newtonInterp (const CanonicalForm & alpha, const CanonicalForm & u,
336              const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
337              const Variable & x)
338{
339  CanonicalForm interPoly;
340
341  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
342                            *newtonPoly;
343  return interPoly;
344}
345
346CanonicalForm
347resultantFp (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
348             bool prob)
349{
350  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
351
352  if (A.isZero() || B.isZero())
353    return 0;
354
355  int degAx= degree (A, x);
356  int degBx= degree (B, x);
357  if (A.level() < x.level())
358    return power (A, degBx);
359  if (B.level() < x.level())
360    return power (B, degAx);
361
362  if (degAx == 0)
363    return power (A, degBx);
364  else if (degBx == 0)
365    return power (B, degAx);
366
367  if (A.isUnivariate() && B.isUnivariate() && A.level() == B.level())
368    return uniResultant (A, B);
369
370  CanonicalForm F= A;
371  CanonicalForm G= B;
372
373  CFMap M, N;
374  myCompress (F, G, M, N, x);
375
376  F= M (F);
377  G= M (G);
378
379  Variable y= Variable (2);
380
381  CanonicalForm GEval, FEval, recResult, H;
382  CanonicalForm newtonPoly= 1;
383  CanonicalForm modResult= 0;
384
385  Variable z= Variable (1);
386  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
387
388  int p= getCharacteristic();
389  CanonicalForm minpoly;
390  Variable alpha= Variable (tmax (F.level(), G.level()) + 1);
391  bool algExt= hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha);
392  CFGenerator * gen;
393  bool extOfExt= false;
394  Variable v= alpha;
395  CanonicalForm primElemAlpha, imPrimElemAlpha;
396  CFList source,dest;
397  if (!algExt && (p < (1 << 28)))
398  {
399    // pass to an extension of size at least 2^29
400    // for very very large input that is maybe too small though
401    int deg= ceil (29.0*((double) log (2)/log (p)))+1;
402    minpoly= randomIrredpoly (deg, z);
403    alpha= rootOf (minpoly);
404    AlgExtGenerator AlgExtGen (alpha);
405    gen= AlgExtGen.clone();
406    for (int i= 0; i < p; i++) // skip values from the prime field
407      (*gen).next();
408  }
409  else if (!algExt)
410  {
411    FFGenerator FFGen;
412    gen= FFGen.clone();
413  }
414  else
415  {
416    int deg= ceil (29.0*((double) log (2)/log (p)));
417    if (degree (getMipo (alpha)) < deg)
418    {
419      mpz_t field_size;
420      mpz_init (field_size);
421      mpz_ui_pow_ui (field_size, p,
422                 deg + degree (getMipo (alpha)) - deg%degree (getMipo (alpha)));
423
424      // field_size needs to fit in an int because of mapUp, mapDown, length of lists etc.
425      if (mpz_fits_sint_p (field_size))
426      {
427        minpoly= randomIrredpoly (deg + degree (getMipo (alpha))
428                                  - deg%degree (getMipo (alpha)), z);
429        v= rootOf (minpoly);
430        Variable V_buf2;
431        bool primFail= false;
432        extOfExt= true;
433        primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
434        ASSERT (!primFail, "failure in integer factorizer");
435        if (primFail)
436          ; //ERROR
437        else
438          imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
439        F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
440        G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
441      }
442      else
443      {
444        deg= deg - deg % degree (getMipo (alpha));
445        mpz_ui_pow_ui (field_size, p, deg);
446        while (deg / degree (getMipo (alpha)) >= 2 && !mpz_fits_sint_p (field_size))
447        {
448          deg -= degree (getMipo (alpha));
449          mpz_ui_pow_ui (field_size, p, deg);
450        }
451        if (deg != degree (getMipo (alpha)))
452        {
453           minpoly= randomIrredpoly (deg, z);
454           v= rootOf (minpoly);
455           Variable V_buf2;
456           bool primFail= false;
457           extOfExt= true;
458           primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
459           ASSERT (!primFail, "failure in integer factorizer");
460           if (primFail)
461             ; //ERROR
462           else
463             imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
464           F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
465           G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
466        }
467      }
468      mpz_clear (field_size);
469    }
470    AlgExtGenerator AlgExtGen (v);
471    gen= AlgExtGen.clone();
472    for (int i= 0; i < p; i++)
473      (*gen).next();
474  }
475  int count= 0;
476  int equalCount= 0;
477  CanonicalForm point;
478  do
479  {
480    evalPoint (F, G, FEval, GEval, *gen);
481
482    recResult= resultantFp (FEval, GEval, z, prob);
483
484    H= newtonInterp ((*gen).item(), recResult, newtonPoly, modResult, y);
485
486    if (H == modResult)
487      equalCount++;
488    else
489      equalCount= 0;
490
491    count++;
492    if (count > bound || (prob && equalCount == 2 && !H.isOne()))
493    {
494      if (!algExt && degree (H, alpha) <= 0)
495        break;
496      else if (algExt)
497      {
498        if (extOfExt && !isInExtension (H, imPrimElemAlpha, 1, primElemAlpha,
499                                        dest, source))
500        {
501          H= mapDown (H, primElemAlpha, imPrimElemAlpha, alpha, dest, source);
502          prune (v);
503          break;
504        }
505        else if (!extOfExt)
506          break;
507      }
508    }
509
510    modResult= H;
511    newtonPoly *= (y - (*gen).item());
512    if ((*gen).hasItems())
513        (*gen).next();
514    else
515      STICKYASSERT (0, "out of evaluation points");
516  } while (1);
517
518  delete gen;
519
520  return N (H);
521}
522
523static inline
524CanonicalForm
525balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
526{
527    Variable x = f.mvar();
528    CanonicalForm result = 0, qh = q / 2;
529    CanonicalForm c;
530    CFIterator i;
531    for ( i = f; i.hasTerms(); i++ ) {
532        c = mod( i.coeff(), q );
533        if ( c > qh )
534            result += power( x, i.exp() ) * (c - q);
535        else
536            result += power( x, i.exp() ) * c;
537    }
538    return result;
539}
540
541static inline
542CanonicalForm
543symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
544{
545  CanonicalForm result= 0;
546  if (f.isUnivariate() || f.inCoeffDomain())
547    return balanceUni (f, q);
548  else
549  {
550    Variable x= f.mvar();
551    for (CFIterator i= f; i.hasTerms(); i++)
552      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
553  }
554  return result;
555}
556
557CanonicalForm
558resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
559            bool prob)
560{
561  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
562#ifndef NOASSERT
563  bool isRat= isOn (SW_RATIONAL);
564  On (SW_RATIONAL);
565  ASSERT (bCommonDen (A).isOne(), "input A is rational");
566  ASSERT (bCommonDen (B).isOne(), "input B is rational");
567  if (!isRat)
568    Off (SW_RATIONAL);
569#endif
570
571  int degAx= degree (A, x);
572  int degBx= degree (B, x);
573  if (A.level() < x.level())
574    return power (A, degBx);
575  if (B.level() < x.level())
576    return power (B, degAx);
577
578  if (degAx == 0)
579    return power (A, degBx);
580  else if (degBx == 0)
581    return power (B, degAx);
582
583  CanonicalForm F= A;
584  CanonicalForm G= B;
585
586  Variable X= x;
587  if (F.level() != x.level() || G.level() != x.level())
588  {
589    if (F.level() > G.level())
590      X= F.mvar();
591    else
592      X= G.mvar();
593    F= swapvar (F, X, x);
594    G= swapvar (G, X, x);
595  }
596
597  // now X is the main variable
598
599  CanonicalForm d= 0;
600  CanonicalForm dd= 0;
601  CanonicalForm buf;
602  for (CFIterator i= F; i.hasTerms(); i++)
603  {
604    buf= oneNorm (i.coeff());
605    d= (buf > d) ? buf : d;
606  }
607  CanonicalForm e= 0, ee= 0;
608  for (CFIterator i= G; i.hasTerms(); i++)
609  {
610    buf= oneNorm (i.coeff());
611    e= (buf > e) ? buf : e;
612  }
613  d= power (d, degBx);
614  e= power (e, degAx);
615  CanonicalForm bound= 1;
616  for (int i= degBx + degAx; i > 1; i--)
617    bound *= i;
618  bound *= d*e;
619  bound *= 2;
620
621  bool onRational= isOn (SW_RATIONAL);
622  if (onRational)
623    Off (SW_RATIONAL);
624  int i = cf_getNumBigPrimes() - 1;
625  int p;
626  CanonicalForm l= lc (F)*lc(G);
627  CanonicalForm resultModP, q (0), newResult, newQ;
628  CanonicalForm result;
629  int equalCount= 0;
630  CanonicalForm test, newTest;
631  int count= 0;
632  do
633  {
634    p = cf_getBigPrime( i );
635    i--;
636    while ( i >= 0 && mod( l, p ) == 0)
637    {
638      p = cf_getBigPrime( i );
639      i--;
640    }
641
642    ASSERT (i >= 0, "ran out of primes"); //sic
643
644    setCharacteristic (p);
645
646    TIMING_START (fac_resultant_p);
647    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
648    TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
649
650    setCharacteristic (0);
651
652    count++;
653    if ( q.isZero() )
654    {
655      result= mapinto(resultModP);
656      q= p;
657    }
658    else
659    {
660      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
661      q= newQ;
662      result= newResult;
663      test= symmetricRemainder (result,q);
664      if (test != newTest)
665      {
666        newTest= test;
667        equalCount= 0;
668      }
669      else
670        equalCount++;
671      if (newQ > bound || (prob && equalCount == 2))
672      {
673        result= test;
674        break;
675      }
676    }
677  } while (1);
678
679  if (onRational)
680    On (SW_RATIONAL);
681  return swapvar (result, X, x);
682}
683#endif
684
Note: See TracBrowser for help on using the repository browser.