source: git/factory/facBivar.cc @ 667ba1

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