source: git/factory/cfModResultant.cc @ 538e06

fieker-DuValspielwiese
Last change on this file since 538e06 was d0dbd7, checked in by Hans Schoenemann <hannes@…>, 6 years ago
chg: more NEW_ARRAY/DELETE_ARRY in factory
  • Property mode set to 100644
File size: 16.6 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_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#ifdef HAVE_NTL
258  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
259  if (F.inCoeffDomain() && G.inCoeffDomain())
260    return 1;
261  if (F.isZero() || G.isZero())
262    return 0;
263  Variable alpha;
264
265#ifdef HAVE_FLINT
266  if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha))
267  {
268    nmod_poly_t FLINTF, FLINTG;
269    convertFacCF2nmod_poly_t (FLINTF, F);
270    convertFacCF2nmod_poly_t (FLINTG, G);
271    mp_limb_t FLINTresult= nmod_poly_resultant (FLINTF, FLINTG);
272    nmod_poly_clear (FLINTF);
273    nmod_poly_clear (FLINTG);
274    return CanonicalForm ((long) FLINTresult);
275  }
276#else
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#endif
292  //at this point F or G has an algebraic var.
293  if (fac_NTL_char != getCharacteristic())
294  {
295    fac_NTL_char= getCharacteristic();
296    zz_p::init (getCharacteristic());
297  }
298  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
299  zz_pE::init (NTLMipo);
300  zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
301  zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
302  zz_pE NTLResult= resultant (NTLF, NTLG);
303
304  return convertNTLzzpE2CF (NTLResult, alpha);
305#else
306  return resultant (F, G, F.mvar());
307#endif
308}
309
310static inline
311void evalPoint (const CanonicalForm& F, const CanonicalForm& G,
312                CanonicalForm& FEval, CanonicalForm& GEval,
313                CFGenerator& evalPoint)
314{
315  int degF, degG;
316  Variable x= Variable (1);
317  degF= degree (F, x);
318  degG= degree (G, x);
319  do
320  {
321    if (!evalPoint.hasItems())
322      break;
323    FEval= F (evalPoint.item(), 2);
324    GEval= G (evalPoint.item(), 2);
325    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
326    {
327      evalPoint.next();
328      continue;
329    }
330    else
331      return;
332  }
333  while (evalPoint.hasItems());
334}
335
336static inline CanonicalForm
337newtonInterp (const CanonicalForm & alpha, const CanonicalForm & u,
338              const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
339              const Variable & x)
340{
341  CanonicalForm interPoly;
342
343  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
344                            *newtonPoly;
345  return interPoly;
346}
347
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
525static inline
526CanonicalForm
527balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
528{
529    Variable x = f.mvar();
530    CanonicalForm result = 0, qh = q / 2;
531    CanonicalForm c;
532    CFIterator i;
533    for ( i = f; i.hasTerms(); i++ ) {
534        c = mod( i.coeff(), q );
535        if ( c > qh )
536            result += power( x, i.exp() ) * (c - q);
537        else
538            result += power( x, i.exp() ) * c;
539    }
540    return result;
541}
542
543static inline
544CanonicalForm
545symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
546{
547  CanonicalForm result= 0;
548  if (f.isUnivariate() || f.inCoeffDomain())
549    return balanceUni (f, q);
550  else
551  {
552    Variable x= f.mvar();
553    for (CFIterator i= f; i.hasTerms(); i++)
554      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
555  }
556  return result;
557}
558
559CanonicalForm
560resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
561            bool prob)
562{
563  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
564#ifndef NOASSERT
565  bool isRat= isOn (SW_RATIONAL);
566  On (SW_RATIONAL);
567  ASSERT (bCommonDen (A).isOne(), "input A is rational");
568  ASSERT (bCommonDen (B).isOne(), "input B is rational");
569  if (!isRat)
570    Off (SW_RATIONAL);
571#endif
572
573  int degAx= degree (A, x);
574  int degBx= degree (B, x);
575  if (A.level() < x.level())
576    return power (A, degBx);
577  if (B.level() < x.level())
578    return power (B, degAx);
579
580  if (degAx == 0)
581    return power (A, degBx);
582  else if (degBx == 0)
583    return power (B, degAx);
584
585  CanonicalForm F= A;
586  CanonicalForm G= B;
587
588  Variable X= x;
589  if (F.level() != x.level() || G.level() != x.level())
590  {
591    if (F.level() > G.level())
592      X= F.mvar();
593    else
594      X= G.mvar();
595    F= swapvar (F, X, x);
596    G= swapvar (G, X, x);
597  }
598
599  // now X is the main variable
600
601  CanonicalForm d= 0;
602  CanonicalForm dd= 0;
603  CanonicalForm buf;
604  for (CFIterator i= F; i.hasTerms(); i++)
605  {
606    buf= oneNorm (i.coeff());
607    d= (buf > d) ? buf : d;
608  }
609  CanonicalForm e= 0, ee= 0;
610  for (CFIterator i= G; i.hasTerms(); i++)
611  {
612    buf= oneNorm (i.coeff());
613    e= (buf > e) ? buf : e;
614  }
615  d= power (d, degBx);
616  e= power (e, degAx);
617  CanonicalForm bound= 1;
618  for (int i= degBx + degAx; i > 1; i--)
619    bound *= i;
620  bound *= d*e;
621  bound *= 2;
622
623  bool onRational= isOn (SW_RATIONAL);
624  if (onRational)
625    Off (SW_RATIONAL);
626  int i = cf_getNumBigPrimes() - 1;
627  int p;
628  CanonicalForm l= lc (F)*lc(G);
629  CanonicalForm resultModP, q (0), newResult, newQ;
630  CanonicalForm result;
631  int equalCount= 0;
632  CanonicalForm test, newTest;
633  int count= 0;
634  do
635  {
636    p = cf_getBigPrime( i );
637    i--;
638    while ( i >= 0 && mod( l, p ) == 0)
639    {
640      p = cf_getBigPrime( i );
641      i--;
642    }
643
644    if (i <= 0)
645      return resultant (A, B, x);
646
647    setCharacteristic (p);
648
649    TIMING_START (fac_resultant_p);
650    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
651    TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
652
653    setCharacteristic (0);
654
655    count++;
656    if ( q.isZero() )
657    {
658      result= mapinto(resultModP);
659      q= p;
660    }
661    else
662    {
663      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
664      q= newQ;
665      result= newResult;
666      test= symmetricRemainder (result,q);
667      if (test != newTest)
668      {
669        newTest= test;
670        equalCount= 0;
671      }
672      else
673        equalCount++;
674      if (newQ > bound || (prob && equalCount == 2))
675      {
676        result= test;
677        break;
678      }
679    }
680  } while (1);
681
682  if (onRational)
683    On (SW_RATIONAL);
684  return swapvar (result, X, x);
685}
686#endif
687
Note: See TracBrowser for help on using the repository browser.