source: git/factory/facAbsFact.cc @ efd410

jengelh-datetimespielwiese
Last change on this file since efd410 was efd410, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: top level abs factorization interface
  • Property mode set to 100644
File size: 6.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facAbsFact.cc
5 *
6 * @author Martin Lee
7 *
8 **/
9/*****************************************************************************/
10
11#include "facAbsFact.h"
12#include "cf_reval.h"
13#include "cf_primes.h"
14#include "cf_algorithm.h"
15#ifdef HAVE_FLINT
16#include "FLINTconvert.h"
17#include <flint/fmpz_poly_factor.h>
18#endif
19#ifdef HAVE_NTL
20#include "NTLconvert.h"
21#include <NTL/LLL.h>
22#endif
23
24#ifdef HAVE_FLINT
25#ifdef HAVE_NTL
26
27//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)
28int choosePoint (const CanonicalForm& F, int tdegF, CFArray& eval)
29{
30  REvaluation E1 (1, 1, IntRandom (25));
31  REvaluation E2 (2, 2, IntRandom (25));
32  //E1.nextpoint();
33  //E2.nextpoint();
34  CanonicalForm f, Fp;
35  int i;
36  eval=CFArray (2);
37  while (1)
38  {
39    f= E1(F);
40    if (!f.isZero() && factorize (f).length() == 2)
41    {
42      Off (SW_RATIONAL);
43      f= E2(f);
44      if (!f.isZero() && f > cf_getSmallPrime (cf_getNumSmallPrimes()))
45      {
46        for (i= cf_getNumPrimes()-1; i > 0; i--)
47        {
48          if (f % CanonicalForm (cf_getPrime (i)) == 0)
49          {
50            Fp= mod (F,cf_getPrime(i));
51            if (totaldegree (Fp) == tdegF)
52            {
53              eval[0]= E1[1];
54              eval[1]= E2[2];
55              return cf_getPrime(i);
56            }
57          }
58        }
59      }
60      else if (!f.isZero())
61      {
62        for (i= cf_getNumSmallPrimes()-1; i > 0; i--)
63        {
64          if (f % CanonicalForm (cf_getSmallPrime (i)) == 0)
65          {
66            Fp= mod (F,cf_getSmallPrime(i));
67            if (totaldegree (Fp) == tdegF)
68            {
69              eval[0]= E1[1];
70              eval[1]= E2[2];
71              return cf_getSmallPrime(i);
72            }
73          }
74        }
75      }
76      E2.nextpoint();
77      On (SW_RATIONAL);
78    }
79    E1.nextpoint();
80  }
81  return 0;
82}
83
84CFAFList absFactorizeMain (const CanonicalForm& G)
85{
86  //F is assumed to be bivariate, irreducible over Q, primitive wrt x and y, compressed
87
88  CanonicalForm F= bCommonDen (G)*G;
89  Off (SW_RATIONAL);
90  F /= icontent (F);
91  On (SW_RATIONAL);
92  CFArray eval;
93  int minTdeg, tdegF= totaldegree (F);
94  CanonicalForm Fp, smallestFactor;
95  int p;
96  while (1)
97  {
98    p= choosePoint (F, tdegF, eval);
99
100    setCharacteristic (p);
101    Fp=F.mapinto();
102    CFFList factors= factorize (Fp);
103
104    factors.removeFirst();
105    CFFListIterator iter= factors;
106    smallestFactor= iter.getItem().factor();
107    minTdeg= totaldegree (smallestFactor);
108    iter++;
109    for (; iter.hasItem(); iter++)
110    {
111      if (totaldegree (iter.getItem().factor()) < minTdeg)
112      {
113        smallestFactor= iter.getItem().factor();
114        minTdeg= totaldegree (smallestFactor);
115      }
116    }
117    if (tdegF % minTdeg == 0)
118      break;
119    //TODO else
120  }
121  CanonicalForm Gp= Fp/smallestFactor;
122  Gp= Gp (eval[0].mapinto(), 1);
123  Fp= Fp (eval[0].mapinto(), 1);
124  CanonicalForm smallestFactorEval= smallestFactor (eval[0].mapinto(),1);
125  setCharacteristic (0);
126  CanonicalForm F1= F(eval[0],1);
127  int s= tdegF/minTdeg + 1;
128  CanonicalForm bound= power (maxNorm (F1), 2*(s-1));
129  bound *= power (CanonicalForm (s),s-1);
130  bound *= power (CanonicalForm (2), ((s-1)*(s-1))/2); //possible int overflow
131
132  CanonicalForm B = p;
133  long k = 1L;
134  while ( B < bound ) {
135    B *= p;
136    k++;
137  }
138
139  //TODO take floor (log2(k))
140  k= k+1;
141  fmpz_poly_t FLINTF1;
142  convertFacCF2Fmpz_poly_t (FLINTF1, F1);
143  setCharacteristic (p);
144  nmod_poly_t FLINTFp, FLINTGp;
145  convertFacCF2nmod_poly_t (FLINTFp, smallestFactorEval/lc (smallestFactorEval));
146  convertFacCF2nmod_poly_t (FLINTGp, Gp/lc (Gp));
147  nmod_poly_factor_t nmodFactors;
148  nmod_poly_factor_init (nmodFactors);
149  nmod_poly_factor_insert (nmodFactors, FLINTFp, 1L);
150  nmod_poly_factor_insert (nmodFactors, FLINTGp, 1L);
151
152  long * link= new long [2];
153  fmpz_poly_t *v= new fmpz_poly_t[2];
154  fmpz_poly_t *w= new fmpz_poly_t[2];
155  fmpz_poly_init(v[0]);
156  fmpz_poly_init(v[1]);
157  fmpz_poly_init(w[0]);
158  fmpz_poly_init(w[1]);
159
160  fmpz_poly_factor_t liftedFactors;
161  fmpz_poly_factor_init (liftedFactors);
162  _fmpz_poly_hensel_start_lift(liftedFactors, link, v, w, FLINTF1, nmodFactors, k); //lift factors up to p^k
163
164  nmod_poly_factor_clear (nmodFactors);
165  nmod_poly_clear (FLINTFp);
166  nmod_poly_clear (FLINTGp);
167
168  setCharacteristic(0);
169  modpk pk= modpk (p,k);
170  CanonicalForm liftedSmallestFactor= convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[0],Variable (2));
171
172  CanonicalForm otherFactor= convertFmpz_poly_t2FacCF ((fmpz_poly_t &)liftedFactors->p[1],Variable (2));
173  CanonicalForm test= pk (otherFactor*liftedSmallestFactor);
174
175  Off (SW_SYMMETRIC_FF);
176  liftedSmallestFactor= pk (liftedSmallestFactor);
177  liftedSmallestFactor= liftedSmallestFactor (eval[1],2);
178  On (SW_SYMMETRIC_FF);
179  CFMatrix M= CFMatrix (s, s);
180  M(s,s)= power (CanonicalForm (p), k);
181  for (int i= 1; i < s; i++)
182  {
183    M (i,i)= 1;
184    M (i+1,i)= -liftedSmallestFactor;
185  }
186
187  mat_ZZ NTLM= *convertFacCFMatrix2NTLmat_ZZ (M);
188
189  ZZ det;
190
191  transpose (NTLM, NTLM);
192  long r=LLL (det, NTLM, 1L, 1L); //use floating point LLL ?
193  transpose (NTLM, NTLM);
194  M= *convertNTLmat_ZZ2FacCFMatrix (NTLM);
195
196  CanonicalForm mipo= 0;
197  Variable x= Variable (1);
198  for (int i= 1; i <= s; i++)
199    mipo += M (i,1)*power (x,s-i);
200
201  CFFList mipoFactors= factorize (mipo);
202  mipoFactors.removeFirst();
203  //TODO check if mipoFactors has length 1 and multiplicity 1 - if not choose a new point!
204  On (SW_RATIONAL);
205  Variable alpha= rootOf (mipo);
206  CFFList QaFactors= factorize (F1, alpha);
207
208  QaFactors.append (CFFactor (mipo, 1)); //last factor is the minimal polynomial that defines the extension
209  fmpz_poly_clear (v[0]);
210  fmpz_poly_clear (v[1]);
211  fmpz_poly_clear (w[0]);
212  fmpz_poly_clear (w[1]);
213  delete [] v;
214  delete [] w;
215  delete [] link;
216  fmpz_poly_factor_clear (liftedFactors);
217  return QaFactors;
218}
219#endif
220#endif
221
222
Note: See TracBrowser for help on using the repository browser.