source: git/factory/facSparseHensel.cc @ fce807

spielwiese
Last change on this file since fce807 was 7b5cb2, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: faster retrieval of terms chg: LucksWangSparseHeuristic now returns also partial factors
  • 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  do
112  {
113    evaluate (A, solution, F.level() + 1);
114    if (isZero (A))
115      break;
116    if (!simplify (A, newSolution, F.level() + 1))
117      break;
118    if (isZero (newSolution))
119      break;
120    if (!merge (solution, newSolution))
121      break;
122  } while (1);
123
124
125  result= CFList();
126  if (isEqual (startingSolution, solution))
127  {
128    delete [] stripped2;
129    return 0;
130  }
131  CanonicalForm factor;
132  num= 0;
133  for (i= 0; i < factors.length(); i++)
134  {
135    k= stripped2[i].size();
136    factor= 0;
137    for (j= 0; j < k; j++, num++)
138    {
139      if (solution [num].isZero())
140        continue;
141      factor += solution [num]*stripped2[i][j];
142    }
143    result.append (factor);
144  }
145
146  delete [] stripped2;
147  if (result.length() > 0)
148    return 1;
149  return 0;
150}
151
152CFList
153sparseHeuristic (const CanonicalForm& A, const CFList& biFactors,
154                 CFList*& moreBiFactors, const CFList& evaluation,
155                 int minFactorsLength)
156{
157  int j= A.level() - 1;
158  int i;
159
160  //initialize storage
161  CFArray *** storeFactors= new CFArray** [j];
162  for (i= 0; i < j; i++)
163    storeFactors [i]= new CFArray* [2];
164
165  CFArray eval= CFArray (j);
166  i= j - 1;
167  for (CFListIterator iter= evaluation; iter.hasItem(); iter++,i--)
168    eval[i]= iter.getItem();
169  storeFactors [0] [0]= new CFArray [minFactorsLength];
170  storeFactors [0] [1]= new CFArray [minFactorsLength];
171  for (i= 1; i < j; i++)
172  {
173    storeFactors[i] [0]= new CFArray [minFactorsLength];
174    storeFactors[i] [1]= new CFArray [minFactorsLength];
175  }
176  //
177
178  CFList * normalizingFactors= new CFList [j];
179  CFList uniFactors;
180  normalizingFactors [0]= findNormalizingFactor1 (biFactors,
181                                              evaluation.getLast(), uniFactors);
182  for (i= j - 1; i > 0; i--)
183  {
184    if (moreBiFactors[i-1].length() != minFactorsLength)
185    {
186      moreBiFactors[i-1]= 
187        recombination (moreBiFactors [i-1], uniFactors, 1,
188                       moreBiFactors[i-1].length()-uniFactors.length()+1,
189                       eval[i], Variable (i + 2)
190                      );
191    }
192    normalizingFactors [i]= findNormalizingFactor2 (moreBiFactors [i - 1],
193                                                    eval[i], uniFactors);
194  }
195
196  CFList tmp;
197  tmp= normalize (biFactors, normalizingFactors[0]);
198  getTerms2 (tmp, storeFactors [0] [0]);
199  storeFactors [0] [1]= evaluate (storeFactors [0] [0], minFactorsLength,
200                                  evaluation.getLast(), Variable (2));
201  for (i= j - 1; i > 0; i--)
202  {
203    tmp= normalize (moreBiFactors [i-1], normalizingFactors [i]);
204    getTerms2 (tmp, storeFactors [i] [0]);
205    storeFactors [i] [1]= evaluate (storeFactors [i] [0], minFactorsLength,
206                                    eval[i], Variable (i + 2));
207  }
208
209
210  int k, l, m, mm, count, sizeOfUniFactors= 0;
211  int*** seperator= new int** [j];
212  Variable x= Variable (1);
213
214  for (i= 0; i < j; i++)
215    seperator [i]= new int* [minFactorsLength];
216  for (k= 0; k < minFactorsLength; k++)
217  {
218    for (i= 0; i < j; i++)
219    {
220      count= 0;
221      for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
222      {
223        if (degree (storeFactors[i][0][k][l], x) <
224            degree (storeFactors[i][0][k][l+1], x))
225          count++;
226      }
227      if (i == 0)
228        sizeOfUniFactors= count;
229      else if (sizeOfUniFactors != count)
230      {
231        for (m= 0; m < j; m++)
232        {
233          delete [] storeFactors [m] [0];
234          delete [] storeFactors [m] [1];
235          delete [] storeFactors [m];
236          for (mm= 0; mm < k; mm++)
237            delete [] seperator [m][mm];
238          delete [] seperator [m];
239        }
240        delete [] storeFactors;
241        delete [] seperator;
242        return CFList();
243      }
244      seperator [i][k]= new int [count + 3];
245      seperator [i][k][0]= count + 1;
246      seperator [i][k][1]= 0;
247      count= 2;
248      for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
249      {
250        if (degree (storeFactors[i][0][k][l], x) <
251            degree (storeFactors[i][0][k][l+1], x))
252        {
253          seperator[i][k][count]=l + 1;
254          count++;
255        }
256      }
257      seperator [i][k][count]= storeFactors[i][0][k].size();
258    }
259  }
260
261  CanonicalForm tmp1, factor, quot;
262  CanonicalForm B= A;
263  CFList result;
264  int maxTerms, n, index1, index2, mmm, found, columns, oneCount;
265  int ** mat;
266
267  for (k= 0; k < minFactorsLength; k++)
268  {
269    factor= 0;
270    sizeOfUniFactors= seperator [0][k][0];
271    for (n= 1; n <= sizeOfUniFactors; n++)
272    {
273      columns= 0;
274      maxTerms= 1;
275      index1= j - 1;
276      for (i= j - 1; i >= 0; i--)
277      {
278        if (maxTerms < seperator[i][k][n+1]-seperator[i][k][n])
279        {
280          maxTerms= seperator[i][k][n + 1]-seperator[i][k][n];
281          index1= i;
282        }
283      }
284      for (i= j - 1; i >= 0; i--)
285      {
286        if (i == index1)
287          continue;
288        columns += seperator [i][k][n+1]-seperator[i][k][n];
289      }
290      mat= new int *[maxTerms];
291      mm= 0;
292      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
293      {
294        tmp1= storeFactors [index1][1][k][m];
295        mat[mm]= new int [columns];
296        for (i= 0; i < columns; i++)
297          mat[mm][i]= 0;
298        index2= 0;
299        for (i= j - 1; i >= 0; i--)
300        {
301          if (i == index1)
302            continue;
303          found= -1;
304          if ((found= search (storeFactors[i][1][k], tmp1, 
305                              seperator[i][k][n], seperator[i][k][n+1])) >= 0)
306            mat[mm][index2 + found - seperator [i][k][n]]= 1;
307          index2 += seperator [i][k][n+1]-seperator[i][k][n];
308        }
309      }
310
311      index2= 0;
312      for (i= j - 1; i >= 0; i--)
313      {
314        if (i == index1)
315          continue;
316        oneCount= 0;
317        for (mm= 0; mm < seperator [i][k][n + 1] - seperator [i][k][n]; mm++)
318        {
319          for (m= 0; m < maxTerms; m++)
320          {
321            if (mat[m][mm+index2] == 1)
322              oneCount++;
323          }
324        }
325        if (oneCount == seperator [i][k][n+1]-seperator[i][k][n] - 1)
326        {
327          for (mm= 0; mm < seperator [i][k][n+1]-seperator[i][k][n]; mm++)
328          {
329            oneCount= 0;
330            for (m= 0; m < maxTerms; m++)
331              if (mat[m][mm+index2] == 1)
332                oneCount++;
333            if (oneCount > 0)
334              continue;
335            for (m= 0; m < maxTerms; m++)
336            {
337              oneCount= 0;
338              for (mmm= 0; mmm < seperator[i][k][n+1]-seperator[i][k][n]; mmm++)
339              {
340                if (mat[m][mmm+index2] == 1)
341                  oneCount++;
342              }
343              if (oneCount > 0)
344                continue;
345              mat[m][mm+index2]= 1;
346            }
347          }
348        }
349        index2 += seperator [i][k][n+1] - seperator [i][k][n];
350      }
351
352      //read off solution
353      mm= 0;
354      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
355      {
356        tmp1= storeFactors [index1][0][k][m];
357        index2= 0;
358        for (i= j - 1; i > -1; i--)
359        {
360          if (i == index1)
361            continue;
362          for (mmm= 0; mmm < seperator [i][k][n+1]-seperator[i][k][n]; mmm++)
363            if (mat[mm][mmm+index2] == 1)
364              tmp1= patch (tmp1, storeFactors[i][0][k][seperator[i][k][n]+mmm],
365                           eval[i]);
366          index2 += seperator [i][k][n+1]-seperator[i][k][n];
367        }
368        factor += tmp1;
369      }
370
371      for (m= 0; m < maxTerms; m++)
372        delete [] mat [m];
373      delete [] mat;
374    }
375
376    if (fdivides (factor, B, quot))
377    {
378      result.append (factor);
379      B= quot;
380      if (result.length() == biFactors.length() - 1)
381      {
382        result.append (quot);
383        break;
384      }
385    }
386  }
387
388  //delete
389  for (i= 0; i < j; i++)
390  {
391    delete [] storeFactors [i] [0];
392    delete [] storeFactors [i] [1];
393    delete [] storeFactors [i];
394    for (k= 0; k < minFactorsLength; k++)
395      delete [] seperator [i][k];
396    delete [] seperator [i];
397  }
398  delete [] seperator;
399  delete [] storeFactors;
400  //
401
402  return result;
403}
Note: See TracBrowser for help on using the repository browser.