source: git/factory/facBivar.cc @ f9da5e

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