source: git/factory/facSparseHensel.h @ 2f6b737

jengelh-datetimespielwiese
Last change on this file since 2f6b737 was 2f6b737, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: getTerms2 and patch for coeffs
  • Property mode set to 100644
File size: 12.8 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqFactorize.h
5 *
6 * This file provides functions for sparse heuristic Hensel lifting
7 *
8 * @author Martin Lee
9 *
10 * @internal @version \$Id$
11 *
12 **/
13/*****************************************************************************/
14
15#ifndef FAC_SPARSE_HENSEL_H
16#define FAC_SPARSE_HENSEL_H
17
18#include "canonicalform.h"
19#include "cf_map_ext.h"
20#include "cf_iter.h"
21#include "templates/ftmpl_functions.h"
22
23/// compare polynomials
24inline
25int comp (const CanonicalForm& A, const CanonicalForm& B)
26{
27  if (A.inCoeffDomain() && !B.inCoeffDomain())
28    return -1;
29  else if (!A.inCoeffDomain() && B.inCoeffDomain())
30    return 1;
31  else if (A.inCoeffDomain() && B.inCoeffDomain())
32    return 0;
33  else if (degree (A, 1) > degree (B, 1))
34    return 1;
35  else if (degree (A, 1) < degree (B, 1))
36    return -1;
37  // here A and B are not in CoeffDomain
38  int n= tmax (A.level(), B.level());
39  for (int i= 2; i <= n; i++)
40  {
41    if (degree (A,i) > degree (B,i))
42      return 1;
43    else if (degree (A,i) < degree (B,i))
44      return -1;
45  }
46  return 0;
47}
48
49/// compare two polynomials up to level @a level
50inline
51int comp (const CanonicalForm& A, const CanonicalForm& B, int level)
52{
53  if (A.inCoeffDomain() && !B.inCoeffDomain() && B.level() <= level)
54    return -1;
55  else if (!A.inCoeffDomain() && A.level() <= level && B.inCoeffDomain())
56    return 1;
57  else if (A.inCoeffDomain() && B.inCoeffDomain())
58    return 0;
59  else if (degree (A, 1) > degree (B, 1))
60    return 1;
61  else if (degree (A, 1) < degree (B, 1))
62    return -1;
63  // here A and B are not in coeffdomain
64  for (int i= 2; i <= level; i++)
65  {
66    if (degree (A,i) > degree (B,i))
67      return 1;
68    else if (degree (A,i) < degree (B,i))
69      return -1;
70  }
71  return 0;
72}
73
74/// swap entry @a i and @a j in @a A
75inline
76void swap (CFArray& A, int i, int j)
77{
78  CanonicalForm tmp= A[i];
79  A[i]= A[j];
80  A[j]= tmp;
81}
82
83/// quick sort helper function
84inline
85void quickSort (int lo, int hi, CFArray& A)
86{
87  int i= lo, j= hi;
88  CanonicalForm tmp= A[(lo+hi)/2];
89  while (i <= j)
90  {
91    while (comp (A [i], tmp) < 0 && i < hi) i++;
92    while (comp (tmp, A[j]) < 0 && j > lo) j--;
93    if (i <= j)
94    {
95      swap (A, i, j);
96      i++;
97      j--;
98    }
99  }
100  if (lo < j) quickSort (lo, j, A);
101  if (i < hi) quickSort (i, hi, A);
102}
103
104/// quick sort @a A
105inline
106void sort (CFArray& A)
107{
108  quickSort (0, A.size() - 1, A);
109}
110
111/// find normalizing factors for @a biFactors and build monic univariate factors
112/// from @a biFactors
113inline CFList
114findNormalizingFactor1 (const CFList& biFactors, const CanonicalForm& evalPoint,
115                        CFList& uniFactors)
116{
117  CFList result;
118  CanonicalForm tmp;
119  for (CFListIterator i= biFactors; i.hasItem(); i++)
120  {
121    tmp= i.getItem() (evalPoint);
122    uniFactors.append (tmp /Lc (tmp));
123    result.append (Lc (tmp));
124  }
125  return result;
126}
127
128/// find normalizing factors for @a biFactors and sort @a biFactors s.t.
129/// the returned @a biFactors evaluated at evalPoint coincide with @a uniFactors
130inline CFList
131findNormalizingFactor2 (CFList& biFactors, const CanonicalForm& evalPoint,
132                        const CFList& uniFactors)
133{
134  CFList result;
135  CFList uniBiFactors= biFactors;
136  CFList newBiFactors;
137  CFList l;
138  int pos;
139  CFListIterator iter;
140  for (iter= uniBiFactors; iter.hasItem(); iter++)
141  {
142    iter.getItem()= iter.getItem() (evalPoint);
143    l.append (Lc (iter.getItem()));
144    iter.getItem() /= Lc (iter.getItem());
145  }
146  for (CFListIterator i= uniFactors; i.hasItem(); i++)
147  {
148    pos= findItem (uniBiFactors, i.getItem());
149    newBiFactors.append (getItem (biFactors, pos));
150    result.append (getItem (l, pos));
151  }
152  biFactors= newBiFactors;
153  return result;
154}
155
156/// get terms of @a F
157inline CFArray
158getTerms (const CanonicalForm& F)
159{
160  if (F.inCoeffDomain())
161  {
162    CFArray result= CFArray (1);
163    result [0]= F;
164    return result;
165  }
166  if (F.isUnivariate())
167  {
168    CFArray result= CFArray (size(F));
169    int j= 0;
170    for (CFIterator i= F; i.hasTerms(); i++, j++)
171      result[j]= i.coeff()*power (F.mvar(), i.exp());
172    return result;
173  }
174  int numMon= size (F);
175  CFArray result= CFArray (numMon);
176  int j= 0;
177  CFArray recResult;
178  Variable x= F.mvar();
179  CanonicalForm powX;
180  for (CFIterator i= F; i.hasTerms(); i++)
181  {
182    powX= power (x, i.exp());
183    recResult= getTerms (i.coeff());
184    for (int k= 0; k < recResult.size(); k++)
185      result[j+k]= powX*recResult[k];
186    j += recResult.size();
187  }
188  return result;
189}
190
191/// build a poly from entries in @a A
192inline CanonicalForm
193buildPolyFromArray (const CFArray& A)
194{
195  CanonicalForm result= 0;
196  for (int i= A.size() - 1; i > -1; i--)
197    result += A[i];
198  return result;
199}
200
201/// group together elements in @a A, where entries in @a A are put put together
202/// if they coincide up to level @a level
203inline void
204groupTogether (CFArray& A, int level)
205{
206  int n= A.size() - 1;
207  int k= A.size();
208  for (int i= 0; i < n; i++)
209  {
210    if (comp (A[i],A[i+1], level) == 0)
211    {
212      A[i+1] += A[i];
213      A[i]= 0;
214      k--;
215    }
216  }
217  CFArray B= CFArray (k);
218  n++;
219  k= 0;
220  for (int i= 0; i < n; i++)
221  {
222    if (!A[i].isZero())
223    {
224      B[k]= A[i];
225      k++;
226    }
227  }
228  A= B;
229}
230
231/// strip off those parts of entries in @a F whose level is less than or equal
232/// than @a level and store the stripped off parts in @a G
233inline void
234strip (CFArray& F, CFArray& G, int level)
235{
236  int n, m, i, j;
237  CanonicalForm g;
238  m= F.size();
239  G= CFArray (m);
240  for (j= 0; j < m; j++)
241  {
242    g= 1;
243    for (i= 1; i <= level; i++)
244    {
245      if ((n= degree (F[j],i)) > 0)
246        g *= power (Variable (i), n);
247    }
248    F[j] /= g;
249    G[j]= g;
250  }
251}
252
253/// s.a. stripped off parts are not returned
254inline void
255strip (CFArray& F, int level)
256{
257  int n, m, i;
258  CanonicalForm g;
259  m= F.size();
260  for (int j= 0; j < m; j++)
261  {
262    g= 1;
263    for (i= 1; i <= level; i++)
264    {
265      if ((n= degree (F[j],i)) > 0)
266        g *= power (Variable (i), n);
267    }
268    F[j] /= g;
269  }
270}
271
272/// get equations for LucksWangSparseHeuristic
273inline
274CFArray getEquations (const CFArray& A, const CFArray& B)
275{
276  ASSERT (A.size() == B.size(), "size of A and B has to coincide");
277  CFArray result= CFArray (A.size());
278  int n= A.size();
279  for (int i= 0; i < n; i++)
280    result[i]= A[i]-B[i];
281  return result;
282}
283
284/// evaluate every entry of @a A at @a B and level @a level
285inline void
286evaluate (CFArray& A, const CanonicalForm& B, int level)
287{
288  int n= A.size();
289  for (int i= 0; i < n; i++)
290  {
291    if (A[i].level() < level)
292      continue;
293    else
294      A[i]= A[i] (B, level);
295  }
296}
297
298/// evaluate every entry of @a A at every entry of @a B starting at level @a
299/// level
300inline void
301evaluate (CFArray& A, const CFArray& B, int level)
302{
303  int n= B.size();
304  for (int i= 0; i < n; i++)
305  {
306    if (!B[i].isZero())
307      evaluate (A, B[i], level + i);
308  }
309}
310
311/// simplify @a A if possible, i.e. @a A consists of 2 terms and contains only
312/// one variable of level greater or equal than @a level
313inline CanonicalForm
314simplify (const CanonicalForm& A, int level)
315{
316  CanonicalForm F= 0;
317  if (size (A, A.level()) == 2)
318  {
319    CanonicalForm C= getVars (A);
320    if ((C/C.mvar()).level() < level)
321    {
322      CanonicalForm B= LC (A);
323      if (B.level() < level)
324        F= -tailcoeff (A/B);
325    }
326  }
327  return F;
328}
329
330///  if possible simplify @a A as described above and store result in @a B
331inline bool
332simplify (CFArray& A, CFArray& B, int level)
333{
334  int n= A.size();
335  CanonicalForm f;
336  int index;
337  for (int i= 0; i < n; i++)
338  {
339    if (!A[i].isZero())
340    {
341      f= simplify (A[i], level);
342      if (!f.isZero())
343      {
344        index= A[i].level() - level;
345        if (index < 0 || index >= B.size())
346          return false;
347        if (!B[index].isZero() && B[index] != f)
348          return false;
349        else if (B[index].isZero())
350        {
351          B[index]= f;
352          A[i]= 0;
353        }
354      }
355    }
356  }
357  return true;
358}
359
360/// merge @a B into @a A if possible, i.e. every non-zero entry in @a A should
361/// be zero in @a B
362inline bool
363merge (CFArray& A, CFArray& B)
364{
365  if (A.size() != B.size())
366    return false;
367  int n= A.size();
368  for (int i= 0; i < n; i++)
369  {
370    if (!B[i].isZero())
371    {
372      if (A[i].isZero())
373      {
374        A[i]= B[i];
375        B[i]= 0;
376      }
377      else if (A[i] == B[i])
378        B[i]= 0;
379      else
380        return false;
381    }
382  }
383  return true;
384}
385
386/// checks if entries of @a A are zero
387inline bool
388isZero (const CFArray& A)
389{
390  int n= A.size();
391  for (int i= 0; i < n; i++)
392    if (!A[i].isZero())
393      return false;
394  return true;
395}
396
397/// checks if @a A equals @a B
398inline bool
399isEqual (const CFArray& A, const CFArray& B)
400{
401  if (A.size() != B.size())
402    return false;
403  int i, n= B.size();
404  for (i= 0; i < n; i++)
405    if (A[i] != B[i])
406      return false;
407  return true;
408}
409
410/// get terms of @a F wrt. Variable (1)
411inline CFArray
412getTerms2 (const CanonicalForm& F)
413{
414  if (F.inCoeffDomain())
415  {
416    CFArray result= CFArray (1);
417    result[0]= F;
418    return result;
419  }
420  CFArray result= CFArray (size (F));
421  int j= 0;
422  Variable x= F.mvar();
423  Variable y= Variable (1);
424  CFIterator k;
425  for (CFIterator i= F; i.hasTerms(); i++)
426  {
427    if (i.coeff().inCoeffDomain())
428    {
429      result[j]= i.coeff()*power (x,i.exp());
430      j++;
431    }
432    else
433    {
434      for (k= i.coeff(); k.hasTerms(); k++, j++)
435        result[j]= k.coeff()*power (x,i.exp())*power (y,k.exp());
436    }
437  }
438  sort (result);
439  return result;
440}
441
442/// get terms of entries in @a F and put them in @a result
443inline
444void getTerms2 (const CFList& F, CFArray* result)
445{
446  int j= 0;
447  for (CFListIterator i= F; i.hasItem(); i++, j++)
448    result[j]= getTerms2 (i.getItem());
449}
450
451/// evaluate entries in @a A at @a eval and @a y
452inline CFArray
453evaluate (const CFArray& A, const CanonicalForm& eval, const Variable& y)
454{
455  int j= A.size();
456  CFArray result= CFArray (j);
457  for (int i= 0; i < j; i++)
458    result [i]= A[i] (eval, y);
459  return result;
460}
461
462/// s.a.
463inline CFArray*
464evaluate (CFArray* const& A, int sizeA, const CanonicalForm& eval,
465          const Variable& y)
466{
467  CFArray* result= new CFArray [sizeA];
468  for (int i= 0; i < sizeA; i++)
469    result[i]= evaluate (A[i], eval, y);
470  return result;
471}
472
473/// normalize entries in @a L with @a normalizingFactor
474inline
475CFList normalize (const CFList& L, const CFList& normalizingFactor)
476{
477  CFList result;
478  CFListIterator j= normalizingFactor;
479  for (CFListIterator i= L; i.hasItem(); i++, j++)
480    result.append (i.getItem() / j.getItem());
481  return result;
482}
483
484/// search for @a F in @a A between index @a i and @a j
485inline
486int search (const CFArray& A, const CanonicalForm& F, int i, int j)
487{
488  for (; i < j; i++)
489    if (A[i] == F)
490      return i;
491  return -1;
492}
493
494/// patch together @a F1 and @a F2 and normalize by a power of @a eval
495/// @a F1 and @a F2 are assumed to be bivariate with one variable having level 1
496inline
497CanonicalForm patch (const CanonicalForm& F1, const CanonicalForm& F2,
498                     const CanonicalForm& eval)
499{
500  CanonicalForm result= F1;
501  if (F2.level() != 1 && !F2.inCoeffDomain())
502  {
503    int d= degree (F2);
504    result *= power (F2.mvar(), d);
505    result /= power (eval, d);
506  }
507  return result;
508}
509
510/// sparse heuristic lifting by Wang and Lucks
511///
512/// @return @a LucksWangSparseHeuristic returns true if it was successful
513bool
514LucksWangSparseHeuristic (const CanonicalForm& F,     ///<[in] polynomial to be
515                                                      ///< factored
516                          const CFList& factors,      ///<[in] factors of F
517                                                      ///< lifted to level
518                          int level,                  ///<[in] level of lifted
519                                                      ///< factors
520                          const CFList& leadingCoeffs,///<[in] leading
521                                                      ///< coefficients of
522                                                      ///< factors
523                          CFList& result              ///<[in,out] result
524                         );
525
526/// sparse heuristic which patches together bivariate factors of @a A wrt.
527/// different second variables by their univariate images
528///
529/// @return @a sparseHeuristic returns a list of found factors of A
530CFList
531sparseHeuristic (const CanonicalForm& A,  ///<[in] polynomial to be factored
532                 const CFList& biFactors, ///<[in] bivariate factors of A where
533                                          ///< the second variable has level 2
534                 CFList*& moreBiFactors,  ///<[in] more bivariate factorizations
535                                          ///< wrt. different second variables
536                 const CFList& evaluation,///<[in] evaluation point
537                 int minFactorsLength     ///<[in] minimal length of bivariate
538                                          ///< factorizations
539                );
540
541#endif
Note: See TracBrowser for help on using the repository browser.