source: git/factory/facSparseHensel.cc @ 5583e1

jengelh-datetimespielwiese
Last change on this file since 5583e1 was 5583e1, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: deleted too much Signed-off-by: Martin Lee <martinlee84@web.de>
  • 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        delete [] seperator;
237        return CFList();
238      }
239      seperator [i][k]= new int [count + 3];
240      seperator [i][k][0]= count + 1;
241      seperator [i][k][1]= 0;
242      count= 2;
243      for (l= 0; l < storeFactors [i][0][k].size() - 1; l++)
244      {
245        if (degree (storeFactors[i][0][k][l], x) <
246            degree (storeFactors[i][0][k][l+1], x))
247        {
248          seperator[i][k][count]=l + 1;
249          count++;
250        }
251      }
252      seperator [i][k][count]= storeFactors[i][0][k].size();
253    }
254  }
255
256  CanonicalForm tmp1, factor, quot;
257  CanonicalForm B= A;
258  CFList result;
259  int maxTerms, n, index1, index2, mmm, found, columns, oneCount;
260  int ** mat;
261
262  for (k= 0; k < minFactorsLength; k++)
263  {
264    factor= 0;
265    sizeOfUniFactors= seperator [0][k][0];
266    for (n= 1; n <= sizeOfUniFactors; n++)
267    {
268      columns= 0;
269      maxTerms= 1;
270      index1= j - 1;
271      for (i= j - 1; i >= 0; i--)
272      {
273        if (maxTerms < seperator[i][k][n+1]-seperator[i][k][n])
274        {
275          maxTerms= seperator[i][k][n + 1]-seperator[i][k][n];
276          index1= i;
277        }
278      }
279      for (i= j - 1; i >= 0; i--)
280      {
281        if (i == index1)
282          continue;
283        columns += seperator [i][k][n+1]-seperator[i][k][n];
284      }
285      mat= new int *[maxTerms];
286      mm= 0;
287      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
288      {
289        tmp1= storeFactors [index1][1][k][m];
290        mat[mm]= new int [columns];
291        for (i= 0; i < columns; i++)
292          mat[mm][i]= 0;
293        index2= 0;
294        for (i= j - 1; i >= 0; i--)
295        {
296          if (i == index1)
297            continue;
298          found= -1;
299          if ((found= search (storeFactors[i][1][k], tmp1, 
300                              seperator[i][k][n], seperator[i][k][n+1])) >= 0)
301            mat[mm][index2 + found - seperator [i][k][n]]= 1;
302          index2 += seperator [i][k][n+1]-seperator[i][k][n];
303        }
304      }
305
306      index2= 0;
307      for (i= j - 1; i >= 0; i--)
308      {
309        if (i == index1)
310          continue;
311        oneCount= 0;
312        for (mm= 0; mm < seperator [i][k][n + 1] - seperator [i][k][n]; mm++)
313        {
314          for (m= 0; m < maxTerms; m++)
315          {
316            if (mat[m][mm+index2] == 1)
317              oneCount++;
318          }
319        }
320        if (oneCount == seperator [i][k][n+1]-seperator[i][k][n] - 1)
321        {
322          for (mm= 0; mm < seperator [i][k][n+1]-seperator[i][k][n]; mm++)
323          {
324            oneCount= 0;
325            for (m= 0; m < maxTerms; m++)
326              if (mat[m][mm+index2] == 1)
327                oneCount++;
328            if (oneCount > 0)
329              continue;
330            for (m= 0; m < maxTerms; m++)
331            {
332              oneCount= 0;
333              for (mmm= 0; mmm < seperator[i][k][n+1]-seperator[i][k][n]; mmm++)
334              {
335                if (mat[m][mmm+index2] == 1)
336                  oneCount++;
337              }
338              if (oneCount > 0)
339                continue;
340              mat[m][mm+index2]= 1;
341            }
342          }
343        }
344        index2 += seperator [i][k][n+1] - seperator [i][k][n];
345      }
346
347      //read off solution
348      mm= 0;
349      for (m= seperator[index1][k][n]; m < seperator[index1][k][n+1]; m++, mm++)
350      {
351        tmp1= storeFactors [index1][0][k][m];
352        index2= 0;
353        for (i= j - 1; i > -1; i--)
354        {
355          if (i == index1)
356            continue;
357          for (mmm= 0; mmm < seperator [i][k][n+1]-seperator[i][k][n]; mmm++)
358            if (mat[mm][mmm+index2] == 1)
359              tmp1= patch (tmp1, storeFactors[i][0][k][seperator[i][k][n]+mmm],
360                           eval[i]);
361          index2 += seperator [i][k][n+1]-seperator[i][k][n];
362        }
363        factor += tmp1;
364      }
365
366      for (m= 0; m < maxTerms; m++)
367        delete [] mat [m];
368      delete [] mat;
369    }
370
371    if (fdivides (factor, B, quot))
372    {
373      result.append (factor);
374      B= quot;
375      if (result.length() == biFactors.length() - 1)
376      {
377        result.append (quot);
378        break;
379      }
380    }
381  }
382
383  //delete
384  for (i= 0; i < j; i++)
385  {
386    delete [] storeFactors [i] [0];
387    delete [] storeFactors [i] [1];
388    delete [] storeFactors [i];
389    for (k= 0; k < minFactorsLength; k++)
390      delete [] seperator [i][k];
391    delete [] seperator [i];
392  }
393  delete [] seperator;
394  delete [] storeFactors;
395  //
396
397  return result;
398}
Note: See TracBrowser for help on using the repository browser.