source: git/factory/facBivar.cc @ 70c40f

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