source: git/factory/facSparseHensel.h @ 7b5cb2

jengelh-datetimespielwiese
Last change on this file since 7b5cb2 was 7b5cb2, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: faster retrieval of terms chg: LucksWangSparseHeuristic now returns also partial factors
  • Property mode set to 100644
File size: 14.4 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 **/
11/*****************************************************************************/
12
13#ifndef FAC_SPARSE_HENSEL_H
14#define FAC_SPARSE_HENSEL_H
15
16#include "canonicalform.h"
17#include "cf_map_ext.h"
18#include "cf_iter.h"
19#include "templates/ftmpl_functions.h"
20#include "cf_algorithm.h"
21#include "cf_map.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, int l)
86{
87  int i= lo, j= hi;
88  CanonicalForm tmp= A[(lo+hi)/2];
89  while (i <= j)
90  {
91    if (l > 0)
92    {
93      while (comp (A [i], tmp, l) < 0 && i < hi) i++;
94      while (comp (tmp, A[j], l) < 0 && j > lo) j--;
95    }
96    else
97    {
98      while (comp (A [i], tmp) < 0 && i < hi) i++;
99      while (comp (tmp, A[j]) < 0 && j > lo) j--;
100    }
101    if (i <= j)
102    {
103      swap (A, i, j);
104      i++;
105      j--;
106    }
107  }
108  if (lo < j) quickSort (lo, j, A, l);
109  if (i < hi) quickSort (i, hi, A, l);
110}
111
112/// quick sort @a A
113inline
114void sort (CFArray& A, int l= 0)
115{
116  quickSort (0, A.size() - 1, A, l);
117}
118
119
120/// find normalizing factors for @a biFactors and build monic univariate factors
121/// from @a biFactors
122inline CFList
123findNormalizingFactor1 (const CFList& biFactors, const CanonicalForm& evalPoint,
124                        CFList& uniFactors)
125{
126  CFList result;
127  CanonicalForm tmp;
128  for (CFListIterator i= biFactors; i.hasItem(); i++)
129  {
130    tmp= i.getItem() (evalPoint);
131    uniFactors.append (tmp /Lc (tmp));
132    result.append (Lc (tmp));
133  }
134  return result;
135}
136
137/// find normalizing factors for @a biFactors and sort @a biFactors s.t.
138/// the returned @a biFactors evaluated at evalPoint coincide with @a uniFactors
139inline CFList
140findNormalizingFactor2 (CFList& biFactors, const CanonicalForm& evalPoint,
141                        const CFList& uniFactors)
142{
143  CFList result;
144  CFList uniBiFactors= biFactors;
145  CFList newBiFactors;
146  CFList l;
147  int pos;
148  CFListIterator iter;
149  for (iter= uniBiFactors; iter.hasItem(); iter++)
150  {
151    iter.getItem()= iter.getItem() (evalPoint);
152    l.append (Lc (iter.getItem()));
153    iter.getItem() /= Lc (iter.getItem());
154  }
155  for (CFListIterator i= uniFactors; i.hasItem(); i++)
156  {
157    pos= findItem (uniBiFactors, i.getItem());
158    newBiFactors.append (getItem (biFactors, pos));
159    result.append (getItem (l, pos));
160  }
161  biFactors= newBiFactors;
162  return result;
163}
164
165/// get terms of @a F
166inline CFArray
167getTerms (const CanonicalForm& F)
168{
169  if (F.inCoeffDomain())
170  {
171    CFArray result= CFArray (1);
172    result [0]= F;
173    return result;
174  }
175  if (F.isUnivariate())
176  {
177    CFArray result= CFArray (size(F));
178    int j= 0;
179    for (CFIterator i= F; i.hasTerms(); i++, j++)
180      result[j]= i.coeff()*power (F.mvar(), i.exp());
181    return result;
182  }
183  int numMon= size (F);
184  CFArray result= CFArray (numMon);
185  int j= 0;
186  CFArray recResult;
187  Variable x= F.mvar();
188  CanonicalForm powX;
189  for (CFIterator i= F; i.hasTerms(); i++)
190  {
191    powX= power (x, i.exp());
192    recResult= getTerms (i.coeff());
193    for (int k= 0; k < recResult.size(); k++)
194      result[j+k]= powX*recResult[k];
195    j += recResult.size();
196  }
197  return result;
198}
199
200/// helper function for getBiTerms
201inline CFArray
202getBiTerms_helper (const CanonicalForm& F, const CFMap& M)
203{
204  CFArray buf= CFArray (size (F));
205  int k= 0, level= F.level() - 1;
206  Variable x= F.mvar();
207  Variable y= Variable (F.level() - 1);
208  Variable one= Variable (1);
209  Variable two= Variable (2);
210  CFIterator j;
211  for (CFIterator i= F; i.hasTerms(); i++)
212  {
213    if (i.coeff().level() < level)
214    {
215      buf[k]= M (i.coeff())*power (one,i.exp());
216      k++;
217      continue;
218    }
219    j= i.coeff();
220    for (;j.hasTerms(); j++, k++)
221      buf[k]= power (one,i.exp())*power (two,j.exp())*M (j.coeff());
222  }
223  CFArray result= CFArray (k);
224  for (int i= 0; i < k; i++)
225    result[i]= buf[i];
226  return result;
227}
228
229/// get terms of @a F where F is considered a bivariate poly in Variable(1),
230/// Variable (2)
231inline CFArray
232getBiTerms (const CanonicalForm& F)
233{
234  if (F.inCoeffDomain())
235  {
236    CFArray result= CFArray (1);
237    result [0]= F;
238    return result;
239  }
240  if (F.isUnivariate())
241  {
242    CFArray result= CFArray (size(F));
243    int j= 0;
244    for (CFIterator i= F; i.hasTerms(); i++, j++)
245      result[j]= i.coeff()*power (F.mvar(), i.exp());
246    return result;
247  }
248
249  CanonicalForm G= F;
250
251  CFMap M;
252  M.newpair (Variable (1), F.mvar());
253  M.newpair (Variable (2), Variable (F.level() - 1));
254  G= swapvar (F, Variable (1), F.mvar());
255  G= swapvar (G, Variable (2), Variable (F.level() - 1));
256
257  CFArray result= getBiTerms_helper (G, M);
258  return result;
259}
260
261/// build a poly from entries in @a A
262inline CanonicalForm
263buildPolyFromArray (const CFArray& A)
264{
265  CanonicalForm result= 0;
266  for (int i= A.size() - 1; i > -1; i--)
267    result += A[i];
268  return result;
269}
270
271/// group together elements in @a A, where entries in @a A are put together
272/// if they coincide up to level @a level
273inline void
274groupTogether (CFArray& A, int level)
275{
276  int n= A.size() - 1;
277  int k= A.size();
278  for (int i= 0; i < n; i++)
279  {
280    if (comp (A[i],A[i+1], level) == 0)
281    {
282      A[i+1] += A[i];
283      A[i]= 0;
284      k--;
285    }
286  }
287  if (A[n].isZero())
288    k--;
289  CFArray B= CFArray (k);
290  n++;
291  k= 0;
292  for (int i= 0; i < n; i++)
293  {
294    if (!A[i].isZero())
295    {
296      B[k]= A[i];
297      k++;
298    }
299  }
300  A= B;
301}
302
303/// strip off those parts of entries in @a F whose level is less than or equal
304/// than @a level and store the stripped off parts in @a G
305inline void
306strip (CFArray& F, CFArray& G, int level)
307{
308  int n, m, i, j;
309  CanonicalForm g;
310  m= F.size();
311  G= CFArray (m);
312  for (j= 0; j < m; j++)
313  {
314    g= 1;
315    for (i= 1; i <= level; i++)
316    {
317      if ((n= degree (F[j],i)) > 0)
318        g *= power (Variable (i), n);
319    }
320    F[j] /= g;
321    G[j]= g;
322  }
323}
324
325/// s.a. stripped off parts are not returned
326inline void
327strip (CFArray& F, int level)
328{
329  int n, m, i;
330  CanonicalForm g;
331  m= F.size();
332  for (int j= 0; j < m; j++)
333  {
334    g= 1;
335    for (i= 1; i <= level; i++)
336    {
337      if ((n= degree (F[j],i)) > 0)
338        g *= power (Variable (i), n);
339    }
340    F[j] /= g;
341  }
342}
343
344/// get equations for LucksWangSparseHeuristic
345inline
346CFArray getEquations (const CFArray& A, const CFArray& B)
347{
348  ASSERT (A.size() == B.size(), "size of A and B has to coincide");
349  CFArray result= CFArray (A.size());
350  int n= A.size();
351  for (int i= 0; i < n; i++)
352    result[i]= A[i]-B[i];
353  return result;
354}
355
356/// evaluate every entry of @a A at @a B and level @a level
357inline void
358evaluate (CFArray& A, const CanonicalForm& B, int level)
359{
360  int n= A.size();
361  for (int i= 0; i < n; i++)
362  {
363    if (A[i].level() < level)
364      continue;
365    else
366      A[i]= A[i] (B, level);
367  }
368}
369
370/// evaluate every entry of @a A at every entry of @a B starting at level @a
371/// level
372inline void
373evaluate (CFArray& A, const CFArray& B, int level)
374{
375  int n= B.size();
376  for (int i= 0; i < n; i++)
377  {
378    if (!B[i].isZero())
379      evaluate (A, B[i], level + i);
380  }
381}
382
383/// simplify @a A if possible, i.e. @a A consists of 2 terms and contains only
384/// one variable of level greater or equal than @a level
385inline CanonicalForm
386simplify (const CanonicalForm& A, int level)
387{
388  CanonicalForm F= 0;
389  if (size (A, A.level()) == 2)
390  {
391    CanonicalForm C= getVars (A);
392    if ((C/C.mvar()).level() < level)
393    {
394      CanonicalForm B= LC (A);
395      if (B.level() < level)
396        F= -tailcoeff (A/B);
397    }
398  }
399  return F;
400}
401
402///  if possible simplify @a A as described above and store result in @a B
403inline bool
404simplify (CFArray& A, CFArray& B, int level)
405{
406  int n= A.size();
407  CanonicalForm f;
408  int index;
409  for (int i= 0; i < n; i++)
410  {
411    if (!A[i].isZero())
412    {
413      f= simplify (A[i], level);
414      if (!f.isZero())
415      {
416        index= A[i].level() - level;
417        if (index < 0 || index >= B.size())
418          return false;
419        if (!B[index].isZero() && B[index] != f)
420          return false;
421        else if (B[index].isZero())
422        {
423          B[index]= f;
424          A[i]= 0;
425        }
426      }
427    }
428  }
429  return true;
430}
431
432/// merge @a B into @a A if possible, i.e. every non-zero entry in @a A should
433/// be zero in @a B
434inline bool
435merge (CFArray& A, CFArray& B)
436{
437  if (A.size() != B.size())
438    return false;
439  int n= A.size();
440  for (int i= 0; i < n; i++)
441  {
442    if (!B[i].isZero())
443    {
444      if (A[i].isZero())
445      {
446        A[i]= B[i];
447        B[i]= 0;
448      }
449      else if (A[i] == B[i])
450        B[i]= 0;
451      else
452        return false;
453    }
454  }
455  return true;
456}
457
458/// checks if entries of @a A are zero
459inline bool
460isZero (const CFArray& A)
461{
462  int n= A.size();
463  for (int i= 0; i < n; i++)
464    if (!A[i].isZero())
465      return false;
466  return true;
467}
468
469/// checks if @a A equals @a B
470inline bool
471isEqual (const CFArray& A, const CFArray& B)
472{
473  if (A.size() != B.size())
474    return false;
475  int i, n= B.size();
476  for (i= 0; i < n; i++)
477    if (A[i] != B[i])
478      return false;
479  return true;
480}
481
482/// get terms of @a F wrt. Variable (1)
483inline CFArray
484getTerms2 (const CanonicalForm& F)
485{
486  if (F.inCoeffDomain())
487  {
488    CFArray result= CFArray (1);
489    result[0]= F;
490    return result;
491  }
492  CFArray result= CFArray (size (F));
493  int j= 0;
494  Variable x= F.mvar();
495  Variable y= Variable (1);
496  CFIterator k;
497  for (CFIterator i= F; i.hasTerms(); i++)
498  {
499    if (i.coeff().inCoeffDomain())
500    {
501      result[j]= i.coeff()*power (x,i.exp());
502      j++;
503    }
504    else
505    {
506      for (k= i.coeff(); k.hasTerms(); k++, j++)
507        result[j]= k.coeff()*power (x,i.exp())*power (y,k.exp());
508    }
509  }
510  sort (result);
511  return result;
512}
513
514/// get terms of entries in @a F and put them in @a result
515inline
516void getTerms2 (const CFList& F, CFArray* result)
517{
518  int j= 0;
519  for (CFListIterator i= F; i.hasItem(); i++, j++)
520    result[j]= getTerms2 (i.getItem());
521}
522
523/// evaluate entries in @a A at @a eval and @a y
524inline CFArray
525evaluate (const CFArray& A, const CanonicalForm& eval, const Variable& y)
526{
527  int j= A.size();
528  CFArray result= CFArray (j);
529  for (int i= 0; i < j; i++)
530    result [i]= A[i] (eval, y);
531  return result;
532}
533
534/// s.a.
535inline CFArray*
536evaluate (CFArray* const& A, int sizeA, const CanonicalForm& eval,
537          const Variable& y)
538{
539  CFArray* result= new CFArray [sizeA];
540  for (int i= 0; i < sizeA; i++)
541    result[i]= evaluate (A[i], eval, y);
542  return result;
543}
544
545/// normalize entries in @a L with @a normalizingFactor
546inline
547CFList normalize (const CFList& L, const CFList& normalizingFactor)
548{
549  CFList result;
550  CFListIterator j= normalizingFactor;
551  for (CFListIterator i= L; i.hasItem(); i++, j++)
552    result.append (i.getItem() / j.getItem());
553  return result;
554}
555
556/// search for @a F in @a A between index @a i and @a j
557inline
558int search (const CFArray& A, const CanonicalForm& F, int i, int j)
559{
560  for (; i < j; i++)
561    if (A[i] == F)
562      return i;
563  return -1;
564}
565
566/// patch together @a F1 and @a F2 and normalize by a power of @a eval
567/// @a F1 and @a F2 are assumed to be bivariate with one variable having level 1
568inline
569CanonicalForm patch (const CanonicalForm& F1, const CanonicalForm& F2,
570                     const CanonicalForm& eval)
571{
572  CanonicalForm result= F1;
573  if (F2.level() != 1 && !F2.inCoeffDomain())
574  {
575    int d= degree (F2);
576    result *= power (F2.mvar(), d);
577    result /= power (eval, d);
578  }
579  return result;
580}
581
582/// sparse heuristic lifting by Wang and Lucks
583///
584/// @return @a LucksWangSparseHeuristic returns true if it was successful
585int
586LucksWangSparseHeuristic (const CanonicalForm& F,     ///<[in] polynomial to be
587                                                      ///< factored
588                          const CFList& factors,      ///<[in] factors of F
589                                                      ///< lifted to level
590                          int level,                  ///<[in] level of lifted
591                                                      ///< factors
592                          const CFList& leadingCoeffs,///<[in] leading
593                                                      ///< coefficients of
594                                                      ///< factors
595                          CFList& result              ///<[in,out] result
596                         );
597
598/// sparse heuristic which patches together bivariate factors of @a A wrt.
599/// different second variables by their univariate images
600///
601/// @return @a sparseHeuristic returns a list of found factors of A
602CFList
603sparseHeuristic (const CanonicalForm& A,  ///<[in] polynomial to be factored
604                 const CFList& biFactors, ///<[in] bivariate factors of A where
605                                          ///< the second variable has level 2
606                 CFList*& moreBiFactors,  ///<[in] more bivariate factorizations
607                                          ///< wrt. different second variables
608                 const CFList& evaluation,///<[in] evaluation point
609                 int minFactorsLength     ///<[in] minimal length of bivariate
610                                          ///< factorizations
611                );
612
613#endif
Note: See TracBrowser for help on using the repository browser.