My Project
facSparseHensel.cc
Go to the documentation of this file.
1 /*****************************************************************************\
2  * Computer Algebra System SINGULAR
3 \*****************************************************************************/
4 /** @file facSparseHensel.cc
5  *
6  * This file implements functions for sparse heuristic Hensel lifting
7  *
8  * ABSTRACT: "A fast implementation of polynomial factorization" by M. Lucks and
9  * "Effective polynomial computation" by R. Zippel
10  *
11  * @author Martin Lee
12  *
13  **/
14 /*****************************************************************************/
15 
16 
17 #include "config.h"
18 
19 #include "cf_assert.h"
20 #include "facSparseHensel.h"
21 #include "cf_algorithm.h"
22 #include "cfModGcd.h"
23 #include "facFqFactorize.h"
24 
25 int
26 LucksWangSparseHeuristic (const CanonicalForm& F, const CFList& factors,
27  int level, const CFList& leadingCoeffs, CFList& result)
28 {
29  int threshold= 450;
30  CFArray termsF= getBiTerms (F, threshold);
31  int si=termsF.size();
32  int fl=factors.length();
33  //printf("size:%d, length=%d, char=%d\n",si,fl,getCharacteristic());
34  if ((si > threshold)
35  || (si <= fl)
36  )
37  return 0;
38  sort (termsF, level);
39 
40  CFArray* monoms= new CFArray [factors.length()];
41  int i= 0;
42  int num= 0;
43  for (CFListIterator iter= factors; iter.hasItem(); iter++, i++)
44  {
45  monoms[i]= getTerms (iter.getItem());
46  num += monoms [i].size();
47  sort (monoms [i]);
48  }
49 
50  i= 0;
51  CFArray* monomsLead= new CFArray [leadingCoeffs.length()];
52  for (CFListIterator iter= leadingCoeffs; iter.hasItem(); iter++, i++)
53  {
54  monomsLead[i]= getTerms (iter.getItem());
55  sort (monomsLead [i]);
56  groupTogether (monomsLead [i], level);
57  strip (monomsLead [i], level);
58  }
59 
60  CFArray solution= CFArray (num);
61  int k, d, count, j= F.level() + 1;
62  num= 0;
63  i= 0;
64  for (CFListIterator iter= factors; iter.hasItem(); i++, iter++)
65  {
66  d= degree (iter.getItem(), 1);
67  count= 0;
68  for (k= 0; k < monoms[i].size(); k++, j++, num++)
69  {
70  monoms [i][k] *= Variable (j);
71  if (degree (monoms[i][k], 1) == d)
72  {
73  solution[num]= monomsLead [i] [count];
74  count++;
75  }
76  }
77  }
78 
79  delete [] monomsLead;
80 
81  CFList tmp;
82  CFArray* stripped2= new CFArray [factors.length()];
83  for (i= factors.length() - 1; i > -1; i--)
84  {
85  tmp.insert (buildPolyFromArray (monoms [i]));
86  strip (monoms[i], stripped2 [i], level);
87  }
88  delete [] monoms;
89 
90  CanonicalForm H= prod (tmp);
91  CFArray monomsH= getMonoms (H);
92  sort (monomsH,F.level());
93 
94  groupTogether (monomsH, F.level());
95 
96  if (monomsH.size() != termsF.size())
97  {
98  delete [] stripped2;
99  return 0;
100  }
101 
102  CFArray strippedH;
103  strip (monomsH, strippedH, level);
104  CFArray strippedF;
105  strip (termsF, strippedF, level);
106 
107  if (!isEqual (strippedH, strippedF))
108  {
109  delete [] stripped2;
110  return 0;
111  }
112 
113  CFArray A= getEquations (monomsH, termsF);
114  CFArray startingSolution= solution;
115  CFArray newSolution= CFArray (solution.size());
116  result= CFList();
117  do
118  {
119  evaluate (A, solution, F.level() + 1);
120  if (isZero (A))
121  break;
122  if (!simplify (A, newSolution, F.level() + 1))
123  {
124  delete [] stripped2;
125  return 0;
126  }
127  if (isZero (newSolution))
128  break;
129  if (!merge (solution, newSolution))
130  break;
131  } while (1);
132 
133  if (isEqual (startingSolution, solution))
134  {
135  delete [] stripped2;
136  return 0;
137  }
139  num= 0;
140  for (i= 0; i < factors.length(); i++)
141  {
142  k= stripped2[i].size();
143  factor= 0;
144  for (j= 0; j < k; j++, num++)
145  {
146  if (solution [num].isZero())
147  continue;
148  factor += solution [num]*stripped2[i][j];
149  }
150  result.append (factor);
151  }
152 
153  delete [] stripped2;
154  if (result.length() > 0)
155  return 1;
156  return 0;
157 }
158 
159 CFList
160 sparseHeuristic (const CanonicalForm& A, const CFList& biFactors,
161  CFList*& moreBiFactors, const CFList& evaluation,
162  int minFactorsLength)
163 {
164  int j= A.level() - 1;
165  int i;
166 
167  //initialize storage
168  CFArray *** storeFactors= new CFArray** [j];
169  for (i= 0; i < j; i++)
170  storeFactors [i]= new CFArray* [2];
171 
172  CFArray eval= CFArray (j);
173  i= j - 1;
175  eval[i]= iter.getItem();
176  storeFactors [0] [0]= new CFArray [minFactorsLength];
177  storeFactors [0] [1]= new CFArray [minFactorsLength];
178  for (i= 1; i < j; i++)
179  {
180  storeFactors[i] [0]= new CFArray [minFactorsLength];
181  storeFactors[i] [1]= new CFArray [minFactorsLength];
182  }
183  //
184 
185  CFList * normalizingFactors= new CFList [j];
186  CFList uniFactors;
187  normalizingFactors [0]= findNormalizingFactor1 (biFactors,
188  evaluation.getLast(), uniFactors);
189  for (i= j - 1; i > 0; i--)
190  {
191  if (moreBiFactors[i-1].length() != minFactorsLength)
192  {
193  moreBiFactors[i-1]=
194  recombination (moreBiFactors [i-1], uniFactors, 1,
195  moreBiFactors[i-1].length()-uniFactors.length()+1,
196  eval[i], Variable (i + 2)
197  );
198  }
199  normalizingFactors [i]= findNormalizingFactor2 (moreBiFactors [i - 1],
200  eval[i], uniFactors);
201  }
202 
203  CFList tmp;
204  tmp= normalize (biFactors, normalizingFactors[0]);
205  getTerms2 (tmp, storeFactors [0] [0]);
206  storeFactors [0] [1]= evaluate (storeFactors [0] [0], minFactorsLength,
207  evaluation.getLast(), Variable (2));
208  for (i= j - 1; i > 0; i--)
209  {
210  tmp= normalize (moreBiFactors [i-1], normalizingFactors [i]);
211  getTerms2 (tmp, storeFactors [i] [0]);
212  storeFactors [i] [1]= evaluate (storeFactors [i] [0], minFactorsLength,
213  eval[i], Variable (i + 2));
214  }
215 
216 
217  int k, l, m, mm, count, sizeOfUniFactors= 0;
218  int*** seperator= new int** [j];
219  Variable x= Variable (1);
220 
221  for (i= 0; i < j; i++)
222  seperator [i]= new int* [minFactorsLength];
223  for (k= 0; k < minFactorsLength; k++)
224  {
225  for (i= 0; i < j; i++)
226  {
227  count= 0;
228  for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
229  {
230  if (degree (storeFactors[i][0][k][l], x) <
231  degree (storeFactors[i][0][k][l+1], x))
232  count++;
233  }
234  if (i == 0)
235  sizeOfUniFactors= count;
236  else if (sizeOfUniFactors != count)
237  {
238  for (m= 0; m < j; m++)
239  {
240  delete [] storeFactors [m] [0];
241  delete [] storeFactors [m] [1];
242  delete [] storeFactors [m];
243  for (mm= 0; mm < k; mm++)
244  delete [] seperator [m][mm];
245  delete [] seperator [m];
246  }
247  delete [] storeFactors;
248  delete [] seperator;
249  delete [] normalizingFactors;
250  return CFList();
251  }
252  seperator [i][k]= new int [count + 3];
253  seperator [i][k][0]= count + 1;
254  seperator [i][k][1]= 0;
255  count= 2;
256  for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
257  {
258  if (degree (storeFactors[i][0][k][l], x) <
259  degree (storeFactors[i][0][k][l+1], x))
260  {
261  seperator[i][k][count]=l + 1;
262  count++;
263  }
264  }
265  seperator [i][k][count]= storeFactors[i][0][k].size();
266  }
267  }
268 
269  CanonicalForm tmp1, factor, quot;
270  CanonicalForm B= A;
271  CFList result;
272  int maxTerms, n, index1, index2, mmm, found, columns, oneCount;
273  int ** mat;
274 
275  for (k= 0; k < minFactorsLength; k++)
276  {
277  factor= 0;
278  sizeOfUniFactors= seperator [0][k][0];
279  for (n= 1; n <= sizeOfUniFactors; n++)
280  {
281  columns= 0;
282  maxTerms= 1;
283  index1= j - 1;
284  for (i= j - 1; i >= 0; i--)
285  {
286  if (maxTerms < seperator[i][k][n+1]-seperator[i][k][n])
287  {
288  maxTerms= seperator[i][k][n + 1]-seperator[i][k][n];
289  index1= i;
290  }
291  }
292  for (i= j - 1; i >= 0; i--)
293  {
294  if (i == index1)
295  continue;
296  columns += seperator [i][k][n+1]-seperator[i][k][n];
297  }
298  mat= new int *[maxTerms];
299  mm= 0;
300  for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
301  {
302  tmp1= storeFactors [index1][1][k][m];
303  mat[mm]= new int [columns];
304  for (i= 0; i < columns; i++)
305  mat[mm][i]= 0;
306  index2= 0;
307  for (i= j - 1; i >= 0; i--)
308  {
309  if (i == index1)
310  continue;
311  found= -1;
312  if ((found= search (storeFactors[i][1][k], tmp1,
313  seperator[i][k][n], seperator[i][k][n+1])) >= 0)
314  mat[mm][index2 + found - seperator [i][k][n]]= 1;
315  index2 += seperator [i][k][n+1]-seperator[i][k][n];
316  }
317  }
318 
319  index2= 0;
320  for (i= j - 1; i >= 0; i--)
321  {
322  if (i == index1)
323  continue;
324  oneCount= 0;
325  for (mm= 0; mm < seperator [i][k][n + 1] - seperator [i][k][n]; mm++)
326  {
327  for (m= 0; m < maxTerms; m++)
328  {
329  if (mat[m][mm+index2] == 1)
330  oneCount++;
331  }
332  }
333  if (oneCount == seperator [i][k][n+1]-seperator[i][k][n] - 1)
334  {
335  for (mm= 0; mm < seperator [i][k][n+1]-seperator[i][k][n]; mm++)
336  {
337  oneCount= 0;
338  for (m= 0; m < maxTerms; m++)
339  if (mat[m][mm+index2] == 1)
340  oneCount++;
341  if (oneCount > 0)
342  continue;
343  for (m= 0; m < maxTerms; m++)
344  {
345  oneCount= 0;
346  for (mmm= 0; mmm < seperator[i][k][n+1]-seperator[i][k][n]; mmm++)
347  {
348  if (mat[m][mmm+index2] == 1)
349  oneCount++;
350  }
351  if (oneCount > 0)
352  continue;
353  mat[m][mm+index2]= 1;
354  }
355  }
356  }
357  index2 += seperator [i][k][n+1] - seperator [i][k][n];
358  }
359 
360  //read off solution
361  mm= 0;
362  for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
363  {
364  tmp1= storeFactors [index1][0][k][m];
365  index2= 0;
366  for (i= j - 1; i > -1; i--)
367  {
368  if (i == index1)
369  continue;
370  for (mmm= 0; mmm < seperator [i][k][n+1]-seperator[i][k][n]; mmm++)
371  if (mat[mm][mmm+index2] == 1)
372  tmp1= patch (tmp1, storeFactors[i][0][k][seperator[i][k][n]+mmm],
373  eval[i]);
374  index2 += seperator [i][k][n+1]-seperator[i][k][n];
375  }
376  factor += tmp1;
377  }
378 
379  for (m= 0; m < maxTerms; m++)
380  delete [] mat [m];
381  delete [] mat;
382  }
383 
384  if (fdivides (factor, B, quot))
385  {
386  result.append (factor);
387  B= quot;
388  if (result.length() == biFactors.length() - 1)
389  {
390  result.append (quot);
391  break;
392  }
393  }
394  }
395 
396  //delete
397  for (i= 0; i < j; i++)
398  {
399  delete [] storeFactors [i] [0];
400  delete [] storeFactors [i] [1];
401  delete [] storeFactors [i];
402  for (k= 0; k < minFactorsLength; k++)
403  delete [] seperator [i][k];
404  delete [] seperator [i];
405  }
406  delete [] seperator;
407  delete [] storeFactors;
408  delete [] normalizingFactors;
409  //
410 
411  return result;
412 }
int degree(const CanonicalForm &f)
Array< CanonicalForm > CFArray
CanonicalForm num(const CanonicalForm &f)
int level(const CanonicalForm &f)
List< CanonicalForm > CFList
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
bool isEqual(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:948
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition: cfModGcd.cc:1959
Variable x
Definition: cfModGcd.cc:4084
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition: cfModGcd.cc:2033
modular and sparse modular GCD algorithms over finite fields and Z.
int ** merge(int **points1, int sizePoints1, int **points2, int sizePoints2, int &sizeResult)
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
declarations of higher level algorithms.
void getTerms(const CanonicalForm &f, const CanonicalForm &t, CFList &result)
get_Terms: Split the polynomial in the containing terms.
Definition: cf_factor.cc:279
assertions for Factory
int size() const
Definition: ftmpl_array.cc:92
factory's main class
Definition: canonicalform.h:86
int level() const
level() returns the level of CO.
T & getItem() const
Definition: ftmpl_list.cc:431
int length() const
Definition: ftmpl_list.cc:273
void insert(const T &)
Definition: ftmpl_list.cc:193
factory's class for variables
Definition: factory.h:134
CFFListIterator iter
Definition: facAbsBiFact.cc:53
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int const CFList & evaluation
Definition: facAbsFact.cc:52
CanonicalForm factor
Definition: facAbsFact.cc:97
CanonicalForm H
Definition: facAbsFact.cc:60
b *CanonicalForm B
Definition: facBivar.cc:52
bool found
Definition: facFactorize.cc:55
CFList & eval
Definition: facFactorize.cc:47
CFList int & minFactorsLength
[in,out] minimal length of bivariate factors
Definition: facFactorize.h:33
CFList tmp1
Definition: facFqBivar.cc:72
CFList recombination(const CFList &factors1, const CFList &factors2, int s, int thres, const CanonicalForm &evalPoint, const Variable &x)
recombination of bivariate factors factors1 s. t. the result evaluated at evalPoint coincides with fa...
This file provides functions for factorizing a multivariate polynomial over , or GF.
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
int LucksWangSparseHeuristic(const CanonicalForm &F, const CFList &factors, int level, const CFList &leadingCoeffs, CFList &result)
sparse heuristic lifting by Wang and Lucks
CFList sparseHeuristic(const CanonicalForm &A, const CFList &biFactors, CFList *&moreBiFactors, const CFList &evaluation, int minFactorsLength)
sparse heuristic which patches together bivariate factors of A wrt. different second variables by the...
This file provides functions for sparse heuristic Hensel lifting.
CanonicalForm buildPolyFromArray(const CFArray &A)
build a poly from entries in A
void groupTogether(CFArray &A, int level)
group together elements in A, where entries in A are put together if they coincide up to level level
CFArray getTerms2(const CanonicalForm &F)
get terms of F wrt. Variable (1)
CFArray getBiTerms(const CanonicalForm &F, int threshold)
get terms of F where F is considered a bivariate poly in Variable(1), Variable (2)
CanonicalForm simplify(const CanonicalForm &A, int level)
simplify A if possible, i.e. A consists of 2 terms and contains only one variable of level greater or...
bool isZero(const CFArray &A)
checks if entries of A are zero
CFList findNormalizingFactor1(const CFList &biFactors, const CanonicalForm &evalPoint, CFList &uniFactors)
find normalizing factors for biFactors and build monic univariate factors from biFactors
CFList findNormalizingFactor2(CFList &biFactors, const CanonicalForm &evalPoint, const CFList &uniFactors)
find normalizing factors for biFactors and sort biFactors s.t. the returned biFactors evaluated at ev...
int search(const CFArray &A, const CanonicalForm &F, int i, int j)
search for F in A between index i and j
CanonicalForm patch(const CanonicalForm &F1, const CanonicalForm &F2, const CanonicalForm &eval)
patch together F1 and F2 and normalize by a power of eval F1 and F2 are assumed to be bivariate with ...
void strip(CFArray &F, CFArray &G, int level)
strip off those parts of entries in F whose level is less than or equal than level and store the stri...
void sort(CFArray &A, int l=0)
quick sort A
CFArray getEquations(const CFArray &A, const CFArray &B)
get equations for LucksWangSparseHeuristic
static BOOLEAN length(leftv result, leftv arg)
Definition: interval.cc:257
int status int void size_t count
Definition: si_signals.h:59
#define A
Definition: sirandom.c:24
static poly normalize(poly next_p, ideal add_generators, syStrategy syzstr, int *g_l, int *p_l, int crit_comp)
Definition: syz3.cc:1026