source: git/factory/cfModResultant.cc @ b7cc2b

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