source: git/factory/facFqFactorize.cc @ 04cdf06

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