source: git/factory/facBivar.cc @ 03c742

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