source: git/factory/facSparseHensel.cc @ ce21bd8

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