source: git/factory/facAbsFact.cc @ 712a5a

jengelh-datetimespielwiese
Last change on this file since 712a5a was 712a5a, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: include issues
  • Property mode set to 100644
File size: 13.3 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_FLINT
32#ifdef HAVE_NTL
33
34TIMING_DEFINE_PRINT(fac_Qa_factorize)
35TIMING_DEFINE_PRINT(fac_evalpoint)
36
37//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)
38int choosePoint (const CanonicalForm& F, int tdegF, CFArray& eval, bool rec)
39{
40  REvaluation E1 (1, 1, IntRandom (25));
41  REvaluation E2 (2, 2, IntRandom (25));
42  if (rec)
43  {
44    E1.nextpoint();
45    E2.nextpoint();
46  }
47  CanonicalForm f, f1, f2, Fp;
48  int i, p;
49  eval=CFArray (2);
50  while (1)
51  {
52    f1= E1(F);
53    if (!f1.isZero() && factorize (f1).length() == 2)
54    {
55      Off (SW_RATIONAL);
56      f= E2(f1);
57      f2= E2 (F);
58      if ((!f.isZero()) && (abs(f)>cf_getSmallPrime (cf_getNumSmallPrimes()-1)))
59      {
60        for (i= cf_getNumPrimes()-1; i >= 0; i--)
61        {
62          if (f % CanonicalForm (cf_getPrime (i)) == 0)
63          {
64            p= cf_getPrime(i);
65            Fp= mod (F,p);
66            if (totaldegree (Fp) == tdegF &&
67                degree (mod (f2,p), 1) == degree (F,1) &&
68                degree (mod (f1, p),2) == degree (F,2))
69            {
70              eval[0]= E1[1];
71              eval[1]= E2[2];
72              return p;
73            }
74          }
75        }
76      }
77      else if (!f.isZero())
78      {
79        for (i= cf_getNumSmallPrimes()-1; i >= 0; i--)
80        {
81          if (f % CanonicalForm (cf_getSmallPrime (i)) == 0)
82          {
83            p= cf_getSmallPrime (i);
84            Fp= mod (F,p);
85            if (totaldegree (Fp) == tdegF &&
86                degree (mod (f2, p),1) == degree (F,1) &&
87                degree (mod (f1,p),2) == degree (F,2))
88            {
89              eval[0]= E1[1];
90              eval[1]= E2[2];
91              return p;
92            }
93          }
94        }
95      }
96      E2.nextpoint();
97      On (SW_RATIONAL);
98    }
99    E1.nextpoint();
100  }
101  return 0;
102}
103
104//G is assumed to be bivariate, irreducible over Q, primitive wrt x and y, compressed
105CFAFList absFactorizeMain (const CanonicalForm& G)
106{
107  CanonicalForm F= bCommonDen (G)*G;
108  Off (SW_RATIONAL);
109  F /= icontent (F);
110  On (SW_RATIONAL);
111  CFArray eval;
112  int minTdeg, tdegF= totaldegree (F);
113  CanonicalForm Fp, smallestFactor;
114  int p;
115  CFFList factors;
116  Variable alpha;
117  bool rec= false;
118  Variable x= Variable (1);
119  Variable y= Variable (2);
120  CanonicalForm bufF= F;
121  CFFListIterator iter;
122differentevalpoint:
123  while (1)
124  {
125    TIMING_START (fac_evalpoint);
126    p= choosePoint (F, tdegF, eval, rec);
127    TIMING_END_AND_PRINT (fac_evalpoint, "time to find eval point: ");
128
129    setCharacteristic (p);
130    Fp=F.mapinto();
131    factors= factorize (Fp);
132
133    if (factors.getFirst().factor().inCoeffDomain())
134      factors.removeFirst();
135
136    if (factors.length() == 1 && factors.getFirst().exp() == 1)
137    {
138      if (absIrredTest (Fp))
139      {
140        setCharacteristic(0);
141        return CFAFList (CFAFactor (G, 1, 1));
142      }
143      else
144      {
145        setCharacteristic (0);
146        if (modularIrredTestWithShift (F))
147        {
148          return CFAFList (CFAFactor (G, 1, 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.