source: git/factory/facAbsFact.cc @ 247b0f

spielwiese
Last change on this file since 247b0f was 386b3d, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: some simplifications
  • 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 "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        return CFAFList (CFAFactor (G, 1, 1));
140      }
141      else
142      {
143        setCharacteristic (0);
144        if (modularIrredTestWithShift (F))
145        {
146          return CFAFList (CFAFactor (G, 1, 1));
147        }
148        rec= true;
149        continue;
150      }
151    }
152    iter= factors;
153    smallestFactor= iter.getItem().factor();
154    while (smallestFactor.isUnivariate() && iter.hasItem())
155    {
156      iter++;
157      if (!iter.hasItem())
158        break;
159      smallestFactor= iter.getItem().factor();
160    }
161
162    minTdeg= totaldegree (smallestFactor);
163    if (iter.hasItem())
164      iter++;
165    for (; iter.hasItem(); iter++)
166    {
167      if (!iter.getItem().factor().isUnivariate() &&
168          (totaldegree (iter.getItem().factor()) < minTdeg))
169      {
170        smallestFactor= iter.getItem().factor();
171        minTdeg= totaldegree (smallestFactor);
172      }
173    }
174    if (tdegF % minTdeg == 0)
175      break;
176    setCharacteristic(0);
177    rec=true;
178  }
179  CanonicalForm Gp= Fp/smallestFactor;
180  Gp= Gp /Lc (Gp);
181
182  CanonicalForm Gpy= Gp (eval[0].mapinto(), 1);
183  CanonicalForm smallestFactorEvaly= smallestFactor (eval[0].mapinto(),1);
184  CanonicalForm Gpx= Gp (eval[1].mapinto(), 2);
185  CanonicalForm smallestFactorEvalx= smallestFactor (eval[1].mapinto(),2);
186
187  bool xValid= !(Gpx.inCoeffDomain() || smallestFactorEvalx.inCoeffDomain() ||
188               !gcd (Gpx, smallestFactorEvalx).inCoeffDomain());
189  bool yValid= !(Gpy.inCoeffDomain() || smallestFactorEvaly.inCoeffDomain() ||
190               !gcd (Gpy, smallestFactorEvaly).inCoeffDomain());
191  if (!xValid && !yValid)
192  {
193    rec= true;
194    setCharacteristic (0);
195    goto differentevalpoint;
196  }
197
198  setCharacteristic (0);
199
200  CanonicalForm mipo;
201
202  int loop, i;
203  if (xValid && yValid)
204  {
205    loop= 3;
206    i=1;
207  }
208  else if (xValid)
209  {
210    loop= 3;
211    i=2;
212  }
213  else
214  {
215    loop= 2;
216    i=1;
217  }
218
219  CFArray mipos= CFArray (loop-i);
220  for (; i < loop; i++)
221  {
222    CanonicalForm Fi= F(eval[i-1],i);
223
224    int s= tdegF/minTdeg + 1;
225    CanonicalForm bound= power (maxNorm (Fi), 2*(s-1));
226    bound *= power (CanonicalForm (s),s-1);
227    bound *= power (CanonicalForm (2), ((s-1)*(s-1))/2); //possible int overflow
228
229    CanonicalForm B = p;
230    long k = 1L;
231    while ( B < bound ) {
232      B *= p;
233      k++;
234    }
235
236    //TODO take floor (log2(k))
237    k= k+1;
238    fmpz_poly_t FLINTFi;
239    convertFacCF2Fmpz_poly_t (FLINTFi, Fi);
240    setCharacteristic (p);
241    nmod_poly_t FLINTFpi, FLINTGpi;
242    if (i == 2)
243    {
244      convertFacCF2nmod_poly_t (FLINTFpi,
245                                smallestFactorEvalx/lc (smallestFactorEvalx));
246      convertFacCF2nmod_poly_t (FLINTGpi, Gpx/lc (Gpx));
247    }
248    else
249    {
250      convertFacCF2nmod_poly_t (FLINTFpi,
251                                smallestFactorEvaly/lc (smallestFactorEvaly));
252      convertFacCF2nmod_poly_t (FLINTGpi, Gpy/lc (Gpy));
253    }
254    nmod_poly_factor_t nmodFactors;
255    nmod_poly_factor_init (nmodFactors);
256    nmod_poly_factor_insert (nmodFactors, FLINTFpi, 1L);
257    nmod_poly_factor_insert (nmodFactors, FLINTGpi, 1L);
258
259    long * link= new long [2];
260    fmpz_poly_t *v= new fmpz_poly_t[2];
261    fmpz_poly_t *w= new fmpz_poly_t[2];
262    fmpz_poly_init(v[0]);
263    fmpz_poly_init(v[1]);
264    fmpz_poly_init(w[0]);
265    fmpz_poly_init(w[1]);
266
267    fmpz_poly_factor_t liftedFactors;
268    fmpz_poly_factor_init (liftedFactors);
269    _fmpz_poly_hensel_start_lift (liftedFactors, link, v, w, FLINTFi,
270                                  nmodFactors, k);
271
272    nmod_poly_factor_clear (nmodFactors);
273    nmod_poly_clear (FLINTFpi);
274    nmod_poly_clear (FLINTGpi);
275
276    setCharacteristic(0);
277    modpk pk= modpk (p,k);
278    CanonicalForm liftedSmallestFactor=
279    convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[0],Variable (1));
280
281    CanonicalForm otherFactor=
282    convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[1],Variable (1));
283    CanonicalForm test= pk (otherFactor*liftedSmallestFactor);
284
285    Off (SW_SYMMETRIC_FF);
286    liftedSmallestFactor= pk (liftedSmallestFactor);
287    if (i == 2)
288      liftedSmallestFactor= liftedSmallestFactor (eval[0],1);
289    else
290      liftedSmallestFactor= liftedSmallestFactor (eval[1],1);
291
292    On (SW_SYMMETRIC_FF);
293    CFMatrix M= CFMatrix (s, s);
294    M(s,s)= power (CanonicalForm (p), k);
295    for (int j= 1; j < s; j++)
296    {
297      M (j,j)= 1;
298      M (j+1,j)= -liftedSmallestFactor;
299    }
300
301    mat_ZZ NTLM= *convertFacCFMatrix2NTLmat_ZZ (M);
302
303    ZZ det;
304
305    transpose (NTLM, NTLM);
306    (void) LLL (det, NTLM, 1L, 1L); //use floating point LLL ?
307    transpose (NTLM, NTLM);
308    M= *convertNTLmat_ZZ2FacCFMatrix (NTLM);
309
310    mipo= 0;
311    for (int j= 1; j <= s; j++)
312      mipo += M (j,1)*power (x,s-j);
313
314    CFFList mipoFactors= factorize (mipo);
315    mipoFactors.removeFirst();
316
317    fmpz_poly_clear (v[0]);
318    fmpz_poly_clear (v[1]);
319    fmpz_poly_clear (w[0]);
320    fmpz_poly_clear (w[1]);
321    delete [] v;
322    delete [] w;
323    delete [] link;
324    fmpz_poly_factor_clear (liftedFactors);
325
326    if (mipoFactors.length() > 1 ||
327        (mipoFactors.length() == 1 && mipoFactors.getFirst().exp() > 1))
328    {
329      if (i+1 >= loop && ((loop-i == 1) || (loop-i==2 && mipos[0].isZero())))
330      {
331        rec=true;
332        goto differentevalpoint;
333      }
334    }
335    else
336      mipos[loop-i-1]= mipo;
337  }
338
339  On (SW_RATIONAL);
340  if (xValid && yValid && !mipos[0].isZero() && !mipos[1].isZero())
341  {
342    if (maxNorm (mipos[0]) < maxNorm (mipos[1]))
343      alpha= rootOf (mipos[0]);
344    else
345      alpha= rootOf (mipos[1]);
346  }
347  else if (xValid && yValid)
348  {
349    if (mipos[0].isZero())
350      alpha= rootOf (mipos[1]);
351    else
352      alpha= rootOf (mipos[0]);
353  }
354  else
355    alpha= rootOf (mipo);
356
357  CanonicalForm F1;
358  CFFList QaF1Factors;
359  int wrongMipo= 0;
360  if (xValid && yValid)
361  {
362    if (degree (F,1) > minTdeg)
363      F1= F (eval[1], 2);
364    else
365      F1= F (eval[0], 1);
366  }
367  else if (xValid)
368    F1= F (eval[1], 2);
369  else
370    F1= F (eval[0], 1);
371
372  bool swap= false;
373  if (F1.level() == 2)
374  {
375    swap= true;
376    F1=swapvar (F1, x, y);
377    F= swapvar (F, x, y);
378  }
379
380  QaF1Factors= factorize (F1, alpha);
381  if (QaF1Factors.getFirst().factor().inCoeffDomain())
382    QaF1Factors.removeFirst();
383  for (iter= QaF1Factors; iter.hasItem(); iter++)
384  {
385    if (degree (iter.getItem().factor()) > minTdeg)
386      wrongMipo++;
387  }
388
389  if (wrongMipo == QaF1Factors.length())
390  {
391    if (xValid && yValid && !mipos[0].isZero() && !mipos[1].isZero())
392    {
393      if (maxNorm (mipos[0]) < maxNorm (mipos[1])) //try the other minpoly
394        alpha= rootOf (mipos[1]);
395      else
396        alpha= rootOf (mipos[0]);
397    }
398    else
399    {
400      rec= true;
401      F= bufF;
402      goto differentevalpoint;
403    }
404
405    wrongMipo= 0;
406    QaF1Factors= factorize (F1, alpha);
407    if (QaF1Factors.getFirst().factor().inCoeffDomain())
408      QaF1Factors.removeFirst();
409    for (iter= QaF1Factors; iter.hasItem(); iter++)
410    {
411      if (degree (iter.getItem().factor()) > minTdeg)
412        wrongMipo++;
413    }
414    if (wrongMipo == QaF1Factors.length())
415    {
416      rec= true;
417      F= bufF;
418      goto differentevalpoint;
419    }
420  }
421
422  CanonicalForm evaluation;
423  if (swap)
424    evaluation= eval[0];
425  else
426    evaluation= eval[1];
427
428  F *= bCommonDen (F);
429  F= F (y + evaluation, y);
430
431  int liftBound= degree (F,y) + 1;
432
433  modpk b= modpk();
434
435  CanonicalForm den= 1;
436
437  mipo= getMipo (alpha);
438
439  CFList uniFactors;
440  for (iter=QaF1Factors; iter.hasItem(); iter++)
441    uniFactors.append (iter.getItem().factor());
442
443  F /= Lc (F1);
444  ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
445  ZZX NTLLcf= convertFacCF2NTLZZX (Lc (F*bCommonDen (F)));
446  ZZ NTLf= resultant (NTLmipo, NTLLcf);
447  ZZ NTLD= discriminant (NTLmipo);
448  den= abs (convertZZ2CF (NTLD*NTLf));
449
450  // make factors elements of Z(a)[x] disable for modularDiophant
451  CanonicalForm multiplier= 1;
452  for (CFListIterator i= uniFactors; i.hasItem(); i++)
453  {
454    multiplier *= bCommonDen (i.getItem());
455    i.getItem()= i.getItem()*bCommonDen(i.getItem());
456  }
457  F *= multiplier;
458  F *= bCommonDen (F);
459
460  Off (SW_RATIONAL);
461  int ii= 0;
462  CanonicalForm discMipo= convertZZ2CF (NTLD);
463  findGoodPrime (bufF*discMipo,ii);
464  findGoodPrime (F1*discMipo,ii);
465  findGoodPrime (F*discMipo,ii);
466
467  int pp=cf_getBigPrime(ii);
468  b = coeffBound( F, pp, mipo );
469  modpk bb= coeffBound (F1, pp, mipo);
470  if (bb.getk() > b.getk() ) b=bb;
471    bb= coeffBound (F, pp, mipo);
472  if (bb.getk() > b.getk() ) b=bb;
473
474  ExtensionInfo dummy= ExtensionInfo (alpha, false);
475  DegreePattern degs= DegreePattern (uniFactors);
476
477  bool earlySuccess= false;
478  CFList earlyFactors;
479  TIMING_START (fac_bi_hensel_lift);
480  uniFactors= henselLiftAndEarly
481              (F, earlySuccess, earlyFactors, degs, liftBound,
482               uniFactors, dummy, evaluation, b, den);
483  TIMING_END_AND_PRINT (fac_bi_hensel_lift,
484                        "time for bivariate hensel lifting over Q: ");
485  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
486
487  CanonicalForm MODl= power (y, liftBound);
488
489  On (SW_RATIONAL);
490  F *= bCommonDen (F);
491  Off (SW_RATIONAL);
492
493  CFList biFactors;
494
495  TIMING_START (fac_bi_factor_recombination);
496  biFactors= factorRecombination (uniFactors, F, MODl, degs, evaluation, 1,
497                                  uniFactors.length()/2, b, den);
498  TIMING_END_AND_PRINT (fac_bi_factor_recombination,
499                        "time for bivariate factor recombination over Q: ");
500
501  On (SW_RATIONAL);
502
503  if (earlySuccess)
504    biFactors= Union (earlyFactors, biFactors);
505  else if (!earlySuccess && degs.getLength() == 1)
506    biFactors= earlyFactors;
507
508  bool swap2= false;
509  appendSwapDecompress (biFactors, CFList(), CFList(), swap, swap2, CFMap());
510  if (isOn (SW_RATIONAL))
511    normalize (biFactors);
512
513  CFAFList result;
514  bool found= false;
515
516  for (CFListIterator i= biFactors; i.hasItem(); i++)
517  {
518    if (totaldegree (i.getItem()) == minTdeg)
519    {
520      result= CFAFList (CFAFactor (i.getItem(), getMipo (alpha), 1));
521      found= true;
522      break;
523    }
524  }
525
526  if (!found)
527  {
528    rec= true;
529    F= bufF;
530    goto differentevalpoint;
531  }
532
533  return result;
534}
535
536#endif
537#endif
538
539
Note: See TracBrowser for help on using the repository browser.