source: git/factory/facFqFactorize.cc @ 11bf82

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