source: git/factory/facBivar.cc @ f5d2647

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