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