source: git/factory/facBivar.cc @ 5f92d8

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