source: git/factory/facBivar.cc @ 4e6d2a

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