source: git/factory/facBivar.cc @ 27ab36

spielwiese
Last change on this file since 27ab36 was 27ab36, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: better debug output
  • Property mode set to 100644
File size: 19.1 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#include "cf_primes.h"
27
28#ifdef HAVE_NTL
29TIMING_DEFINE_PRINT(uni_factorize)
30TIMING_DEFINE_PRINT(hensel_lift12)
31
32// bound on coeffs of f (cf. Musser: Multivariate Polynomial Factorization,
33//                          Gelfond: Transcendental and Algebraic Numbers)
34modpk
35coeffBound ( const CanonicalForm & f, int p )
36{
37    int * degs = degrees( f );
38    int M = 0, i, k = f.level();
39    CanonicalForm b= 1;
40    for ( i = 1; i <= k; i++ )
41    {
42      M += degs[i];
43      b *= degs[i] + 1;
44    }
45    b /= power (CanonicalForm (2), k);
46    b= b.sqrt() + 1;
47    b *= 2 * maxNorm( f ) * power( CanonicalForm( 2 ), M );
48    CanonicalForm B = p;
49    k = 1;
50    while ( B < b ) {
51        B *= p;
52        k++;
53    }
54    return modpk( p, k );
55}
56
57void findGoodPrime(const CanonicalForm &f, int &start)
58{
59  if (! f.inBaseDomain() )
60  {
61    CFIterator i = f;
62    for(;;)
63    {
64      if  ( i.hasTerms() )
65      {
66        findGoodPrime(i.coeff(),start);
67        if (0==cf_getBigPrime(start)) return;
68        if((i.exp()!=0) && ((i.exp() % cf_getBigPrime(start))==0))
69        {
70          start++;
71          i=f;
72        }
73        else  i++;
74      }
75      else break;
76    }
77  }
78  else
79  {
80    if (f.inZ())
81    {
82      if (0==cf_getBigPrime(start)) return;
83      while((!f.isZero()) && (mod(f,cf_getBigPrime(start))==0))
84      {
85        start++;
86        if (0==cf_getBigPrime(start)) return;
87      }
88    }
89  }
90}
91
92modpk
93coeffBound ( const CanonicalForm & f, int p, const CanonicalForm& mipo )
94{
95    int * degs = degrees( f );
96    int M = 0, i, k = f.level();
97    CanonicalForm K= 1;
98    for ( i = 1; i <= k; i++ )
99    {
100        M += degs[i];
101        K *= degs[i] + 1;
102    }
103    K /= power (CanonicalForm (2), k);
104    K= K.sqrt()+1;
105    K *= power (CanonicalForm (2), M);
106    int N= degree (mipo);
107    CanonicalForm b;
108    b= 2*power (maxNorm (f), N)*power (maxNorm (mipo), 4*N)*K*
109       power (CanonicalForm (2), N)*(CanonicalForm (M+1).sqrt()+1)*
110       power (CanonicalForm (N+1).sqrt()+1, 7*N);
111    b /= power (abs (lc (mipo)), N);
112
113    ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
114    ZZX NTLLcf= convertFacCF2NTLZZX (Lc (f));
115    ZZ NTLf= resultant (NTLmipo, NTLLcf);
116    ZZ NTLD= discriminant (NTLmipo);
117    b /= abs (convertZZ2CF (NTLf))*abs (convertZZ2CF (NTLD));
118
119    CanonicalForm B = p;
120    k = 1;
121    while ( B < b ) {
122        B *= p;
123        k++;
124    }
125    return modpk( p, k );
126}
127
128CFList conv (const CFFList& L)
129{
130  CFList result;
131  for (CFFListIterator i= L; i.hasItem(); i++)
132    result.append (i.getItem().factor());
133  return result;
134}
135
136bool testPoint (const CanonicalForm& F, CanonicalForm& G, int i)
137{
138  G= F (i, 2);
139  if (G.inCoeffDomain() || degree (F, 1) > degree (G, 1))
140    return false;
141
142  if (degree (gcd (deriv (G, G.mvar()), G)) > 0)
143    return false;
144  return true;
145}
146
147CanonicalForm evalPoint (const CanonicalForm& F, int& i)
148{
149  Variable x= Variable (1);
150  Variable y= Variable (2);
151  CanonicalForm result;
152
153  int k;
154
155  if (i == 0)
156  {
157    if (testPoint (F, result, i))
158      return result;
159  }
160  do
161  {
162    if (i > 0)
163      k= 1;
164    else
165      k= 2;
166    while (k < 3)
167    {
168      if (k == 1)
169      {
170        if (testPoint (F, result, i))
171          return result;
172      }
173      else
174      {
175        if (testPoint (F, result, -i))
176        {
177          i= -i;
178          return result;
179        }
180        else if (i < 0)
181          i= -i;
182      }
183      k++;
184    }
185    i++;
186  } while (1);
187}
188
189CFList
190earlyFactorDetection0 (CanonicalForm& F, CFList& factors,int& adaptedLiftBound,
191                      DegreePattern& degs, bool& success, int deg)
192{
193  DegreePattern bufDegs1= degs;
194  DegreePattern bufDegs2;
195  CFList result;
196  CFList T= factors;
197  CanonicalForm buf= F;
198  CanonicalForm LCBuf= LC (buf, Variable (1));
199  CanonicalForm g, quot;
200  CanonicalForm M= power (F.mvar(), deg);
201  adaptedLiftBound= 0;
202  int d= degree (F) + degree (LCBuf, F.mvar());
203  for (CFListIterator i= factors; i.hasItem(); i++)
204  {
205    if (!bufDegs1.find (degree (i.getItem(), 1)))
206      continue;
207    else
208    {
209      g= i.getItem() (0, 1);
210      g *= LCBuf;
211      g= mod (g, M);
212      if (fdivides (LC (g), LCBuf))
213      {
214        g= mulMod2 (i.getItem(), LCBuf, M);
215        g /= content (g, Variable (1));
216        if (fdivides (g, buf, quot))
217        {
218          result.append (g);
219          buf= quot;
220          d -= degree (g) + degree (LC (g, Variable (1)), F.mvar());
221          LCBuf= LC (buf, Variable (1));
222          T= Difference (T, CFList (i.getItem()));
223
224          // compute new possible degree pattern
225          bufDegs2= DegreePattern (T);
226          bufDegs1.intersect (bufDegs2);
227          bufDegs1.refine ();
228          if (bufDegs1.getLength() <= 1)
229          {
230            result.append (buf);
231            break;
232          }
233        }
234      }
235    }
236  }
237  adaptedLiftBound= d + 1;
238  if (d < deg)
239  {
240    factors= T;
241    degs= bufDegs1;
242    F= buf;
243    success= true;
244  }
245  if (bufDegs1.getLength() <= 1)
246    degs= bufDegs1;
247  return result;
248}
249
250
251CFList
252henselLiftAndEarly0 (CanonicalForm& A, bool& earlySuccess, CFList&
253                    earlyFactors, DegreePattern& degs, int& liftBound,
254                    const CFList& uniFactors, const CanonicalForm& eval)
255{
256  int sizeOfLiftPre;
257  int * liftPre= getLiftPrecisions (A, sizeOfLiftPre, degree (LC (A, 1), 2));
258
259  Variable x= Variable (1);
260  Variable y= Variable (2);
261  CFArray Pi;
262  CFList diophant;
263  CFList bufUniFactors= uniFactors;
264  bufUniFactors.insert (LC (A, x));
265  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
266  earlySuccess= false;
267  int newLiftBound= 0;
268  int smallFactorDeg= tmin (11, liftPre [sizeOfLiftPre- 1] + 1);//this is a tunable parameter
269  if (smallFactorDeg >= liftBound || degree (A,y) <= 4)
270    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
271  else if (sizeOfLiftPre > 1 && sizeOfLiftPre < 30)
272  {
273    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
274    earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
275                                        degs, earlySuccess,
276                                        smallFactorDeg);
277    if (degs.getLength() > 1 && !earlySuccess &&
278        smallFactorDeg != liftPre [sizeOfLiftPre-1] + 1)
279    {
280      if (newLiftBound > liftPre[sizeOfLiftPre-1]+1)
281      {
282        bufUniFactors.insert (LC (A, x));
283        henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
284                            liftPre[sizeOfLiftPre-1] + 1, Pi, diophant, M);
285        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
286                      degs, earlySuccess, liftPre[sizeOfLiftPre-1] + 1);
287      }
288    }
289    else if (earlySuccess)
290      liftBound= newLiftBound;
291    int i= sizeOfLiftPre - 1;
292    while (degs.getLength() > 1 && !earlySuccess && i - 1 >= 0)
293    {
294      if (newLiftBound > liftPre[i] + 1)
295      {
296        bufUniFactors.insert (LC (A, x));
297        henselLiftResume12 (A, bufUniFactors, liftPre[i] + 1,
298                            liftPre[i-1] + 1, Pi, diophant, M);
299        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
300                      degs, earlySuccess, liftPre[i-1] + 1);
301      }
302      else
303      {
304        liftBound= newLiftBound;
305        break;
306      }
307      i--;
308    }
309    if (earlySuccess)
310      liftBound= newLiftBound;
311    //after here all factors are lifted to liftPre[sizeOfLiftPre-1]
312  }
313  else
314  {
315    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
316    earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
317                                        degs, earlySuccess,
318                                        smallFactorDeg);
319    int i= 1;
320    while ((degree (A,y)/4)*i + 4 <= smallFactorDeg)
321      i++;
322    if (degs.getLength() > 1 && !earlySuccess)
323    {
324      bufUniFactors.insert (LC (A, x));
325      henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
326                          (degree (A, y)/4)*i + 4, Pi, diophant, M);
327      earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
328                    degs, earlySuccess, (degree (A, y)/4)*i + 4);
329    }
330    while (degs.getLength() > 1 && !earlySuccess && i < 4)
331    {
332      if (newLiftBound > (degree (A, y)/4)*i + 4)
333      {
334        bufUniFactors.insert (LC (A, x));
335        henselLiftResume12 (A, bufUniFactors, (degree (A,y)/4)*i + 4,
336                            (degree (A, y)/4)*(i+1) + 4, Pi, diophant, M);
337        earlyFactors= earlyFactorDetection0 (A, bufUniFactors, newLiftBound,
338                      degs, earlySuccess, (degree (A, y)/4)*(i+1) + 4);
339      }
340      else
341      {
342        liftBound= newLiftBound;
343        break;
344      }
345      i++;
346    }
347    if (earlySuccess)
348      liftBound= newLiftBound;
349  }
350
351  return bufUniFactors;
352}
353
354CFList biFactorize (const CanonicalForm& F, const Variable& v)
355{
356  if (F.inCoeffDomain())
357    return CFList(F);
358
359  bool extension= (v.level() != 1);
360  CanonicalForm A;
361  if (isOn (SW_RATIONAL))
362    A= F*bCommonDen (F);
363  else
364    A= F;
365
366  CanonicalForm mipo;
367  if (extension)
368    mipo= getMipo (v);
369
370  if (A.isUnivariate())
371  {
372    CFFList buf;
373    if (extension)
374      buf= factorize (A, v);
375    else
376      buf= factorize (A, true);
377    CFList result= conv (buf);
378    if (result.getFirst().inCoeffDomain())
379      result.removeFirst();
380    return result;
381  }
382
383  CFMap N;
384  A= compress (A, N);
385  Variable y= A.mvar();
386
387  if (y.level() > 2) return CFList (F);
388  Variable x= Variable (1);
389
390  CanonicalForm contentAx= content (A, x);
391  CanonicalForm contentAy= content (A);
392
393  A= A/(contentAx*contentAy);
394  CFFList contentAxFactors, contentAyFactors;
395
396  if (extension)
397  {
398    if (!contentAx.inCoeffDomain())
399    {
400      contentAxFactors= factorize (contentAx, v);
401      if (contentAxFactors.getFirst().factor().inCoeffDomain())
402        contentAxFactors.removeFirst();
403    }
404    if (!contentAy.inCoeffDomain())
405    {
406      contentAyFactors= factorize (contentAy, v);
407      if (contentAyFactors.getFirst().factor().inCoeffDomain())
408        contentAyFactors.removeFirst();
409    }
410  }
411  else
412  {
413    if (!contentAx.inCoeffDomain())
414    {
415      contentAxFactors= factorize (contentAx, true);
416      if (contentAxFactors.getFirst().factor().inCoeffDomain())
417        contentAxFactors.removeFirst();
418    }
419    if (!contentAy.inCoeffDomain())
420    {
421      contentAyFactors= factorize (contentAy, true);
422      if (contentAyFactors.getFirst().factor().inCoeffDomain())
423        contentAyFactors.removeFirst();
424    }
425  }
426
427  //check trivial case
428  if (degree (A) == 1 || degree (A, 1) == 1)
429  {
430    CFList factors;
431    factors.append (A);
432
433    appendSwapDecompress (factors, conv (contentAxFactors),
434                          conv (contentAyFactors), false, false, N);
435
436    normalize (factors);
437    return factors;
438  }
439
440  //trivial case
441  CFList factors;
442  if (A.inCoeffDomain())
443  {
444    append (factors, conv (contentAxFactors));
445    append (factors, conv (contentAyFactors));
446    decompress (factors, N);
447    return factors;
448  }
449  else if (A.isUnivariate())
450  {
451    if (extension)
452      factors= conv (factorize (A, v));
453    else
454      factors= conv (factorize (A, true));
455    append (factors, conv (contentAxFactors));
456    append (factors, conv (contentAyFactors));
457    decompress (factors, N);
458    return factors;
459  }
460
461  bool swap= false;
462  if (degree (A) > degree (A, x))
463  {
464    A= swapvar (A, y, x);
465    swap= true;
466  }
467
468  CanonicalForm Aeval, bufAeval, buf;
469  CFList uniFactors, list, bufUniFactors;
470  DegreePattern degs;
471  DegreePattern bufDegs;
472
473  CanonicalForm Aeval2, bufAeval2;
474  CFList bufUniFactors2, list2, uniFactors2;
475  DegreePattern degs2;
476  DegreePattern bufDegs2;
477  bool swap2= false;
478
479  // several univariate factorizations to obtain more information about the
480  // degree pattern therefore usually less combinations have to be tried during
481  // the recombination process
482  int factorNums= 2;
483  int subCheck1= substituteCheck (A, x);
484  int subCheck2= substituteCheck (A, y);
485  buf= swapvar (A,x,y);
486  int evaluation, evaluation2, bufEvaluation= 0, bufEvaluation2= 0;
487  for (int i= 0; i < factorNums; i++)
488  {
489    bufAeval= A;
490    bufAeval= evalPoint (A, bufEvaluation);
491
492    bufAeval2= buf;
493    bufAeval2= evalPoint (buf, bufEvaluation2);
494
495    // univariate factorization
496    TIMING_START (uni_factorize);
497
498    if (extension)
499      bufUniFactors= conv (factorize (bufAeval, v));
500    else
501      bufUniFactors= conv (factorize (bufAeval, true));
502    TIMING_END_AND_PRINT (uni_factorize,
503                          "time for univariate factorization: ");
504    DEBOUTLN (cerr, "prod (bufUniFactors)== bufAeval " <<
505              (prod (bufUniFactors) == bufAeval));
506
507    TIMING_START (uni_factorize);
508    if (extension)
509      bufUniFactors2= conv (factorize (bufAeval2, v));
510    else
511      bufUniFactors2= conv (factorize (bufAeval2, true));
512    TIMING_END_AND_PRINT (uni_factorize,
513                          "time for univariate factorization in y: ");
514    DEBOUTLN (cerr, "prod (bufuniFactors2)== bufAeval2 " <<
515              (prod (bufUniFactors2) == bufAeval2));
516
517    if (bufUniFactors.getFirst().inCoeffDomain())
518      bufUniFactors.removeFirst();
519    if (bufUniFactors2.getFirst().inCoeffDomain())
520      bufUniFactors2.removeFirst();
521    if (bufUniFactors.length() == 1 || bufUniFactors2.length() == 1)
522    {
523      factors.append (A);
524
525      appendSwapDecompress (factors, conv (contentAxFactors),
526                            conv (contentAyFactors), swap, swap2, N);
527
528      if (isOn (SW_RATIONAL))
529        normalize (factors);
530      return factors;
531    }
532
533    if (i == 0)
534    {
535      if (subCheck1 > 0)
536      {
537        int subCheck= substituteCheck (bufUniFactors);
538
539        if (subCheck > 1 && (subCheck1%subCheck == 0))
540        {
541          CanonicalForm bufA= A;
542          subst (bufA, bufA, subCheck, x);
543          factors= biFactorize (bufA, v);
544          reverseSubst (factors, subCheck, x);
545          appendSwapDecompress (factors, conv (contentAxFactors),
546                                conv (contentAyFactors), swap, swap2, N);
547          if (isOn (SW_RATIONAL))
548            normalize (factors);
549          return factors;
550        }
551      }
552
553      if (subCheck2 > 0)
554      {
555        int subCheck= substituteCheck (bufUniFactors2);
556
557        if (subCheck > 1 && (subCheck2%subCheck == 0))
558        {
559          CanonicalForm bufA= A;
560          subst (bufA, bufA, subCheck, y);
561          factors= biFactorize (bufA, v);
562          reverseSubst (factors, subCheck, y);
563          appendSwapDecompress (factors, conv (contentAxFactors),
564                                conv (contentAyFactors), swap, swap2, N);
565          if (isOn (SW_RATIONAL))
566            normalize (factors);
567          return factors;
568        }
569      }
570    }
571
572    // degree analysis
573    bufDegs = DegreePattern (bufUniFactors);
574    bufDegs2= DegreePattern (bufUniFactors2);
575
576    if (i == 0)
577    {
578      Aeval= bufAeval;
579      evaluation= bufEvaluation;
580      uniFactors= bufUniFactors;
581      degs= bufDegs;
582      Aeval2= bufAeval2;
583      evaluation2= bufEvaluation2;
584      uniFactors2= bufUniFactors2;
585      degs2= bufDegs2;
586    }
587    else
588    {
589      degs.intersect (bufDegs);
590      degs2.intersect (bufDegs2);
591      if (bufUniFactors2.length() < uniFactors2.length())
592      {
593        uniFactors2= bufUniFactors2;
594        Aeval2= bufAeval2;
595        evaluation2= bufEvaluation2;
596      }
597      if (bufUniFactors.length() < uniFactors.length())
598      {
599        uniFactors= bufUniFactors;
600        Aeval= bufAeval;
601        evaluation= bufEvaluation;
602      }
603    }
604    if (bufEvaluation > 0)
605      bufEvaluation++;
606    else
607      bufEvaluation= -bufEvaluation + 1;
608    if (bufEvaluation > 0)
609      bufEvaluation2++;
610    else
611      bufEvaluation2= -bufEvaluation2 + 1;
612  }
613
614  if (uniFactors.length() > uniFactors2.length() ||
615      (uniFactors.length() == uniFactors2.length()
616       && degs.getLength() > degs2.getLength()))
617  {
618    degs= degs2;
619    uniFactors= uniFactors2;
620    evaluation= evaluation2;
621    Aeval= Aeval2;
622    A= buf;
623    swap2= true;
624  }
625
626  if (degs.getLength() == 1) // A is irreducible
627  {
628    factors.append (A);
629    appendSwapDecompress (factors, conv (contentAxFactors),
630                          conv (contentAyFactors), swap, swap2, N);
631    if (isOn (SW_RATIONAL))
632      normalize (factors);
633    return factors;
634  }
635
636  A= A (y + evaluation, y);
637
638  int liftBound= degree (A, y) + 1;
639
640  modpk b= modpk();
641  bool mipoHasDen= false;
642  if (!extension)
643  {
644    Off (SW_RATIONAL);
645    int i= 0;
646    findGoodPrime(F,i);
647    findGoodPrime(Aeval,i);
648    findGoodPrime(A,i);
649    if (i >= cf_getNumBigPrimes())
650      printf ("out of primes\n"); //TODO exit
651
652    int p=cf_getBigPrime(i);
653    b = coeffBound( A, p );
654    modpk bb= coeffBound (Aeval, p);
655    if (bb.getk() > b.getk() ) b=bb;
656      bb= coeffBound (F, p);
657    if (bb.getk() > b.getk() ) b=bb;
658  }
659  else
660  {
661    A /= Lc (Aeval);
662    // make factors elements of Z(a)[x] disable for modularDiophant
663    CanonicalForm multiplier= 1;
664    for (CFListIterator i= uniFactors; i.hasItem(); i++)
665    {
666      multiplier *= bCommonDen (i.getItem());
667      i.getItem()= i.getItem()*bCommonDen(i.getItem());
668    }
669    A *= multiplier;
670    A *= bCommonDen (A);
671
672    mipoHasDen= !bCommonDen(mipo).isOne();
673    mipo *= bCommonDen (mipo);
674    Off (SW_RATIONAL);
675    int i= 0;
676    ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
677    CanonicalForm discMipo= convertZZ2CF (discriminant (NTLmipo));
678    findGoodPrime (F*discMipo,i);
679    findGoodPrime (Aeval*discMipo,i);
680    findGoodPrime (A*discMipo,i);
681
682    int p=cf_getBigPrime(i);
683    b = coeffBound( A, p, mipo );
684    modpk bb= coeffBound (Aeval, p, mipo);
685    if (bb.getk() > b.getk() ) b=bb;
686      bb= coeffBound (F, p, mipo);
687    if (bb.getk() > b.getk() ) b=bb;
688  }
689
690  ExtensionInfo dummy= ExtensionInfo (false);
691  if (extension)
692    dummy= ExtensionInfo (v, false);
693  bool earlySuccess= false;
694  CFList earlyFactors;
695  TIMING_START (fac_hensel_lift);
696  uniFactors= henselLiftAndEarly
697              (A, earlySuccess, earlyFactors, degs, liftBound,
698               uniFactors, dummy, evaluation, b);
699  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
700  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
701
702  CanonicalForm MODl= power (y, liftBound);
703
704  if (mipoHasDen)
705  {
706    Variable vv;
707    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
708      if (hasFirstAlgVar (iter.getItem(), vv))
709        break;
710    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
711      iter.getItem()= replacevar (iter.getItem(), vv, v);
712  }
713
714  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
715                                uniFactors.length()/2, b);
716
717  On (SW_RATIONAL);
718
719  if (earlySuccess)
720    factors= Union (earlyFactors, factors);
721  else if (!earlySuccess && degs.getLength() == 1)
722    factors= earlyFactors;
723
724  for (CFListIterator i= factors; i.hasItem(); i++)
725    i.getItem()= i.getItem() (y - evaluation, y);
726
727  appendSwapDecompress (factors, conv (contentAxFactors),
728                        conv (contentAyFactors), swap, swap2, N);
729  if (isOn (SW_RATIONAL))
730    normalize (factors);
731
732  return factors;
733}
734
735#endif
Note: See TracBrowser for help on using the repository browser.