source: git/factory/facSparseHensel.cc @ aba90e8

spielwiese
Last change on this file since aba90e8 was ee4a7d, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: changed behaviour of Lucks Wang sparse heuristic
  • Property mode set to 100644
File size: 10.6 KB
Line 
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#include "config.h"
17#include "cf_assert.h"
18#include "facSparseHensel.h"
19#include "cf_algorithm.h"
20#include "cf_gcd_smallp.h"
21#include "facFqFactorize.h"
22
23int
24LucksWangSparseHeuristic (const CanonicalForm& F, const CFList& factors,
25                          int level, const CFList& leadingCoeffs, CFList& result)
26{
27  CFArray* monoms= new CFArray [factors.length()];
28  int i= 0;
29  int num= 0;
30  for (CFListIterator iter= factors; iter.hasItem(); iter++, i++)
31  {
32    monoms[i]= getTerms (iter.getItem());
33    num += monoms [i].size();
34    sort (monoms [i]);
35  }
36
37  i= 0;
38  CFArray* monomsLead= new CFArray [leadingCoeffs.length()];
39  for (CFListIterator iter= leadingCoeffs; iter.hasItem(); iter++, i++)
40  {
41    monomsLead[i]= getTerms (iter.getItem());
42    sort (monomsLead [i]);
43    groupTogether (monomsLead [i], level);
44    strip (monomsLead [i], level);
45  }
46
47  CFArray solution= CFArray (num);
48  int k, d, count, j= F.level() + 1;
49  num= 0;
50  i= 0;
51  for (CFListIterator iter= factors; iter.hasItem(); i++, iter++)
52  {
53    d= degree (iter.getItem(), 1);
54    count= 0;
55    for (k= 0; k < monoms[i].size(); k++, j++, num++)
56    {
57      monoms [i][k] *= Variable (j);
58      if (degree (monoms[i][k], 1) == d)
59      {
60        solution[num]= monomsLead [i] [count];
61        count++;
62      }
63    }
64  }
65
66  delete [] monomsLead;
67
68  CFArray termsF= getBiTerms (F);
69  if (termsF.size() > 450)
70  {
71    delete [] monoms;
72    return 0;
73  }
74  sort (termsF, level);
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  }
133  CanonicalForm factor;
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
154CFList
155sparseHeuristic (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}
Note: See TracBrowser for help on using the repository browser.