source: git/factory/facSparseHensel.cc @ 72bfc8

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