source:git/factory/facBivar.cc@fd4146

spielwiese
Last change on this file since fd4146 was fd4146, checked in by Martin Lee <martinlee84@…>, 12 years ago
fix: bug in evaluation point search
• 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);
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  }
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  CanonicalForm multiplier= 1;
362  if (isOn (SW_RATIONAL))
363    A= F*bCommonDen (F);
364  else
365    A= F;
366
367  CanonicalForm mipo;
368  if (extension)
369    mipo= getMipo (v);
370
371  if (A.isUnivariate())
372  {
373    CFFList buf;
374    if (extension)
375      buf= factorize (A, v);
376    else
377      buf= factorize (A, true);
378    CFList result= conv (buf);
379    if (result.getFirst().inCoeffDomain())
380      result.removeFirst();
381    return result;
382  }
383
384  CFMap N;
385  A= compress (A, N);
386  Variable y= A.mvar();
387
388  if (y.level() > 2) return CFList (F);
389  Variable x= Variable (1);
390
391  CanonicalForm contentAx= content (A, x);
392  CanonicalForm contentAy= content (A);
393
394  A= A/(contentAx*contentAy);
395  CFFList contentAxFactors, contentAyFactors;
396
397  if (extension)
398  {
399    if (!contentAx.inCoeffDomain())
400    {
401      contentAxFactors= factorize (contentAx, v);
402      if (contentAxFactors.getFirst().factor().inCoeffDomain())
403        contentAxFactors.removeFirst();
404    }
405    if (!contentAy.inCoeffDomain())
406    {
407      contentAyFactors= factorize (contentAy, v);
408      if (contentAyFactors.getFirst().factor().inCoeffDomain())
409        contentAyFactors.removeFirst();
410    }
411  }
412  else
413  {
414    if (!contentAx.inCoeffDomain())
415    {
416      contentAxFactors= factorize (contentAx, true);
417      if (contentAxFactors.getFirst().factor().inCoeffDomain())
418        contentAxFactors.removeFirst();
419    }
420    if (!contentAy.inCoeffDomain())
421    {
422      contentAyFactors= factorize (contentAy, true);
423      if (contentAyFactors.getFirst().factor().inCoeffDomain())
424        contentAyFactors.removeFirst();
425    }
426  }
427
428  //check trivial case
429  if (degree (A) == 1 || degree (A, 1) == 1)
430  {
431    CFList factors;
432    factors.append (A);
433
434    appendSwapDecompress (factors, conv (contentAxFactors),
435                          conv (contentAyFactors), false, false, N);
436
437    normalize (factors);
438    return factors;
439  }
440
441  //trivial case
442  CFList factors;
443  if (A.inCoeffDomain())
444  {
445    append (factors, conv (contentAxFactors));
446    append (factors, conv (contentAyFactors));
447    decompress (factors, N);
448    return factors;
449  }
450  else if (A.isUnivariate())
451  {
452    if (extension)
453      factors= conv (factorize (A, v));
454    else
455      factors= conv (factorize (A, true));
456    append (factors, conv (contentAxFactors));
457    append (factors, conv (contentAyFactors));
458    decompress (factors, N);
459    return factors;
460  }
461
462  bool swap= false;
463  if (degree (A) > degree (A, x))
464  {
465    A= swapvar (A, y, x);
466    swap= true;
467  }
468
469  CanonicalForm Aeval, bufAeval, buf;
470  CFList uniFactors, list, bufUniFactors;
471  DegreePattern degs;
472  DegreePattern bufDegs;
473
474  CanonicalForm Aeval2, bufAeval2;
475  CFList bufUniFactors2, list2, uniFactors2;
476  DegreePattern degs2;
477  DegreePattern bufDegs2;
478  bool swap2= false;
479
481  // degree pattern therefore usually less combinations have to be tried during
482  // the recombination process
483  int factorNums= 2;
484  int subCheck1= substituteCheck (A, x);
485  int subCheck2= substituteCheck (A, y);
486  buf= swapvar (A,x,y);
487  int evaluation, evaluation2, bufEvaluation= 0, bufEvaluation2= 0;
488  for (int i= 0; i < factorNums; i++)
489  {
490    bufAeval= A;
491    bufAeval= evalPoint (A, bufEvaluation);
492
493    bufAeval2= buf;
494    bufAeval2= evalPoint (buf, bufEvaluation2);
495
496    // univariate factorization
497    TIMING_START (uni_factorize);
498
499    if (extension)
500      bufUniFactors= conv (factorize (bufAeval, v));
501    else
502      bufUniFactors= conv (factorize (bufAeval, true));
503    TIMING_END_AND_PRINT (uni_factorize,
504                          "time for univariate factorization: ");
505    DEBOUTLN (cerr, "Lc (bufAeval)*prod (bufUniFactors)== bufAeval " <<
506              (prod (bufUniFactors)*Lc (bufAeval) == bufAeval));
507
508    TIMING_START (uni_factorize);
509    if (extension)
510      bufUniFactors2= conv (factorize (bufAeval2, v));
511    else
512      bufUniFactors2= conv (factorize (bufAeval2, true));
513    TIMING_END_AND_PRINT (uni_factorize,
514                          "time for univariate factorization in y: ");
515    DEBOUTLN (cerr, "Lc (Aeval2)*prod (uniFactors2)== Aeval2 " <<
516              (prod (bufUniFactors2)*Lc (bufAeval2) == bufAeval2));
517
518    if (bufUniFactors.getFirst().inCoeffDomain())
519      bufUniFactors.removeFirst();
520    if (bufUniFactors2.getFirst().inCoeffDomain())
521      bufUniFactors2.removeFirst();
522    if (bufUniFactors.length() == 1 || bufUniFactors2.length() == 1)
523    {
524      factors.append (A);
525
526      appendSwapDecompress (factors, conv (contentAxFactors),
527                            conv (contentAyFactors), swap, swap2, N);
528
529      if (isOn (SW_RATIONAL))
530        normalize (factors);
531      return factors;
532    }
533
534    if (i == 0)
535    {
536      if (subCheck1 > 0)
537      {
538        int subCheck= substituteCheck (bufUniFactors);
539
540        if (subCheck > 1 && (subCheck1%subCheck == 0))
541        {
542          CanonicalForm bufA= A;
543          subst (bufA, bufA, subCheck, x);
544          factors= biFactorize (bufA, v);
545          reverseSubst (factors, subCheck, x);
546          appendSwapDecompress (factors, conv (contentAxFactors),
547                                conv (contentAyFactors), swap, swap2, N);
548          if (isOn (SW_RATIONAL))
549            normalize (factors);
550          return factors;
551        }
552      }
553
554      if (subCheck2 > 0)
555      {
556        int subCheck= substituteCheck (bufUniFactors2);
557
558        if (subCheck > 1 && (subCheck2%subCheck == 0))
559        {
560          CanonicalForm bufA= A;
561          subst (bufA, bufA, subCheck, y);
562          factors= biFactorize (bufA, v);
563          reverseSubst (factors, subCheck, y);
564          appendSwapDecompress (factors, conv (contentAxFactors),
565                                conv (contentAyFactors), swap, swap2, N);
566          if (isOn (SW_RATIONAL))
567            normalize (factors);
568          return factors;
569        }
570      }
571    }
572
573    // degree analysis
574    bufDegs = DegreePattern (bufUniFactors);
575    bufDegs2= DegreePattern (bufUniFactors2);
576
577    if (i == 0)
578    {
579      Aeval= bufAeval;
580      evaluation= bufEvaluation;
581      uniFactors= bufUniFactors;
582      degs= bufDegs;
583      Aeval2= bufAeval2;
584      evaluation2= bufEvaluation2;
585      uniFactors2= bufUniFactors2;
586      degs2= bufDegs2;
587    }
588    else
589    {
590      degs.intersect (bufDegs);
591      degs2.intersect (bufDegs2);
592      if (bufUniFactors2.length() < uniFactors2.length())
593      {
594        uniFactors2= bufUniFactors2;
595        Aeval2= bufAeval2;
596        evaluation2= bufEvaluation2;
597      }
598      if (bufUniFactors.length() < uniFactors.length())
599      {
600        uniFactors= bufUniFactors;
601        Aeval= bufAeval;
602        evaluation= bufEvaluation;
603      }
604    }
605    if (bufEvaluation > 0)
606      bufEvaluation++;
607    else
608      bufEvaluation= -bufEvaluation + 1;
609    if (bufEvaluation > 0)
610      bufEvaluation2++;
611    else
612      bufEvaluation2= -bufEvaluation2 + 1;
613  }
614
615  if (uniFactors.length() > uniFactors2.length() ||
616      (uniFactors.length() == uniFactors2.length()
617       && degs.getLength() > degs2.getLength()))
618  {
619    degs= degs2;
620    uniFactors= uniFactors2;
621    evaluation= evaluation2;
622    Aeval= Aeval2;
623    A= buf;
624    swap2= true;
625  }
626
627  if (degs.getLength() == 1) // A is irreducible
628  {
629    factors.append (A);
630    appendSwapDecompress (factors, conv (contentAxFactors),
631                          conv (contentAyFactors), swap, swap2, N);
632    if (isOn (SW_RATIONAL))
633      normalize (factors);
634    return factors;
635  }
636
637  A= A (y + evaluation, y);
638
639  int liftBound= degree (A, y) + 1;
640
641  modpk b= modpk();
642  bool mipoHasDen= false;
643  if (!extension)
644  {
645    Off (SW_RATIONAL);
646    int i= 0;
647    findGoodPrime(F,i);
648    findGoodPrime(Aeval,i);
649    findGoodPrime(A,i);
650    if (i >= cf_getNumBigPrimes())
651      printf ("out of primes\n"); //TODO exit
652
653    int p=cf_getBigPrime(i);
654    b = coeffBound( A, p );
655    modpk bb= coeffBound (Aeval, p);
656    if (bb.getk() > b.getk() ) b=bb;
657      bb= coeffBound (F, p);
658    if (bb.getk() > b.getk() ) b=bb;
659  }
660  else
661  {
662    A /= Lc (Aeval);
663    A *= bCommonDen (A);
664    // make factors elements of Z(a)[x] disable for modularDiophant
665    for (CFListIterator i= uniFactors; i.hasItem(); i++)
666      i.getItem()= i.getItem()*bCommonDen(i.getItem());
667    mipoHasDen= !bCommonDen(mipo).isOne();
668    mipo *= bCommonDen (mipo);
669    Off (SW_RATIONAL);
670    int i= 0;
671    ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
672    CanonicalForm discMipo= convertZZ2CF (discriminant (NTLmipo));
673    findGoodPrime (F*discMipo,i);
674    findGoodPrime (Aeval*discMipo,i);
675    findGoodPrime (A*discMipo,i);
676
677    int p=cf_getBigPrime(i);
678    b = coeffBound( A, p, mipo );
679    modpk bb= coeffBound (Aeval, p, mipo);
680    if (bb.getk() > b.getk() ) b=bb;
681      bb= coeffBound (F, p, mipo);
682    if (bb.getk() > b.getk() ) b=bb;
683  }
684
685  ExtensionInfo dummy= ExtensionInfo (false);
686  if (extension)
687    dummy= ExtensionInfo (v, false);
688  bool earlySuccess= false;
689  CFList earlyFactors;
690  TIMING_START (fac_hensel_lift);
691  uniFactors= henselLiftAndEarly
692              (A, earlySuccess, earlyFactors, degs, liftBound,
693               uniFactors, dummy, evaluation, b);
694  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
695  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
696
697  CanonicalForm MODl= power (y, liftBound);
698
699  if (mipoHasDen)
700  {
701    Variable vv;
702    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
703      if (hasFirstAlgVar (iter.getItem(), vv))
704        break;
705    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
706      iter.getItem()= replacevar (iter.getItem(), vv, v);
707  }
708
709  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
710                                uniFactors.length()/2, b);
711
712  On (SW_RATIONAL);
713
714  if (earlySuccess)
715    factors= Union (earlyFactors, factors);
716  else if (!earlySuccess && degs.getLength() == 1)
717    factors= earlyFactors;
718
719  for (CFListIterator i= factors; i.hasItem(); i++)
720    i.getItem()= i.getItem() (y - evaluation, y);
721
722  appendSwapDecompress (factors, conv (contentAxFactors),
723                        conv (contentAyFactors), swap, swap2, N);
724  if (isOn (SW_RATIONAL))
725    normalize (factors);
726
727  return factors;
728}
729
730#endif
Note: See TracBrowser for help on using the repository browser.