source: git/factory/cfModResultant.cc @ c1f4d51

spielwiese
Last change on this file since c1f4d51 was c1f4d51, checked in by Hans Schoenemann <hannes@…>, 3 years ago
factory: implement findMinPoly/mapPrimElem and related
  • 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
347CanonicalForm
348resultantFp (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
349             bool prob)
350{
351  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
352
353  if (A.isZero() || B.isZero())
354    return 0;
355
356  int degAx= degree (A, x);
357  int degBx= degree (B, x);
358  if (A.level() < x.level())
359    return power (A, degBx);
360  if (B.level() < x.level())
361    return power (B, degAx);
362
363  if (degAx == 0)
364    return power (A, degBx);
365  else if (degBx == 0)
366    return power (B, degAx);
367
368  if (A.isUnivariate() && B.isUnivariate() && A.level() == B.level())
369    return uniResultant (A, B);
370
371  CanonicalForm F= A;
372  CanonicalForm G= B;
373
374  CFMap M, N;
375  myCompress (F, G, M, N, x);
376
377  F= M (F);
378  G= M (G);
379
380  Variable y= Variable (2);
381
382  CanonicalForm GEval, FEval, recResult, H;
383  CanonicalForm newtonPoly= 1;
384  CanonicalForm modResult= 0;
385
386  Variable z= Variable (1);
387  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
388
389  int p= getCharacteristic();
390  CanonicalForm minpoly;
391  Variable alpha= Variable (tmax (F.level(), G.level()) + 1);
392  bool algExt= hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha);
393  CFGenerator * gen;
394  bool extOfExt= false;
395  Variable v= alpha;
396  CanonicalForm primElemAlpha, imPrimElemAlpha;
397  CFList source,dest;
398  if (!algExt && (p < (1 << 28)))
399  {
400    // pass to an extension of size at least 2^29
401    // for very very large input that is maybe too small though
402    int deg= ceil (29.0*((double) log (2)/log (p)))+1;
403    minpoly= randomIrredpoly (deg, z);
404    alpha= rootOf (minpoly);
405    AlgExtGenerator AlgExtGen (alpha);
406    gen= AlgExtGen.clone();
407    for (int i= 0; i < p; i++) // skip values from the prime field
408      (*gen).next();
409  }
410  else if (!algExt)
411  {
412    FFGenerator FFGen;
413    gen= FFGen.clone();
414  }
415  else
416  {
417    int deg= ceil (29.0*((double) log (2)/log (p)));
418    if (degree (getMipo (alpha)) < deg)
419    {
420      mpz_t field_size;
421      mpz_init (field_size);
422      mpz_ui_pow_ui (field_size, p,
423                 deg + degree (getMipo (alpha)) - deg%degree (getMipo (alpha)));
424
425      // field_size needs to fit in an int because of mapUp, mapDown, length of lists etc.
426      if (mpz_fits_sint_p (field_size))
427      {
428        minpoly= randomIrredpoly (deg + degree (getMipo (alpha))
429                                  - deg%degree (getMipo (alpha)), z);
430        v= rootOf (minpoly);
431        Variable V_buf2;
432        bool primFail= false;
433        extOfExt= true;
434        primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
435        ASSERT (!primFail, "failure in integer factorizer");
436        if (primFail)
437          ; //ERROR
438        else
439          imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
440        F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
441        G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
442      }
443      else
444      {
445        deg= deg - deg % degree (getMipo (alpha));
446        mpz_ui_pow_ui (field_size, p, deg);
447        while (deg / degree (getMipo (alpha)) >= 2 && !mpz_fits_sint_p (field_size))
448        {
449          deg -= degree (getMipo (alpha));
450          mpz_ui_pow_ui (field_size, p, deg);
451        }
452        if (deg != degree (getMipo (alpha)))
453        {
454           minpoly= randomIrredpoly (deg, z);
455           v= rootOf (minpoly);
456           Variable V_buf2;
457           bool primFail= false;
458           extOfExt= true;
459           primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
460           ASSERT (!primFail, "failure in integer factorizer");
461           if (primFail)
462             ; //ERROR
463           else
464             imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
465           F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
466           G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
467        }
468      }
469      mpz_clear (field_size);
470    }
471    AlgExtGenerator AlgExtGen (v);
472    gen= AlgExtGen.clone();
473    for (int i= 0; i < p; i++)
474      (*gen).next();
475  }
476  int count= 0;
477  int equalCount= 0;
478  CanonicalForm point;
479  do
480  {
481    evalPoint (F, G, FEval, GEval, *gen);
482
483    recResult= resultantFp (FEval, GEval, z, prob);
484
485    H= newtonInterp ((*gen).item(), recResult, newtonPoly, modResult, y);
486
487    if (H == modResult)
488      equalCount++;
489    else
490      equalCount= 0;
491
492    count++;
493    if (count > bound || (prob && equalCount == 2 && !H.inCoeffDomain()))
494    {
495      if (!algExt && degree (H, alpha) <= 0)
496        break;
497      else if (algExt)
498      {
499        if (extOfExt && !isInExtension (H, imPrimElemAlpha, 1, primElemAlpha,
500                                        dest, source))
501        {
502          H= mapDown (H, primElemAlpha, imPrimElemAlpha, alpha, dest, source);
503          prune (v);
504          break;
505        }
506        else if (!extOfExt)
507          break;
508      }
509    }
510
511    modResult= H;
512    newtonPoly *= (y - (*gen).item());
513    if ((*gen).hasItems())
514        (*gen).next();
515    else
516      STICKYASSERT (0, "out of evaluation points");
517  } while (1);
518
519  delete gen;
520
521  return N (H);
522}
523
524static inline
525CanonicalForm
526balanceUni ( const CanonicalForm & f, const CanonicalForm & q )
527{
528    Variable x = f.mvar();
529    CanonicalForm result = 0, qh = q / 2;
530    CanonicalForm c;
531    CFIterator i;
532    for ( i = f; i.hasTerms(); i++ ) {
533        c = mod( i.coeff(), q );
534        if ( c > qh )
535            result += power( x, i.exp() ) * (c - q);
536        else
537            result += power( x, i.exp() ) * c;
538    }
539    return result;
540}
541
542static inline
543CanonicalForm
544symmetricRemainder (const CanonicalForm& f, const CanonicalForm& q)
545{
546  CanonicalForm result= 0;
547  if (f.isUnivariate() || f.inCoeffDomain())
548    return balanceUni (f, q);
549  else
550  {
551    Variable x= f.mvar();
552    for (CFIterator i= f; i.hasTerms(); i++)
553      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
554  }
555  return result;
556}
557#ifdef HAVE_NTL // mapPrimitiveElem
558CanonicalForm
559resultantZ (const CanonicalForm& A, const CanonicalForm& B, const Variable& x,
560            bool prob)
561{
562  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
563#ifndef NOASSERT
564  bool isRat= isOn (SW_RATIONAL);
565  On (SW_RATIONAL);
566  ASSERT (bCommonDen (A).isOne(), "input A is rational");
567  ASSERT (bCommonDen (B).isOne(), "input B is rational");
568  if (!isRat)
569    Off (SW_RATIONAL);
570#endif
571
572  int degAx= degree (A, x);
573  int degBx= degree (B, x);
574  if (A.level() < x.level())
575    return power (A, degBx);
576  if (B.level() < x.level())
577    return power (B, degAx);
578
579  if (degAx == 0)
580    return power (A, degBx);
581  else if (degBx == 0)
582    return power (B, degAx);
583
584  CanonicalForm F= A;
585  CanonicalForm G= B;
586
587  Variable X= x;
588  if (F.level() != x.level() || G.level() != x.level())
589  {
590    if (F.level() > G.level())
591      X= F.mvar();
592    else
593      X= G.mvar();
594    F= swapvar (F, X, x);
595    G= swapvar (G, X, x);
596  }
597
598  // now X is the main variable
599
600  CanonicalForm d= 0;
601  CanonicalForm dd= 0;
602  CanonicalForm buf;
603  for (CFIterator i= F; i.hasTerms(); i++)
604  {
605    buf= oneNorm (i.coeff());
606    d= (buf > d) ? buf : d;
607  }
608  CanonicalForm e= 0, ee= 0;
609  for (CFIterator i= G; i.hasTerms(); i++)
610  {
611    buf= oneNorm (i.coeff());
612    e= (buf > e) ? buf : e;
613  }
614  d= power (d, degBx);
615  e= power (e, degAx);
616  CanonicalForm bound= 1;
617  for (int i= degBx + degAx; i > 1; i--)
618    bound *= i;
619  bound *= d*e;
620  bound *= 2;
621
622  bool onRational= isOn (SW_RATIONAL);
623  if (onRational)
624    Off (SW_RATIONAL);
625  int i = cf_getNumBigPrimes() - 1;
626  int p;
627  CanonicalForm l= lc (F)*lc(G);
628  CanonicalForm resultModP, q (0), newResult, newQ;
629  CanonicalForm result;
630  int equalCount= 0;
631  CanonicalForm test, newTest;
632  int count= 0;
633  do
634  {
635    p = cf_getBigPrime( i );
636    i--;
637    while ( i >= 0 && mod( l, p ) == 0)
638    {
639      p = cf_getBigPrime( i );
640      i--;
641    }
642
643    if (i <= 0)
644      return resultant (A, B, x);
645
646    setCharacteristic (p);
647
648    TIMING_START (fac_resultant_p);
649    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
650    TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
651
652    setCharacteristic (0);
653
654    count++;
655    if ( q.isZero() )
656    {
657      result= mapinto(resultModP);
658      q= p;
659    }
660    else
661    {
662      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
663      q= newQ;
664      result= newResult;
665      test= symmetricRemainder (result,q);
666      if (test != newTest)
667      {
668        newTest= test;
669        equalCount= 0;
670      }
671      else
672        equalCount++;
673      if (newQ > bound || (prob && equalCount == 2))
674      {
675        result= test;
676        break;
677      }
678    }
679  } while (1);
680
681  if (onRational)
682    On (SW_RATIONAL);
683  return swapvar (result, X, x);
684}
685#endif
686#endif
687
Note: See TracBrowser for help on using the repository browser.