source: git/factory/facBivar.cc @ 3ef2d6

fieker-DuValspielwiese
Last change on this file since 3ef2d6 was 0e6668, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: more changes in order to replace old factorization over integers
  • Property mode set to 100644
File size: 15.5 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facBivar.cc
5 *
6 * bivariate factorization over Q(a)
7 *
8 * @author Martin Lee
9 *
10 * @internal @version \$Id$
11 *
12 **/
13/*****************************************************************************/
14
15#include "config.h"
16
17#include "assert.h"
18#include "debug.h"
19#include "timing.h"
20
21#include "cf_algorithm.h"
22#include "facFqBivar.h"
23#include "facBivar.h"
24#include "facHensel.h"
25
26#ifdef HAVE_NTL
27TIMING_DEFINE_PRINT(uni_factorize)
28TIMING_DEFINE_PRINT(hensel_lift12)
29
30CFList conv (const CFFList& L)
31{
32  CFList result;
33  for (CFFListIterator i= L; i.hasItem(); i++)
34    result.append (i.getItem().factor());
35  return result;
36}
37
38bool testPoint (const CanonicalForm& F, CanonicalForm& G, int i)
39{
40  G= F (i, 2);
41  if (G.inCoeffDomain() || degree (F, 1) > degree (G, 1))
42    return false;
43
44  if (degree (gcd (deriv (G, G.mvar()), G)) > 0)
45    return false;
46  return true;
47}
48
49CanonicalForm evalPoint (const CanonicalForm& F, int& i)
50{
51  Variable x= Variable (1);
52  Variable y= Variable (2);
53  CanonicalForm result;
54
55  int k;
56
57  if (i == 0)
58  {
59    if (testPoint (F, result, i))
60      return result;
61  }
62  do
63  {
64    if (i > 0)
65      k= 1;
66    else
67      k= 2;
68    while (k < 3)
69    {
70      if (k == 1)
71      {
72        if (testPoint (F, result, i))
73          return result;
74      }
75      else
76      {
77        if (testPoint (F, result, -i))
78          return result;
79      }
80      k++;
81    }
82    i++;
83  } while (1);
84}
85
86CFList
87earlyFactorDetection0 (CanonicalForm& F, CFList& factors,int& adaptedLiftBound,
88                      DegreePattern& degs, bool& success, int deg)
89{
90  DegreePattern bufDegs1= degs;
91  DegreePattern bufDegs2;
92  CFList result;
93  CFList T= factors;
94  CanonicalForm buf= F;
95  CanonicalForm LCBuf= LC (buf, Variable (1));
96  CanonicalForm g, quot;
97  CanonicalForm M= power (F.mvar(), deg);
98  adaptedLiftBound= 0;
99  int d= degree (F) + degree (LCBuf, F.mvar());
100  for (CFListIterator i= factors; i.hasItem(); i++)
101  {
102    if (!bufDegs1.find (degree (i.getItem(), 1)))
103      continue;
104    else
105    {
106      g= i.getItem() (0, 1);
107      g *= LCBuf;
108      g= mod (g, M);
109      if (fdivides (LC (g), LCBuf))
110      {
111        g= mulMod2 (i.getItem(), LCBuf, M);
112        g /= content (g, Variable (1));
113        if (fdivides (g, buf, quot))
114        {
115          result.append (g);
116          buf= quot;
117          d -= degree (g) + degree (LC (g, Variable (1)), F.mvar());
118          LCBuf= LC (buf, Variable (1));
119          T= Difference (T, CFList (i.getItem()));
120
121          // compute new possible degree pattern
122          bufDegs2= DegreePattern (T);
123          bufDegs1.intersect (bufDegs2);
124          bufDegs1.refine ();
125          if (bufDegs1.getLength() <= 1)
126          {
127            result.append (buf);
128            break;
129          }
130        }
131      }
132    }
133  }
134  adaptedLiftBound= d + 1;
135  if (d < deg)
136  {
137    factors= T;
138    degs= bufDegs1;
139    F= buf;
140    success= true;
141  }
142  if (bufDegs1.getLength() <= 1)
143    degs= bufDegs1;
144  return result;
145}
146
147
148CFList
149henselLiftAndEarly0 (CanonicalForm& A, bool& earlySuccess, CFList&
150                    earlyFactors, DegreePattern& degs, int& liftBound,
151                    const CFList& uniFactors, const CanonicalForm& eval)
152{
153  int sizeOfLiftPre;
154  int * liftPre= getLiftPrecisions (A, sizeOfLiftPre, degree (LC (A, 1), 2));
155
156  Variable x= Variable (1);
157  Variable y= Variable (2);
158  CFArray Pi;
159  CFList diophant;
160  CFList bufUniFactors= uniFactors;
161  bufUniFactors.insert (LC (A, x));
162  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
163  earlySuccess= false;
164  int newLiftBound= 0;
165  int smallFactorDeg= tmin (11, liftPre [sizeOfLiftPre- 1] + 1);//this is a tunable parameter
166  if (smallFactorDeg >= liftBound || degree (A,y) <= 4)
167    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
168  else if (sizeOfLiftPre > 1 && sizeOfLiftPre < 30)
169  {
170    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
171    earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
172                                        degs, earlySuccess,
173                                        smallFactorDeg);
174    if (degs.getLength() > 1 && !earlySuccess &&
175        smallFactorDeg != liftPre [sizeOfLiftPre-1] + 1)
176    {
177      if (newLiftBound > liftPre[sizeOfLiftPre-1]+1)
178      {
179        bufUniFactors.insert (LC (A, x));
180        henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
181                            liftPre[sizeOfLiftPre-1] + 1, Pi, diophant, M);
182        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
183                      degs, earlySuccess, liftPre[sizeOfLiftPre-1] + 1);
184      }
185    }
186    else if (earlySuccess)
187      liftBound= newLiftBound;
188    int i= sizeOfLiftPre - 1;
189    while (degs.getLength() > 1 && !earlySuccess && i - 1 >= 0)
190    {
191      if (newLiftBound > liftPre[i] + 1)
192      {
193        bufUniFactors.insert (LC (A, x));
194        henselLiftResume12 (A, bufUniFactors, liftPre[i] + 1,
195                            liftPre[i-1] + 1, Pi, diophant, M);
196        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
197                      degs, earlySuccess, liftPre[i-1] + 1);
198      }
199      else
200      {
201        liftBound= newLiftBound;
202        break;
203      }
204      i--;
205    }
206    if (earlySuccess)
207      liftBound= newLiftBound;
208    //after here all factors are lifted to liftPre[sizeOfLiftPre-1]
209  }
210  else
211  {
212    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
213    earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
214                                        degs, earlySuccess,
215                                        smallFactorDeg);
216    int i= 1;
217    while ((degree (A,y)/4)*i + 4 <= smallFactorDeg)
218      i++;
219    if (degs.getLength() > 1 && !earlySuccess)
220    {
221      bufUniFactors.insert (LC (A, x));
222      henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
223                          (degree (A, y)/4)*i + 4, Pi, diophant, M);
224      earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
225                    degs, earlySuccess, (degree (A, y)/4)*i + 4);
226    }
227    while (degs.getLength() > 1 && !earlySuccess && i < 4)
228    {
229      if (newLiftBound > (degree (A, y)/4)*i + 4)
230      {
231        bufUniFactors.insert (LC (A, x));
232        henselLiftResume12 (A, bufUniFactors, (degree (A,y)/4)*i + 4,
233                            (degree (A, y)/4)*(i+1) + 4, Pi, diophant, M);
234        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
235                      degs, earlySuccess, (degree (A, y)/4)*(i+1) + 4);
236      }
237      else
238      {
239        liftBound= newLiftBound;
240        break;
241      }
242      i++;
243    }
244    if (earlySuccess)
245      liftBound= newLiftBound;
246  }
247
248  return bufUniFactors;
249}
250
251CFList biFactorize (const CanonicalForm& F, const Variable& v)
252{
253  if (F.inCoeffDomain())
254    return CFList(F);
255
256  bool extension= (v.level() != 1);
257  CanonicalForm A;
258  CanonicalForm multiplier= 1;
259  if (isOn (SW_RATIONAL))
260    A= F*bCommonDen (F);
261  else
262    A= F;
263
264  if (A.isUnivariate())
265  {
266    CFFList buf;
267    if (extension)
268      buf= factorize (A, v);
269    else
270      buf= factorize (A, true);
271    CFList result= conv (buf);
272    if (result.getFirst().inCoeffDomain())
273      result.removeFirst();
274    return result;
275  }
276
277  CFMap N;
278  A= compress (A, N);
279  Variable y= A.mvar();
280
281  if (y.level() > 2) return CFList (F);
282  Variable x= Variable (1);
283
284  CanonicalForm contentAx= content (A, x);
285  CanonicalForm contentAy= content (A);
286
287  A= A/(contentAx*contentAy);
288  CFFList contentAxFactors, contentAyFactors;
289
290  if (extension)
291  {
292    if (!contentAx.inCoeffDomain())
293    {
294      contentAxFactors= factorize (contentAx, v);
295      if (contentAxFactors.getFirst().factor().inCoeffDomain())
296        contentAxFactors.removeFirst();
297    }
298    if (!contentAy.inCoeffDomain())
299    {
300      contentAyFactors= factorize (contentAy, v);
301      if (contentAyFactors.getFirst().factor().inCoeffDomain())
302        contentAyFactors.removeFirst();
303    }
304  }
305  else
306  {
307    if (!contentAx.inCoeffDomain())
308    {
309      contentAxFactors= factorize (contentAx, true);
310      if (contentAxFactors.getFirst().factor().inCoeffDomain())
311        contentAxFactors.removeFirst();
312    }
313    if (!contentAy.inCoeffDomain())
314    {
315      contentAyFactors= factorize (contentAy, true);
316      if (contentAyFactors.getFirst().factor().inCoeffDomain())
317        contentAyFactors.removeFirst();
318    }
319  }
320
321  //check trivial case
322  if (degree (A) == 1 || degree (A, 1) == 1)
323  {
324    CFList factors;
325    factors.append (A);
326
327    appendSwapDecompress (factors, conv (contentAxFactors),
328                          conv (contentAyFactors), false, false, N);
329
330    normalize (factors);
331    return factors;
332  }
333
334  //trivial case
335  CFList factors;
336  if (A.inCoeffDomain())
337  {
338    append (factors, conv (contentAxFactors));
339    append (factors, conv (contentAyFactors));
340    decompress (factors, N);
341    return factors;
342  }
343  else if (A.isUnivariate())
344  {
345    if (extension)
346      factors= conv (factorize (A, v));
347    else
348      factors= conv (factorize (A, true));
349    append (factors, conv (contentAxFactors));
350    append (factors, conv (contentAyFactors));
351    decompress (factors, N);
352    return factors;
353  }
354
355  bool swap= false;
356  if (degree (A) > degree (A, x))
357  {
358    A= swapvar (A, y, x);
359    swap= true;
360  }
361
362  CanonicalForm Aeval, bufAeval, buf;
363  CFList uniFactors, list, bufUniFactors;
364  DegreePattern degs;
365  DegreePattern bufDegs;
366
367  CanonicalForm Aeval2, bufAeval2;
368  CFList bufUniFactors2, list2, uniFactors2;
369  DegreePattern degs2;
370  DegreePattern bufDegs2;
371  bool swap2= false;
372
373  // several univariate factorizations to obtain more information about the
374  // degree pattern therefore usually less combinations have to be tried during
375  // the recombination process
376  int factorNums= 2;
377  int subCheck1= substituteCheck (A, x);
378  int subCheck2= substituteCheck (A, y);
379  buf= swapvar (A,x,y);
380  int evaluation, evaluation2, bufEvaluation= 0, bufEvaluation2= 0;
381  for (int i= 0; i < factorNums; i++)
382  {
383    bufAeval= A;
384    bufAeval= evalPoint (A, bufEvaluation);
385
386    bufAeval2= buf;
387    bufAeval2= evalPoint (buf, bufEvaluation2);
388
389
390    // univariate factorization
391    TIMING_START (uni_factorize);
392
393    if (extension)
394      bufUniFactors= conv (factorize (bufAeval, v));
395    else
396      bufUniFactors= conv (factorize (bufAeval, true));
397    TIMING_END_AND_PRINT (uni_factorize,
398                          "time for univariate factorization: ");
399    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
400              (prod (bufUniFactors)*Lc (bufAeval) == bufAeval));
401
402    TIMING_START (uni_factorize);
403    if (extension)
404      bufUniFactors2= conv (factorize (bufAeval2, v));
405    else
406      bufUniFactors2= conv (factorize (bufAeval2, true));
407    TIMING_END_AND_PRINT (uni_factorize,
408                          "time for univariate factorization in y: ");
409    DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
410              (prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2));
411
412    if (bufUniFactors.getFirst().inCoeffDomain())
413      bufUniFactors.removeFirst();
414    if (bufUniFactors2.getFirst().inCoeffDomain())
415      bufUniFactors2.removeFirst();
416    if (bufUniFactors.length() == 1 || bufUniFactors2.length() == 1)
417    {
418      factors.append (A);
419
420      appendSwapDecompress (factors, conv (contentAxFactors),
421                            conv (contentAyFactors), swap, swap2, N);
422
423      if (isOn (SW_RATIONAL))
424        normalize (factors);
425      return factors;
426    }
427
428    if (i == 0)
429    {
430      if (subCheck1 > 0)
431      {
432        int subCheck= substituteCheck (bufUniFactors);
433
434        if (subCheck > 1 && (subCheck1%subCheck == 0))
435        {
436          CanonicalForm bufA= A;
437          subst (bufA, bufA, subCheck, x);
438          factors= biFactorize (bufA, v);
439          reverseSubst (factors, subCheck, x);
440          appendSwapDecompress (factors, conv (contentAxFactors),
441                                conv (contentAyFactors), swap, swap2, N);
442          if (isOn (SW_RATIONAL))
443            normalize (factors);
444          return factors;
445        }
446      }
447
448      if (subCheck2 > 0)
449      {
450        int subCheck= substituteCheck (bufUniFactors2);
451
452        if (subCheck > 1 && (subCheck2%subCheck == 0))
453        {
454          CanonicalForm bufA= A;
455          subst (bufA, bufA, subCheck, y);
456          factors= biFactorize (bufA, v);
457          reverseSubst (factors, subCheck, y);
458          appendSwapDecompress (factors, conv (contentAxFactors),
459                                conv (contentAyFactors), swap, swap2, N);
460          if (isOn (SW_RATIONAL))
461            normalize (factors);
462          return factors;
463        }
464      }
465    }
466
467    // degree analysis
468    bufDegs = DegreePattern (bufUniFactors);
469    bufDegs2= DegreePattern (bufUniFactors2);
470
471    if (i == 0)
472    {
473      Aeval= bufAeval;
474      evaluation= bufEvaluation;
475      uniFactors= bufUniFactors;
476      degs= bufDegs;
477      Aeval2= bufAeval2;
478      evaluation2= bufEvaluation2;
479      uniFactors2= bufUniFactors2;
480      degs2= bufDegs2;
481    }
482    else
483    {
484      degs.intersect (bufDegs);
485      degs2.intersect (bufDegs2);
486      if (bufUniFactors2.length() < uniFactors2.length())
487      {
488        uniFactors2= bufUniFactors2;
489        Aeval2= bufAeval2;
490        evaluation2= bufEvaluation2;
491      }
492      if (bufUniFactors.length() < uniFactors.length())
493      {
494        if (evaluation != 0)
495        {
496          uniFactors= bufUniFactors;
497          Aeval= bufAeval;
498          evaluation= bufEvaluation;
499        }
500      }
501    }
502    bufEvaluation++;
503    bufEvaluation2++;
504  }
505
506  if (evaluation != 0 && (uniFactors.length() > uniFactors2.length() ||
507      (uniFactors.length() == uniFactors2.length()
508       && degs.getLength() > degs2.getLength())))
509  {
510    degs= degs2;
511    uniFactors= uniFactors2;
512    evaluation= evaluation2;
513    Aeval= Aeval2;
514    A= buf;
515    swap2= true;
516  }
517
518
519  if (degs.getLength() == 1) // A is irreducible
520  {
521    factors.append (A);
522    appendSwapDecompress (factors, conv (contentAxFactors),
523                          conv (contentAyFactors), swap, swap2, N);
524    if (isOn (SW_RATIONAL))
525      normalize (factors);
526    return factors;
527  }
528
529  if (!extension)
530  {
531    for (CFListIterator i= uniFactors; i.hasItem(); i++)
532      i.getItem() /= lc (i.getItem());
533  }
534
535  A= A (y + evaluation, y);
536
537  int liftBound= degree (A, y) + 1;
538
539  bool earlySuccess= false;
540  CFList earlyFactors;
541  TIMING_START (fac_hensel_lift);
542  //out_cf ("A before= ",A, "\n");
543  uniFactors= henselLiftAndEarly0
544             (A, earlySuccess, earlyFactors, degs, liftBound,
545              uniFactors, evaluation);
546  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
547  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
548
549  //printf ("earlyFactors.length()= %d\n", earlyFactors.length());
550  //printf ("liftBound after= %d\n", liftBound);
551  //printf ("earlySuccess= %d\n", earlySuccess);
552  CanonicalForm MODl= power (y, liftBound);
553
554  CanonicalForm test= prod (uniFactors);
555  test *= LC (A,1);
556  test= mod (test, MODl);
557  //printf ("test == A %d\n", test == A);
558  //out_cf ("test= ", test, "\n");
559  //out_cf ("A= ", A, "\n");
560  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
561                                uniFactors.length()/2);
562
563  if (earlySuccess)
564    factors= Union (earlyFactors, factors);
565  else if (!earlySuccess && degs.getLength() == 1)
566    factors= earlyFactors;
567
568  for (CFListIterator i= factors; i.hasItem(); i++)
569    i.getItem()= i.getItem() (y - evaluation, y);
570
571  appendSwapDecompress (factors, conv (contentAxFactors),
572                        conv (contentAyFactors), swap, swap2, N);
573  if (isOn (SW_RATIONAL))
574    normalize (factors);
575
576  return factors;
577}
578
579#endif
Note: See TracBrowser for help on using the repository browser.