source: git/factory/facBivar.cc @ 7e3d56

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