source: git/factory/facSparseHensel.cc @ 16f511

spielwiese
Last change on this file since 16f511 was 16f511, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fixed the usage of "config.h" (if defined HAVE_CONFIG_H)
  • Property mode set to 100644
File size: 10.7 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#ifdef HAVE_CONFIG_H
17#include "config.h"
18#endif /* HAVE_CONFIG_H */
19#include "cf_assert.h"
20#include "facSparseHensel.h"
21#include "cf_algorithm.h"
22#include "cf_gcd_smallp.h"
23#include "facFqFactorize.h"
24
25int
26LucksWangSparseHeuristic (const CanonicalForm& F, const CFList& factors,
27                          int level, const CFList& leadingCoeffs, CFList& result)
28{
29  CFArray* monoms= new CFArray [factors.length()];
30  int i= 0;
31  int num= 0;
32  for (CFListIterator iter= factors; iter.hasItem(); iter++, i++)
33  {
34    monoms[i]= getTerms (iter.getItem());
35    num += monoms [i].size();
36    sort (monoms [i]);
37  }
38
39  i= 0;
40  CFArray* monomsLead= new CFArray [leadingCoeffs.length()];
41  for (CFListIterator iter= leadingCoeffs; iter.hasItem(); iter++, i++)
42  {
43    monomsLead[i]= getTerms (iter.getItem());
44    sort (monomsLead [i]);
45    groupTogether (monomsLead [i], level);
46    strip (monomsLead [i], level);
47  }
48
49  CFArray solution= CFArray (num);
50  int k, d, count, j= F.level() + 1;
51  num= 0;
52  i= 0;
53  for (CFListIterator iter= factors; iter.hasItem(); i++, iter++)
54  {
55    d= degree (iter.getItem(), 1);
56    count= 0;
57    for (k= 0; k < monoms[i].size(); k++, j++, num++)
58    {
59      monoms [i][k] *= Variable (j);
60      if (degree (monoms[i][k], 1) == d)
61      {
62        solution[num]= monomsLead [i] [count];
63        count++;
64      }
65    }
66  }
67
68  delete [] monomsLead;
69
70  CFArray termsF= getBiTerms (F);
71  if (termsF.size() > 450)
72  {
73    delete [] monoms;
74    return 0;
75  }
76  sort (termsF, level);
77
78  CFList tmp;
79  CFArray* stripped2= new CFArray [factors.length()];
80  for (i= factors.length() - 1; i > -1; i--)
81  {
82    tmp.insert (buildPolyFromArray (monoms [i]));
83    strip (monoms[i], stripped2 [i], level);
84  }
85  delete [] monoms;
86
87  CanonicalForm H= prod (tmp);
88  CFArray monomsH= getMonoms (H);
89  sort (monomsH,F.level());
90
91  groupTogether (monomsH, F.level());
92
93  if (monomsH.size() != termsF.size())
94  {
95    delete [] stripped2;
96    return 0;
97  }
98
99  CFArray strippedH;
100  strip (monomsH, strippedH, level);
101  CFArray strippedF;
102  strip (termsF, strippedF, level);
103
104  if (!isEqual (strippedH, strippedF))
105  {
106    delete [] stripped2;
107    return 0;
108  }
109
110  CFArray A= getEquations (monomsH, termsF);
111  CFArray startingSolution= solution;
112  CFArray newSolution= CFArray (solution.size());
113  result= CFList();
114  do
115  {
116    evaluate (A, solution, F.level() + 1);
117    if (isZero (A))
118      break;
119    if (!simplify (A, newSolution, F.level() + 1))
120    {
121      delete [] stripped2;
122      return 0;
123    }
124    if (isZero (newSolution))
125      break;
126    if (!merge (solution, newSolution))
127      break;
128  } while (1);
129
130  if (isEqual (startingSolution, solution))
131  {
132    delete [] stripped2;
133    return 0;
134  }
135  CanonicalForm factor;
136  num= 0;
137  for (i= 0; i < factors.length(); i++)
138  {
139    k= stripped2[i].size();
140    factor= 0;
141    for (j= 0; j < k; j++, num++)
142    {
143      if (solution [num].isZero())
144        continue;
145      factor += solution [num]*stripped2[i][j];
146    }
147    result.append (factor);
148  }
149
150  delete [] stripped2;
151  if (result.length() > 0)
152    return 1;
153  return 0;
154}
155
156CFList
157sparseHeuristic (const CanonicalForm& A, const CFList& biFactors,
158                 CFList*& moreBiFactors, const CFList& evaluation,
159                 int minFactorsLength)
160{
161  int j= A.level() - 1;
162  int i;
163
164  //initialize storage
165  CFArray *** storeFactors= new CFArray** [j];
166  for (i= 0; i < j; i++)
167    storeFactors [i]= new CFArray* [2];
168
169  CFArray eval= CFArray (j);
170  i= j - 1;
171  for (CFListIterator iter= evaluation; iter.hasItem(); iter++,i--)
172    eval[i]= iter.getItem();
173  storeFactors [0] [0]= new CFArray [minFactorsLength];
174  storeFactors [0] [1]= new CFArray [minFactorsLength];
175  for (i= 1; i < j; i++)
176  {
177    storeFactors[i] [0]= new CFArray [minFactorsLength];
178    storeFactors[i] [1]= new CFArray [minFactorsLength];
179  }
180  //
181
182  CFList * normalizingFactors= new CFList [j];
183  CFList uniFactors;
184  normalizingFactors [0]= findNormalizingFactor1 (biFactors,
185                                              evaluation.getLast(), uniFactors);
186  for (i= j - 1; i > 0; i--)
187  {
188    if (moreBiFactors[i-1].length() != minFactorsLength)
189    {
190      moreBiFactors[i-1]= 
191        recombination (moreBiFactors [i-1], uniFactors, 1,
192                       moreBiFactors[i-1].length()-uniFactors.length()+1,
193                       eval[i], Variable (i + 2)
194                      );
195    }
196    normalizingFactors [i]= findNormalizingFactor2 (moreBiFactors [i - 1],
197                                                    eval[i], uniFactors);
198  }
199
200  CFList tmp;
201  tmp= normalize (biFactors, normalizingFactors[0]);
202  getTerms2 (tmp, storeFactors [0] [0]);
203  storeFactors [0] [1]= evaluate (storeFactors [0] [0], minFactorsLength,
204                                  evaluation.getLast(), Variable (2));
205  for (i= j - 1; i > 0; i--)
206  {
207    tmp= normalize (moreBiFactors [i-1], normalizingFactors [i]);
208    getTerms2 (tmp, storeFactors [i] [0]);
209    storeFactors [i] [1]= evaluate (storeFactors [i] [0], minFactorsLength,
210                                    eval[i], Variable (i + 2));
211  }
212
213
214  int k, l, m, mm, count, sizeOfUniFactors= 0;
215  int*** seperator= new int** [j];
216  Variable x= Variable (1);
217
218  for (i= 0; i < j; i++)
219    seperator [i]= new int* [minFactorsLength];
220  for (k= 0; k < minFactorsLength; k++)
221  {
222    for (i= 0; i < j; i++)
223    {
224      count= 0;
225      for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
226      {
227        if (degree (storeFactors[i][0][k][l], x) <
228            degree (storeFactors[i][0][k][l+1], x))
229          count++;
230      }
231      if (i == 0)
232        sizeOfUniFactors= count;
233      else if (sizeOfUniFactors != count)
234      {
235        for (m= 0; m < j; m++)
236        {
237          delete [] storeFactors [m] [0];
238          delete [] storeFactors [m] [1];
239          delete [] storeFactors [m];
240          for (mm= 0; mm < k; mm++)
241            delete [] seperator [m][mm];
242          delete [] seperator [m];
243        }
244        delete [] storeFactors;
245        delete [] seperator;
246        return CFList();
247      }
248      seperator [i][k]= new int [count + 3];
249      seperator [i][k][0]= count + 1;
250      seperator [i][k][1]= 0;
251      count= 2;
252      for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
253      {
254        if (degree (storeFactors[i][0][k][l], x) <
255            degree (storeFactors[i][0][k][l+1], x))
256        {
257          seperator[i][k][count]=l + 1;
258          count++;
259        }
260      }
261      seperator [i][k][count]= storeFactors[i][0][k].size();
262    }
263  }
264
265  CanonicalForm tmp1, factor, quot;
266  CanonicalForm B= A;
267  CFList result;
268  int maxTerms, n, index1, index2, mmm, found, columns, oneCount;
269  int ** mat;
270
271  for (k= 0; k < minFactorsLength; k++)
272  {
273    factor= 0;
274    sizeOfUniFactors= seperator [0][k][0];
275    for (n= 1; n <= sizeOfUniFactors; n++)
276    {
277      columns= 0;
278      maxTerms= 1;
279      index1= j - 1;
280      for (i= j - 1; i >= 0; i--)
281      {
282        if (maxTerms < seperator[i][k][n+1]-seperator[i][k][n])
283        {
284          maxTerms= seperator[i][k][n + 1]-seperator[i][k][n];
285          index1= i;
286        }
287      }
288      for (i= j - 1; i >= 0; i--)
289      {
290        if (i == index1)
291          continue;
292        columns += seperator [i][k][n+1]-seperator[i][k][n];
293      }
294      mat= new int *[maxTerms];
295      mm= 0;
296      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
297      {
298        tmp1= storeFactors [index1][1][k][m];
299        mat[mm]= new int [columns];
300        for (i= 0; i < columns; i++)
301          mat[mm][i]= 0;
302        index2= 0;
303        for (i= j - 1; i >= 0; i--)
304        {
305          if (i == index1)
306            continue;
307          found= -1;
308          if ((found= search (storeFactors[i][1][k], tmp1, 
309                              seperator[i][k][n], seperator[i][k][n+1])) >= 0)
310            mat[mm][index2 + found - seperator [i][k][n]]= 1;
311          index2 += seperator [i][k][n+1]-seperator[i][k][n];
312        }
313      }
314
315      index2= 0;
316      for (i= j - 1; i >= 0; i--)
317      {
318        if (i == index1)
319          continue;
320        oneCount= 0;
321        for (mm= 0; mm < seperator [i][k][n + 1] - seperator [i][k][n]; mm++)
322        {
323          for (m= 0; m < maxTerms; m++)
324          {
325            if (mat[m][mm+index2] == 1)
326              oneCount++;
327          }
328        }
329        if (oneCount == seperator [i][k][n+1]-seperator[i][k][n] - 1)
330        {
331          for (mm= 0; mm < seperator [i][k][n+1]-seperator[i][k][n]; mm++)
332          {
333            oneCount= 0;
334            for (m= 0; m < maxTerms; m++)
335              if (mat[m][mm+index2] == 1)
336                oneCount++;
337            if (oneCount > 0)
338              continue;
339            for (m= 0; m < maxTerms; m++)
340            {
341              oneCount= 0;
342              for (mmm= 0; mmm < seperator[i][k][n+1]-seperator[i][k][n]; mmm++)
343              {
344                if (mat[m][mmm+index2] == 1)
345                  oneCount++;
346              }
347              if (oneCount > 0)
348                continue;
349              mat[m][mm+index2]= 1;
350            }
351          }
352        }
353        index2 += seperator [i][k][n+1] - seperator [i][k][n];
354      }
355
356      //read off solution
357      mm= 0;
358      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
359      {
360        tmp1= storeFactors [index1][0][k][m];
361        index2= 0;
362        for (i= j - 1; i > -1; i--)
363        {
364          if (i == index1)
365            continue;
366          for (mmm= 0; mmm < seperator [i][k][n+1]-seperator[i][k][n]; mmm++)
367            if (mat[mm][mmm+index2] == 1)
368              tmp1= patch (tmp1, storeFactors[i][0][k][seperator[i][k][n]+mmm],
369                           eval[i]);
370          index2 += seperator [i][k][n+1]-seperator[i][k][n];
371        }
372        factor += tmp1;
373      }
374
375      for (m= 0; m < maxTerms; m++)
376        delete [] mat [m];
377      delete [] mat;
378    }
379
380    if (fdivides (factor, B, quot))
381    {
382      result.append (factor);
383      B= quot;
384      if (result.length() == biFactors.length() - 1)
385      {
386        result.append (quot);
387        break;
388      }
389    }
390  }
391
392  //delete
393  for (i= 0; i < j; i++)
394  {
395    delete [] storeFactors [i] [0];
396    delete [] storeFactors [i] [1];
397    delete [] storeFactors [i];
398    for (k= 0; k < minFactorsLength; k++)
399      delete [] seperator [i][k];
400    delete [] seperator [i];
401  }
402  delete [] seperator;
403  delete [] storeFactors;
404  //
405
406  return result;
407}
Note: See TracBrowser for help on using the repository browser.