source: git/factory/facSparseHensel.cc @ 767da0

spielwiese
Last change on this file since 767da0 was 8c73aaf, checked in by Hans Schoenemann <hannes@…>, 4 years ago
use flint, comment unused stuff
  • Property mode set to 100644
File size: 10.8 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
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
25int
26LucksWangSparseHeuristic (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  }
138  CanonicalForm factor;
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
159CFList
160sparseHeuristic (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;
174  for (CFListIterator iter= evaluation; iter.hasItem(); iter++,i--)
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}
Note: See TracBrowser for help on using the repository browser.