source: git/factory/facBivar.cc @ 2a95b2

fieker-DuValspielwiese
Last change on this file since 2a95b2 was afbebe, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: deleted some unused code chg: joined precomputeLeadingCoeff from facFactorize.cc and facFqFactorize.cc
  • Property mode set to 100644
File size: 14.4 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 "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)
30
31// bound on coeffs of f (cf. Musser: Multivariate Polynomial Factorization,
32//                          Gelfond: Transcendental and Algebraic Numbers)
33modpk
34coeffBound ( const CanonicalForm & f, int p )
35{
36    int * degs = degrees( f );
37    int M = 0, i, k = f.level();
38    CanonicalForm b= 1;
39    for ( i = 1; i <= k; i++ )
40    {
41      M += degs[i];
42      b *= degs[i] + 1;
43    }
44    b /= power (CanonicalForm (2), k);
45    b= b.sqrt() + 1;
46    b *= 2 * maxNorm( f ) * power( CanonicalForm( 2 ), M );
47    CanonicalForm B = p;
48    k = 1;
49    while ( B < b ) {
50        B *= p;
51        k++;
52    }
53    return modpk( p, k );
54}
55
56void findGoodPrime(const CanonicalForm &f, int &start)
57{
58  if (! f.inBaseDomain() )
59  {
60    CFIterator i = f;
61    for(;;)
62    {
63      if  ( i.hasTerms() )
64      {
65        findGoodPrime(i.coeff(),start);
66        if (0==cf_getBigPrime(start)) return;
67        if((i.exp()!=0) && ((i.exp() % cf_getBigPrime(start))==0))
68        {
69          start++;
70          i=f;
71        }
72        else  i++;
73      }
74      else break;
75    }
76  }
77  else
78  {
79    if (f.inZ())
80    {
81      if (0==cf_getBigPrime(start)) return;
82      while((!f.isZero()) && (mod(f,cf_getBigPrime(start))==0))
83      {
84        start++;
85        if (0==cf_getBigPrime(start)) return;
86      }
87    }
88  }
89}
90
91modpk
92coeffBound ( const CanonicalForm & f, int p, const CanonicalForm& mipo )
93{
94    int * degs = degrees( f );
95    int M = 0, i, k = f.level();
96    CanonicalForm K= 1;
97    for ( i = 1; i <= k; i++ )
98    {
99        M += degs[i];
100        K *= degs[i] + 1;
101    }
102    K /= power (CanonicalForm (2), k);
103    K= K.sqrt()+1;
104    K *= power (CanonicalForm (2), M);
105    int N= degree (mipo);
106    CanonicalForm b;
107    b= 2*power (maxNorm (f), N)*power (maxNorm (mipo), 4*N)*K*
108       power (CanonicalForm (2), N)*(CanonicalForm (M+1).sqrt()+1)*
109       power (CanonicalForm (N+1).sqrt()+1, 7*N);
110    b /= power (abs (lc (mipo)), N);
111
112    ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
113    ZZX NTLLcf= convertFacCF2NTLZZX (Lc (f));
114    ZZ NTLf= resultant (NTLmipo, NTLLcf);
115    ZZ NTLD= discriminant (NTLmipo);
116    b /= abs (convertZZ2CF (NTLf))*abs (convertZZ2CF (NTLD));
117
118    CanonicalForm B = p;
119    k = 1;
120    while ( B < b ) {
121        B *= p;
122        k++;
123    }
124    return modpk( p, k );
125}
126
127CFList conv (const CFFList& L)
128{
129  CFList result;
130  for (CFFListIterator i= L; i.hasItem(); i++)
131    result.append (i.getItem().factor());
132  return result;
133}
134
135bool testPoint (const CanonicalForm& F, CanonicalForm& G, int i)
136{
137  G= F (i, 2);
138  if (G.inCoeffDomain() || degree (F, 1) > degree (G, 1))
139    return false;
140
141  if (degree (gcd (deriv (G, G.mvar()), G)) > 0)
142    return false;
143  return true;
144}
145
146CanonicalForm evalPoint (const CanonicalForm& F, int& i)
147{
148  Variable x= Variable (1);
149  Variable y= Variable (2);
150  CanonicalForm result;
151
152  int k;
153
154  if (i == 0)
155  {
156    if (testPoint (F, result, i))
157      return result;
158  }
159  do
160  {
161    if (i > 0)
162      k= 1;
163    else
164      k= 2;
165    while (k < 3)
166    {
167      if (k == 1)
168      {
169        if (testPoint (F, result, i))
170          return result;
171      }
172      else
173      {
174        if (testPoint (F, result, -i))
175        {
176          i= -i;
177          return result;
178        }
179        else if (i < 0)
180          i= -i;
181      }
182      k++;
183    }
184    i++;
185  } while (1);
186}
187
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 && gcd (degree (A), degree (A,1)).isOne()))
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    bufAeval= evalPoint (A, bufEvaluation);
337
338    bufAeval2= buf;
339    bufAeval2= evalPoint (buf, bufEvaluation2);
340
341    // univariate factorization
342    TIMING_START (uni_factorize);
343
344    if (extension)
345      bufUniFactors= conv (factorize (bufAeval, v));
346    else
347      bufUniFactors= conv (factorize (bufAeval, true));
348    TIMING_END_AND_PRINT (uni_factorize,
349                          "time for univariate factorization: ");
350    DEBOUTLN (cerr, "prod (bufUniFactors)== bufAeval " <<
351              (prod (bufUniFactors) == bufAeval));
352
353    TIMING_START (uni_factorize);
354    if (extension)
355      bufUniFactors2= conv (factorize (bufAeval2, v));
356    else
357      bufUniFactors2= conv (factorize (bufAeval2, true));
358    TIMING_END_AND_PRINT (uni_factorize,
359                          "time for univariate factorization in y: ");
360    DEBOUTLN (cerr, "prod (bufuniFactors2)== bufAeval2 " <<
361              (prod (bufUniFactors2) == bufAeval2));
362
363    if (bufUniFactors.getFirst().inCoeffDomain())
364      bufUniFactors.removeFirst();
365    if (bufUniFactors2.getFirst().inCoeffDomain())
366      bufUniFactors2.removeFirst();
367    if (bufUniFactors.length() == 1 || bufUniFactors2.length() == 1)
368    {
369      factors.append (A);
370
371      appendSwapDecompress (factors, conv (contentAxFactors),
372                            conv (contentAyFactors), swap, swap2, N);
373
374      if (isOn (SW_RATIONAL))
375        normalize (factors);
376      return factors;
377    }
378
379    if (i == 0)
380    {
381      if (subCheck1 > 0)
382      {
383        int subCheck= substituteCheck (bufUniFactors);
384
385        if (subCheck > 1 && (subCheck1%subCheck == 0))
386        {
387          CanonicalForm bufA= A;
388          subst (bufA, bufA, subCheck, x);
389          factors= biFactorize (bufA, v);
390          reverseSubst (factors, subCheck, x);
391          appendSwapDecompress (factors, conv (contentAxFactors),
392                                conv (contentAyFactors), swap, swap2, N);
393          if (isOn (SW_RATIONAL))
394            normalize (factors);
395          return factors;
396        }
397      }
398
399      if (subCheck2 > 0)
400      {
401        int subCheck= substituteCheck (bufUniFactors2);
402
403        if (subCheck > 1 && (subCheck2%subCheck == 0))
404        {
405          CanonicalForm bufA= A;
406          subst (bufA, bufA, subCheck, y);
407          factors= biFactorize (bufA, v);
408          reverseSubst (factors, subCheck, y);
409          appendSwapDecompress (factors, conv (contentAxFactors),
410                                conv (contentAyFactors), swap, swap2, N);
411          if (isOn (SW_RATIONAL))
412            normalize (factors);
413          return factors;
414        }
415      }
416    }
417
418    // degree analysis
419    bufDegs = DegreePattern (bufUniFactors);
420    bufDegs2= DegreePattern (bufUniFactors2);
421
422    if (i == 0)
423    {
424      Aeval= bufAeval;
425      evaluation= bufEvaluation;
426      uniFactors= bufUniFactors;
427      degs= bufDegs;
428      Aeval2= bufAeval2;
429      evaluation2= bufEvaluation2;
430      uniFactors2= bufUniFactors2;
431      degs2= bufDegs2;
432    }
433    else
434    {
435      degs.intersect (bufDegs);
436      degs2.intersect (bufDegs2);
437      if (bufUniFactors2.length() < uniFactors2.length())
438      {
439        uniFactors2= bufUniFactors2;
440        Aeval2= bufAeval2;
441        evaluation2= bufEvaluation2;
442      }
443      if (bufUniFactors.length() < uniFactors.length())
444      {
445        uniFactors= bufUniFactors;
446        Aeval= bufAeval;
447        evaluation= bufEvaluation;
448      }
449    }
450    if (bufEvaluation > 0)
451      bufEvaluation++;
452    else
453      bufEvaluation= -bufEvaluation + 1;
454    if (bufEvaluation > 0)
455      bufEvaluation2++;
456    else
457      bufEvaluation2= -bufEvaluation2 + 1;
458  }
459
460  if (uniFactors.length() > uniFactors2.length() ||
461      (uniFactors.length() == uniFactors2.length()
462       && degs.getLength() > degs2.getLength()))
463  {
464    degs= degs2;
465    uniFactors= uniFactors2;
466    evaluation= evaluation2;
467    Aeval= Aeval2;
468    A= buf;
469    swap2= true;
470  }
471
472  if (degs.getLength() == 1) // A is irreducible
473  {
474    factors.append (A);
475    appendSwapDecompress (factors, conv (contentAxFactors),
476                          conv (contentAyFactors), swap, swap2, N);
477    if (isOn (SW_RATIONAL))
478      normalize (factors);
479    return factors;
480  }
481
482  A *= bCommonDen (A);
483  A= A (y + evaluation, y);
484
485  int liftBound= degree (A, y) + 1;
486
487  modpk b= modpk();
488  bool mipoHasDen= false;
489  if (!extension)
490  {
491    Off (SW_RATIONAL);
492    int i= 0;
493    findGoodPrime(F,i);
494    findGoodPrime(Aeval,i);
495    findGoodPrime(A,i);
496    if (i >= cf_getNumBigPrimes())
497      printf ("out of primes\n"); //TODO exit
498
499    int p=cf_getBigPrime(i);
500    b = coeffBound( A, p );
501    modpk bb= coeffBound (Aeval, p);
502    if (bb.getk() > b.getk() ) b=bb;
503      bb= coeffBound (F, p);
504    if (bb.getk() > b.getk() ) b=bb;
505  }
506  else
507  {
508    A /= Lc (Aeval);
509    // make factors elements of Z(a)[x] disable for modularDiophant
510    CanonicalForm multiplier= 1;
511    for (CFListIterator i= uniFactors; i.hasItem(); i++)
512    {
513      multiplier *= bCommonDen (i.getItem());
514      i.getItem()= i.getItem()*bCommonDen(i.getItem());
515    }
516    A *= multiplier;
517    A *= bCommonDen (A);
518
519    mipoHasDen= !bCommonDen(mipo).isOne();
520    mipo *= bCommonDen (mipo);
521    Off (SW_RATIONAL);
522    int i= 0;
523    ZZX NTLmipo= convertFacCF2NTLZZX (mipo);
524    CanonicalForm discMipo= convertZZ2CF (discriminant (NTLmipo));
525    findGoodPrime (F*discMipo,i);
526    findGoodPrime (Aeval*discMipo,i);
527    findGoodPrime (A*discMipo,i);
528
529    int p=cf_getBigPrime(i);
530    b = coeffBound( A, p, mipo );
531    modpk bb= coeffBound (Aeval, p, mipo);
532    if (bb.getk() > b.getk() ) b=bb;
533      bb= coeffBound (F, p, mipo);
534    if (bb.getk() > b.getk() ) b=bb;
535  }
536
537  ExtensionInfo dummy= ExtensionInfo (false);
538  if (extension)
539    dummy= ExtensionInfo (v, false);
540  bool earlySuccess= false;
541  CFList earlyFactors;
542  TIMING_START (fac_bi_hensel_lift);
543  uniFactors= henselLiftAndEarly
544              (A, earlySuccess, earlyFactors, degs, liftBound,
545               uniFactors, dummy, evaluation, b);
546  TIMING_END_AND_PRINT (fac_bi_hensel_lift, "time for hensel lifting: ");
547  DEBOUTLN (cerr, "lifted factors= " << uniFactors);
548
549  CanonicalForm MODl= power (y, liftBound);
550
551  if (mipoHasDen)
552  {
553    Variable vv;
554    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
555      if (hasFirstAlgVar (iter.getItem(), vv))
556        break;
557    for (CFListIterator iter= uniFactors; iter.hasItem(); iter++)
558      iter.getItem()= replacevar (iter.getItem(), vv, v);
559  }
560
561  On (SW_RATIONAL);
562  A *= bCommonDen (A);
563  Off (SW_RATIONAL);
564
565  factors= factorRecombination (uniFactors, A, MODl, degs, 1,
566                                uniFactors.length()/2, b);
567
568  On (SW_RATIONAL);
569
570  if (earlySuccess)
571    factors= Union (earlyFactors, factors);
572  else if (!earlySuccess && degs.getLength() == 1)
573    factors= earlyFactors;
574
575  for (CFListIterator i= factors; i.hasItem(); i++)
576    i.getItem()= i.getItem() (y - evaluation, y);
577
578  appendSwapDecompress (factors, conv (contentAxFactors),
579                        conv (contentAyFactors), swap, swap2, N);
580  if (isOn (SW_RATIONAL))
581    normalize (factors);
582
583  return factors;
584}
585
586#endif
Note: See TracBrowser for help on using the repository browser.