source: git/factory/facFactorize.cc @ 5295f9

spielwiese
Last change on this file since 5295f9 was 5295f9, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: more heuristics to detect right leading coeff
  • Property mode set to 100644
File size: 35.9 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqFactorize.cc
5 *
6 * multivariate 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 "facFqFactorizeUtil.h"
21#include "facFactorize.h"
22#include "facFqFactorize.h"
23#include "cf_random.h"
24#include "facHensel.h"
25#include "cf_gcd_smallp.h"
26#include "cf_map_ext.h"
27#include "algext.h"
28#include "cf_reval.h"
29#include "facSparseHensel.h"
30
31TIMING_DEFINE_PRINT(fac_bi_factorize)
32TIMING_DEFINE_PRINT(fac_hensel_lift)
33TIMING_DEFINE_PRINT(fac_factor_recombination)
34
35#ifdef HAVE_NTL
36CFList evalPoints (const CanonicalForm& F, CFList& eval, Evaluation& E)
37{
38  CFList result;
39  Variable x= Variable (1);
40
41  bool found= false;
42  bool allZero= true;
43  bool foundZero= false;
44  CanonicalForm deriv_x, gcd_deriv;
45  CFListIterator iter;
46  do
47  {
48    eval.insert (F);
49    bool bad= false;
50    for (int i= E.max(); i >= E.min(); i--)
51    {
52      eval.insert (eval.getFirst()( E [i], i));
53      result.append (E[i]);
54      if (!E[i].isZero())
55        allZero= false;
56      else
57        foundZero= true;
58      if (!allZero && foundZero)
59      {
60        result= CFList();
61        eval= CFList();
62        bad= true;
63        foundZero= false;
64        break;
65      }
66      if (degree (eval.getFirst(), i - 1) != degree (F, i - 1))
67      {
68        result= CFList();
69        eval= CFList();
70        bad= true;
71        break;
72      }
73    }
74
75    if (bad)
76    {
77      E.nextpoint();
78      continue;
79    }
80
81    if (degree (eval.getFirst()) != degree (F, 1))
82    {
83      result= CFList();
84      eval= CFList();
85      E.nextpoint();
86      continue;
87    }
88
89    deriv_x= deriv (eval.getFirst(), x);
90    gcd_deriv= gcd (eval.getFirst(), deriv_x);
91    if (degree (gcd_deriv) > 0)
92    {
93      result= CFList();
94      eval= CFList();
95      E.nextpoint();
96      continue;
97    }
98    iter= eval;
99    iter++;
100    CanonicalForm contentx= content (iter.getItem(), x);
101    if (degree (contentx) > 0)
102    {
103      result= CFList();
104      eval= CFList();
105      E.nextpoint();
106      continue;
107    }
108    found= true;
109  }
110  while (!found);
111
112  if (!eval.isEmpty())
113    eval.removeFirst();
114  return result;
115}
116
117void
118factorizationWRTDifferentSecondVars (const CanonicalForm& A, CFList*& Aeval,
119                                     int& minFactorsLength, bool& irred,
120                                     const Variable& w)
121{
122  Variable x= Variable (1);
123  minFactorsLength= 0;
124  irred= false;
125  Variable v;
126  CFList factors;
127  CanonicalForm LCA= LC (A,1);
128  if (!LCA.inCoeffDomain())
129  {
130    for (int j= 0; j < A.level() - 2; j++)
131    {
132      if (!Aeval[j].isEmpty() && (degree (LCA, j+3) > 0))
133      {
134        v= Variable (Aeval[j].getFirst().level());
135
136        factors= ratBiSqrfFactorize (Aeval[j].getFirst(), w);
137
138        if (factors.getFirst().inCoeffDomain())
139          factors.removeFirst();
140
141        if (minFactorsLength == 0)
142          minFactorsLength= factors.length();
143        else
144          minFactorsLength= tmin (minFactorsLength, factors.length());
145
146        if (factors.length() == 1)
147        {
148          irred= true;
149          return;
150        }
151        sortList (factors, x);
152        Aeval [j]= factors;
153      }
154      else if (!Aeval[j].isEmpty())
155      {
156        Aeval[j]=CFList();
157      }
158    }
159  }
160  else
161  {
162    for (int j= 0; j < A.level() - 2; j++)
163      Aeval[j]= CFList();
164  }
165}
166
167int
168testFactors (const CanonicalForm& G, const CFList& uniFactors,
169             CanonicalForm& sqrfPartF, CFList& factors,
170             CFFList*& bufSqrfFactors, CFList& evalSqrfPartF,
171             const CFArray& evalPoint)
172{
173  CanonicalForm tmp;
174  CFListIterator j;
175  for (CFListIterator i= uniFactors; i.hasItem(); i++)
176  {
177    tmp= i.getItem();
178    if (i.hasItem())
179      i++;
180    else
181      break;
182    for (j= i; j.hasItem(); j++)
183    {
184      if (tmp == j.getItem())
185        return 0;
186    }
187  }
188
189  CanonicalForm F= G;
190  CFFList sqrfFactorization= sqrFree (F);
191
192  sqrfPartF= 1;
193  for (CFFListIterator i= sqrfFactorization; i.hasItem(); i++)
194    sqrfPartF *= i.getItem().factor();
195
196  evalSqrfPartF= evaluateAtEval (sqrfPartF, evalPoint);
197
198  CanonicalForm test= evalSqrfPartF.getFirst() (evalPoint[0], 2);
199
200  if (degree (test) != degree (sqrfPartF, 1) || test.inCoeffDomain())
201    return 0;
202
203  CFFList sqrfFactors;
204  CFList tmp2;
205  int k= 0;
206  factors= uniFactors;
207  CFFListIterator iter;
208  for (CFListIterator i= factors; i.hasItem(); i++, k++)
209  {
210    tmp= 1;
211    sqrfFactors= sqrFree (i.getItem());
212
213    for (iter= sqrfFactors; iter.hasItem(); iter++)
214    {
215      tmp2.append (iter.getItem().factor());
216      tmp *= iter.getItem().factor();
217    }
218    i.getItem()= tmp/Lc(tmp);
219    bufSqrfFactors [k]= sqrfFactors;
220  }
221
222  for (int i= 0; i < factors.length() - 1; i++)
223  {
224    for (int k= i + 1; k < factors.length(); k++)
225    {
226      gcdFreeBasis (bufSqrfFactors [i], bufSqrfFactors[k]);
227    }
228  }
229
230  factors= CFList();
231  for (int i= 0; i < uniFactors.length(); i++)
232  {
233    if (i == 0)
234    {
235      for (iter= bufSqrfFactors [i]; iter.hasItem(); iter++)
236      {
237        if (iter.getItem().factor().inCoeffDomain())
238          continue;
239        iter.getItem()= CFFactor (iter.getItem().factor()/
240                                  Lc (iter.getItem().factor()),
241                                  iter.getItem().exp());
242        factors.append (iter.getItem().factor());
243      }
244    }
245    else
246    {
247      for (iter= bufSqrfFactors [i]; iter.hasItem(); iter++)
248      {
249        if (iter.getItem().factor().inCoeffDomain())
250          continue;
251        iter.getItem()= CFFactor (iter.getItem().factor()/
252                               Lc (iter.getItem().factor()),
253                               iter.getItem().exp());
254        if (!find (factors, iter.getItem().factor()))
255          factors.append (iter.getItem().factor());
256      }
257    }
258  }
259
260  test= prod (factors);
261  tmp= evalSqrfPartF.getFirst() (evalPoint[0],2);
262  if (test/Lc (test) != tmp/Lc (tmp))
263    return 0;
264  else
265    return 1;
266}
267
268CFList
269precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
270                        const CFList& evaluation,CFList*& differentSecondVarLCs,
271                        int length, Variable& y
272                       )
273{
274  y= Variable (1);
275  if (LCF.inCoeffDomain())
276  {
277    CFList result;
278    for (int i= 1; i <= LCFFactors.length() + 1; i++)
279      result.append (1);
280    return result;
281  }
282
283  CFMap N, M;
284  CFArray dummy= CFArray (2);
285  dummy [0]= LCF;
286  dummy [1]= Variable (2);
287  compress (dummy, M, N);
288  CanonicalForm F= M (LCF);
289  if (LCF.isUnivariate())
290  {
291    CFList result;
292    int LCFLevel= LCF.level();
293    bool found= false;
294    if (LCFLevel == 2)
295    {
296    //bivariate leading coefficients are already the true leading coefficients
297      result= LCFFactors;
298      found= true;
299    }
300    else
301    {
302      CFListIterator j;
303      for (int i= 0; i < length; i++)
304      {
305        for (j= differentSecondVarLCs[i]; j.hasItem(); j++)
306        {
307          if (j.getItem().level() == LCFLevel)
308          {
309            found= true;
310            break;
311          }
312        }
313        if (found)
314        {
315          result= differentSecondVarLCs [i];
316          break;
317        }
318      }
319      if (!found)
320        result= LCFFactors;
321    }
322    if (found)
323      result.insert (Lc (LCF));
324    else
325      result.append (LCF);
326    return result;
327  }
328
329  CFList factors= LCFFactors;
330
331  for (CFListIterator i= factors; i.hasItem(); i++)
332    i.getItem()= M (i.getItem());
333
334  CanonicalForm sqrfPartF;
335  CFFList * bufSqrfFactors= new CFFList [factors.length()];
336  CFList evalSqrfPartF, bufFactors;
337  CFArray evalPoint= CFArray (evaluation.length() - 1);
338  CFArray buf= CFArray (evaluation.length());
339  CFArray swap= CFArray (evaluation.length());
340  CFListIterator iter= evaluation;
341  CanonicalForm vars=getVars (LCF)*Variable (2);
342  for (int i= evaluation.length() +1; i > 1; i--, iter++)
343  {
344    buf[i-2]=iter.getItem();
345    if (degree (vars, i) > 0)
346      swap[M(Variable (i)).level()-1]=buf[i-2];
347  }
348  buf= swap;
349  for (int i= 0; i < evaluation.length() - 1; i++)
350    evalPoint[i]= buf[i+1];
351
352  //TODO sqrfPartF einmal berechnen nicht stÀndig
353  int pass= testFactors (F, factors, sqrfPartF,
354                         bufFactors, bufSqrfFactors, evalSqrfPartF, evalPoint);
355
356  bool foundDifferent= false;
357  Variable z;
358  Variable x= y;
359  int j= 0;
360  if (!pass)
361  {
362    int lev= 0;
363    CanonicalForm bufF;
364    CFList bufBufFactors;
365    for (int i= 0; i < length; i++)
366    {
367      if (!differentSecondVarLCs [i].isEmpty())
368      {
369        bool allConstant= true;
370        for (iter= differentSecondVarLCs[i]; iter.hasItem(); iter++)
371        {
372          if (!iter.getItem().inCoeffDomain())
373          {
374            allConstant= false;
375            y= Variable (iter.getItem().level());
376            lev= M(y).level();
377          }
378        }
379        if (allConstant)
380          continue;
381
382        bufFactors= differentSecondVarLCs [i];
383        for (iter= bufFactors; iter.hasItem(); iter++)
384          iter.getItem()= swapvar (iter.getItem(), x, y);
385        bufF= F;
386        z= Variable (lev);
387        bufF= swapvar (bufF, x, z);
388        bufBufFactors= bufFactors;
389        evalPoint= CFArray (evaluation.length() - 1);
390        for (int k= 0; k < evaluation.length()-1; k++)
391        {
392          if (N (Variable (k+1)).level() != y.level())
393            evalPoint[k]= buf[k+1];
394          else
395            evalPoint[k]= buf[0];
396        }
397        pass= testFactors (bufF, bufBufFactors, sqrfPartF, bufFactors,
398                           bufSqrfFactors, evalSqrfPartF, evalPoint);
399        if (pass)
400        {
401          foundDifferent= true;
402          F= bufF;
403          CFList l= factors;
404          for (iter= l; iter.hasItem(); iter++)
405            iter.getItem()= swapvar (iter.getItem(), x, y);
406          differentSecondVarLCs [i]= l;
407          j= i;
408          break;
409        }
410        if (!pass && i == length - 1)
411        {
412          CFList result;
413          result.append (LCF);
414          for (int j= 1; j <= factors.length(); j++)
415            result.append (LCF);
416          y= Variable (1);
417          delete [] bufSqrfFactors;
418          return result;
419        }
420      }
421    }
422  }
423  if (!pass)
424  {
425    CFList result;
426    result.append (LCF);
427    for (int j= 1; j <= factors.length(); j++)
428      result.append (LCF);
429    y= Variable (1);
430    delete [] bufSqrfFactors;
431    return result;
432  }
433  else
434    factors= bufFactors;
435
436
437  bufFactors= factors;
438  CFList evaluation2;
439  for (int i= 0; i < F.level()-1; i++)
440    evaluation2.insert (evalPoint[i]);
441
442  CFList interMedResult;
443  CanonicalForm oldSqrfPartF= sqrfPartF;
444  sqrfPartF= shift2Zero (sqrfPartF, evalSqrfPartF, evaluation2);
445  if (factors.length() > 1)
446  {
447    CanonicalForm LC1= LC (oldSqrfPartF, 1);
448    CFList leadingCoeffs;
449    for (int i= 0; i < factors.length(); i++)
450      leadingCoeffs.append (LC1);
451
452    CFList LC1eval= evaluateAtEval (LC1, evaluation2,2);
453    CFList oldFactors= factors;
454    for (CFListIterator i= oldFactors; i.hasItem(); i++)
455      i.getItem() *= LC1eval.getFirst()/Lc (i.getItem());
456
457    bool success= false;
458    CanonicalForm oldSqrfPartFPowLC= oldSqrfPartF*power(LC1,factors.length()-1);
459    if (size (oldSqrfPartFPowLC)/getNumVars (oldSqrfPartFPowLC) < 500 &&
460        LucksWangSparseHeuristic (oldSqrfPartFPowLC,
461                                  oldFactors, 2, leadingCoeffs, factors))
462    {
463      interMedResult= recoverFactors (oldSqrfPartF, factors);
464      if (oldFactors.length() == interMedResult.length())
465        success= true;
466    }
467    if (!success)
468    {
469      LC1= LC (evalSqrfPartF.getFirst(), 1);
470
471      CFArray leadingCoeffs= CFArray (factors.length());
472      for (int i= 0; i < factors.length(); i++)
473        leadingCoeffs[i]= LC1;
474
475      for (CFListIterator i= factors; i.hasItem(); i++)
476        i.getItem() *= LC1 (0,2)/Lc (i.getItem());
477      factors.insert (1);
478
479      CanonicalForm
480      newSqrfPartF= evalSqrfPartF.getFirst()*power (LC1, factors.length() - 2);
481
482      int liftBound= degree (newSqrfPartF,2) + 1;
483
484      CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
485      CFArray Pi;
486      CFList diophant;
487      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
488                            leadingCoeffs, false);
489
490      if (sqrfPartF.level() > 2)
491      {
492        int* liftBounds= new int [sqrfPartF.level() - 1];
493        liftBounds [0]= liftBound;
494        bool noOneToOne= false;
495        CFList *leadingCoeffs2= new CFList [sqrfPartF.level()-2];
496        LC1= LC (evalSqrfPartF.getLast(), 1);
497        CFList LCs;
498        for (int i= 0; i < factors.length(); i++)
499          LCs.append (LC1);
500        leadingCoeffs2 [sqrfPartF.level() - 3]= LCs;
501        for (int i= sqrfPartF.level() - 1; i > 2; i--)
502        {
503          for (CFListIterator j= LCs; j.hasItem(); j++)
504            j.getItem()= j.getItem() (0, i + 1);
505          leadingCoeffs2 [i - 3]= LCs;
506        }
507        sqrfPartF *= power (LC1, factors.length()-1);
508
509        int liftBoundsLength= sqrfPartF.level() - 1;
510        for (int i= 1; i < liftBoundsLength; i++)
511          liftBounds [i]= degree (sqrfPartF, i + 2) + 1;
512        evalSqrfPartF= evaluateAtZero (sqrfPartF);
513        evalSqrfPartF.removeFirst();
514        factors= nonMonicHenselLift (evalSqrfPartF, factors, leadingCoeffs2,
515                 diophant, Pi, liftBounds, sqrfPartF.level() - 1, noOneToOne);
516        delete [] leadingCoeffs2;
517        delete [] liftBounds;
518      }
519      for (CFListIterator iter= factors; iter.hasItem(); iter++)
520        iter.getItem()= reverseShift (iter.getItem(), evaluation2);
521
522      interMedResult=
523      recoverFactors (reverseShift(evalSqrfPartF.getLast(),evaluation2),
524                      factors);
525    }
526  }
527  else
528  {
529    CanonicalForm contF=content (oldSqrfPartF,1);
530    factors= CFList (oldSqrfPartF/contF);
531    interMedResult= recoverFactors (oldSqrfPartF, factors);
532  }
533
534  CFList result;
535  CFFListIterator k;
536  CanonicalForm tmp;
537  for (int i= 0; i < LCFFactors.length(); i++)
538  {
539    tmp= 1;
540    for (k= bufSqrfFactors[i]; k.hasItem(); k++)
541    {
542      int pos= findItem (bufFactors, k.getItem().factor());
543      if (pos)
544        tmp *= power (getItem (interMedResult, pos), k.getItem().exp());
545    }
546    result.append (tmp);
547  }
548
549  for (CFListIterator i= result; i.hasItem(); i++)
550  {
551    F /= i.getItem();
552    if (foundDifferent)
553      i.getItem()= swapvar (i.getItem(), x, z);
554    i.getItem()= N (i.getItem());
555  }
556
557  if (foundDifferent)
558  {
559    CFList l= differentSecondVarLCs [j];
560    for (CFListIterator i= l; i.hasItem(); i++)
561      i.getItem()= swapvar (i.getItem(), y, z);
562    differentSecondVarLCs [j]= l;
563    F= swapvar (F, x, z);
564  }
565
566  result.insert (N (F));
567
568  result= distributeContent (result, differentSecondVarLCs, length);
569
570  if (!result.getFirst().inCoeffDomain())
571  {
572    CFListIterator i= result;
573    CanonicalForm tmp;
574    if (foundDifferent)
575      i.getItem()= swapvar (i.getItem(), Variable (2), y);
576
577    tmp= i.getItem();
578
579    i++;
580    for (; i.hasItem(); i++)
581    {
582      if (foundDifferent)
583        i.getItem()= swapvar (i.getItem(), Variable (2), y)*tmp;
584      else
585        i.getItem() *= tmp;
586    }
587  }
588  else
589    y= Variable (1);
590
591  delete [] bufSqrfFactors;
592
593  return result;
594}
595
596
597CFList
598multiFactorize (const CanonicalForm& F, const Variable& v)
599{
600  if (F.inCoeffDomain())
601    return CFList (F);
602
603  // compress and find main Variable
604  CFMap N;
605  CanonicalForm A= myCompress (F, N);
606
607  //univariate case
608  if (F.isUnivariate())
609  {
610    CFList result;
611    if (v.level() != 1)
612      result= conv (factorize (F, v));
613    else
614      result= conv (factorize (F,true));
615    if (result.getFirst().inCoeffDomain())
616      result.removeFirst();
617    return result;
618  }
619
620  //bivariate case
621  if (A.level() == 2)
622  {
623    CFList buf= biFactorize (F, v);
624    return buf;
625  }
626
627  Variable x= Variable (1);
628  Variable y= Variable (2);
629
630  // remove content
631  CFList contentAi;
632  CanonicalForm lcmCont= lcmContent (A, contentAi);
633  A /= lcmCont;
634
635  // trivial after content removal
636  CFList contentAFactors;
637  if (A.inCoeffDomain())
638  {
639    for (CFListIterator i= contentAi; i.hasItem(); i++)
640    {
641      if (i.getItem().inCoeffDomain())
642        continue;
643      else
644      {
645        lcmCont /= i.getItem();
646        contentAFactors=
647        Union (multiFactorize (lcmCont, v),
648               multiFactorize (i.getItem(), v));
649        break;
650      }
651    }
652    decompress (contentAFactors, N);
653    if (isOn (SW_RATIONAL))
654      normalize (contentAFactors);
655    return contentAFactors;
656  }
657
658  // factorize content
659  contentAFactors= multiFactorize (lcmCont, v);
660
661  // univariate after content removal
662  CFList factors;
663  if (A.isUnivariate ())
664  {
665    if (v.level() != 1)
666      factors= conv (factorize (A, v));
667    else
668      factors= conv (factorize (A,true));
669    append (factors, contentAFactors);
670    decompress (factors, N);
671    return factors;
672  }
673
674  A *= bCommonDen (A);
675  // check main variable
676  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
677  int factorNums= 3;
678  CanonicalForm bivarEval;
679  CFList biFactors, bufBiFactors;
680  CanonicalForm evalPoly;
681  int lift, bufLift;
682  CFList* bufAeval2= new CFList [A.level() - 2];
683  CFList* Aeval2= new CFList [A.level() - 2];
684  int counter;
685  int differentSecondVar= 0;
686  CanonicalForm bufA;
687  // several bivariate factorizations
688  REvaluation E (2, A.level(), IntRandom (25));
689  for (int i= 0; i < factorNums; i++)
690  {
691    counter= 0;
692    bufA= A;
693    bufAeval= CFList();
694    bufEvaluation= evalPoints (bufA, bufAeval, E);
695    E.nextpoint();
696    evalPoly= 0;
697
698    bivarEval= bufEvaluation.getLast();
699
700    evaluationWRTDifferentSecondVars (bufAeval2, bufEvaluation, A);
701
702    for (int j= 0; j < A.level() - 1; j++)
703    {
704      if (!bufAeval2[j].isEmpty())
705        counter++;
706    }
707
708    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
709
710    TIMING_START (fac_bi_factorize);
711    bufBiFactors= ratBiSqrfFactorize (bufAeval.getFirst(), v);
712    TIMING_END_AND_PRINT (fac_bi_factorize,
713                          "time for bivariate factorization: ");
714    bufBiFactors.removeFirst();
715
716    if (bufBiFactors.length() == 1)
717    {
718      factors.append (A);
719      appendSwapDecompress (factors, contentAFactors, N, 0, 0, x);
720      if (isOn (SW_RATIONAL))
721        normalize (factors);
722      delete [] bufAeval2;
723      delete [] Aeval2;
724      return factors;
725    }
726
727    if (i == 0)
728    {
729      Aeval= bufAeval;
730      evaluation= bufEvaluation;
731      biFactors= bufBiFactors;
732      lift= bufLift;
733      for (int j= 0; j < A.level() - 2; j++)
734        Aeval2 [j]= bufAeval2 [j];
735      differentSecondVar= counter;
736    }
737    else
738    {
739      if (bufBiFactors.length() < biFactors.length() ||
740          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())) ||
741          counter > differentSecondVar)
742      {
743        Aeval= bufAeval;
744        evaluation= bufEvaluation;
745        biFactors= bufBiFactors;
746        lift= bufLift;
747        for (int j= 0; j < A.level() - 2; j++)
748          Aeval2 [j]= bufAeval2 [j];
749        differentSecondVar= counter;
750      }
751    }
752    int k= 0;
753    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
754      evalPoly += j.getItem()*power (x, k);
755    list.append (evalPoly);
756  }
757
758  delete [] bufAeval2;
759
760  sortList (biFactors, x);
761
762  int minFactorsLength;
763  bool irred= false;
764  factorizationWRTDifferentSecondVars (A, Aeval2, minFactorsLength, irred, v);
765
766  if (irred)
767  {
768    factors.append (A);
769    appendSwapDecompress (factors, contentAFactors, N, 0, 0, x);
770    if (isOn (SW_RATIONAL))
771      normalize (factors);
772    delete [] Aeval2;
773    return factors;
774  }
775
776  if (minFactorsLength == 0)
777    minFactorsLength= biFactors.length();
778  else if (biFactors.length() > minFactorsLength)
779    refineBiFactors (A, biFactors, Aeval2, evaluation, minFactorsLength);
780
781  if (differentSecondVar == A.level() - 2 && getNumVars(LC(A,1)) == A.level()-1)
782  {
783    bool zeroOccured= false;
784    for (CFListIterator iter= evaluation; iter.hasItem(); iter++)
785    {
786      if (iter.getItem().isZero())
787      {
788        zeroOccured= true;
789        break;
790      }
791    }
792    if (!zeroOccured)
793    {
794      factors= sparseHeuristic (A, biFactors, Aeval2, evaluation,
795                                minFactorsLength);
796      if (factors.length() == biFactors.length())
797      {
798        appendSwapDecompress (factors, contentAFactors, N, 0, 0, x);
799        normalize (factors);
800        delete [] Aeval2;
801        return factors;
802      }
803      else
804        factors= CFList();
805      //TODO case where factors.length() > 0
806    }
807  }
808
809  CFList uniFactors= buildUniFactors (biFactors, evaluation.getLast(), y);
810
811  sortByUniFactors (Aeval2, A.level() - 2, uniFactors, evaluation);
812
813  CFList * oldAeval= new CFList [A.level() - 2];
814  for (int i= 0; i < A.level() - 2; i++)
815    oldAeval[i]= Aeval2[i];
816
817  getLeadingCoeffs (A, Aeval2, uniFactors, evaluation);
818
819  CFList biFactorsLCs;
820  for (CFListIterator i= biFactors; i.hasItem(); i++)
821    biFactorsLCs.append (LC (i.getItem(), 1));
822
823  Variable w;
824  CFList leadingCoeffs= precomputeLeadingCoeff (LC (A, 1), biFactorsLCs,
825                                          evaluation, Aeval2, A.level() - 2, w);
826
827  if (w.level() != 1)
828  {
829    A= swapvar (A, y, w);
830    int i= A.level();
831    CanonicalForm evalPoint;
832    for (CFListIterator iter= evaluation; iter.hasItem(); iter++, i--)
833    {
834      if (i == w.level())
835      {
836        evalPoint= iter.getItem();
837        iter.getItem()= evaluation.getLast();
838        evaluation.removeLast();
839        evaluation.append (evalPoint);
840        break;
841      }
842    }
843    for (i= 0; i < A.level() - 2; i++)
844    {
845      if (oldAeval[i].isEmpty())
846        continue;
847      if (oldAeval[i].getFirst().level() == w.level())
848      {
849        CFArray tmp= copy (oldAeval[i]);
850        oldAeval[i]= biFactors; 
851        for (CFListIterator iter= oldAeval[i]; iter.hasItem(); iter++)
852          iter.getItem()= swapvar (iter.getItem(), w, y);
853        for (int ii= 0; ii < tmp.size(); ii++)
854          tmp[ii]= swapvar (tmp[ii], w, y);
855        CFArray tmp2= CFArray (tmp.size());
856        CanonicalForm buf;
857        for (int ii= 0; ii < tmp.size(); ii++)
858        {
859          buf= tmp[ii] (evaluation.getLast(),y);
860          buf /= Lc (buf);
861          tmp2[findItem (uniFactors, buf)-1]=tmp[ii];
862        }
863        biFactors= CFList();
864        for (int j= 0; j < tmp2.size(); j++)
865          biFactors.append (tmp2[j]);
866      }
867    }
868  }
869
870  CFListIterator iter;
871  CanonicalForm oldA= A;
872  CFList oldBiFactors= biFactors;
873  if (!leadingCoeffs.getFirst().inCoeffDomain())
874  {
875    CanonicalForm tmp= power (leadingCoeffs.getFirst(), biFactors.length() - 1);
876    A *= tmp;
877    tmp= leadingCoeffs.getFirst();
878    iter= evaluation;
879    for (int i= A.level(); i > 2; i--, iter++)
880      tmp= tmp (iter.getItem(), i);
881    if (!tmp.inCoeffDomain())
882    {
883      for (CFListIterator i= biFactors; i.hasItem(); i++)
884      {
885        i.getItem() *= tmp/LC (i.getItem(), 1);
886        i.getItem() /= Lc (i.getItem());
887      }
888    }
889  }
890
891  CanonicalForm LCmultiplier= leadingCoeffs.getFirst();
892  leadingCoeffs.removeFirst();
893
894  //prepare leading coefficients
895  CFList* leadingCoeffs2= new CFList [A.level() - 2];
896  prepareLeadingCoeffs (leadingCoeffs2, A.level(), leadingCoeffs, biFactors,
897                        evaluation);
898
899
900  Aeval= evaluateAtEval (A, evaluation, 2);
901
902  CanonicalForm hh= Lc (Aeval.getFirst());
903
904  for (iter= Aeval; iter.hasItem(); iter++)
905    iter.getItem() /= hh;
906
907  A /= hh;
908
909  CFList bufFactors= CFList();
910  if (LucksWangSparseHeuristic (A, biFactors, 2, leadingCoeffs2 [A.level() - 3],
911                                factors))
912  {
913    int check= biFactors.length();
914    int * index= new int [factors.length()];
915    CFList oldFactors= factors;
916    factors= recoverFactors (A, factors, index);
917
918    if (check == factors.length())
919    {
920      if (w.level() != 1)
921      {
922        for (iter= factors; iter.hasItem(); iter++)
923          iter.getItem()= swapvar (iter.getItem(), w, y);
924      }
925
926      appendSwapDecompress (factors, contentAFactors, N, 0, 0, x);
927      normalize (factors);
928      delete [] index;
929      delete [] Aeval2;
930      return factors;
931    }
932    else if (factors.length() > 0)
933    {
934      int oneCount= 0;
935      CFList l;
936      for (int i= 0; i < check; i++)
937      {
938        if (index[i] == 1)
939        {
940          iter=biFactors;
941          for (int j=1; j <= i-oneCount; j++)
942            iter++;
943          iter.remove (1);
944          for (int j= 0; j < A.level() -2; j++)
945          {
946            l= leadingCoeffs2[j];
947            iter= l;
948            for (int k=1; k <= i-oneCount; k++)
949              iter++;
950            iter.remove (1);
951            leadingCoeffs2[j]=l;
952          }
953          oneCount++;
954        }
955      }
956      bufFactors= factors;
957      factors= CFList();
958    }
959    else if (!LCmultiplier.inCoeffDomain() && factors.length() == 0)
960    {
961      factors= oldFactors;
962      CanonicalForm cont;
963      CFList contents, LCs;
964      CFListIterator iter2;
965      int index=1;
966      bool foundTrueMultiplier= false;
967      for (iter= factors; iter.hasItem(); iter++, index++)
968      {
969        cont= content (iter.getItem(), 1);
970        cont= gcd (cont , LCmultiplier);
971        contents.append (cont);
972        if (cont.inCoeffDomain()) // trivial content->LCmultiplier needs to go there
973        {
974          foundTrueMultiplier= true;
975          int index2= 1;
976          for (iter2= leadingCoeffs2[A.level()-3]; iter2.hasItem(); iter2++,
977                                                                    index2++)
978          {
979            if (index2 == index)
980              continue;
981            iter2.getItem() /= LCmultiplier;
982          }
983          A /= power (LCmultiplier, biFactors.length() -1);
984          leadingCoeffs= leadingCoeffs2[A.level()-3];
985          for (int i= A.level()-3; i > -1; i--)
986            leadingCoeffs2[i]= CFList();
987          prepareLeadingCoeffs (leadingCoeffs2, A.level(), leadingCoeffs,
988                                biFactors, evaluation );
989          Aeval= evaluateAtEval (A, evaluation, 2);
990
991          hh= Lc (Aeval.getFirst());
992
993          for (iter2= Aeval; iter2.hasItem(); iter2++)
994            iter2.getItem() /= hh;
995
996          A /= hh;
997          break;
998        }
999        else
1000          LCs.append (LC (iter.getItem()/cont, 1));
1001      }
1002      if (!foundTrueMultiplier)
1003      {
1004        index= 1;
1005        iter2= factors;
1006        bool foundMultiplier= false;
1007        for (iter= contents; iter.hasItem(); iter++, iter2++, index++)
1008        {
1009          if (fdivides (iter.getItem(), LCmultiplier))
1010          {
1011            if ((LCmultiplier/iter.getItem()).inCoeffDomain() &&
1012                !isOnlyLeadingCoeff(iter2.getItem())) //content divides LCmultiplier completely and factor consists of more terms than just the leading coeff
1013            {
1014              int index2= 1;
1015              for (CFListIterator iter3= leadingCoeffs2[A.level()-3];
1016                   iter3.hasItem(); iter3++, index2++)
1017              {
1018                if (index2 == index)
1019                {
1020                  iter3.getItem() /= LCmultiplier;
1021                  break;
1022                }
1023              }
1024              A /= LCmultiplier;
1025              foundMultiplier= true;
1026              iter.getItem()= 1;
1027              break;
1028            }
1029          }
1030        }
1031        // coming from above: divide out more LCmultiplier if possible
1032        if (foundMultiplier)
1033        {
1034          foundMultiplier= false;
1035          index=1;
1036          iter2= factors;
1037          for (iter= contents; iter.hasItem(); iter++, iter2++, index++)
1038          {
1039            if (!(iter.getItem().isOne()) &&
1040                fdivides (iter.getItem(), LCmultiplier))
1041            {
1042              if (!isOnlyLeadingCoeff (iter2.getItem())) // factor is more than just leading coeff
1043              {
1044                int index2= 1;
1045                for (iter2= leadingCoeffs2[A.level()-3]; iter2.hasItem();
1046                     iter2++, index2++)
1047                {
1048                  if (index2 == index)
1049                  {
1050                    iter2.getItem() /= iter.getItem();
1051                    foundMultiplier= true;
1052                    break;
1053                  }
1054                }
1055                A /= iter.getItem();
1056                LCmultiplier /= iter.getItem();
1057                iter.getItem()= 1;
1058                break;
1059              }
1060              else //factor is consists of just leading coeff
1061              {
1062                CanonicalForm vars=getVars (iter.getItem());
1063                CanonicalForm factor;
1064                Variable xx;
1065                bool oneVariableNotInCommon= false;
1066                for (int i= 0; i < A.level()-2; i++)
1067                {
1068                  if (oldAeval[i].isEmpty())
1069                    continue;
1070                  xx= oldAeval[i].getFirst().mvar();
1071                  factor= LC (getItem (oldAeval[i], index),1);
1072                  if ((factor.inCoeffDomain() && degree (vars,xx) > 0) ||
1073                      (degree (factor,xx) > 0 && degree (vars,xx) < 0)) //scan for bivariate factors with leading coeff that does not contain variables which occur in LCmultiplier
1074                  {
1075                    oneVariableNotInCommon= true;
1076                    break;
1077                  }
1078                }
1079                if (oneVariableNotInCommon)
1080                {
1081                  int index2= 1;
1082                  for (iter2= leadingCoeffs2[A.level()-3]; iter2.hasItem();
1083                       iter2++, index2++)
1084                  {
1085                    if (index2 == index)
1086                    {
1087                      iter2.getItem() /= iter.getItem();
1088                      foundMultiplier= true;
1089                      break;
1090                    }
1091                  }
1092                  A /= iter.getItem();
1093                  LCmultiplier /= iter.getItem();
1094                  iter.getItem()= 1;
1095                  break;
1096                }
1097              }
1098            }
1099          }
1100          // wipe out the last LCmultiplier
1101          if (foundMultiplier)
1102          {
1103            index= 1;
1104            for (iter= contents; iter.hasItem(); iter++, index++)
1105            {
1106              if (!iter.getItem().isOne() &&
1107                  fdivides (LCmultiplier, iter.getItem()))
1108              {
1109                int index2= 1;
1110                for (iter2= leadingCoeffs2[A.level()-3]; iter2.hasItem();
1111                     iter2++, index2++)
1112                {
1113                  if (index2 == index)
1114                  {
1115                    iter2.getItem() /= LCmultiplier;
1116                    A /= LCmultiplier;
1117                    iter.getItem() /= LCmultiplier;
1118                  }
1119                }
1120              }
1121            }
1122          }
1123        }
1124        else
1125        {
1126          CanonicalForm pLCs= prod (LCs);
1127          if (fdivides (pLCs, LC (oldA,1)) && (LC(oldA,1)/pLCs).inCoeffDomain()) // check if the product of the lead coeffs of the primitive factors equals the lead coeff of the old A
1128          {
1129            A= oldA;
1130            iter2= leadingCoeffs2[A.level()-3];
1131            for (iter= contents; iter.hasItem(); iter++, iter2++)
1132              iter2.getItem() /= iter.getItem();
1133          }
1134        }
1135
1136        // patch everything together again
1137        leadingCoeffs= leadingCoeffs2[A.level()-3];
1138        for (int i= A.level()-3; i > -1; i--)
1139          leadingCoeffs2[i]= CFList();
1140        prepareLeadingCoeffs (leadingCoeffs2,A.level(),leadingCoeffs, biFactors,
1141                              evaluation);
1142        Aeval= evaluateAtEval (A, evaluation, 2);
1143
1144        hh= Lc (Aeval.getFirst());
1145
1146        for (CFListIterator i= Aeval; i.hasItem(); i++)
1147          i.getItem() /= hh;
1148
1149        A /= hh;
1150      }
1151      factors= CFList();
1152    }
1153    else
1154      factors= CFList();
1155    delete [] index;
1156  }
1157
1158  //shifting to zero
1159  A= shift2Zero (A, Aeval, evaluation);
1160
1161  for (iter= biFactors; iter.hasItem(); iter++)
1162    iter.getItem()= iter.getItem () (y + evaluation.getLast(), y);
1163
1164  for (int i= 0; i < A.level() - 3; i++)
1165    leadingCoeffs2[i]= CFList();
1166  for (iter= leadingCoeffs2[A.level() - 3]; iter.hasItem(); iter++)
1167  {
1168    iter.getItem()= shift2Zero (iter.getItem(), list, evaluation);
1169    for (int i= A.level() - 4; i > -1; i--)
1170    {
1171      if (i + 1 == A.level() - 3)
1172        leadingCoeffs2[i].append (iter.getItem() (0, i + 4));
1173      else
1174        leadingCoeffs2[i].append (leadingCoeffs2[i+1].getLast() (0, i + 4));
1175    }
1176  }
1177
1178  CFArray Pi;
1179  CFList diophant;
1180  int* liftBounds= new int [A.level() - 1];
1181  int liftBoundsLength= A.level() - 1;
1182  for (int i= 0; i < liftBoundsLength; i++)
1183    liftBounds [i]= degree (A, i + 2) + 1;
1184
1185  Aeval.removeFirst();
1186  bool noOneToOne= false;
1187
1188  CFList commonDenominators;
1189  for (CFListIterator iter=biFactors; iter.hasItem(); iter++)
1190    commonDenominators.append (bCommonDen (iter.getItem()));
1191  CanonicalForm tmp1, tmp2, tmp3=1;
1192  CFListIterator iter1, iter2;
1193  for (int i= 0; i < A.level() - 2; i++)
1194  {
1195    iter2= commonDenominators;
1196    for (iter1= leadingCoeffs2[i]; iter1.hasItem(); iter1++, iter2++)
1197    {
1198      tmp1= bCommonDen (iter1.getItem());
1199      Off (SW_RATIONAL);
1200      iter2.getItem()= lcm (iter2.getItem(), tmp1);
1201      On (SW_RATIONAL);
1202    }
1203  }
1204  tmp1= prod (commonDenominators);
1205  for (iter1= Aeval; iter1.hasItem(); iter1++)
1206  {
1207    tmp2= bCommonDen (iter1.getItem());
1208    Off (SW_RATIONAL);
1209    tmp3= lcm (tmp2,tmp3);
1210    On (SW_RATIONAL);
1211  }
1212  CanonicalForm multiplier;
1213  multiplier= tmp3/tmp1;
1214  iter2= commonDenominators;
1215  for (iter1=biFactors; iter1.hasItem(); iter1++, iter2++)
1216    iter1.getItem() *= iter2.getItem()*multiplier;
1217
1218  for (iter1= Aeval; iter1.hasItem(); iter1++)
1219    iter1.getItem() *= tmp3*power (multiplier, biFactors.length() - 1);
1220
1221  for (int i= 0; i < A.level() - 2; i++)
1222  {
1223    iter2= commonDenominators;
1224    for (iter1= leadingCoeffs2[i]; iter1.hasItem(); iter1++, iter2++)
1225      iter1.getItem() *= iter2.getItem()*multiplier;
1226  }
1227
1228
1229  factors= nonMonicHenselLift (Aeval, biFactors, leadingCoeffs2, diophant,
1230                               Pi, liftBounds, liftBoundsLength, noOneToOne);
1231
1232  if (!noOneToOne)
1233  {
1234    int check= factors.length();
1235    A= oldA;
1236    factors= recoverFactors (A, factors, evaluation);
1237    if (check != factors.length())
1238      noOneToOne= true;
1239    else
1240      factors= Union (factors, bufFactors);
1241  }
1242  if (noOneToOne)
1243  {
1244    A= shift2Zero (oldA, Aeval, evaluation);
1245    biFactors= oldBiFactors;
1246    for (iter= biFactors; iter.hasItem(); iter++)
1247      iter.getItem()= iter.getItem () (y + evaluation.getLast(), y);
1248    CanonicalForm LCA= LC (Aeval.getFirst(), 1);
1249    CanonicalForm yToLift= power (y, lift);
1250    CFListIterator i= biFactors;
1251    lift= degree (i.getItem(), 2) + degree (LC (i.getItem(), 1)) + 1;
1252    i++;
1253
1254    for (; i.hasItem(); i++)
1255      lift= tmax (lift, degree (i.getItem(), 2)+degree (LC (i.getItem(), 1))+1);
1256
1257    lift= tmax (degree (Aeval.getFirst() , 2) + 1, lift);
1258
1259    i= biFactors;
1260    yToLift= power (y, lift);
1261    CanonicalForm dummy;
1262    for (; i.hasItem(); i++)
1263    {
1264      LCA= LC (i.getItem(), 1);
1265      extgcd (LCA, yToLift, LCA, dummy);
1266      i.getItem()= mod (i.getItem()*LCA, yToLift);
1267    }
1268
1269    liftBoundsLength= F.level() - 1;
1270    liftBounds= liftingBounds (A, lift);
1271
1272    CFList MOD;
1273    bool earlySuccess;
1274    CFList earlyFactors;
1275    ExtensionInfo info= ExtensionInfo (false);
1276    CFList liftedFactors;
1277    TIMING_START (fac_hensel_lift);
1278    liftedFactors= henselLiftAndEarly
1279                   (A, MOD, liftBounds, earlySuccess, earlyFactors,
1280                    Aeval, biFactors, evaluation, info);
1281    TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1282
1283    TIMING_START (fac_factor_recombination);
1284    factors= factorRecombination (A, liftedFactors, MOD);
1285    TIMING_END_AND_PRINT (fac_factor_recombination,
1286                          "time for factor recombination: ");
1287
1288    if (earlySuccess)
1289      factors= Union (factors, earlyFactors);
1290
1291    for (CFListIterator i= factors; i.hasItem(); i++)
1292    {
1293      int kk= Aeval.getLast().level();
1294      for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
1295      {
1296        if (i.getItem().level() < kk)
1297          continue;
1298       i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
1299      }
1300    }
1301  }
1302
1303  if (w.level() != 1)
1304  {
1305    for (CFListIterator iter= factors; iter.hasItem(); iter++)
1306      iter.getItem()= swapvar (iter.getItem(), w, y);
1307  }
1308
1309  swap (factors, 0, 0, x);
1310  append (factors, contentAFactors);
1311  decompress (factors, N);
1312  if (isOn (SW_RATIONAL))
1313    normalize (factors);
1314
1315  delete [] leadingCoeffs2;
1316  delete [] oldAeval;
1317  delete [] Aeval2;
1318  delete[] liftBounds;
1319
1320  return factors;
1321}
1322
1323#endif
Note: See TracBrowser for help on using the repository browser.