source: git/factory/facFqFactorize.cc @ 1372ae

spielwiese
Last change on this file since 1372ae was 1372ae, checked in by Hans Schoenemann <hannes@…>, 13 years ago
Variable(i) ->v.level() git-svn-id: file:///usr/local/Singular/svn/trunk@13666 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 54.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqFactorize.cc
5 *
6 * This file implements functions for factoring a multivariate polynomial over
7 * a finite field.
8 *
9 * ABSTRACT: "Efficient Multivariate Factorization over Finite Fields" by
10 * L. Bernardin & M. Monagon.
11 *
12 * @author Martin Lee
13 *
14 * @internal @version \$Id$
15 *
16 **/
17/*****************************************************************************/
18
19#include <config.h>
20
21#include "assert.h"
22#include "debug.h"
23#include "timing.h"
24
25#include "facFqFactorizeUtil.h"
26#include "facFqFactorize.h"
27#include "cf_random.h"
28#include "facHensel.h"
29#include "cf_gcd_smallp.h"
30#include "cf_map_ext.h"
31#include "algext.h"
32
33#ifdef HAVE_NTL
34#include <NTL/ZZ_pEX.h>
35#include "NTLconvert.h"
36
37TIMING_DEFINE_PRINT(fac_bi_factorizer);
38TIMING_DEFINE_PRINT(fac_hensel_lift);
39TIMING_DEFINE_PRINT(fac_factor_recombination);
40
41static inline
42CanonicalForm
43listGCD (const CFList& L);
44
45static inline
46CanonicalForm
47myContent (const CanonicalForm& F)
48{
49  CanonicalForm G= swapvar (F, F.mvar(), Variable (1));
50  CFList L;
51  Variable alpha;
52  bool algExt= hasFirstAlgVar (G, alpha);
53  for (CFIterator i= G; i.hasTerms(); i++)
54    L.append (i.coeff());
55  if (L.length() == 2)
56  {
57    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
58    if (GF)
59    {
60      return swapvar (GCD_GF (L.getFirst(), L.getLast()), F.mvar(),
61                      Variable (1));
62    }
63    else if (!GF && algExt)
64    {
65      return swapvar (GCD_Fp_extension (L.getFirst(), L.getLast(), alpha),
66                                       F.mvar(), Variable (1));
67    }
68    else
69      return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(),
70                     Variable (1));
71  }
72  if (L.length() == 1)
73  {
74    return LC (F, Variable (1));
75  }
76  return swapvar (listGCD (L), F.mvar(), Variable (1));
77}
78
79static inline
80CanonicalForm
81listGCD (const CFList& L)
82{
83  if (L.length() == 1)
84    return L.getFirst();
85  Variable alpha;
86  if (L.length() == 2)
87  {
88    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
89    if (GF)
90    {
91      return GCD_GF (L.getFirst(), L.getLast());
92    }
93    else if (!GF && (hasFirstAlgVar (L.getFirst(), alpha) ||
94                     hasFirstAlgVar (L.getLast(), alpha)))
95    {
96      return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
97    }
98    else
99      return GCD_small_p (L.getFirst(), L.getLast());
100  }
101  else
102  {
103    CFList lHi, lLo;
104    CanonicalForm resultHi, resultLo;
105    int length= L.length()/2;
106    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
107    int j= 0;
108    for (CFListIterator i= L; j < length; i++, j++)
109      lHi.append (i.getItem());
110    lLo= Difference (L, lHi);
111    resultHi= listGCD (lHi);
112    resultLo= listGCD (lLo);
113    if (resultHi.isOne() || resultLo.isOne())
114      return 1;
115    if (GF)
116    {
117      return GCD_GF (resultHi, resultLo);
118    }
119    else if (!GF && (hasFirstAlgVar (resultHi, alpha) ||
120                     hasFirstAlgVar (resultLo, alpha)))
121    {
122      return GCD_Fp_extension (resultHi, resultLo, alpha);
123    }
124    else
125      return GCD_small_p (resultHi, resultLo);
126  }
127}
128
129static inline
130CanonicalForm
131myContent (const CanonicalForm& F, const Variable& x)
132{
133  if (degree (F, x) <= 0)
134    return 1;
135  CanonicalForm G= F;
136  bool swap= false;
137  if (x != F.mvar())
138  {
139    swap= true;
140    G= swapvar (F, x, F.mvar());
141  }
142  CFList L;
143  Variable alpha;
144  bool algExt= hasFirstAlgVar (G, alpha);
145  for (CFIterator i= G; i.hasTerms(); i++)
146    L.append (i.coeff());
147  if (L.length() == 2)
148  {
149    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
150    if (GF)
151    {
152      if (swap)
153        return swapvar(GCD_GF (L.getFirst(), L.getLast()), F.mvar(), x);
154      else
155        return GCD_GF (L.getFirst(), L.getLast());
156    }
157    else if (!GF && algExt)
158    {
159      if (swap)
160        return swapvar(GCD_Fp_extension (L.getFirst(), L.getLast(), alpha),
161                                         F.mvar(), x);
162      else
163        return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
164    }
165    else
166    {
167      if (swap)
168        return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(),
169                        x);
170      else
171        return GCD_small_p (L.getFirst(), L.getLast());
172    }
173  }
174  if (L.length() == 1)
175  {
176    return LC (F, x);
177  }
178  if (swap)
179    return swapvar (listGCD (L), F.mvar(), x);
180  else
181    return listGCD (L);
182}
183
184static inline
185CanonicalForm
186myLcm (const CanonicalForm& F, const CanonicalForm& G)
187{
188  if ( F.isZero() || G.isZero() )
189    return 0;
190  else
191  {
192    Variable alpha;
193    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
194    if (GF)
195      return (F/GCD_GF (F, G))*G;
196    else if (!GF && (hasFirstAlgVar (F, alpha) ||
197                     hasFirstAlgVar (G, alpha)))
198      return (F/GCD_Fp_extension (F, G, alpha))*G;
199    else
200      return (F/GCD_small_p (F, G))*G;
201  }
202}
203
204
205static inline
206CanonicalForm myCompress (const CanonicalForm& F, CFMap& N)
207{
208  int n= F.level();
209  int * degsf= new int [n + 1];
210  int ** swap;
211  swap= new int* [n + 1];
212  for (int i= 0; i <= n; i++)
213  {
214    degsf[i]= 0;
215    swap [i]= new int [2];
216    swap [i] [0]= 0;
217    swap [i] [1]= 0;
218  }
219  int i= 1;
220  n= 1;
221  degsf= degrees (F, degsf);
222
223  CanonicalForm result= F;
224  while ( i <= F.level() )
225  {
226    while( degsf[i] == 0 ) i++;
227    swap[n][0]= i;
228    swap[n][1]= degsf[i];
229    if (i != n)
230      result= swapvar (result, Variable (n), Variable(i));
231    n++; i++;
232  }
233
234  int buf1, buf2;
235  n--;
236
237  for (i= 1; i < n; i++)
238  {
239    for (int j= 1; j < n - i + 1; j++)
240    {
241      if (swap[j][1] < swap[j + 1][1])
242      {
243        buf1= swap [j + 1] [0];
244        buf2= swap [j + 1] [1];
245        swap[j + 1] [0]= swap[j] [0];
246        swap[j + 1] [1]= swap[j] [1];
247        swap[j][0]= buf1;
248        swap[j][1]= buf2;
249        result= swapvar (result, Variable (j + 1), Variable (j));
250      }
251    }
252  }
253
254  for (i= n; i > 0; i--)
255  {
256    if (i != swap[i] [0])
257      N.newpair (Variable (i), Variable (swap[i] [0]));
258  }
259
260  for (i= 0; i <= n; i++)
261    delete [] swap[i];
262  delete [] swap;
263
264  delete [] degsf;
265
266  return result;
267}
268
269CFList
270monicFactorRecombi (const CFList& factors,const CanonicalForm& F, const
271                    CanonicalForm& M, const DegreePattern& degs)
272{
273  if (degs.getLength() == 1)
274    return CFList(F);
275
276  CFList T, S;
277
278  T= factors;
279
280  int s= 1;
281  CFList result;
282  CanonicalForm buf= F;
283  CanonicalForm LCBuf= LC (buf, Variable (1));
284  CanonicalForm g, gg;
285  int v [T.length()];
286  for (int i= 0; i < T.length(); i++)
287    v[i]= 0;
288  bool noSubset= false;
289  CFArray TT;
290  int subsetDeg;
291  DegreePattern bufDegs1= degs, bufDegs2;
292  TT= copy (factors);
293  while (T.length() >= 2*s)
294  {
295    while (noSubset == false)
296    {
297      if (T.length() == s)
298      {
299        result.append (prodMod (T, M));
300        return result;
301      }
302      S= subset (v, s, TT, noSubset);
303      if (noSubset) break;
304      subsetDeg= subsetDegree (S);
305      if (!degs.find (subsetDeg))
306        continue;
307      else
308      {
309        g= prodMod0 (S, M);
310        gg= mod (g*LCBuf, M); //univariate polys
311        gg /= content (gg);
312        if (fdivides (LC (gg), LCBuf))
313        {
314          g= prodMod (S, M);
315          gg= mulMod2 (LCBuf, g, M);
316          gg /= content (gg, Variable (1));
317          if (fdivides (gg, buf))
318          {
319            result.append (g);
320            buf /= gg;
321            LCBuf= LC (buf, Variable(1));
322            T= Difference (T, S);
323            // compute new possible degree pattern
324            bufDegs2= DegreePattern (T);
325            bufDegs1.intersect (bufDegs2);
326            bufDegs1.refine ();
327            if (T.length() < 2*s || T.length() == s ||
328                bufDegs1.getLength() == 1)
329            {
330              result.append (prodMod (T, M));
331              return result;
332            }
333            TT= copy (T);
334            indexUpdate (v, s, T.length(), noSubset);
335            if (noSubset) break;
336          }
337        }
338      }
339    }
340    s++;
341    if (T.length() < 2*s || T.length() == s)
342    {
343      result.append (prodMod (T, M));
344      return result;
345    }
346    int v [T.length()];
347    for (int i= 0; i < T.length(); i++)
348      v[i]= 0;
349    noSubset= false;
350  }
351  if (T.length() < 2*s)
352    result.append (prodMod (T, M));
353
354  return result;
355}
356
357CFList
358earlyMonicFactorDetect (CanonicalForm& F, CFList& factors,
359                        int& adaptedLiftBound, DegreePattern& degs,
360                        bool& success, int deg, const int bound)
361{
362  DegreePattern bufDegs1= degs;
363  DegreePattern bufDegs2;
364  CFList result;
365  CFList T= factors;
366  CanonicalForm buf= F;
367  CanonicalForm LCBuf= LC (buf, Variable (1));
368  CanonicalForm g, gg;
369  CanonicalForm M= power (F.mvar(), deg);
370  int d= bound;
371  int e= 0;
372  adaptedLiftBound= 0;
373  for (CFListIterator i= factors; i.hasItem(); i++)
374  {
375    if (!bufDegs1.find (degree (i.getItem(), 1)))
376      continue;
377    else
378    {
379      gg= i.getItem() (0, 1);
380      gg *= LCBuf;
381      gg= mod (gg, M);
382      if (fdivides (LC (gg), LCBuf))
383      {
384        g= i.getItem();
385        gg= mulMod2 (g, LCBuf, M);
386        gg /= content (gg, Variable (1));
387        if (fdivides (gg, buf))
388        {
389          result.append (g);
390          CanonicalForm bufbuf= buf;
391          buf /= gg;
392          d -= degree (gg) + degree (LC (gg, 1));
393          e= tmax (e, degree (gg) + degree (LC (gg, 1)));
394          LCBuf= LC (buf, Variable (1));
395          T= Difference (T, CFList (i.getItem()));
396
397          // compute new possible degree pattern
398          bufDegs2= DegreePattern (T);
399          bufDegs1.intersect (bufDegs2);
400          bufDegs1.refine ();
401          if (bufDegs1.getLength() <= 1)
402          {
403            result.append (prodMod (T, M));
404            break;
405          }
406        }
407      }
408    }
409  }
410  adaptedLiftBound= d;
411
412  if (adaptedLiftBound < deg)
413  {
414    factors= T;
415    degs= bufDegs1;
416    if (adaptedLiftBound < degree (F) + 1)
417    {
418      if (d == 1)
419      {
420        if (e + 1 > deg)
421        {
422          success= false;
423          adaptedLiftBound= deg;
424        }
425        else
426        {
427          F= buf;
428          success= true;
429          if (e + 1 < degree (F) + 1)
430            adaptedLiftBound= deg;
431          else
432            adaptedLiftBound= e + 1;
433        }
434      }
435      else
436      {
437        success= true;
438        F= buf;
439        adaptedLiftBound= deg;
440      }
441    }
442    else
443    {
444      F= buf;
445      success= true;
446    }
447  }
448  if (bufDegs1.getLength() <= 1)
449  {
450    degs= bufDegs1;
451    if (!success)
452      adaptedLiftBound= deg;
453  }
454
455  return result;
456}
457
458CFList biFactorizer (const CanonicalForm& F, const Variable& alpha,
459                     CanonicalForm& bivarEval, int& liftBound)
460{
461  CanonicalForm A= F;
462  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
463
464  if (A.isUnivariate())
465    return uniFactorizer (F, alpha, GF);
466
467
468  CFMap N;
469  A= compress (A, N);
470  Variable y= A.mvar();
471  A /= Lc(A);
472
473  if (y.level() > 2) return CFList (F);
474  Variable x= Variable (1);
475
476  //remove content and factorize content
477  CanonicalForm contentAy= content (A, y);
478  CanonicalForm contentAx= content (A, x);
479  A= A/(contentAy*contentAx);
480  CFList contentAyFactors= uniFactorizer (contentAy, alpha, GF);
481
482  //trivial case
483  CFList factors;
484  if (A.inCoeffDomain())
485  {
486    append (factors, contentAyFactors);
487    decompress (factors, N);
488    bivarEval= N (y - bivarEval);
489    return factors;
490  }
491  else if (A.isUnivariate())
492  {
493    if (A.mvar() == x)
494      factors= uniFactorizer (A, alpha, GF);
495    append (factors, contentAyFactors);
496    decompress (factors, N);
497    bivarEval= N (y - bivarEval);
498    return factors;
499  }
500  bool fail;
501  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
502  CFList uniFactors, list, bufUniFactors;
503  DegreePattern degs;
504  DegreePattern bufDegs;
505
506  int factorNums= 5;
507  double logarithm= (double) ilog2 (totaldegree (A));
508  logarithm /= log2exp;
509  logarithm= ceil (logarithm);
510  if (factorNums < (int) logarithm)
511    factorNums= (int) logarithm;
512  for (int i= 0; i < factorNums; i++)
513  {
514    if (i == 0)
515      Aeval= A (bivarEval, y);
516    else if (i > 0 && contentAx.inCoeffDomain())
517    {
518      Aeval= A;
519      evaluation= evalPoint (A, Aeval, alpha, list, GF, fail);
520    }
521
522    if (fail && (i != 0))
523      break;
524
525    // univariate factorization
526    uniFactors= uniFactorizer (Aeval, alpha, GF);
527
528    if (uniFactors.length() == 1)
529    {
530      append (factors, contentAyFactors);
531      if (contentAyFactors.isEmpty())
532        factors.append (F/lc(F));
533      else
534      {
535        A= A (y - bivarEval, y);
536        A= A/lc (A);
537        if (!LC(A, 1).isOne())
538        {
539          CanonicalForm s,t;
540          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
541          A *= s;
542          A= mod (A, power (y, liftBound));
543        }
544        factors.append (A);
545      }
546      decompress (factors, N);
547      bivarEval= N (bivarEval);
548      return factors;
549    }
550
551    // degree analysis
552    degs= DegreePattern (uniFactors);
553
554    if (i == 0)
555    {
556      bufAeval= Aeval;
557      bufEvaluation= bivarEval;
558      bufUniFactors= uniFactors;
559      bufDegs= degs;
560      if (!contentAx.inCoeffDomain())
561        break;
562    }
563    else
564    {
565      bufDegs.intersect (degs);
566      if (uniFactors.length() < bufUniFactors.length())
567      {
568        bufUniFactors= uniFactors;
569        bufAeval= Aeval;
570        bufEvaluation= evaluation;
571      }
572    }
573
574    if (bufDegs.getLength() == 1)
575    {
576      append (factors, contentAyFactors);
577      if (contentAyFactors.isEmpty())
578        factors.append (F/lc(F));
579      else
580      {
581        A= A (y - bivarEval, y);
582        A= A/lc (A);
583        if (!LC(A, 1).isOne())
584        {
585          CanonicalForm s,t;
586          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
587          A *= s;
588          A= mod (A, power (y, liftBound));
589        }
590        factors.append (A);
591      }
592      decompress (factors, N);
593      bivarEval= N (bivarEval);
594      return factors;
595    }
596    list.append (evaluation);
597  }
598
599  bivarEval= y - bufEvaluation;
600  A= A (y + bufEvaluation, y);
601  bufUniFactors.insert (LC (A, x));
602
603  // Hensel lifting
604  CFList diophant;
605  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
606  CFArray Pi;
607  bool earlySuccess= false;
608  int newLiftBound= 0;
609  CFList earlyFactors;
610  int smallFactorDeg= 11; //tunable parameter
611  if (smallFactorDeg >= liftBound)
612    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
613  else if (smallFactorDeg >= degree (A, y) + 1)
614  {
615    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
616    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
617                                           bufDegs, earlySuccess,
618                                           degree (A, y) + 1, liftBound);
619    if (bufDegs.getLength() > 1 && !earlySuccess)
620    {
621      if (newLiftBound > degree (A, y) + 1)
622      {
623        liftBound= newLiftBound;
624        bufUniFactors.insert (LC(A, x));
625        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
626                            Pi, diophant, M);
627      }
628    }
629    else if (earlySuccess)
630      liftBound= newLiftBound;
631  }
632  else if (smallFactorDeg < degree (A, y) + 1)
633  {
634    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
635    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
636                                           bufDegs, earlySuccess,
637                                           smallFactorDeg, liftBound);
638    if (bufDegs.getLength() > 1 && !earlySuccess)
639    {
640      bufUniFactors.insert (LC (A, x));
641      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) +
642                          1, Pi, diophant, M);
643      earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
644                                             bufDegs, earlySuccess,
645                                             degree (A, y) + 1, liftBound);
646      if (bufDegs.getLength() > 1 && !earlySuccess)
647      {
648        if (newLiftBound > degree (A, y) + 1)
649        {
650          if (newLiftBound < newLiftBound)
651            liftBound= newLiftBound;
652          bufUniFactors.insert (LC(A, x));
653          henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
654                              Pi, diophant, M);
655        }
656      }
657      else if (earlySuccess)
658        liftBound= newLiftBound;
659    }
660    else if (earlySuccess)
661      liftBound= newLiftBound;
662  }
663
664  if (newLiftBound > 0)
665    liftBound= newLiftBound;
666
667  CanonicalForm MODl= power (y, liftBound);
668
669  if (bufDegs.getLength() > 1)
670    factors= monicFactorRecombi (bufUniFactors, A, MODl, bufDegs);
671
672  if (earlySuccess)
673    factors= Union (earlyFactors, factors);
674  else if (!earlySuccess && bufDegs.getLength() == 1)
675    factors= earlyFactors;
676
677  decompressAppend (factors, contentAyFactors, N);
678
679  bivarEval= N (bivarEval);
680
681  return factors;
682}
683
684CFList
685extFactorRecombination (const CFList& factors, const CanonicalForm& F,
686                        const CFList& M, const ExtensionInfo& info,
687                        const CFList& evaluation)
688{
689  Variable alpha= info.getAlpha();
690  Variable beta= info.getBeta();
691  CanonicalForm gamma= info.getGamma();
692  CanonicalForm delta= info.getDelta();
693  int k= info.getGFDegree();
694  CFList source, dest;
695  if (factors.length() == 1)
696  {
697    CanonicalForm buf= reverseShift (F, evaluation);
698    return CFList (mapDown (buf, info, source, dest));
699  }
700  if (factors.length() < 1)
701    return CFList();
702
703  int degMipoBeta;
704  if (!k && beta.level() == 1)
705    degMipoBeta= 1;
706  else if (!k && beta.level() != 1)
707    degMipoBeta= degree (getMipo (beta));
708
709  CFList T, S;
710  T= factors;
711
712  int s= 1;
713  CFList result;
714  CanonicalForm buf;
715  if (beta != Variable (1))
716    buf= mapDown (F, gamma, delta, alpha, source, dest);
717  else
718    buf= F;
719
720  CanonicalForm g, LCBuf= LC (buf, Variable (1));
721  CanonicalForm buf2;
722  int v [T.length()];
723  for (int i= 0; i < T.length(); i++)
724    v[i]= 0;
725  bool noSubset= false;
726  CFArray TT;
727  int subsetDeg;
728  TT= copy (factors);
729  bool recombination= false;
730  while (T.length() >= 2*s)
731  {
732    while (noSubset == false)
733    {
734      if (T.length() == s)
735      {
736        if (recombination)
737        {
738          T.insert (LCBuf);
739          g= prodMod (T, M);
740          T.removeFirst();
741          result.append (g/myContent (g));
742          g= reverseShift (g, evaluation);
743          g /= Lc (g);
744          appendTestMapDown (result, g, info, source, dest);
745          return result;
746        }
747        else
748        {
749          buf= reverseShift (buf, evaluation);
750          return CFList (buf);
751        }
752        return result;
753      }
754
755      S= subset (v, s, TT, noSubset);
756      if (noSubset) break;
757
758      S.insert (LCBuf);
759      g= prodMod (S, M);
760      S.removeFirst();
761      g /= myContent (g);
762      if (fdivides (g, buf))
763      {
764        buf2= reverseShift (g, evaluation);
765        buf2 /= Lc (buf2);
766        if (!k && beta == Variable (1))
767        {
768          if (degree (buf2, alpha) < degMipoBeta)
769          {
770            appendTestMapDown (result, buf2, info, source, dest);
771            buf /= g;
772            LCBuf= LC (buf, Variable (1));
773            recombination= true;
774          }
775        }
776        else
777        {
778          if (!isInExtension (buf2, delta, k))
779          {
780            appendTestMapDown (result, buf2, info, source, dest);
781            buf /= g;
782            LCBuf= LC (buf, Variable (1));
783            recombination= true;
784          }
785        }
786        T= Difference (T, S);
787
788        if (T.length() < 2*s || T.length() == s)
789        {
790          buf= reverseShift (buf, evaluation);
791          buf /= Lc (buf);
792          appendTestMapDown (result, buf, info, source, dest);
793          return result;
794        }
795        TT= copy (T);
796        indexUpdate (v, s, T.length(), noSubset);
797        if (noSubset) break;
798      }
799    }
800    s++;
801    if (T.length() < 2*s || T.length() == s)
802    {
803      buf= reverseShift (buf, evaluation);
804      appendTestMapDown (result, buf, info, source, dest);
805      return result;
806    }
807    int v [T.length()];
808    for (int i= 0; i < T.length(); i++)
809      v[i]= 0;
810    noSubset= false;
811  }
812  if (T.length() < 2*s)
813  {
814    buf= reverseShift (F, evaluation);
815    appendMapDown (result, buf, info, source, dest);
816  }
817
818  return result;
819}
820
821CFList
822factorRecombination (const CanonicalForm& F, const CFList& factors,
823                     const CFList& M)
824{
825  if (factors.length() == 1)
826    return CFList(F);
827  if (factors.length() < 1)
828    return CFList();
829
830  CFList T, S;
831
832  T= factors;
833
834  int s= 1;
835  CFList result;
836  CanonicalForm LCBuf= LC (F, Variable (1));
837  CanonicalForm g, buf= F;
838  int v [T.length()];
839  for (int i= 0; i < T.length(); i++)
840    v[i]= 0;
841  bool noSubset= false;
842  CFArray TT;
843  int subsetDeg;
844  TT= copy (factors);
845  Variable y= F.level() - 1;
846  bool recombination= false;
847  CanonicalForm h;
848  while (T.length() >= 2*s)
849  {
850    while (noSubset == false)
851    {
852      if (T.length() == s)
853      {
854        if (recombination)
855        {
856          T.insert (LC (buf));
857          g= prodMod (T, M);
858          result.append (g/myContent (g));
859          return result;
860        }
861        else
862          return CFList (F);
863      }
864      S= subset (v, s, TT, noSubset);
865      if (noSubset) break;
866      S.insert (LCBuf);
867      g= prodMod (S, M);
868      S.removeFirst();
869      g /= myContent (g);
870      if (fdivides (g, buf))
871      {
872        recombination= true;
873        result.append (g);
874        buf /= g;
875        LCBuf= LC (buf, Variable(1));
876        T= Difference (T, S);
877        if (T.length() < 2*s || T.length() == s)
878        {
879          result.append (buf);
880          return result;
881        }
882        TT= copy (T);
883        indexUpdate (v, s, T.length(), noSubset);
884        if (noSubset) break;
885      }
886    }
887    s++;
888    if (T.length() < 2*s || T.length() == s)
889    {
890      result.append (buf);
891      return result;
892    }
893    int v [T.length()];
894    for (int i= 0; i < T.length(); i++)
895      v[i]= 0;
896    noSubset= false;
897  }
898  if (T.length() < 2*s)
899    result.append (F);
900
901  return result;
902}
903
904int
905liftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
906                   success, const int deg, const CFList& MOD, const int bound)
907{
908  int adaptedLiftBound= 0;
909  CanonicalForm buf= F;
910  Variable y= F.mvar();
911  CanonicalForm LCBuf= LC (buf, Variable (1));
912  CanonicalForm g;
913  CFList M= MOD;
914  M.append (power (y, deg));
915  int d= bound;
916  int e= 0;
917  int nBuf;
918  for (CFListIterator i= factors; i.hasItem(); i++)
919  {
920    g= mulMod (i.getItem(), LCBuf, M);
921    g /= myContent (g);
922    if (fdivides (g, buf))
923    {
924      nBuf= degree (g, y) + degree (LC (g, 1), y);
925      d -= nBuf;
926      e= tmax (e, nBuf);
927      buf /= g;
928      LCBuf= LC (buf, Variable (1));
929    }
930  }
931  adaptedLiftBound= d;
932
933  if (adaptedLiftBound < deg)
934  {
935    if (adaptedLiftBound < degree (F) + 1)
936    {
937      if (d == 1)
938      {
939        if (e + 1 > deg)
940        {
941          adaptedLiftBound= deg;
942          success= false;
943        }
944        else
945        {
946          success= true;
947          if (e + 1 < degree (F) + 1)
948            adaptedLiftBound= deg;
949          else
950            adaptedLiftBound= e + 1;
951        }
952      }
953      else
954      {
955        success= true;
956        adaptedLiftBound= deg;
957      }
958    }
959    else
960    {
961      success= true;
962    }
963  }
964  return adaptedLiftBound;
965}
966
967int
968extLiftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
969                      success, const ExtensionInfo& info, const CFList& eval,
970                      const int deg, const CFList& MOD, const int bound)
971{
972  Variable alpha= info.getAlpha();
973  Variable beta= info.getBeta();
974  CanonicalForm gamma= info.getGamma();
975  CanonicalForm delta= info.getDelta();
976  int k= info.getGFDegree();
977  int adaptedLiftBound= 0;
978  CanonicalForm buf= F;
979  Variable y= F.mvar();
980  CanonicalForm LCBuf= LC (buf, Variable (1));
981  CanonicalForm g, gg;
982  CFList M= MOD;
983  M.append (power (y, deg));
984  adaptedLiftBound= 0;
985  int d= bound;
986  int e= 0;
987  int nBuf;
988  int degMipoBeta;
989  if (!k && beta.level() == 1)
990    degMipoBeta= 1;
991  else if (!k && beta.level() != 1)
992    degMipoBeta= degree (getMipo (beta));
993
994  for (CFListIterator i= factors; i.hasItem(); i++)
995  {
996    g= mulMod (i.getItem(), LCBuf, M);
997    g /= myContent (g);
998    if (fdivides (g, buf))
999    {
1000      gg= reverseShift (g, eval);
1001      gg /= Lc (gg);
1002      if (!k && beta == Variable (1))
1003      {
1004        if (degree (gg, alpha) < degMipoBeta)
1005        {
1006          buf /= g;
1007          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1008          d -= nBuf;
1009          e= tmax (e, nBuf);
1010          LCBuf= LC (buf, Variable (1));
1011        }
1012      }
1013      else
1014      {
1015        if (!isInExtension (gg, delta, k))
1016        {
1017          buf /= g;
1018          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1019          d -= nBuf;
1020          e= tmax (e, nBuf);
1021          LCBuf= LC (buf, Variable (1));
1022        }
1023      }
1024    }
1025  }
1026  adaptedLiftBound= d;
1027
1028  if (adaptedLiftBound < deg)
1029  {
1030    if (adaptedLiftBound < degree (F) + 1)
1031    {
1032      if (d == 1)
1033      {
1034        if (e + 1 > deg)
1035        {
1036          adaptedLiftBound= deg;
1037          success= false;
1038        }
1039        else
1040        {
1041          success= true;
1042          if (e + 1 < degree (F) + 1)
1043            adaptedLiftBound= deg;
1044          else
1045            adaptedLiftBound= e + 1;
1046        }
1047      }
1048      else
1049      {
1050        success= true;
1051        adaptedLiftBound= deg;
1052      }
1053    }
1054    else
1055    {
1056      success= true;
1057    }
1058  }
1059
1060  return adaptedLiftBound;
1061}
1062
1063CFList
1064earlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1065                   bool& success, const int deg, const CFList& MOD,
1066                   const int bound)
1067{
1068  CFList result;
1069  CFList T= factors;
1070  CanonicalForm buf= F;
1071  Variable y= F.mvar();
1072  CanonicalForm LCBuf= LC (buf, Variable (1));
1073  CanonicalForm g;
1074  CFList M= MOD;
1075  M.append (power (y, deg));
1076  adaptedLiftBound= 0;
1077  int d= bound;
1078  int e= 0;
1079  int nBuf;
1080  for (CFListIterator i= factors; i.hasItem(); i++)
1081  {
1082    g= mulMod (i.getItem(), LCBuf, M);
1083    g /= myContent (g);
1084    if (fdivides (g, buf))
1085    {
1086      result.append (g);
1087      nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1088      d -= nBuf;
1089      e= tmax (e, nBuf);
1090      buf /= g;
1091      LCBuf= LC (buf, Variable (1));
1092      T= Difference (T, CFList (i.getItem()));
1093    }
1094  }
1095  adaptedLiftBound= d;
1096
1097  if (adaptedLiftBound < deg)
1098  {
1099    if (adaptedLiftBound < degree (F) + 1)
1100    {
1101      if (d == 1)
1102        adaptedLiftBound= tmin (e + 1, deg);
1103      else
1104        adaptedLiftBound= deg;
1105    }
1106    factors= T;
1107    F= buf;
1108    success= true;
1109  }
1110  return result;
1111}
1112
1113CFList
1114extEarlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1115                      bool& success, const ExtensionInfo& info, const CFList&
1116                      eval, const int deg, const CFList& MOD, const int bound)
1117{
1118  Variable alpha= info.getAlpha();
1119  Variable beta= info.getBeta();
1120  CanonicalForm gamma= info.getGamma();
1121  CanonicalForm delta= info.getDelta();
1122  int k= info.getGFDegree();
1123  CFList result;
1124  CFList T= factors;
1125  CanonicalForm buf= F;
1126  Variable y= F.mvar();
1127  CanonicalForm LCBuf= LC (buf, Variable (1));
1128  CanonicalForm g, gg;
1129  CFList M= MOD;
1130  M.append (power (y, deg));
1131  adaptedLiftBound= 0;
1132  int d= bound;
1133  int e= 0;
1134  int nBuf;
1135  CFList source, dest;
1136
1137  int degMipoBeta;
1138  if (!k && beta.level() == 1)
1139    degMipoBeta= 1;
1140  else if (!k && beta.level() != 1)
1141    degMipoBeta= degree (getMipo (beta));
1142
1143  for (CFListIterator i= factors; i.hasItem(); i++)
1144  {
1145    g= mulMod (i.getItem(), LCBuf, M);
1146    g /= myContent (g);
1147    if (fdivides (g, buf))
1148    {
1149      gg= reverseShift (g, eval);
1150      gg /= Lc (gg);
1151          if (!k && beta == Variable (1))
1152          {
1153            if (degree (gg, alpha) < degMipoBeta)
1154            {
1155              appendTestMapDown (result, gg, info, source, dest);
1156              buf /= g;
1157              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1158              d -= nBuf;
1159              e= tmax (e, nBuf);
1160              LCBuf= LC (buf, Variable (1));
1161            }
1162          }
1163          else
1164          {
1165            if (!isInExtension (gg, delta, k))
1166            {
1167              appendTestMapDown (result, gg, info, source, dest);
1168              buf /= g;
1169              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1170              d -= nBuf;
1171              e= tmax (e, nBuf);
1172              LCBuf= LC (buf, Variable (1));
1173            }
1174          }
1175      T= Difference (T, CFList (i.getItem()));
1176    }
1177  }
1178  adaptedLiftBound= d;
1179
1180  if (adaptedLiftBound < deg)
1181  {
1182    if (adaptedLiftBound < degree (F) + 1)
1183    {
1184      if (d == 1)
1185        adaptedLiftBound= tmin (e + 1, deg);
1186      else
1187        adaptedLiftBound= deg;
1188    }
1189    success= true;
1190    factors= T;
1191    F= buf;
1192  }
1193  return result;
1194}
1195
1196CFList
1197evalPoints (const CanonicalForm& F, CFList & eval, const Variable& alpha,
1198            CFList& list, const bool& GF, bool& fail)
1199{
1200  int k= F.level() - 1;
1201  Variable x= Variable (1);
1202  CFList result;
1203  FFRandom genFF;
1204  GFRandom genGF;
1205  int p= getCharacteristic ();
1206  double bound;
1207  if (alpha != Variable (1))
1208  {
1209    bound= pow ((double) p, (double) degree (getMipo(alpha)));
1210    bound= pow ((double) bound, (double) k);
1211  }
1212  else if (GF)
1213  {
1214    bound= pow ((double) p, (double) getGFDegree());
1215    bound= pow ((double) bound, (double) k);
1216  }
1217  else
1218    bound= pow ((double) p, (double) k);
1219
1220  CanonicalForm random;
1221  CanonicalForm deriv_x, gcd_deriv;
1222  do
1223  {
1224    random= 0;
1225    // possible overflow if list.length() does not fit into a int
1226    if (list.length() >= bound)
1227    {
1228      fail= true;
1229      break;
1230    }
1231    for (int i= 0; i < k; i++)
1232    {
1233      if (list.isEmpty())
1234        result.append (0);
1235      else if (GF)
1236      {
1237        result.append (genGF.generate());
1238        random += result.getLast()*power (x, i);
1239      }
1240      else if (alpha.level() != 1)
1241      {
1242        AlgExtRandomF genAlgExt (alpha);
1243        result.append (genAlgExt.generate());
1244        random += result.getLast()*power (x, i);
1245      }
1246      else
1247      {
1248        result.append (genFF.generate());
1249        random += result.getLast()*power (x, i);
1250      }
1251    }
1252    if (find (list, random))
1253    {
1254      result= CFList();
1255      continue;
1256    }
1257    int l= F.level();
1258    eval.insert (F);
1259    for (CFListIterator i= result; i.hasItem(); i++, l--)
1260      eval.insert (eval.getFirst()(i.getItem(), l));
1261
1262    if (degree (eval.getFirst()) != degree (F, 1))
1263    {
1264      if (!find (list, random))
1265        list.append (random);
1266      result= CFList();
1267      eval= CFList();
1268      continue;
1269    }
1270
1271    deriv_x= deriv (eval.getFirst(), x);
1272    gcd_deriv= gcd (eval.getFirst(), deriv_x);
1273    if (degree (gcd_deriv) > 0)
1274    {
1275      if (!find (list, random))
1276        list.append (random);
1277      result= CFList();
1278      eval= CFList();
1279      continue;
1280    }
1281    CFListIterator i= eval;
1282    i++;
1283    CanonicalForm contentx= content (i.getItem(), x);
1284    if (degree (contentx) > 0)
1285    {
1286      if (!find (list, random))
1287        list.append (random);
1288      result= CFList();
1289      eval= CFList();
1290      continue;
1291    }
1292
1293    if (list.length() >= bound)
1294    {
1295      fail= true;
1296      break;
1297    }
1298  } while (find (list, random));
1299
1300  if (!eval.isEmpty())
1301    eval.removeFirst();
1302
1303  return result;
1304}
1305
1306static inline
1307int newMainVariableSearch (CanonicalForm& A, CFList& Aeval, CFList&
1308                           evaluation, const Variable& alpha, const int lev)
1309{
1310  Variable x= Variable (1);
1311  CanonicalForm derivI, buf;
1312  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1313  int swapLevel= 0;
1314  CFList list;
1315  bool fail= false;
1316  buf= A;
1317  Aeval= CFList();
1318  evaluation= CFList();
1319  for (int i= lev; i <= A.level(); i++)
1320  {
1321    derivI= deriv (buf, Variable (i));
1322    if (!derivI.isZero())
1323    {
1324      buf= swapvar (buf, x, Variable (i));
1325      Aeval= CFList();
1326      evaluation= CFList();
1327      fail= false;
1328      evaluation= evalPoints (buf, Aeval, alpha, list, GF, fail);
1329      if (!fail)
1330      {
1331        A= buf;
1332        swapLevel= i;
1333        break;
1334      }
1335      else
1336        buf= A;
1337    }
1338  }
1339  return swapLevel;
1340}
1341
1342static inline
1343CanonicalForm lcmContent (const CanonicalForm& A, CFList& contentAi)
1344{
1345  int i= A.level();
1346  contentAi.append (myContent (A, i));
1347  contentAi.append (myContent (A, i - 1));
1348  CanonicalForm result= myLcm (contentAi.getFirst(), contentAi.getLast());
1349  for (i= i - 2; i > 0; i--)
1350  {
1351    contentAi.append (content (A, i));
1352    result= myLcm (result, contentAi.getLast());
1353  }
1354  return result;
1355}
1356
1357CFList
1358henselLiftAndEarly (CanonicalForm& A, CFList& MOD, int*& liftBounds, bool&
1359                    earlySuccess, CFList& earlyFactors, const CFList& Aeval,
1360                    const CFList& biFactors, const CFList& evaluation,
1361                    const ExtensionInfo& info)
1362{
1363  bool extension= info.isInExtension();
1364  CFList bufFactors= biFactors;
1365  bufFactors.insert (LC (Aeval.getFirst(), 1));
1366
1367  sortList (bufFactors, Variable (1));
1368
1369  CFList diophant;
1370  CFArray Pi;
1371  int smallFactorDeg= 11; //tunable parameter
1372  CFList result;
1373  int newLiftBound= 0;
1374  int adaptedLiftBound= 0;
1375  int liftBound= liftBounds[1];
1376
1377  earlySuccess= false;
1378  bool earlyReconst= false;
1379  CFList earlyReconstFactors;
1380  CFListIterator j= Aeval;
1381  j++;
1382  CanonicalForm buf= j.getItem();
1383  CFMatrix Mat= CFMatrix (liftBound, bufFactors.length() - 1);
1384  MOD= CFList (power (Variable (2), liftBounds[0]));
1385  if (smallFactorDeg >= liftBound)
1386  {
1387    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1388  }
1389  else if (smallFactorDeg >= degree (buf) + 1)
1390  {
1391    liftBounds[1]= degree (buf) + 1;
1392    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1393    if (Aeval.length() == 2)
1394    {
1395      if (!extension)
1396        earlyFactors= earlyFactorDetect
1397                       (buf, result, adaptedLiftBound, earlySuccess,
1398                        degree (buf) + 1, MOD, liftBound);
1399      else
1400        earlyFactors= extEarlyFactorDetect
1401                       (buf, result, adaptedLiftBound, earlySuccess,
1402                        info, evaluation, degree
1403                        (buf) + 1, MOD, liftBound);
1404    }
1405    else
1406    {
1407      if (!extension)
1408        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1409                                             degree (buf) + 1, MOD, liftBound);
1410      else
1411        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1412                                                evaluation, degree (buf) + 1,
1413                                                MOD, liftBound);
1414    }
1415    if (!earlySuccess)
1416    {
1417      result.insert (LC (buf, 1));
1418      liftBounds[1]= adaptedLiftBound;
1419      liftBound= adaptedLiftBound;
1420      henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1421                        Pi, diophant, Mat, MOD);
1422    }
1423    else
1424      liftBounds[1]= adaptedLiftBound;
1425  }
1426  else if (smallFactorDeg < degree (buf) + 1)
1427  {
1428    liftBounds[1]= smallFactorDeg;
1429    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1430    if (Aeval.length() == 2)
1431    {
1432      if (!extension)
1433        earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1434                                         earlySuccess, smallFactorDeg, MOD,
1435                                         liftBound);
1436      else
1437        earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1438                                            earlySuccess, info, evaluation,
1439                                            smallFactorDeg, MOD, liftBound);
1440    }
1441    else
1442    {
1443      if (!extension)
1444        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1445                                             smallFactorDeg, MOD, liftBound);
1446      else
1447        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1448                                                evaluation, smallFactorDeg, MOD,
1449                                                liftBound);
1450    }
1451
1452    if (!earlySuccess)
1453    {
1454      result.insert (LC (buf, 1));
1455      henselLiftResume (buf, result, smallFactorDeg, degree (buf) + 1,
1456                        Pi, diophant, Mat, MOD);
1457      if (Aeval.length() == 2)
1458      {
1459         if (!extension)
1460           earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1461                                            earlySuccess, degree (buf) + 1,
1462                                            MOD, liftBound);
1463         else
1464           earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1465                                               earlySuccess, info, evaluation,
1466                                               degree (buf) + 1, MOD,
1467                                               liftBound);
1468      }
1469      else
1470      {
1471        if (!extension)
1472          adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1473                                               degree (buf) + 1, MOD,liftBound);
1474        else
1475          adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess,
1476                                                  info, evaluation,
1477                                                  degree (buf) + 1, MOD,
1478                                                  liftBound);
1479      }
1480      if (!earlySuccess)
1481      {
1482        result.insert (LC (buf, 1));
1483        liftBounds[1]= adaptedLiftBound;
1484        liftBound= adaptedLiftBound;
1485        henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1486                          Pi, diophant, Mat, MOD);
1487      }
1488      else
1489        liftBounds[1]= adaptedLiftBound;
1490    }
1491    else
1492      liftBounds[1]= adaptedLiftBound;
1493  }
1494
1495  MOD.append (power (Variable (3), liftBounds[1]));
1496
1497  if (Aeval.length() > 2)
1498  {
1499    CFListIterator j= Aeval;
1500    j++;
1501    CFList bufEval;
1502    bufEval.append (j.getItem());
1503    j++;
1504    int liftBoundsLength= Aeval.getLast().level() - 1;
1505    for (int i= 2; i <= liftBoundsLength && j.hasItem(); i++, j++)
1506    {
1507      earlySuccess= false;
1508      result.insert (LC (bufEval.getFirst(), 1));
1509      bufEval.append (j.getItem());
1510      liftBound= liftBounds[i];
1511      Mat= CFMatrix (liftBounds[i], result.length() - 1);
1512
1513      buf= j.getItem();
1514      if (smallFactorDeg >= liftBound)
1515        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1516                            liftBounds[i -  1], liftBounds[i]);
1517      else if (smallFactorDeg >= degree (buf) + 1)
1518      {
1519        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1520                            liftBounds[i -  1], degree (buf) + 1);
1521
1522        if (Aeval.length() == i + 1)
1523        {
1524          if (!extension)
1525            earlyFactors= earlyFactorDetect
1526                           (buf, result, adaptedLiftBound, earlySuccess,
1527                            degree (buf) + 1, MOD, liftBound);
1528          else
1529            earlyFactors= extEarlyFactorDetect
1530                           (buf, result, adaptedLiftBound, earlySuccess,
1531                            info, evaluation, degree (buf) + 1, MOD, liftBound);
1532        }
1533        else
1534        {
1535          if (!extension)
1536            adaptedLiftBound= liftBoundAdaption
1537                                (buf, result, earlySuccess, degree (buf)
1538                                 + 1,  MOD, liftBound);
1539          else
1540            adaptedLiftBound= extLiftBoundAdaption
1541                                (buf, result, earlySuccess, info, evaluation,
1542                                 degree (buf) + 1, MOD, liftBound);
1543        }
1544
1545        if (!earlySuccess)
1546        {
1547          result.insert (LC (buf, 1));
1548          liftBounds[i]= adaptedLiftBound;
1549          liftBound= adaptedLiftBound;
1550          henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1551                            Pi, diophant, Mat, MOD);
1552        }
1553        else
1554        {
1555          liftBounds[i]= adaptedLiftBound;
1556        }
1557      }
1558      else if (smallFactorDeg < degree (buf) + 1)
1559      {
1560        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1561                            liftBounds[i -  1], smallFactorDeg);
1562
1563        if (Aeval.length() == i + 1)
1564        {
1565          if (!extension)
1566            earlyFactors= earlyFactorDetect
1567                           (buf, result, adaptedLiftBound, earlySuccess,
1568                            smallFactorDeg, MOD, liftBound);
1569          else
1570            earlyFactors= extEarlyFactorDetect
1571                           (buf, result, adaptedLiftBound, earlySuccess,
1572                            info, evaluation, smallFactorDeg, MOD, liftBound);
1573        }
1574        else
1575        {
1576          if (!extension)
1577            adaptedLiftBound= liftBoundAdaption
1578                                (buf, result, earlySuccess,
1579                                 smallFactorDeg, MOD, liftBound);
1580          else
1581            adaptedLiftBound= extLiftBoundAdaption
1582                                (buf, result, earlySuccess, info, evaluation,
1583                                 smallFactorDeg, MOD, liftBound);
1584        }
1585
1586        if (!earlySuccess)
1587        {
1588          result.insert (LC (buf, 1));
1589          henselLiftResume (buf, result, smallFactorDeg,
1590                            degree (buf) + 1, Pi, diophant, Mat, MOD);
1591          if (Aeval.length() == i + 1)
1592          {
1593            if (!extension)
1594              earlyFactors= earlyFactorDetect
1595                             (buf, result, adaptedLiftBound, earlySuccess,
1596                              degree (buf) +  1,  MOD, liftBound);
1597            else
1598              earlyFactors= extEarlyFactorDetect
1599                             (buf, result, adaptedLiftBound, earlySuccess,
1600                              info, evaluation, degree (buf) + 1, MOD,
1601                              liftBound);
1602          }
1603          else
1604          {
1605            if (!extension)
1606              adaptedLiftBound= liftBoundAdaption
1607                                  (buf, result, earlySuccess, degree
1608                                   (buf) +  1,  MOD, liftBound);
1609            else
1610              adaptedLiftBound= extLiftBoundAdaption
1611                                  (buf, result, earlySuccess, info, evaluation,
1612                                   degree (buf) + 1,  MOD, liftBound);
1613          }
1614
1615          if (!earlySuccess)
1616          {
1617            result.insert (LC (buf, 1));
1618            liftBounds[i]= adaptedLiftBound;
1619            liftBound= adaptedLiftBound;
1620            henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1621                              Pi, diophant, Mat, MOD);
1622          }
1623          else
1624            liftBounds[i]= adaptedLiftBound;
1625        }
1626        else
1627          liftBounds[i]= adaptedLiftBound;
1628      }
1629      MOD.append (power (Variable (i + 2), liftBounds[i]));
1630      bufEval.removeFirst();
1631    }
1632    bufFactors= result;
1633  }
1634  else
1635    bufFactors= result;
1636
1637  if (earlySuccess)
1638    A= buf;
1639  return result;
1640}
1641
1642CFList
1643extFactorize (const CanonicalForm& F, const ExtensionInfo& info);
1644
1645CFList
1646multiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1647{
1648
1649  if (F.inCoeffDomain())
1650    return CFList (F);
1651
1652  // compress and find main Variable
1653  CFMap N;
1654  CanonicalForm A= myCompress (F, N);
1655
1656  A /= Lc (A); // make monic
1657
1658  Variable alpha= info.getAlpha();
1659  Variable beta= info.getBeta();
1660  CanonicalForm gamma= info.getGamma();
1661  CanonicalForm delta= info.getDelta();
1662  int k= info.getGFDegree();
1663  bool extension= info.isInExtension();
1664  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1665  //univariate case
1666  if (F.isUnivariate())
1667  {
1668    if (extension == false)
1669      return uniFactorizer (F, alpha, GF);
1670    else
1671    {
1672      CFList source, dest;
1673      A= mapDown (F, info, source, dest);
1674      return uniFactorizer (A, beta, GF);
1675    }
1676  }
1677
1678  //bivariate case
1679  if (A.level() == 2)
1680  {
1681    CFList buf= biFactorize (F, info);
1682    return buf;
1683  }
1684
1685  Variable x= Variable (1);
1686  Variable y= Variable (2);
1687
1688  // remove content
1689  CFList contentAi;
1690  CanonicalForm lcmCont= lcmContent (A, contentAi);
1691  A /= lcmCont;
1692
1693  // trivial after content removal
1694  CFList contentAFactors;
1695  if (A.inCoeffDomain())
1696  {
1697    for (CFListIterator i= contentAi; i.hasItem(); i++)
1698    {
1699      if (i.getItem().inCoeffDomain())
1700        continue;
1701      else
1702      {
1703        lcmCont /= i.getItem();
1704        contentAFactors=
1705        Union (multiFactorize (lcmCont, info),
1706               multiFactorize (i.getItem(), info));
1707        break;
1708      }
1709    }
1710    decompress (contentAFactors, N);
1711    normalize (contentAFactors);
1712    return contentAFactors;
1713  }
1714
1715  // factorize content
1716  contentAFactors= multiFactorize (lcmCont, info);
1717
1718  // univariate after content removal
1719  CFList factors;
1720  if (A.isUnivariate ())
1721  {
1722    factors= uniFactorizer (A, alpha, GF);
1723    append (factors, contentAFactors);
1724    decompress (factors, N);
1725    return factors;
1726  }
1727
1728  // check main variable
1729  int swapLevel= 0;
1730  CanonicalForm derivZ;
1731  CanonicalForm gcdDerivZ;
1732  CanonicalForm bufA= A;
1733  Variable z;
1734  for (int i= 1; i <= A.level(); i++)
1735  {
1736    z= Variable (i);
1737    derivZ= deriv (bufA, z);
1738    if (derivZ.isZero())
1739    {
1740      if (i == 1)
1741        swapLevel= 1;
1742      else
1743        continue;
1744    }
1745    else
1746    {
1747      if (swapLevel == 1)
1748      {
1749        swapLevel= i;
1750        A= swapvar (A, x, z);
1751      }
1752      if (GF)
1753        gcdDerivZ= GCD_GF (bufA, derivZ);
1754      else if (alpha == Variable (1))
1755        gcdDerivZ= GCD_small_p (bufA, derivZ);
1756      else
1757        gcdDerivZ= GCD_Fp_extension (bufA, derivZ, alpha);
1758      if (degree (gcdDerivZ) > 0 && !derivZ.isZero())
1759      {
1760        CanonicalForm g= bufA/gcdDerivZ;
1761        CFList factorsG=
1762        Union (multiFactorize (g, info),
1763               multiFactorize (gcdDerivZ, info));
1764        appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
1765        normalize (factorsG);
1766        return factorsG;
1767      }
1768    }
1769  }
1770
1771
1772  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
1773  bool fail= false;
1774  int swapLevel2= 0;
1775  int level;
1776  int factorNums= 3;
1777  CanonicalForm bivarEval;
1778  CFList biFactors, bufBiFactors;
1779  CanonicalForm evalPoly;
1780  int lift, bufLift;
1781  double logarithm= (double) ilog2 (totaldegree (A));
1782  logarithm /= log2exp;
1783  logarithm= ceil (logarithm);
1784  if (factorNums < (int) logarithm)
1785    factorNums= (int) logarithm;
1786  // several bivariate factorizations
1787  for (int i= 0; i < factorNums; i++)
1788  {
1789    bufA= A;
1790    bufAeval= CFList();
1791    bufEvaluation= evalPoints (bufA, bufAeval, alpha, list, GF, fail);
1792    evalPoly= 0;
1793
1794    if (fail && (i == 0))
1795    {
1796      if (!swapLevel)
1797        level= 2;
1798      else
1799        level= swapLevel + 1;
1800
1801      swapLevel2= newMainVariableSearch (A, Aeval, evaluation, alpha, level);
1802
1803      if (!swapLevel2) // need to pass to an extension
1804      {
1805        factors= extFactorize (A, info);
1806        appendSwapDecompress (factors, contentAFactors, N, swapLevel, x);
1807        normalize (factors);
1808        return factors;
1809      }
1810      else
1811      {
1812        fail= false;
1813        bufAeval= Aeval;
1814        bufA= A;
1815        bufEvaluation= evaluation;
1816      }
1817    }
1818    else if (fail && (i > 0))
1819      break;
1820
1821    bivarEval= bufEvaluation.getLast();
1822    bufEvaluation.removeLast();
1823
1824    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
1825
1826    TIMING_START (fac_bi_factorizer);
1827    bufBiFactors= biFactorizer (bufAeval.getFirst(), alpha, bivarEval, bufLift);
1828    TIMING_END_AND_PRINT (fac_bi_factorizer,
1829                          "time for bivariate factorization: ");
1830
1831    if (bufBiFactors.length() == 1)
1832    {
1833      if (extension)
1834      {
1835        CFList source, dest;
1836        A= mapDown (A, info, source, dest);
1837      }
1838      factors.append (A);
1839      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
1840                            swapLevel2, x);
1841      normalize (factors);
1842      return factors;
1843    }
1844
1845    bufEvaluation.append (-bivarEval[0]);
1846    if (i == 0)
1847    {
1848      Aeval= bufAeval;
1849      evaluation= bufEvaluation;
1850      biFactors= bufBiFactors;
1851      lift= bufLift;
1852    }
1853    else
1854    {
1855      if (bufBiFactors.length() < biFactors.length() ||
1856          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())))
1857      {
1858        Aeval= bufAeval;
1859        evaluation= bufEvaluation;
1860        biFactors= bufBiFactors;
1861        lift= bufLift;
1862      }
1863    }
1864    int k= 0;
1865    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
1866      evalPoly += j.getItem()*power (x, k);
1867    list.append (evalPoly);
1868  }
1869
1870  //shifting to zero
1871  A= shift2Zero (A, Aeval, evaluation);
1872
1873  int* liftBounds;
1874  int liftBoundsLength= F.level() - 1;
1875  liftBounds= liftingBounds (A, lift);
1876
1877  CFList MOD;
1878  bool earlySuccess;
1879  CFList earlyFactors;
1880  TIMING_START (fac_hensel_lift);
1881  CFList liftedFactors= henselLiftAndEarly
1882                        (A, MOD, liftBounds, earlySuccess, earlyFactors,
1883                         Aeval, biFactors, evaluation, info);
1884  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1885
1886  if (!extension)
1887  {
1888    TIMING_START (fac_factor_recombination);
1889    factors= factorRecombination (A, liftedFactors, MOD);
1890    TIMING_END_AND_PRINT (fac_factor_recombination,
1891                          "time for factor recombination: ");
1892  }
1893  else
1894  {
1895    TIMING_START (fac_factor_recombination);
1896    factors= extFactorRecombination (liftedFactors, A, MOD, info, evaluation);
1897    TIMING_END_AND_PRINT (fac_factor_recombination,
1898                          "time for factor recombination: ");
1899  }
1900
1901  if (earlySuccess)
1902    factors= Union (factors, earlyFactors);
1903
1904  if (!extension)
1905  {
1906    for (CFListIterator i= factors; i.hasItem(); i++)
1907    {
1908      int kk= Aeval.getLast().level();
1909      for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
1910      {
1911        if (i.getItem().level() < kk)
1912          continue;
1913        i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
1914      }
1915    }
1916  }
1917
1918  swap (factors, swapLevel, swapLevel2, x);
1919  append (factors, contentAFactors);
1920  decompress (factors, N);
1921  normalize (factors);
1922
1923  delete[] liftBounds;
1924
1925  return factors;
1926}
1927
1928/// multivariate factorization over an extension of the initial field
1929CFList
1930extFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1931{
1932  CanonicalForm A= F;
1933
1934  Variable alpha= info.getAlpha();
1935  Variable beta= info.getBeta();
1936  int k= info.getGFDegree();
1937  char cGFName= info.getGFName();
1938  CanonicalForm delta= info.getDelta();
1939  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1940  Variable w= Variable (1);
1941
1942  CFList factors;
1943  if (!GF && alpha == w)  // we are in F_p
1944  {
1945    CFList factors;
1946    bool extension= true;
1947    int p= getCharacteristic();
1948    if (p*p < (1<<16)) // pass to GF if possible
1949    {
1950      setCharacteristic (getCharacteristic(), 2, 'Z');
1951      ExtensionInfo info= ExtensionInfo (extension);
1952      A= A.mapinto();
1953      factors= multiFactorize (A, info);
1954
1955      Variable vBuf= rootOf (gf_mipo);
1956      setCharacteristic (getCharacteristic());
1957      for (CFListIterator j= factors; j.hasItem(); j++)
1958        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1959    }
1960    else  // not able to pass to GF, pass to F_p(\alpha)
1961    {
1962      Variable v= chooseExtension (A, beta);
1963      ExtensionInfo info= ExtensionInfo (v, extension);
1964      factors= multiFactorize (A, info);
1965    }
1966    return factors;
1967  }
1968  else if (!GF && (alpha != w)) // we are in F_p(\alpha)
1969  {
1970    if (k == 1) // need factorization over F_p
1971    {
1972      int extDeg= degree (getMipo (alpha));
1973      extDeg++;
1974      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1975      Variable v= rootOf (mipo);
1976      ExtensionInfo info= ExtensionInfo (v);
1977      factors= biFactorize (A, info);
1978    }
1979    else
1980    {
1981      if (beta == Variable (1))
1982      {
1983        Variable v= chooseExtension (A, alpha);
1984        CanonicalForm primElem, imPrimElem;
1985        bool primFail= false;
1986        Variable vBuf;
1987        primElem= primitiveElement (alpha, vBuf, primFail);
1988        ASSERT (!primFail, "failure in integer factorizer");
1989        if (primFail)
1990          ; //ERROR
1991        else
1992          imPrimElem= mapPrimElem (primElem, vBuf, v);
1993
1994        CFList source, dest;
1995        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1996                                   source, dest);
1997        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1998        factors= biFactorize (bufA, info);
1999      }
2000      else
2001      {
2002        Variable v= chooseExtension (A, alpha);
2003        CanonicalForm primElem, imPrimElem;
2004        bool primFail= false;
2005        Variable vBuf;
2006        ASSERT (!primFail, "failure in integer factorizer");
2007        if (primFail)
2008          ; //ERROR
2009        else
2010          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
2011
2012        CFList source, dest;
2013        CanonicalForm bufA= mapDown (A, info, source, dest);
2014        source= CFList();
2015        dest= CFList();
2016        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
2017        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
2018        factors= biFactorize (bufA, info);
2019      }
2020    }
2021    return factors;
2022  }
2023  else // we are in GF (p^k)
2024  {
2025    int p= getCharacteristic();
2026    int extensionDeg= getGFDegree();
2027    bool extension= true;
2028    if (k == 1) // need factorization over F_p
2029    {
2030      extensionDeg++;
2031      if (pow ((double) p, (double) extensionDeg) < (1<<16))
2032      // pass to GF(p^k+1)
2033      {
2034        setCharacteristic (p);
2035        Variable vBuf= rootOf (gf_mipo);
2036        A= GF2FalphaRep (A, vBuf);
2037        setCharacteristic (p, extensionDeg, 'Z');
2038        ExtensionInfo info= ExtensionInfo (extension);
2039        factors= multiFactorize (A.mapinto(), info);
2040      }
2041      else // not able to pass to another GF, pass to F_p(\alpha)
2042      {
2043        setCharacteristic (p);
2044        Variable vBuf= rootOf (gf_mipo);
2045        A= GF2FalphaRep (A, vBuf);
2046        Variable v= chooseExtension (A, beta);
2047        ExtensionInfo info= ExtensionInfo (v, extension);
2048        factors= multiFactorize (A, info);
2049      }
2050    }
2051    else // need factorization over GF (p^k)
2052    {
2053      if (pow ((double) p, (double) 2*extensionDeg) < (1<<16))
2054      // pass to GF(p^2k)
2055      {
2056        setCharacteristic (p, 2*extensionDeg, 'Z');
2057        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
2058        factors= multiFactorize (GFMapUp (A, extensionDeg), info);
2059        setCharacteristic (p, extensionDeg, cGFName);
2060      }
2061      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
2062      {
2063        setCharacteristic (p);
2064        Variable v1= rootOf (gf_mipo);
2065        A= GF2FalphaRep (A, v1);
2066        Variable v2= chooseExtension (A, v1);
2067        CanonicalForm primElem, imPrimElem;
2068        bool primFail= false;
2069        Variable vBuf;
2070        primElem= primitiveElement (v1, vBuf, primFail);
2071        if (primFail)
2072        {
2073          ; //ERROR
2074        }
2075        else
2076        {
2077          imPrimElem= mapPrimElem (primElem, vBuf, v2);
2078        }
2079        CFList source, dest;
2080        CanonicalForm bufA= mapUp (A, alpha, beta, primElem, imPrimElem,
2081                                     source, dest);
2082        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
2083        factors= multiFactorize (bufA, info);
2084        setCharacteristic (p, k, cGFName);
2085        for (CFListIterator i= factors; i.hasItem(); i++)
2086          i.getItem()= Falpha2GFRep (i.getItem());
2087      }
2088    }
2089    return factors;
2090  }
2091}
2092
2093#endif
2094/* HAVE_NTL */
2095
Note: See TracBrowser for help on using the repository browser.