source: git/factory/facAbsFact.cc @ e1a221

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