source: git/factory/facAbsFact.cc @ e4e36c

jengelh-datetimespielwiese
Last change on this file since e4e36c was e4e36c, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: in case there is no FLINT use NTL
  • Property mode set to 100644
File size: 14.2 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facAbsFact.cc
5 *
6 * @author Martin Lee
7 *
8 **/
9/*****************************************************************************/
10
11#include "config.h"
12
13#include "timing.h"
14#include "debug.h"
15
16#include "facAbsFact.h"
17#include "facBivar.h"
18#include "facFqBivar.h"
19#include "cf_reval.h"
20#include "cf_primes.h"
21#include "cf_algorithm.h"
22#ifdef HAVE_FLINT
23#include "FLINTconvert.h"
24#include <flint/fmpz_poly_factor.h>
25#endif
26#ifdef HAVE_NTL
27#include "NTLconvert.h"
28#include <NTL/LLL.h>
29#endif
30
31#ifdef HAVE_NTL
32
33TIMING_DEFINE_PRINT(fac_Qa_factorize)
34TIMING_DEFINE_PRINT(fac_evalpoint)
35
36//TODO optimize choice of p -> choose p as large as possible (better than small p since factorization mod p does not require field extension, also less lifting)
37int choosePoint (const CanonicalForm& F, int tdegF, CFArray& eval, bool rec)
38{
39  REvaluation E1 (1, 1, IntRandom (25));
40  REvaluation E2 (2, 2, IntRandom (25));
41  if (rec)
42  {
43    E1.nextpoint();
44    E2.nextpoint();
45  }
46  CanonicalForm f, f1, f2, Fp;
47  int i, p;
48  eval=CFArray (2);
49  while (1)
50  {
51    f1= E1(F);
52    if (!f1.isZero() && factorize (f1).length() == 2)
53    {
54      Off (SW_RATIONAL);
55      f= E2(f1);
56      f2= E2 (F);
57      if ((!f.isZero()) && (abs(f)>cf_getSmallPrime (cf_getNumSmallPrimes()-1)))
58      {
59        for (i= cf_getNumPrimes()-1; i >= 0; i--)
60        {
61          if (f % CanonicalForm (cf_getPrime (i)) == 0)
62          {
63            p= cf_getPrime(i);
64            Fp= mod (F,p);
65            if (totaldegree (Fp) == tdegF &&
66                degree (mod (f2,p), 1) == degree (F,1) &&
67                degree (mod (f1, p),2) == degree (F,2))
68            {
69              eval[0]= E1[1];
70              eval[1]= E2[2];
71              return p;
72            }
73          }
74        }
75      }
76      else if (!f.isZero())
77      {
78        for (i= cf_getNumSmallPrimes()-1; i >= 0; i--)
79        {
80          if (f % CanonicalForm (cf_getSmallPrime (i)) == 0)
81          {
82            p= cf_getSmallPrime (i);
83            Fp= mod (F,p);
84            if (totaldegree (Fp) == tdegF &&
85                degree (mod (f2, p),1) == degree (F,1) &&
86                degree (mod (f1,p),2) == degree (F,2))
87            {
88              eval[0]= E1[1];
89              eval[1]= E2[2];
90              return p;
91            }
92          }
93        }
94      }
95      E2.nextpoint();
96      On (SW_RATIONAL);
97    }
98    E1.nextpoint();
99  }
100  return 0;
101}
102
103//G is assumed to be bivariate, irreducible over Q, primitive wrt x and y, compressed
104CFAFList absFactorizeMain (const CanonicalForm& G)
105{
106  CanonicalForm F= bCommonDen (G)*G;
107  Off (SW_RATIONAL);
108  F /= icontent (F);
109  On (SW_RATIONAL);
110  CFArray eval;
111  int minTdeg, tdegF= totaldegree (F);
112  CanonicalForm Fp, smallestFactor;
113  int p;
114  CFFList factors;
115  Variable alpha;
116  bool rec= false;
117  Variable x= Variable (1);
118  Variable y= Variable (2);
119  CanonicalForm bufF= F;
120  CFFListIterator iter;
121differentevalpoint:
122  while (1)
123  {
124    TIMING_START (fac_evalpoint);
125    p= choosePoint (F, tdegF, eval, rec);
126    TIMING_END_AND_PRINT (fac_evalpoint, "time to find eval point: ");
127
128    setCharacteristic (p);
129    Fp=F.mapinto();
130    factors= factorize (Fp);
131
132    if (factors.getFirst().factor().inCoeffDomain())
133      factors.removeFirst();
134
135    if (factors.length() == 1 && factors.getFirst().exp() == 1)
136    {
137      if (absIrredTest (Fp))
138      {
139        setCharacteristic(0);
140        return CFAFList (CFAFactor (G, 1, 1));
141      }
142      else
143      {
144        setCharacteristic (0);
145        if (modularIrredTestWithShift (F))
146        {
147          return CFAFList (CFAFactor (G, 1, 1));
148        }
149        rec= true;
150        continue;
151      }
152    }
153    iter= factors;
154    smallestFactor= iter.getItem().factor();
155    while (smallestFactor.isUnivariate() && iter.hasItem())
156    {
157      iter++;
158      if (!iter.hasItem())
159        break;
160      smallestFactor= iter.getItem().factor();
161    }
162
163    minTdeg= totaldegree (smallestFactor);
164    if (iter.hasItem())
165      iter++;
166    for (; iter.hasItem(); iter++)
167    {
168      if (!iter.getItem().factor().isUnivariate() &&
169          (totaldegree (iter.getItem().factor()) < minTdeg))
170      {
171        smallestFactor= iter.getItem().factor();
172        minTdeg= totaldegree (smallestFactor);
173      }
174    }
175    if (tdegF % minTdeg == 0)
176      break;
177    setCharacteristic(0);
178    rec=true;
179  }
180  CanonicalForm Gp= Fp/smallestFactor;
181  Gp= Gp /Lc (Gp);
182
183  CanonicalForm Gpy= Gp (eval[0].mapinto(), 1);
184  CanonicalForm smallestFactorEvaly= smallestFactor (eval[0].mapinto(),1);
185  CanonicalForm Gpx= Gp (eval[1].mapinto(), 2);
186  CanonicalForm smallestFactorEvalx= smallestFactor (eval[1].mapinto(),2);
187
188  bool xValid= !(Gpx.inCoeffDomain() || smallestFactorEvalx.inCoeffDomain() ||
189               !gcd (Gpx, smallestFactorEvalx).inCoeffDomain());
190  bool yValid= !(Gpy.inCoeffDomain() || smallestFactorEvaly.inCoeffDomain() ||
191               !gcd (Gpy, smallestFactorEvaly).inCoeffDomain());
192  if (!xValid && !yValid)
193  {
194    rec= true;
195    setCharacteristic (0);
196    goto differentevalpoint;
197  }
198
199  setCharacteristic (0);
200
201  CanonicalForm mipo;
202
203  int loop, i;
204  if (xValid && yValid)
205  {
206    loop= 3;
207    i=1;
208  }
209  else if (xValid)
210  {
211    loop= 3;
212    i=2;
213  }
214  else
215  {
216    loop= 2;
217    i=1;
218  }
219
220  CFArray mipos= CFArray (loop-i);
221  for (; i < loop; i++)
222  {
223    CanonicalForm Fi= F(eval[i-1],i);
224
225    int s= tdegF/minTdeg + 1;
226    CanonicalForm bound= power (maxNorm (Fi), 2*(s-1));
227    bound *= power (CanonicalForm (s),s-1);
228    bound *= power (CanonicalForm (2), ((s-1)*(s-1))/2); //possible int overflow
229
230    CanonicalForm B = p;
231    long k = 1L;
232    while ( B < bound ) {
233      B *= p;
234      k++;
235    }
236
237    //TODO take floor (log2(k))
238    k= k+1;
239#ifdef HAVE_FLINT
240    fmpz_poly_t FLINTFi;
241    convertFacCF2Fmpz_poly_t (FLINTFi, Fi);
242    setCharacteristic (p);
243    nmod_poly_t FLINTFpi, FLINTGpi;
244    if (i == 2)
245    {
246      convertFacCF2nmod_poly_t (FLINTFpi,
247                                smallestFactorEvalx/lc (smallestFactorEvalx));
248      convertFacCF2nmod_poly_t (FLINTGpi, Gpx/lc (Gpx));
249    }
250    else
251    {
252      convertFacCF2nmod_poly_t (FLINTFpi,
253                                smallestFactorEvaly/lc (smallestFactorEvaly));
254      convertFacCF2nmod_poly_t (FLINTGpi, Gpy/lc (Gpy));
255    }
256    nmod_poly_factor_t nmodFactors;
257    nmod_poly_factor_init (nmodFactors);
258    nmod_poly_factor_insert (nmodFactors, FLINTFpi, 1L);
259    nmod_poly_factor_insert (nmodFactors, FLINTGpi, 1L);
260
261    long * link= new long [2];
262    fmpz_poly_t *v= new fmpz_poly_t[2];
263    fmpz_poly_t *w= new fmpz_poly_t[2];
264    fmpz_poly_init(v[0]);
265    fmpz_poly_init(v[1]);
266    fmpz_poly_init(w[0]);
267    fmpz_poly_init(w[1]);
268
269    fmpz_poly_factor_t liftedFactors;
270    fmpz_poly_factor_init (liftedFactors);
271    _fmpz_poly_hensel_start_lift (liftedFactors, link, v, w, FLINTFi,
272                                  nmodFactors, k);
273
274    nmod_poly_factor_clear (nmodFactors);
275    nmod_poly_clear (FLINTFpi);
276    nmod_poly_clear (FLINTGpi);
277
278    setCharacteristic(0);
279    CanonicalForm liftedSmallestFactor=
280    convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[0],Variable (1));
281
282    CanonicalForm otherFactor=
283    convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[1],Variable (1));
284    modpk pk= modpk (p, k);
285#else
286    modpk pk= modpk (p, k);
287    ZZX NTLFi=convertFacCF2NTLZZX (pk (Fi*pk.inverse (lc(Fi))));
288    setCharacteristic (p);
289
290    if (fac_NTL_char != p)
291    {
292      fac_NTL_char= p;
293      zz_p::init (p);
294    }
295    zz_pX NTLFpi, NTLGpi;
296    if (i == 2)
297    {
298      NTLFpi= convertFacCF2NTLzzpX (smallestFactorEvalx/lc (smallestFactorEvalx));
299      NTLGpi= convertFacCF2NTLzzpX (Gpx/lc (Gpx));
300    }
301    else
302    {
303      NTLFpi= convertFacCF2NTLzzpX (smallestFactorEvaly/lc (smallestFactorEvaly));
304      NTLGpi= convertFacCF2NTLzzpX (Gpy/lc (Gpy));
305    }
306    vec_zz_pX modFactors;
307    modFactors.SetLength(2);
308    modFactors[0]= NTLFpi;
309    modFactors[1]= NTLGpi;
310    vec_ZZX liftedFactors;
311    MultiLift (liftedFactors, modFactors, NTLFi, (long) k);
312    setCharacteristic(0);
313    CanonicalForm liftedSmallestFactor=
314                  convertNTLZZX2CF (liftedFactors[0],Variable (1));
315
316    CanonicalForm otherFactor=
317                  convertNTLZZX2CF (liftedFactors[1], Variable (1));
318#endif
319
320    Off (SW_SYMMETRIC_FF);
321    liftedSmallestFactor= pk (liftedSmallestFactor);
322    if (i == 2)
323      liftedSmallestFactor= liftedSmallestFactor (eval[0],1);
324    else
325      liftedSmallestFactor= liftedSmallestFactor (eval[1],1);
326
327    On (SW_SYMMETRIC_FF);
328    CFMatrix M= CFMatrix (s, s);
329    M(s,s)= power (CanonicalForm (p), k);
330    for (int j= 1; j < s; j++)
331    {
332      M (j,j)= 1;
333      M (j+1,j)= -liftedSmallestFactor;
334    }
335
336    mat_ZZ NTLM= *convertFacCFMatrix2NTLmat_ZZ (M);
337
338    ZZ det;
339
340    transpose (NTLM, NTLM);
341    (void) LLL (det, NTLM, 1L, 1L); //use floating point LLL ?
342    transpose (NTLM, NTLM);
343    M= *convertNTLmat_ZZ2FacCFMatrix (NTLM);
344
345    mipo= 0;
346    for (int j= 1; j <= s; j++)
347      mipo += M (j,1)*power (x,s-j);
348
349    CFFList mipoFactors= factorize (mipo);
350    mipoFactors.removeFirst();
351
352#ifdef HAVE_FLINT
353    fmpz_poly_clear (v[0]);
354    fmpz_poly_clear (v[1]);
355    fmpz_poly_clear (w[0]);
356    fmpz_poly_clear (w[1]);
357    delete [] v;
358    delete [] w;
359    delete [] link;
360    fmpz_poly_factor_clear (liftedFactors);
361#endif
362
363    if (mipoFactors.length() > 1 ||
364        (mipoFactors.length() == 1 && mipoFactors.getFirst().exp() > 1))
365    {
366      if (i+1 >= loop && ((loop-i == 1) || (loop-i==2 && mipos[0].isZero())))
367      {
368        rec=true;
369        goto differentevalpoint;
370      }
371    }
372    else
373      mipos[loop-i-1]= mipo;
374  }
375
376  On (SW_RATIONAL);
377  if (xValid && yValid && !mipos[0].isZero() && !mipos[1].isZero())
378  {
379    if (maxNorm (mipos[0]) < maxNorm (mipos[1]))
380      alpha= rootOf (mipos[0]);
381    else
382      alpha= rootOf (mipos[1]);
383  }
384  else if (xValid && yValid)
385  {
386    if (mipos[0].isZero())
387      alpha= rootOf (mipos[1]);
388    else
389      alpha= rootOf (mipos[0]);
390  }
391  else
392    alpha= rootOf (mipo);
393
394  CanonicalForm F1;
395  CFFList QaF1Factors;
396  int wrongMipo= 0;
397  if (xValid && yValid)
398  {
399    if (degree (F,1) > minTdeg)
400      F1= F (eval[1], 2);
401    else
402      F1= F (eval[0], 1);
403  }
404  else if (xValid)
405    F1= F (eval[1], 2);
406  else
407    F1= F (eval[0], 1);
408
409  bool swap= false;
410  if (F1.level() == 2)
411  {
412    swap= true;
413    F1=swapvar (F1, x, y);
414    F= swapvar (F, x, y);
415  }
416
417  QaF1Factors= factorize (F1, alpha);
418  if (QaF1Factors.getFirst().factor().inCoeffDomain())
419    QaF1Factors.removeFirst();
420  for (iter= QaF1Factors; iter.hasItem(); iter++)
421  {
422    if (degree (iter.getItem().factor()) > minTdeg)
423      wrongMipo++;
424  }
425
426  if (wrongMipo == QaF1Factors.length())
427  {
428    if (xValid && yValid && !mipos[0].isZero() && !mipos[1].isZero())
429    {
430      if (maxNorm (mipos[0]) < maxNorm (mipos[1])) //try the other minpoly
431        alpha= rootOf (mipos[1]);
432      else
433        alpha= rootOf (mipos[0]);
434    }
435    else
436    {
437      rec= true;
438      F= bufF;
439      goto differentevalpoint;
440    }
441
442    wrongMipo= 0;
443    QaF1Factors= factorize (F1, alpha);
444    if (QaF1Factors.getFirst().factor().inCoeffDomain())
445      QaF1Factors.removeFirst();
446    for (iter= QaF1Factors; iter.hasItem(); iter++)
447    {
448      if (degree (iter.getItem().factor()) > minTdeg)
449        wrongMipo++;
450    }
451    if (wrongMipo == QaF1Factors.length())
452    {
453      rec= true;
454      F= bufF;
455      goto differentevalpoint;
456    }
457  }
458
459  CanonicalForm evaluation;
460  if (swap)
461    evaluation= eval[0];
462  else
463    evaluation= eval[1];
464
465  F *= bCommonDen (F);
466  F= F (y + evaluation, y);
467
468  int liftBound= degree (F,y) + 1;
469
470  modpk b= modpk();
471
472  CanonicalForm den= 1;
473
474  mipo= getMipo (alpha);
475
476  CFList uniFactors;
477  for (iter=QaF1Factors; iter.hasItem(); iter++)
478    uniFactors.append (iter.getItem().factor());
479
480  F /= Lc (F1);
481  ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
482  ZZX NTLLcf= convertFacCF2NTLZZX (Lc (F*bCommonDen (F)));
483  ZZ NTLf= resultant (NTLmipo, NTLLcf);
484  ZZ NTLD= discriminant (NTLmipo);
485  den= abs (convertZZ2CF (NTLD*NTLf));
486
487  // make factors elements of Z(a)[x] disable for modularDiophant
488  CanonicalForm multiplier= 1;
489  for (CFListIterator i= uniFactors; i.hasItem(); i++)
490  {
491    multiplier *= bCommonDen (i.getItem());
492    i.getItem()= i.getItem()*bCommonDen(i.getItem());
493  }
494  F *= multiplier;
495  F *= bCommonDen (F);
496
497  Off (SW_RATIONAL);
498  int ii= 0;
499  CanonicalForm discMipo= convertZZ2CF (NTLD);
500  findGoodPrime (bufF*discMipo,ii);
501  findGoodPrime (F1*discMipo,ii);
502  findGoodPrime (F*discMipo,ii);
503
504  int pp=cf_getBigPrime(ii);
505  b = coeffBound( F, pp, mipo );
506  modpk bb= coeffBound (F1, pp, mipo);
507  if (bb.getk() > b.getk() ) b=bb;
508    bb= coeffBound (F, pp, mipo);
509  if (bb.getk() > b.getk() ) b=bb;
510
511  ExtensionInfo dummy= ExtensionInfo (alpha, false);
512  DegreePattern degs= DegreePattern (uniFactors);
513
514  bool earlySuccess= false;
515  CFList earlyFactors;
516  TIMING_START (fac_bi_hensel_lift);
517  uniFactors= henselLiftAndEarly
518              (F, earlySuccess, earlyFactors, degs, liftBound,
519               uniFactors, dummy, evaluation, b, den);
520  TIMING_END_AND_PRINT (fac_bi_hensel_lift,
521                        "time for bivariate hensel lifting over Q: ");
522  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
523
524  CanonicalForm MODl= power (y, liftBound);
525
526  On (SW_RATIONAL);
527  F *= bCommonDen (F);
528  Off (SW_RATIONAL);
529
530  CFList biFactors;
531
532  TIMING_START (fac_bi_factor_recombination);
533  biFactors= factorRecombination (uniFactors, F, MODl, degs, evaluation, 1,
534                                  uniFactors.length()/2, b, den);
535  TIMING_END_AND_PRINT (fac_bi_factor_recombination,
536                        "time for bivariate factor recombination over Q: ");
537
538  On (SW_RATIONAL);
539
540  if (earlySuccess)
541    biFactors= Union (earlyFactors, biFactors);
542  else if (!earlySuccess && degs.getLength() == 1)
543    biFactors= earlyFactors;
544
545  bool swap2= false;
546  appendSwapDecompress (biFactors, CFList(), CFList(), swap, swap2, CFMap());
547  if (isOn (SW_RATIONAL))
548    normalize (biFactors);
549
550  CFAFList result;
551  bool found= false;
552
553  for (CFListIterator i= biFactors; i.hasItem(); i++)
554  {
555    if (totaldegree (i.getItem()) == minTdeg)
556    {
557      result= CFAFList (CFAFactor (i.getItem(), getMipo (alpha), 1));
558      found= true;
559      break;
560    }
561  }
562
563  if (!found)
564  {
565    rec= true;
566    F= bufF;
567    goto differentevalpoint;
568  }
569
570  return result;
571}
572
573#endif
574
575
Note: See TracBrowser for help on using the repository browser.