source: git/factory/facBivar.cc @ 767da0

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