source: git/factory/facFqFactorize.cc @ 5dad7c7

fieker-DuValspielwiese
Last change on this file since 5dad7c7 was 589ef64, checked in by Martin Lee <martinlee84@…>, 12 years ago
fix: need to buffer factors
  • Property mode set to 100644
File size: 76.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. Precomputation of leading coefficients is
11 * described in "Sparse Hensel lifting" by E. Kaltofen
12 *
13 * @author Martin Lee
14 *
15 **/
16/*****************************************************************************/
17
18#include "config.h"
19
20#include "cf_assert.h"
21#include "debug.h"
22#include "timing.h"
23
24#include "facFqFactorizeUtil.h"
25#include "facFqFactorize.h"
26#include "cf_random.h"
27#include "facHensel.h"
28#include "cf_gcd_smallp.h"
29#include "cf_map_ext.h"
30#include "algext.h"
31#include "facSparseHensel.h"
32#include "facMul.h"
33
34#ifdef HAVE_NTL
35#include "NTLconvert.h"
36
37TIMING_DEFINE_PRINT(fac_fq_bi_factorizer)
38TIMING_DEFINE_PRINT(fac_fq_hensel_lift)
39TIMING_DEFINE_PRINT(fac_fq_factor_recombination)
40
41static inline
42CanonicalForm
43listGCD (const CFList& L);
44
45static inline
46CanonicalForm
47myContent (const CanonicalForm& F)
48{
49  Variable x= Variable (1);
50  CanonicalForm G= swapvar (F, F.mvar(), x);
51  CFList L;
52  for (CFIterator i= G; i.hasTerms(); i++)
53    L.append (i.coeff());
54  if (L.length() == 2)
55    return swapvar (gcd (L.getFirst(), L.getLast()), F.mvar(), x);
56  if (L.length() == 1)
57    return LC (F, x);
58  return swapvar (listGCD (L), F.mvar(), x);
59}
60
61static inline
62CanonicalForm
63listGCD (const CFList& L)
64{
65  if (L.length() == 0)
66    return 0;
67  if (L.length() == 1)
68    return L.getFirst();
69  if (L.length() == 2)
70    return gcd (L.getFirst(), L.getLast());
71  else
72  {
73    CFList lHi, lLo;
74    CanonicalForm resultHi, resultLo;
75    int length= L.length()/2;
76    int j= 0;
77    for (CFListIterator i= L; j < length; i++, j++)
78      lHi.append (i.getItem());
79    lLo= Difference (L, lHi);
80    resultHi= listGCD (lHi);
81    resultLo= listGCD (lLo);
82    if (resultHi.isOne() || resultLo.isOne())
83      return 1;
84    return gcd (resultHi, resultLo);
85  }
86}
87
88static inline
89CanonicalForm
90myContent (const CanonicalForm& F, const Variable& x)
91{
92  if (degree (F, x) <= 0)
93    return 1;
94  CanonicalForm G= F;
95  bool swap= false;
96  if (x != F.mvar())
97  {
98    swap= true;
99    G= swapvar (F, x, F.mvar());
100  }
101  CFList L;
102  Variable alpha;
103  for (CFIterator i= G; i.hasTerms(); i++)
104    L.append (i.coeff());
105  if (L.length() == 2)
106  {
107    if (swap)
108      return swapvar (gcd (L.getFirst(), L.getLast()), F.mvar(), x);
109    else
110      return gcd (L.getFirst(), L.getLast());
111  }
112  if (L.length() == 1)
113  {
114    return LC (F, x);
115  }
116  if (swap)
117    return swapvar (listGCD (L), F.mvar(), x);
118  else
119    return listGCD (L);
120}
121
122CanonicalForm myCompress (const CanonicalForm& F, CFMap& N)
123{
124  int n= F.level();
125  int * degsf= new int [n + 1];
126  int ** swap;
127  swap= new int* [n + 1];
128  for (int i= 0; i <= n; i++)
129  {
130    degsf[i]= 0;
131    swap [i]= new int [2];
132    swap [i] [0]= 0;
133    swap [i] [1]= 0;
134  }
135  int i= 1;
136  n= 1;
137  degsf= degrees (F, degsf);
138
139  CanonicalForm result= F;
140  while ( i <= F.level() )
141  {
142    while( degsf[i] == 0 ) i++;
143    swap[n][0]= i;
144    swap[n][1]= size (LC (F,i));
145    if (i != n)
146      result= swapvar (result, Variable (n), Variable(i));
147    n++; i++;
148  }
149
150  int buf1, buf2;
151  n--;
152
153  for (i= 1; i < n; i++)
154  {
155    for (int j= 1; j < n - i + 1; j++)
156    {
157      if (swap[j][1] > swap[j + 1][1])
158      {
159        buf1= swap [j + 1] [0];
160        buf2= swap [j + 1] [1];
161        swap[j + 1] [0]= swap[j] [0];
162        swap[j + 1] [1]= swap[j] [1];
163        swap[j][0]= buf1;
164        swap[j][1]= buf2;
165        result= swapvar (result, Variable (j + 1), Variable (j));
166      }
167    }
168  }
169
170  for (i= n; i > 0; i--)
171  {
172    if (i != swap[i] [0])
173      N.newpair (Variable (i), Variable (swap[i] [0]));
174  }
175
176  for (i= 0; i <= n; i++)
177    delete [] swap[i];
178  delete [] swap;
179
180  delete [] degsf;
181
182  return result;
183}
184
185CFList
186extFactorRecombination (const CFList& factors, const CanonicalForm& F,
187                        const CFList& M, const ExtensionInfo& info,
188                        const CFList& evaluation)
189{
190  Variable alpha= info.getAlpha();
191  Variable beta= info.getBeta();
192  CanonicalForm gamma= info.getGamma();
193  CanonicalForm delta= info.getDelta();
194  int k= info.getGFDegree();
195  CFList source, dest;
196  if (factors.length() == 1)
197  {
198    CanonicalForm buf= reverseShift (F, evaluation);
199    return CFList (mapDown (buf, info, source, dest));
200  }
201  if (factors.length() < 1)
202    return CFList();
203
204  int degMipoBeta= 1;
205  if (!k && beta.level() != 1)
206    degMipoBeta= degree (getMipo (beta));
207
208  CFList T, S;
209  T= factors;
210
211  int s= 1;
212  CFList result;
213  CanonicalForm buf;
214
215  buf= F;
216
217  Variable x= Variable (1);
218  CanonicalForm g, LCBuf= LC (buf, x);
219  CanonicalForm buf2, quot;
220  int * v= new int [T.length()];
221  for (int i= 0; i < T.length(); i++)
222    v[i]= 0;
223  bool noSubset= false;
224  CFArray TT;
225  TT= copy (factors);
226  bool recombination= false;
227  bool trueFactor= false;
228  while (T.length() >= 2*s)
229  {
230    while (noSubset == false)
231    {
232      if (T.length() == s)
233      {
234        delete [] v;
235        if (recombination)
236        {
237          T.insert (LCBuf);
238          g= prodMod (T, M);
239          T.removeFirst();
240          result.append (g/myContent (g));
241          g= reverseShift (g, evaluation);
242          g /= Lc (g);
243          appendTestMapDown (result, g, info, source, dest);
244          return result;
245        }
246        else
247        {
248          buf= reverseShift (buf, evaluation);
249          return CFList (buf);
250        }
251      }
252
253      S= subset (v, s, TT, noSubset);
254      if (noSubset) break;
255
256      S.insert (LCBuf);
257      g= prodMod (S, M);
258      S.removeFirst();
259      g /= myContent (g);
260      if (fdivides (g, buf, quot))
261      {
262        buf2= reverseShift (g, evaluation);
263        buf2 /= Lc (buf2);
264        if (!k && beta == x)
265        {
266          if (degree (buf2, alpha) < degMipoBeta)
267          {
268            appendTestMapDown (result, buf2, info, source, dest);
269            buf= quot;
270            LCBuf= LC (buf, x);
271            recombination= true;
272            trueFactor= true;
273          }
274        }
275        else
276        {
277          if (!isInExtension (buf2, gamma, k, delta, source, dest))
278          {
279            appendTestMapDown (result, buf2, info, source, dest);
280            buf /= g;
281            LCBuf= LC (buf, x);
282            recombination= true;
283            trueFactor= true;
284          }
285        }
286
287        if (trueFactor)
288        {
289          T= Difference (T, S);
290
291          if (T.length() < 2*s || T.length() == s)
292          {
293            buf= reverseShift (buf, evaluation);
294            buf /= Lc (buf);
295            appendTestMapDown (result, buf, info, source, dest);
296            delete [] v;
297            return result;
298          }
299          trueFactor= false;
300          TT= copy (T);
301          indexUpdate (v, s, T.length(), noSubset);
302          if (noSubset) break;
303        }
304      }
305    }
306    s++;
307    if (T.length() < 2*s || T.length() == s)
308    {
309      buf= reverseShift (buf, evaluation);
310      appendTestMapDown (result, buf, info, source, dest);
311      delete [] v;
312      return result;
313    }
314    for (int i= 0; i < T.length(); i++)
315      v[i]= 0;
316    noSubset= false;
317  }
318  if (T.length() < 2*s)
319  {
320    buf= reverseShift (F, evaluation);
321    appendMapDown (result, buf, info, source, dest);
322  }
323
324  delete [] v;
325  return result;
326}
327
328CFList
329factorRecombination (const CanonicalForm& F, const CFList& factors,
330                     const CFList& M)
331{
332  if (factors.length() == 1)
333    return CFList(F);
334  if (factors.length() < 1)
335    return CFList();
336
337  CFList T, S;
338
339  T= factors;
340
341  int s= 1;
342  CFList result;
343  CanonicalForm LCBuf= LC (F, Variable (1));
344  CanonicalForm g, buf= F;
345  int * v= new int [T.length()];
346  for (int i= 0; i < T.length(); i++)
347    v[i]= 0;
348  bool noSubset= false;
349  CFArray TT;
350  TT= copy (factors);
351  Variable y= F.level() - 1;
352  bool recombination= false;
353  CanonicalForm h, quot;
354  while (T.length() >= 2*s)
355  {
356    while (noSubset == false)
357    {
358      if (T.length() == s)
359      {
360        delete [] v;
361        if (recombination)
362        {
363          T.insert (LC (buf));
364          g= prodMod (T, M);
365          result.append (g/myContent (g));
366          return result;
367        }
368        else
369          return CFList (F);
370      }
371      S= subset (v, s, TT, noSubset);
372      if (noSubset) break;
373      S.insert (LCBuf);
374      g= prodMod (S, M);
375      S.removeFirst();
376      g /= myContent (g);
377      if (fdivides (g, buf, quot))
378      {
379        recombination= true;
380        result.append (g);
381        buf= quot;
382        LCBuf= LC (buf, Variable(1));
383        T= Difference (T, S);
384        if (T.length() < 2*s || T.length() == s)
385        {
386          result.append (buf);
387          delete [] v;
388          return result;
389        }
390        TT= copy (T);
391        indexUpdate (v, s, T.length(), noSubset);
392        if (noSubset) break;
393      }
394    }
395    s++;
396    if (T.length() < 2*s || T.length() == s)
397    {
398      result.append (buf);
399      delete [] v;
400      return result;
401    }
402    for (int i= 0; i < T.length(); i++)
403      v[i]= 0;
404    noSubset= false;
405  }
406  if (T.length() < 2*s)
407    result.append (F);
408
409  delete [] v;
410  return result;
411}
412
413int
414liftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
415                   success, const int deg, const CFList& MOD, const int bound)
416{
417  int adaptedLiftBound= 0;
418  CanonicalForm buf= F;
419  Variable y= F.mvar();
420  Variable x= Variable (1);
421  CanonicalForm LCBuf= LC (buf, x);
422  CanonicalForm g, quot;
423  CFList M= MOD;
424  M.append (power (y, deg));
425  int d= bound;
426  int e= 0;
427  int nBuf;
428  for (CFListIterator i= factors; i.hasItem(); i++)
429  {
430    g= mulMod (i.getItem(), LCBuf, M);
431    g /= myContent (g);
432    if (fdivides (g, buf, quot))
433    {
434      nBuf= degree (g, y) + degree (LC (g, 1), y);
435      d -= nBuf;
436      e= tmax (e, nBuf);
437      buf= quot;
438      LCBuf= LC (buf, x);
439    }
440  }
441  adaptedLiftBound= d;
442
443  if (adaptedLiftBound < deg)
444  {
445    if (adaptedLiftBound < degree (F) + 1)
446    {
447      if (d == 1)
448      {
449        if (e + 1 > deg)
450        {
451          adaptedLiftBound= deg;
452          success= false;
453        }
454        else
455        {
456          success= true;
457          if (e + 1 < degree (F) + 1)
458            adaptedLiftBound= deg;
459          else
460            adaptedLiftBound= e + 1;
461        }
462      }
463      else
464      {
465        success= true;
466        adaptedLiftBound= deg;
467      }
468    }
469    else
470    {
471      success= true;
472    }
473  }
474  return adaptedLiftBound;
475}
476
477int
478extLiftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
479                      success, const ExtensionInfo& info, const CFList& eval,
480                      const int deg, const CFList& MOD, const int bound)
481{
482  Variable alpha= info.getAlpha();
483  Variable beta= info.getBeta();
484  CanonicalForm gamma= info.getGamma();
485  CanonicalForm delta= info.getDelta();
486  int k= info.getGFDegree();
487  int adaptedLiftBound= 0;
488  CanonicalForm buf= F;
489  Variable y= F.mvar();
490  Variable x= Variable (1);
491  CanonicalForm LCBuf= LC (buf, x);
492  CanonicalForm g, gg, quot;
493  CFList M= MOD;
494  M.append (power (y, deg));
495  adaptedLiftBound= 0;
496  int d= bound;
497  int e= 0;
498  int nBuf;
499  int degMipoBeta= 1;
500  if (!k && beta.level() != 1)
501    degMipoBeta= degree (getMipo (beta));
502
503  CFList source, dest;
504  for (CFListIterator i= factors; i.hasItem(); i++)
505  {
506    g= mulMod (i.getItem(), LCBuf, M);
507    g /= myContent (g);
508    if (fdivides (g, buf, quot))
509    {
510      gg= reverseShift (g, eval);
511      gg /= Lc (gg);
512      if (!k && beta == x)
513      {
514        if (degree (gg, alpha) < degMipoBeta)
515        {
516          buf= quot;
517          nBuf= degree (g, y) + degree (LC (g, x), y);
518          d -= nBuf;
519          e= tmax (e, nBuf);
520          LCBuf= LC (buf, x);
521        }
522      }
523      else
524      {
525        if (!isInExtension (gg, gamma, k, delta, source, dest))
526        {
527          buf= quot;
528          nBuf= degree (g, y) + degree (LC (g, x), y);
529          d -= nBuf;
530          e= tmax (e, nBuf);
531          LCBuf= LC (buf, x);
532        }
533      }
534    }
535  }
536  adaptedLiftBound= d;
537
538  if (adaptedLiftBound < deg)
539  {
540    if (adaptedLiftBound < degree (F) + 1)
541    {
542      if (d == 1)
543      {
544        if (e + 1 > deg)
545        {
546          adaptedLiftBound= deg;
547          success= false;
548        }
549        else
550        {
551          success= true;
552          if (e + 1 < degree (F) + 1)
553            adaptedLiftBound= deg;
554          else
555            adaptedLiftBound= e + 1;
556        }
557      }
558      else
559      {
560        success= true;
561        adaptedLiftBound= deg;
562      }
563    }
564    else
565    {
566      success= true;
567    }
568  }
569
570  return adaptedLiftBound;
571}
572
573CFList
574earlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
575                   bool& success, const int deg, const CFList& MOD,
576                   const int bound)
577{
578  CFList result;
579  CFList T= factors;
580  CanonicalForm buf= F;
581  Variable y= F.mvar();
582  Variable x= Variable (1);
583  CanonicalForm LCBuf= LC (buf, x);
584  CanonicalForm g, quot;
585  CFList M= MOD;
586  M.append (power (y, deg));
587  adaptedLiftBound= 0;
588  int d= bound;
589  int e= 0;
590  int nBuf;
591  for (CFListIterator i= factors; i.hasItem(); i++)
592  {
593    g= mulMod (i.getItem(), LCBuf, M);
594    g /= myContent (g);
595    if (fdivides (g, buf, quot))
596    {
597      result.append (g);
598      nBuf= degree (g, y) + degree (LC (g, x), y);
599      d -= nBuf;
600      e= tmax (e, nBuf);
601      buf= quot;
602      LCBuf= LC (buf, x);
603      T= Difference (T, CFList (i.getItem()));
604    }
605  }
606  adaptedLiftBound= d;
607
608  if (adaptedLiftBound < deg)
609  {
610    if (adaptedLiftBound < degree (F) + 1)
611    {
612      if (d == 1)
613        adaptedLiftBound= tmin (e + 1, deg);
614      else
615        adaptedLiftBound= deg;
616    }
617    factors= T;
618    F= buf;
619    success= true;
620  }
621  return result;
622}
623
624CFList
625extEarlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
626                      bool& success, const ExtensionInfo& info, const CFList&
627                      eval, const int deg, const CFList& MOD, const int bound)
628{
629  Variable alpha= info.getAlpha();
630  Variable beta= info.getBeta();
631  CanonicalForm gamma= info.getGamma();
632  CanonicalForm delta= info.getDelta();
633  int k= info.getGFDegree();
634  CFList result;
635  CFList T= factors;
636  CanonicalForm buf= F;
637  Variable y= F.mvar();
638  Variable x= Variable (1);
639  CanonicalForm LCBuf= LC (buf, x);
640  CanonicalForm g, gg, quot;
641  CFList M= MOD;
642  M.append (power (y, deg));
643  adaptedLiftBound= 0;
644  int d= bound;
645  int e= 0;
646  int nBuf;
647  CFList source, dest;
648
649  int degMipoBeta= 1;
650  if (!k && beta.level() != 1)
651    degMipoBeta= degree (getMipo (beta));
652
653  for (CFListIterator i= factors; i.hasItem(); i++)
654  {
655    g= mulMod (i.getItem(), LCBuf, M);
656    g /= myContent (g);
657    if (fdivides (g, buf, quot))
658    {
659      gg= reverseShift (g, eval);
660      gg /= Lc (gg);
661      if (!k && beta == x)
662      {
663        if (degree (gg, alpha) < degMipoBeta)
664        {
665          appendTestMapDown (result, gg, info, source, dest);
666          buf= quot;
667          nBuf= degree (g, y) + degree (LC (g, x), y);
668          d -= nBuf;
669          e= tmax (e, nBuf);
670          LCBuf= LC (buf, x);
671          T= Difference (T, CFList (i.getItem()));
672        }
673      }
674      else
675      {
676        if (!isInExtension (gg, gamma, k, delta, source, dest))
677        {
678          appendTestMapDown (result, gg, info, source, dest);
679          buf= quot;
680          nBuf= degree (g, y) + degree (LC (g, x), y);
681          d -= nBuf;
682          e= tmax (e, nBuf);
683          LCBuf= LC (buf, x);
684          T= Difference (T, CFList (i.getItem()));
685         }
686      }
687    }
688  }
689  adaptedLiftBound= d;
690
691  if (adaptedLiftBound < deg)
692  {
693    if (adaptedLiftBound < degree (F) + 1)
694    {
695      if (d == 1)
696        adaptedLiftBound= tmin (e + 1, deg);
697      else
698        adaptedLiftBound= deg;
699    }
700    success= true;
701    factors= T;
702    F= buf;
703  }
704  return result;
705}
706
707CFList
708evalPoints (const CanonicalForm& F, CFList & eval, const Variable& alpha,
709            CFList& list, const bool& GF, bool& fail)
710{
711  int k= F.level() - 1;
712  Variable x= Variable (1);
713  CFList result;
714  FFRandom genFF;
715  GFRandom genGF;
716  int p= getCharacteristic ();
717  double bound;
718  if (alpha != x)
719  {
720    bound= pow ((double) p, (double) degree (getMipo(alpha)));
721    bound= pow ((double) bound, (double) k);
722  }
723  else if (GF)
724  {
725    bound= pow ((double) p, (double) getGFDegree());
726    bound= pow ((double) bound, (double) k);
727  }
728  else
729    bound= pow ((double) p, (double) k);
730
731  CanonicalForm random;
732  CanonicalForm deriv_x, gcd_deriv;
733  do
734  {
735    random= 0;
736    // possible overflow if list.length() does not fit into a int
737    if (list.length() >= bound)
738    {
739      fail= true;
740      break;
741    }
742    for (int i= 0; i < k; i++)
743    {
744      if (list.isEmpty())
745        result.append (0);
746      else if (GF)
747      {
748        result.append (genGF.generate());
749        random += result.getLast()*power (x, i);
750      }
751      else if (alpha.level() != 1)
752      {
753        AlgExtRandomF genAlgExt (alpha);
754        result.append (genAlgExt.generate());
755        random += result.getLast()*power (x, i);
756      }
757      else
758      {
759        result.append (genFF.generate());
760        random += result.getLast()*power (x, i);
761      }
762    }
763    if (find (list, random))
764    {
765      result= CFList();
766      continue;
767    }
768    int l= F.level();
769    eval.insert (F);
770    bool bad= false;
771    for (CFListIterator i= result; i.hasItem(); i++, l--)
772    {
773      eval.insert (eval.getFirst()(i.getItem(), l));
774      if (degree (eval.getFirst(), l - 1) != degree (F, l - 1))
775      {
776        if (!find (list, random))
777          list.append (random);
778        result= CFList();
779        eval= CFList();
780        bad= true;
781        break;
782      }
783    }
784
785    if (bad)
786      continue;
787
788    if (degree (eval.getFirst()) != degree (F, 1))
789    {
790      if (!find (list, random))
791        list.append (random);
792      result= CFList();
793      eval= CFList();
794      continue;
795    }
796
797    deriv_x= deriv (eval.getFirst(), x);
798    gcd_deriv= gcd (eval.getFirst(), deriv_x);
799    if (degree (gcd_deriv) > 0)
800    {
801      if (!find (list, random))
802        list.append (random);
803      result= CFList();
804      eval= CFList();
805      continue;
806    }
807    CFListIterator i= eval;
808    i++;
809    CanonicalForm contentx= content (i.getItem(), x);
810    if (degree (contentx) > 0)
811    {
812      if (!find (list, random))
813        list.append (random);
814      result= CFList();
815      eval= CFList();
816      continue;
817    }
818
819    if (list.length() >= bound)
820    {
821      fail= true;
822      break;
823    }
824  } while (find (list, random));
825
826  if (!eval.isEmpty())
827    eval.removeFirst();
828
829  return result;
830}
831
832static inline
833int newMainVariableSearch (CanonicalForm& A, CFList& Aeval, CFList&
834                           evaluation, const Variable& alpha, const int lev,
835                           CanonicalForm& g
836                          )
837{
838  Variable x= Variable (1);
839  CanonicalForm derivI, buf;
840  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
841  int swapLevel= 0;
842  CFList list;
843  bool fail= false;
844  buf= A;
845  Aeval= CFList();
846  evaluation= CFList();
847  for (int i= lev; i <= A.level(); i++)
848  {
849    derivI= deriv (buf, Variable (i));
850    if (!derivI.isZero())
851    {
852      g= gcd (buf, derivI);
853      if (degree (g) > 0)
854        return -1;
855
856      buf= swapvar (buf, x, Variable (i));
857      Aeval= CFList();
858      evaluation= CFList();
859      fail= false;
860      evaluation= evalPoints (buf, Aeval, alpha, list, GF, fail);
861      if (!fail)
862      {
863        A= buf;
864        swapLevel= i;
865        break;
866      }
867      else
868        buf= A;
869    }
870  }
871  return swapLevel;
872}
873
874CanonicalForm lcmContent (const CanonicalForm& A, CFList& contentAi)
875{
876  int i= A.level();
877  contentAi.append (myContent (A, i));
878  contentAi.append (myContent (A, i - 1));
879  CanonicalForm result= lcm (contentAi.getFirst(), contentAi.getLast());
880  for (i= i - 2; i > 0; i--)
881  {
882    contentAi.append (content (A, i));
883    result= lcm (result, contentAi.getLast());
884  }
885  return result;
886}
887
888CFList
889henselLiftAndEarly (CanonicalForm& A, CFList& MOD, int*& liftBounds, bool&
890                    earlySuccess, CFList& earlyFactors, const CFList& Aeval,
891                    const CFList& biFactors, const CFList& evaluation,
892                    const ExtensionInfo& info)
893{
894  bool extension= info.isInExtension();
895  CFList bufFactors= biFactors;
896  bufFactors.insert (LC (Aeval.getFirst(), 1));
897
898  sortList (bufFactors, Variable (1));
899
900  CFList diophant;
901  CFArray Pi;
902  int smallFactorDeg= 11; //tunable parameter
903  CFList result;
904  int adaptedLiftBound= 0;
905  int liftBound= liftBounds[1];
906
907  earlySuccess= false;
908  CFList earlyReconstFactors;
909  CFListIterator j= Aeval;
910  j++;
911  CanonicalForm buf= j.getItem();
912  CFMatrix Mat= CFMatrix (liftBound, bufFactors.length() - 1);
913  MOD= CFList (power (Variable (2), liftBounds[0]));
914  if (smallFactorDeg >= liftBound)
915  {
916    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
917  }
918  else if (smallFactorDeg >= degree (buf) + 1)
919  {
920    liftBounds[1]= degree (buf) + 1;
921    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
922    if (Aeval.length() == 2)
923    {
924      if (!extension)
925        earlyFactors= earlyFactorDetect
926                       (buf, result, adaptedLiftBound, earlySuccess,
927                        degree (buf) + 1, MOD, liftBound);
928      else
929        earlyFactors= extEarlyFactorDetect
930                       (buf, result, adaptedLiftBound, earlySuccess,
931                        info, evaluation, degree
932                        (buf) + 1, MOD, liftBound);
933    }
934    else
935    {
936      if (!extension)
937        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
938                                             degree (buf) + 1, MOD, liftBound);
939      else
940        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
941                                                evaluation, degree (buf) + 1,
942                                                MOD, liftBound);
943    }
944    if (!earlySuccess)
945    {
946      result.insert (LC (buf, 1));
947      liftBounds[1]= adaptedLiftBound;
948      liftBound= adaptedLiftBound;
949      henselLiftResume (buf, result, degree (buf) + 1, liftBound,
950                        Pi, diophant, Mat, MOD);
951    }
952    else
953      liftBounds[1]= adaptedLiftBound;
954  }
955  else if (smallFactorDeg < degree (buf) + 1)
956  {
957    liftBounds[1]= smallFactorDeg;
958    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
959    if (Aeval.length() == 2)
960    {
961      if (!extension)
962        earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
963                                         earlySuccess, smallFactorDeg, MOD,
964                                         liftBound);
965      else
966        earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
967                                            earlySuccess, info, evaluation,
968                                            smallFactorDeg, MOD, liftBound);
969    }
970    else
971    {
972      if (!extension)
973        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
974                                             smallFactorDeg, MOD, liftBound);
975      else
976        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
977                                                evaluation, smallFactorDeg, MOD,
978                                                liftBound);
979    }
980
981    if (!earlySuccess)
982    {
983      result.insert (LC (buf, 1));
984      henselLiftResume (buf, result, smallFactorDeg, degree (buf) + 1,
985                        Pi, diophant, Mat, MOD);
986      if (Aeval.length() == 2)
987      {
988         if (!extension)
989           earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
990                                            earlySuccess, degree (buf) + 1,
991                                            MOD, liftBound);
992         else
993           earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
994                                               earlySuccess, info, evaluation,
995                                               degree (buf) + 1, MOD,
996                                               liftBound);
997      }
998      else
999      {
1000        if (!extension)
1001          adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1002                                               degree (buf) + 1, MOD,liftBound);
1003        else
1004          adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess,
1005                                                  info, evaluation,
1006                                                  degree (buf) + 1, MOD,
1007                                                  liftBound);
1008      }
1009      if (!earlySuccess)
1010      {
1011        result.insert (LC (buf, 1));
1012        liftBounds[1]= adaptedLiftBound;
1013        liftBound= adaptedLiftBound;
1014        henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1015                          Pi, diophant, Mat, MOD);
1016      }
1017      else
1018        liftBounds[1]= adaptedLiftBound;
1019    }
1020    else
1021      liftBounds[1]= adaptedLiftBound;
1022  }
1023
1024  MOD.append (power (Variable (3), liftBounds[1]));
1025
1026  if (Aeval.length() > 2)
1027  {
1028    CFListIterator j= Aeval;
1029    j++;
1030    CFList bufEval;
1031    bufEval.append (j.getItem());
1032    j++;
1033    int liftBoundsLength= Aeval.getLast().level() - 1;
1034    for (int i= 2; i <= liftBoundsLength && j.hasItem(); i++, j++)
1035    {
1036      earlySuccess= false;
1037      result.insert (LC (bufEval.getFirst(), 1));
1038      bufEval.append (j.getItem());
1039      liftBound= liftBounds[i];
1040      Mat= CFMatrix (liftBounds[i], result.length() - 1);
1041
1042      buf= j.getItem();
1043      if (smallFactorDeg >= liftBound)
1044        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1045                            liftBounds[i -  1], liftBounds[i]);
1046      else if (smallFactorDeg >= degree (buf) + 1)
1047      {
1048        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1049                            liftBounds[i -  1], degree (buf) + 1);
1050
1051        if (Aeval.length() == i + 1)
1052        {
1053          if (!extension)
1054            earlyFactors= earlyFactorDetect
1055                           (buf, result, adaptedLiftBound, earlySuccess,
1056                            degree (buf) + 1, MOD, liftBound);
1057          else
1058            earlyFactors= extEarlyFactorDetect
1059                           (buf, result, adaptedLiftBound, earlySuccess,
1060                            info, evaluation, degree (buf) + 1, MOD, liftBound);
1061        }
1062        else
1063        {
1064          if (!extension)
1065            adaptedLiftBound= liftBoundAdaption
1066                                (buf, result, earlySuccess, degree (buf)
1067                                 + 1,  MOD, liftBound);
1068          else
1069            adaptedLiftBound= extLiftBoundAdaption
1070                                (buf, result, earlySuccess, info, evaluation,
1071                                 degree (buf) + 1, MOD, liftBound);
1072        }
1073
1074        if (!earlySuccess)
1075        {
1076          result.insert (LC (buf, 1));
1077          liftBounds[i]= adaptedLiftBound;
1078          liftBound= adaptedLiftBound;
1079          henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1080                            Pi, diophant, Mat, MOD);
1081        }
1082        else
1083        {
1084          liftBounds[i]= adaptedLiftBound;
1085        }
1086      }
1087      else if (smallFactorDeg < degree (buf) + 1)
1088      {
1089        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1090                            liftBounds[i -  1], smallFactorDeg);
1091
1092        if (Aeval.length() == i + 1)
1093        {
1094          if (!extension)
1095            earlyFactors= earlyFactorDetect
1096                           (buf, result, adaptedLiftBound, earlySuccess,
1097                            smallFactorDeg, MOD, liftBound);
1098          else
1099            earlyFactors= extEarlyFactorDetect
1100                           (buf, result, adaptedLiftBound, earlySuccess,
1101                            info, evaluation, smallFactorDeg, MOD, liftBound);
1102        }
1103        else
1104        {
1105          if (!extension)
1106            adaptedLiftBound= liftBoundAdaption
1107                                (buf, result, earlySuccess,
1108                                 smallFactorDeg, MOD, liftBound);
1109          else
1110            adaptedLiftBound= extLiftBoundAdaption
1111                                (buf, result, earlySuccess, info, evaluation,
1112                                 smallFactorDeg, MOD, liftBound);
1113        }
1114
1115        if (!earlySuccess)
1116        {
1117          result.insert (LC (buf, 1));
1118          henselLiftResume (buf, result, smallFactorDeg,
1119                            degree (buf) + 1, Pi, diophant, Mat, MOD);
1120          if (Aeval.length() == i + 1)
1121          {
1122            if (!extension)
1123              earlyFactors= earlyFactorDetect
1124                             (buf, result, adaptedLiftBound, earlySuccess,
1125                              degree (buf) +  1,  MOD, liftBound);
1126            else
1127              earlyFactors= extEarlyFactorDetect
1128                             (buf, result, adaptedLiftBound, earlySuccess,
1129                              info, evaluation, degree (buf) + 1, MOD,
1130                              liftBound);
1131          }
1132          else
1133          {
1134            if (!extension)
1135              adaptedLiftBound= liftBoundAdaption
1136                                  (buf, result, earlySuccess, degree
1137                                   (buf) +  1,  MOD, liftBound);
1138            else
1139              adaptedLiftBound= extLiftBoundAdaption
1140                                  (buf, result, earlySuccess, info, evaluation,
1141                                   degree (buf) + 1,  MOD, liftBound);
1142          }
1143
1144          if (!earlySuccess)
1145          {
1146            result.insert (LC (buf, 1));
1147            liftBounds[i]= adaptedLiftBound;
1148            liftBound= adaptedLiftBound;
1149            henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1150                              Pi, diophant, Mat, MOD);
1151          }
1152          else
1153            liftBounds[i]= adaptedLiftBound;
1154        }
1155        else
1156          liftBounds[i]= adaptedLiftBound;
1157      }
1158      MOD.append (power (Variable (i + 2), liftBounds[i]));
1159      bufEval.removeFirst();
1160    }
1161    bufFactors= result;
1162  }
1163  else
1164    bufFactors= result;
1165
1166  if (earlySuccess)
1167    A= buf;
1168  return result;
1169}
1170
1171void
1172gcdFreeBasis (CFFList& factors1, CFFList& factors2)
1173{
1174  CanonicalForm g;
1175  int k= factors1.length();
1176  int l= factors2.length();
1177  int n= 1;
1178  int m;
1179  CFFListIterator j;
1180  for (CFFListIterator i= factors1; (n < k && i.hasItem()); i++, n++)
1181  {
1182    m= 1;
1183    for (j= factors2; (m < l && j.hasItem()); j++, m++)
1184    {
1185      g= gcd (i.getItem().factor(), j.getItem().factor());
1186      if (degree (g,1) > 0)
1187      {
1188        j.getItem()= CFFactor (j.getItem().factor()/g, j.getItem().exp());
1189        i.getItem()= CFFactor (i.getItem().factor()/g, i.getItem().exp());
1190        factors1.append (CFFactor (g, i.getItem().exp()));
1191        factors2.append (CFFactor (g, j.getItem().exp()));
1192      }
1193    }
1194  }
1195}
1196
1197CFList
1198distributeContent (const CFList& L, const CFList* differentSecondVarFactors,
1199                   int length
1200                  )
1201{
1202  CFList l= L;
1203  CanonicalForm content= l.getFirst();
1204
1205  if (content.inCoeffDomain())
1206    return l;
1207
1208  if (l.length() == 1)
1209  {
1210    CFList result;
1211    for (int i= 0; i < length; i++)
1212    {
1213      if (differentSecondVarFactors[i].isEmpty())
1214        continue;
1215      if (result.isEmpty())
1216      {
1217        result= differentSecondVarFactors[i];
1218        for (CFListIterator iter= result; iter.hasItem(); iter++)
1219          content /= iter.getItem();
1220      }
1221      else
1222      {
1223        CFListIterator iter1= result;
1224        for (CFListIterator iter2= differentSecondVarFactors[i];iter2.hasItem();
1225             iter2++, iter1++)
1226        {
1227          iter1.getItem() *= iter2.getItem();
1228          content /= iter2.getItem();
1229        }
1230      }
1231    }
1232    result.insert (content);
1233    return result;
1234  }
1235
1236  Variable v;
1237  CFListIterator iter1, iter2;
1238  CanonicalForm tmp, g;
1239  CFList multiplier;
1240  for (int i= 0; i < length; i++)
1241  {
1242    if (differentSecondVarFactors[i].isEmpty())
1243      continue;
1244    iter1= l;
1245    iter1++;
1246
1247    v= Variable (i + 3);
1248    tmp= 1;
1249    for (iter2= differentSecondVarFactors[i]; iter2.hasItem();
1250         iter2++, iter1++)
1251    {
1252      if (degree (iter2.getItem(),v) == degree (iter1.getItem(),v))
1253      {
1254        multiplier.append (1);
1255        continue;
1256      }
1257      g= gcd (iter2.getItem(), content);
1258      if (!g.inCoeffDomain())
1259      {
1260        tmp *= g;
1261        multiplier.append (g);
1262      }
1263      else
1264        multiplier.append (1);
1265    }
1266    if (!tmp.isOne() && fdivides (tmp, content))
1267    {
1268      iter1= l;
1269      iter1++;
1270      content /= tmp;
1271      for (iter2= multiplier; iter2.hasItem(); iter1++, iter2++)
1272        iter1.getItem() *= iter2.getItem();
1273    }
1274    multiplier= CFList();
1275  }
1276
1277  l.removeFirst();
1278  l.insert (content);
1279  return l;
1280}
1281
1282int
1283testFactors (const CanonicalForm& G, const CFList& uniFactors,
1284             const Variable& alpha, CanonicalForm& sqrfPartF, CFList& factors,
1285             CFFList*& bufSqrfFactors, CFList& evalSqrfPartF,
1286             const CFArray& evalPoint)
1287{
1288  CanonicalForm tmp;
1289  CFListIterator j;
1290  for (CFListIterator i= uniFactors; i.hasItem(); i++)
1291  {
1292    tmp= i.getItem();
1293    if (i.hasItem())
1294      i++;
1295    else
1296      break;
1297    for (j= i; j.hasItem(); j++)
1298    {
1299      if (tmp == j.getItem())
1300        return 0;
1301    }
1302  }
1303
1304  CanonicalForm F= G;
1305  CFFList sqrfFactorization= squarefreeFactorization (F, alpha);
1306
1307  sqrfPartF= 1;
1308  for (CFFListIterator i= sqrfFactorization; i.hasItem(); i++)
1309    sqrfPartF *= i.getItem().factor();
1310
1311  evalSqrfPartF= evaluateAtEval (sqrfPartF, evalPoint);
1312
1313  CanonicalForm test= evalSqrfPartF.getFirst() (evalPoint[0], 2);
1314
1315  if (degree (test) != degree (sqrfPartF, 1) || test.inCoeffDomain())
1316    return 0;
1317
1318  CFFList sqrfFactors;
1319  CFList tmp2;
1320  int k= 0;
1321  factors= uniFactors;
1322  CFFListIterator iter;
1323  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1324  {
1325    tmp= 1;
1326    sqrfFactors= squarefreeFactorization (i.getItem(), alpha);
1327
1328    for (iter= sqrfFactors; iter.hasItem(); iter++)
1329    {
1330      tmp2.append (iter.getItem().factor());
1331      tmp *= iter.getItem().factor();
1332    }
1333    i.getItem()= tmp/Lc(tmp);
1334    bufSqrfFactors [k]= sqrfFactors;
1335  }
1336
1337  for (int i= 0; i < factors.length() - 1; i++)
1338  {
1339    for (k= i + 1; k < factors.length(); k++)
1340    {
1341      gcdFreeBasis (bufSqrfFactors [i], bufSqrfFactors[k]);
1342    }
1343  }
1344
1345  factors= CFList();
1346  for (int i= 0; i < uniFactors.length(); i++)
1347  {
1348    if (i == 0)
1349    {
1350      for (iter= bufSqrfFactors [i]; iter.hasItem(); iter++)
1351      {
1352        if (iter.getItem().factor().inCoeffDomain())
1353          continue;
1354        iter.getItem()= CFFactor (iter.getItem().factor()/
1355                                  Lc (iter.getItem().factor()),
1356                                  iter.getItem().exp());
1357        factors.append (iter.getItem().factor());
1358      }
1359    }
1360    else
1361    {
1362      for (iter= bufSqrfFactors [i]; iter.hasItem(); iter++)
1363      {
1364        if (iter.getItem().factor().inCoeffDomain())
1365          continue;
1366        iter.getItem()= CFFactor (iter.getItem().factor()/
1367                                  Lc (iter.getItem().factor()),
1368                                  iter.getItem().exp());
1369        if (!find (factors, iter.getItem().factor()))
1370          factors.append (iter.getItem().factor());
1371      }
1372    }
1373  }
1374
1375  test= prod (factors);
1376  tmp= evalSqrfPartF.getFirst() (evalPoint[0],2);
1377  if (test/Lc (test) != tmp/Lc (tmp))
1378    return 0;
1379  else
1380    return 1;
1381}
1382
1383CFList
1384precomputeLeadingCoeff (const CanonicalForm& LCF, const CFList& LCFFactors,
1385                        const Variable& alpha, const CFList& evaluation,
1386                        CFList* & differentSecondVarLCs, int lSecondVarLCs,
1387                        Variable& y
1388                       )
1389{
1390  y= Variable (1);
1391  if (LCF.inCoeffDomain())
1392  {
1393    CFList result;
1394    for (int i= 1; i <= LCFFactors.length() + 1; i++)
1395      result.append (1);
1396    return result;
1397  }
1398
1399  CFMap N, M;
1400  CFArray dummy= CFArray (2);
1401  dummy [0]= LCF;
1402  dummy [1]= Variable (2);
1403  compress (dummy, M, N);
1404  CanonicalForm F= M (LCF);
1405  if (LCF.isUnivariate())
1406  {
1407    CFList result;
1408    int LCFLevel= LCF.level();
1409    bool found= false;
1410    if (LCFLevel == 2)
1411    {
1412    //bivariate leading coefficients are already the true leading coefficients
1413      result= LCFFactors;
1414      found= true;
1415    }
1416    else
1417    {
1418      CFListIterator j;
1419      for (int i= 0; i < lSecondVarLCs; i++)
1420      {
1421        for (j= differentSecondVarLCs[i]; j.hasItem(); j++)
1422        {
1423          if (j.getItem().level() == LCFLevel)
1424          {
1425            found= true;
1426            break;
1427          }
1428        }
1429        if (found)
1430        {
1431          result= differentSecondVarLCs [i];
1432          break;
1433        }
1434      }
1435      if (!found)
1436        result= LCFFactors;
1437    }
1438    if (found)
1439      result.insert (Lc (LCF));
1440    else
1441      result.append (LCF);
1442    return result;
1443  }
1444
1445  CFList factors= LCFFactors;
1446
1447  for (CFListIterator i= factors; i.hasItem(); i++)
1448    i.getItem()= M (i.getItem());
1449
1450  CanonicalForm sqrfPartF;
1451  CFFList * bufSqrfFactors= new CFFList [factors.length()];
1452  CFList evalSqrfPartF, bufFactors;
1453  CFArray evalPoint= CFArray (evaluation.length() - 1);
1454  CFArray buf= CFArray (evaluation.length());
1455  CFArray swap= CFArray (evaluation.length());
1456  CFListIterator iter= evaluation;
1457  CanonicalForm vars=getVars (LCF)*Variable (2);
1458  for (int i= evaluation.length() +1; i > 1; i--, iter++)
1459  {
1460    buf[i-2]=iter.getItem();
1461    if (degree (vars, i) > 0)
1462      swap[M(Variable (i)).level()-1]=buf[i-2];
1463  }
1464  buf= swap;
1465  for (int i= 0; i < evaluation.length() - 1; i++)
1466    evalPoint[i]= buf[i+1];
1467
1468  int pass= testFactors (F, factors, alpha, sqrfPartF,
1469                         bufFactors, bufSqrfFactors, evalSqrfPartF, evalPoint);
1470
1471  bool foundDifferent= false;
1472  Variable z, x= y;
1473  int j= 0;
1474  if (!pass)
1475  {
1476    int lev= 0;
1477    // LCF is non-constant here
1478    CFList bufBufFactors;
1479    CanonicalForm bufF;
1480    for (int i= 0; i < lSecondVarLCs; i++)
1481    {
1482      if (!differentSecondVarLCs [i].isEmpty())
1483      {
1484        bool allConstant= true;
1485        for (iter= differentSecondVarLCs[i]; iter.hasItem(); iter++)
1486        {
1487          if (!iter.getItem().inCoeffDomain())
1488          {
1489            allConstant= false;
1490            y= Variable (iter.getItem().level());
1491            lev= M(y).level();
1492          }
1493        }
1494        if (allConstant)
1495          continue;
1496
1497        bufFactors= differentSecondVarLCs [i];
1498        for (iter= bufFactors; iter.hasItem(); iter++)
1499          iter.getItem()= swapvar (iter.getItem(), x, y);
1500        bufF= F;
1501        z= Variable (lev);
1502        bufF= swapvar (bufF, x, z);
1503        bufBufFactors= bufFactors;
1504        evalPoint= CFArray (evaluation.length() - 1);
1505        for (int k= 0; k < evaluation.length()-1; k++)
1506        {
1507          if (N (Variable (k+1)).level() != y.level())
1508            evalPoint[k]= buf[k+1];
1509          else
1510            evalPoint[k]= buf[0];
1511        }
1512        pass= testFactors (bufF, bufBufFactors, alpha, sqrfPartF, bufFactors,
1513                           bufSqrfFactors, evalSqrfPartF, evalPoint);
1514        if (pass)
1515        {
1516          foundDifferent= true;
1517          F= bufF;
1518          CFList l= factors;
1519          for (iter= l; iter.hasItem(); iter++)
1520            iter.getItem()= swapvar (iter.getItem(), x, y);
1521          differentSecondVarLCs [i]= l;
1522          j= i;
1523          break;
1524        }
1525        if (!pass && i == lSecondVarLCs - 1)
1526        {
1527          CFList result;
1528          result.append (LCF);
1529          for (int k= 1; k <= factors.length(); k++)
1530            result.append (LCF);
1531          y= Variable (1);
1532          delete [] bufSqrfFactors;
1533          return result;
1534        }
1535      }
1536    }
1537  }
1538  if (!pass)
1539  {
1540    CFList result;
1541    result.append (LCF);
1542    for (int k= 1; k <= factors.length(); k++)
1543      result.append (LCF);
1544    y= Variable (1);
1545    delete [] bufSqrfFactors;
1546    return result;
1547  }
1548  else
1549    factors= bufFactors;
1550
1551  bufFactors= factors;
1552
1553  CFMap MM, NN;
1554  dummy [0]= sqrfPartF;
1555  dummy [1]= 1;
1556  compress (dummy, MM, NN);
1557  sqrfPartF= MM (sqrfPartF);
1558  CanonicalForm varsSqrfPartF= getVars (sqrfPartF);
1559  for (CFListIterator iter= factors; iter.hasItem(); iter++)
1560    iter.getItem()= MM (iter.getItem());
1561
1562  CFList evaluation2;
1563  for (int i= 2; i <= varsSqrfPartF.level(); i++)
1564    evaluation2.insert (evalPoint[NN (Variable (i)).level()-2]);
1565
1566  CFList interMedResult;
1567  CanonicalForm oldSqrfPartF= sqrfPartF;
1568  sqrfPartF= shift2Zero (sqrfPartF, evalSqrfPartF, evaluation2);
1569  if (factors.length() > 1)
1570  {
1571    CanonicalForm LC1= LC (oldSqrfPartF, 1);
1572    CFList leadingCoeffs;
1573    for (int i= 0; i < factors.length(); i++)
1574      leadingCoeffs.append (LC1);
1575
1576    CFList LC1eval= evaluateAtEval (LC1, evaluation2, 2);
1577    CFList oldFactors= factors;
1578    for (CFListIterator i= oldFactors; i.hasItem(); i++)
1579      i.getItem() *= LC1eval.getFirst()/Lc (i.getItem());
1580
1581    bool success= false;
1582    CanonicalForm oldSqrfPartFPowLC= oldSqrfPartF*power(LC1,factors.length()-1);
1583    CFList heuResult;
1584    if (size (oldSqrfPartFPowLC)/getNumVars (oldSqrfPartFPowLC) < 500 &&
1585        LucksWangSparseHeuristic (oldSqrfPartFPowLC,
1586                                  oldFactors, 2, leadingCoeffs, heuResult))
1587    {
1588      interMedResult= recoverFactors (oldSqrfPartF, heuResult);
1589      if (oldFactors.length() == interMedResult.length())
1590        success= true;
1591    }
1592    if (!success)
1593    {
1594      LC1= LC (evalSqrfPartF.getFirst(), 1);
1595
1596      CFArray leadingCoeffs= CFArray (factors.length());
1597      for (int i= 0; i < factors.length(); i++)
1598        leadingCoeffs[i]= LC1;
1599
1600      for (CFListIterator i= factors; i.hasItem(); i++)
1601        i.getItem() *= LC1 (0,2)/Lc (i.getItem());
1602      factors.insert (1);
1603
1604      CanonicalForm
1605      newSqrfPartF= evalSqrfPartF.getFirst()*power (LC1, factors.length() - 2);
1606
1607      int liftBound= degree (newSqrfPartF,2) + 1;
1608
1609      CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
1610      CFArray Pi;
1611      CFList diophant;
1612      nonMonicHenselLift12 (newSqrfPartF, factors, liftBound, Pi, diophant, M,
1613                            leadingCoeffs, false);
1614
1615      if (sqrfPartF.level() > 2)
1616      {
1617        int* liftBounds= new int [sqrfPartF.level() - 1];
1618        liftBounds [0]= liftBound;
1619        bool noOneToOne= false;
1620        CFList *leadingCoeffs2= new CFList [sqrfPartF.level()-2];
1621        LC1= LC (evalSqrfPartF.getLast(), 1);
1622        CFList LCs;
1623        for (int i= 0; i < factors.length(); i++)
1624          LCs.append (LC1);
1625        leadingCoeffs2 [sqrfPartF.level() - 3]= LCs;
1626        for (int i= sqrfPartF.level() - 1; i > 2; i--)
1627        {
1628          for (CFListIterator j= LCs; j.hasItem(); j++)
1629            j.getItem()= j.getItem() (0, i + 1);
1630          leadingCoeffs2 [i - 3]= LCs;
1631        }
1632        sqrfPartF *= power (LC1, factors.length()-1);
1633
1634        int liftBoundsLength= sqrfPartF.level() - 1;
1635        for (int i= 1; i < liftBoundsLength; i++)
1636          liftBounds [i]= degree (sqrfPartF, i + 2) + 1;
1637        evalSqrfPartF= evaluateAtZero (sqrfPartF);
1638        evalSqrfPartF.removeFirst();
1639        factors= nonMonicHenselLift (evalSqrfPartF, factors, leadingCoeffs2,
1640                 diophant, Pi, liftBounds, sqrfPartF.level() - 1, noOneToOne);
1641        delete [] leadingCoeffs2;
1642        delete [] liftBounds;
1643      }
1644      for (CFListIterator iter= factors; iter.hasItem(); iter++)
1645        iter.getItem()= reverseShift (iter.getItem(), evaluation2);
1646
1647      interMedResult=
1648      recoverFactors (reverseShift(evalSqrfPartF.getLast(),evaluation2),
1649                      factors);
1650    }
1651  }
1652  else
1653  {
1654    CanonicalForm contF=content (oldSqrfPartF,1);
1655    factors= CFList (oldSqrfPartF/contF);
1656    interMedResult= recoverFactors (oldSqrfPartF, factors);
1657  }
1658
1659  for (CFListIterator iter= interMedResult; iter.hasItem(); iter++)
1660    iter.getItem()= NN (iter.getItem());
1661
1662  CFList result;
1663  CFFListIterator k;
1664  for (int i= 0; i < LCFFactors.length(); i++)
1665  {
1666    CanonicalForm tmp= 1;
1667    for (k= bufSqrfFactors[i]; k.hasItem(); k++)
1668    {
1669      int pos= findItem (bufFactors, k.getItem().factor());
1670      if (pos)
1671        tmp *= power (getItem (interMedResult, pos), k.getItem().exp());
1672    }
1673    result.append (tmp);
1674  }
1675
1676  for (CFListIterator i= result; i.hasItem(); i++)
1677  {
1678    F /= i.getItem();
1679    if (foundDifferent)
1680      i.getItem()= swapvar (i.getItem(), x, z);
1681    i.getItem()= N (i.getItem());
1682  }
1683
1684  if (foundDifferent)
1685  {
1686    CFList l= differentSecondVarLCs [j];
1687    for (CFListIterator i= l; i.hasItem(); i++)
1688      i.getItem()= swapvar (i.getItem(), y, z);
1689    differentSecondVarLCs [j]= l;
1690    F= swapvar (F, x, z);
1691  }
1692
1693  result.insert (N (F));
1694
1695  result= distributeContent (result, differentSecondVarLCs, lSecondVarLCs);
1696
1697  if (!result.getFirst().inCoeffDomain())
1698  {
1699    CFListIterator i= result;
1700    CanonicalForm tmp;
1701    if (foundDifferent)
1702      i.getItem()= swapvar (i.getItem(), Variable (2), y);
1703
1704    tmp= i.getItem();
1705
1706    i++;
1707    for (; i.hasItem(); i++)
1708    {
1709      if (foundDifferent)
1710        i.getItem()= swapvar (i.getItem(), Variable (2), y)*tmp;
1711      else
1712        i.getItem() *= tmp;
1713    }
1714  }
1715  else
1716    y= Variable (1);
1717
1718  delete [] bufSqrfFactors;
1719
1720  return result;
1721}
1722
1723void
1724evaluationWRTDifferentSecondVars (CFList*& Aeval, const CFList& evaluation,
1725                                  const CanonicalForm& A)
1726{
1727  CanonicalForm tmp;
1728  CFList tmp2;
1729  CFListIterator iter;
1730  bool preserveDegree= true;
1731  int j;
1732  for (int i= A.level(); i > 2; i--)
1733  {
1734    tmp= A;
1735    tmp2= CFList();
1736    iter= evaluation;
1737    preserveDegree= true;
1738    for (j= A.level(); j > 1; j--, iter++)
1739    {
1740      if (j == i)
1741        continue;
1742      else
1743      {
1744        tmp= tmp (iter.getItem(), j);
1745        tmp2.insert (tmp);
1746        if ((degree (tmp, i) != degree (A, i)) ||
1747            (degree (tmp, 1) != degree (A, 1)))
1748        {
1749          preserveDegree= false;
1750          break;
1751        }
1752        if (!content(tmp).inCoeffDomain() || !content(tmp,1).inCoeffDomain())
1753        {
1754          preserveDegree= false;
1755          break;
1756        }
1757      }
1758    }
1759    if (preserveDegree)
1760      Aeval [i - 3]= tmp2;
1761    else
1762      Aeval [i - 3]= CFList();
1763  }
1764}
1765
1766#endif
1767
1768static inline
1769CanonicalForm prodEval (const CFList& l, const CanonicalForm& evalPoint,
1770                        const Variable& v)
1771{
1772  CanonicalForm result= 1;
1773  for (CFListIterator i= l; i.hasItem(); i++)
1774    result *= i.getItem() (evalPoint, v);
1775  return result;
1776}
1777
1778//recombine bivariate factors in case one bivariate factorization yields less
1779// factors than the other
1780CFList
1781recombination (const CFList& factors1, const CFList& factors2, int s, int thres,
1782               const CanonicalForm& evalPoint, const Variable& x)
1783{
1784  CFList T, S;
1785
1786  T= factors1;
1787  CFList result;
1788  CanonicalForm buf;
1789  int * v= new int [T.length()];
1790  for (int i= 0; i < T.length(); i++)
1791    v[i]= 0;
1792  bool nosubset= false;
1793  CFArray TT;
1794  TT= copy (factors1);
1795  while (T.length() >= 2*s && s <= thres)
1796  {
1797    while (nosubset == false)
1798    {
1799      if (T.length() == s)
1800      {
1801        delete [] v;
1802        result.append (prod (T));
1803        return result;
1804      }
1805      S= subset (v, s, TT, nosubset);
1806      if (nosubset) break;
1807      buf= prodEval (S, evalPoint, x);
1808      buf /= Lc (buf);
1809      if (find (factors2, buf))
1810      {
1811        T= Difference (T, S);
1812        result.append (prod (S));
1813        TT= copy (T);
1814        indexUpdate (v, s, T.length(), nosubset);
1815        if (nosubset) break;
1816      }
1817    }
1818    s++;
1819    if (T.length() < 2*s || T.length() == s) 
1820    {
1821      delete [] v;
1822      result.append (prod (T));
1823      return result;
1824    }
1825    for (int i= 0; i < T.length(); i++)
1826      v[i]= 0;
1827    nosubset= false;
1828  }
1829
1830  delete [] v;
1831  if (T.length() < 2*s)
1832  {
1833    result.append (prod (T));
1834    return result;
1835  }
1836
1837  return result;
1838}
1839
1840#ifdef HAVE_NTL
1841void
1842factorizationWRTDifferentSecondVars (const CanonicalForm& A, CFList*& Aeval,
1843                                     const ExtensionInfo& info,
1844                                     int& minFactorsLength, bool& irred)
1845{
1846  Variable x= Variable (1);
1847  minFactorsLength= 0;
1848  irred= false;
1849  CFList factors;
1850  Variable v;
1851  for (int j= 0; j < A.level() - 2; j++)
1852  {
1853    if (!Aeval[j].isEmpty())
1854    {
1855      v= Variable (Aeval[j].getFirst().level());
1856      if (CFFactory::gettype() == GaloisFieldDomain)
1857        factors= GFBiSqrfFactorize (Aeval[j].getFirst());
1858      else if (info.getAlpha().level() == 1)
1859        factors= FpBiSqrfFactorize (Aeval[j].getFirst());
1860      else
1861        factors= FqBiSqrfFactorize (Aeval[j].getFirst(), info.getAlpha());
1862
1863      factors.removeFirst();
1864      if (minFactorsLength == 0)
1865        minFactorsLength= factors.length();
1866      else
1867        minFactorsLength= tmin (minFactorsLength, factors.length());
1868
1869      if (factors.length() == 1)
1870      {
1871        irred= true;
1872        return;
1873      }
1874      sortList (factors, x);
1875      Aeval [j]= factors;
1876    }
1877  }
1878}
1879
1880CFList conv (const CFArray & A)
1881{
1882  CFList result;
1883  for (int i= A.max(); i >= A.min(); i--)
1884    result.insert (A[i]);
1885  return result;
1886}
1887
1888
1889void getLeadingCoeffs (const CanonicalForm& A, CFList*& Aeval,
1890                       const CFList& uniFactors, const CFList& evaluation
1891                      )
1892{
1893  CFListIterator iter;
1894  CFList LCs;
1895  for (int j= 0; j < A.level() - 2; j++)
1896  {
1897    if (!Aeval[j].isEmpty())
1898    {
1899      LCs= CFList();
1900      for (iter= Aeval[j]; iter.hasItem(); iter++)
1901        LCs.append (LC (iter.getItem(), 1));
1902      //normalize (LCs);
1903      Aeval[j]= LCs;
1904    }
1905  }
1906}
1907
1908void sortByUniFactors (CFList*& Aeval, int AevalLength,
1909                       const CFList& uniFactors, const CFList& evaluation
1910                      )
1911{
1912  CanonicalForm evalPoint;
1913  int i;
1914  CFListIterator iter, iter2;
1915  Variable v;
1916  CFList LCs, buf;
1917  CFArray l;
1918  int pos, index;
1919  for (int j= 0; j < AevalLength; j++)
1920  {
1921    if (!Aeval[j].isEmpty())
1922    {
1923      i= evaluation.length() + 1;
1924      for (iter= evaluation; iter.hasItem(); iter++, i--)
1925      {
1926        if (i == Aeval[j].getFirst().level())
1927        {
1928          evalPoint= iter.getItem();
1929          break;
1930        }
1931      }
1932
1933      v= Variable (i);
1934      if (Aeval[j].length() > uniFactors.length())
1935        Aeval[j]= recombination (Aeval[j], uniFactors, 1,
1936                                 Aeval[j].length() - uniFactors.length() + 1,
1937                                 evalPoint, v);
1938
1939      buf= buildUniFactors (Aeval[j], evalPoint, v);
1940      l= CFArray (uniFactors.length());
1941      index= 1;
1942      for (iter= buf; iter.hasItem(); iter++, index++)
1943      {
1944        pos= findItem (uniFactors, iter.getItem());
1945        if (pos)
1946          l[pos-1]= getItem (Aeval[j], index);
1947      }
1948      buf= conv (l);
1949      Aeval [j]= buf;
1950
1951      buf= buildUniFactors (Aeval[j], evalPoint, v);
1952    }
1953  }
1954}
1955
1956CFList
1957buildUniFactors (const CFList& biFactors, const CanonicalForm& evalPoint,
1958                 const Variable& y)
1959{
1960  CFList result;
1961  CanonicalForm tmp;
1962  for (CFListIterator i= biFactors; i.hasItem(); i++)
1963  {
1964    tmp= mod (i.getItem(), y - evalPoint);
1965    tmp /= Lc (tmp);
1966    result.append (tmp);
1967  }
1968  return result;
1969}
1970
1971void refineBiFactors (const CanonicalForm& A, CFList& biFactors,
1972                      CFList* const& Aeval, const CFList& evaluation,
1973                      int minFactorsLength)
1974{
1975  CFListIterator iter;
1976  CanonicalForm evalPoint;
1977  int i;
1978  Variable v;
1979  Variable y= Variable (2);
1980  CFList list;
1981  for (int j= 0; j < A.level() - 2; j++)
1982  {
1983    if (Aeval[j].length() == minFactorsLength)
1984    {
1985      i= A.level();
1986
1987      for (iter= evaluation; iter.hasItem(); iter++, i--)
1988      {
1989        if (i == Aeval[j].getFirst().level())
1990        {
1991          evalPoint= iter.getItem();
1992          break;
1993        }
1994      }
1995
1996      v= Variable (i);
1997      list= buildUniFactors (Aeval[j], evalPoint, v);
1998
1999      biFactors= recombination (biFactors, list, 1,
2000                                biFactors.length() - list.length() + 1,
2001                                evaluation.getLast(), y);
2002      return;
2003    }
2004  }
2005}
2006
2007void prepareLeadingCoeffs (CFList*& LCs, int n, const CFList& leadingCoeffs,
2008                           const CFList& biFactors, const CFList& evaluation)
2009{
2010  CFList l= leadingCoeffs;
2011  LCs [n-3]= l;
2012  CFListIterator j;
2013  CFListIterator iter= evaluation;
2014  for (int i= n - 1; i > 2; i--, iter++)
2015  {
2016    for (j= l; j.hasItem(); j++)
2017      j.getItem()= j.getItem() (iter.getItem(), i + 1);
2018    LCs [i - 3]= l;
2019  }
2020  l= LCs [0];
2021  for (CFListIterator i= l; i.hasItem(); i++)
2022    i.getItem()= i.getItem() (iter.getItem(), 3);
2023  CFListIterator ii= biFactors;
2024  CFList normalizeFactor;
2025  for (CFListIterator i= l; i.hasItem(); i++, ii++)
2026    normalizeFactor.append (Lc (LC (ii.getItem(), 1))/Lc (i.getItem()));
2027  for (int i= 0; i < n-2; i++)
2028  {
2029    ii= normalizeFactor;
2030    for (j= LCs [i]; j.hasItem(); j++, ii++)
2031      j.getItem() *= ii.getItem();
2032  }
2033}
2034
2035CFList
2036extNonMonicFactorRecombination (const CFList& factors, const CanonicalForm& F,
2037                                const ExtensionInfo& info)
2038{
2039  Variable alpha= info.getAlpha();
2040  Variable beta= info.getBeta();
2041  CanonicalForm gamma= info.getGamma();
2042  CanonicalForm delta= info.getDelta();
2043  int k= info.getGFDegree();
2044  CFList source, dest;
2045
2046  int degMipoBeta= 1;
2047  if (!k && beta != Variable(1))
2048    degMipoBeta= degree (getMipo (beta));
2049
2050  CFList T, S;
2051  T= factors;
2052  int s= 1;
2053  CFList result;
2054  CanonicalForm quot, buf= F;
2055
2056  CanonicalForm g;
2057  CanonicalForm buf2;
2058  int * v= new int [T.length()];
2059  for (int i= 0; i < T.length(); i++)
2060    v[i]= 0;
2061  bool noSubset= false;
2062  CFArray TT;
2063  TT= copy (factors);
2064  bool recombination= false;
2065  bool trueFactor= false;
2066  while (T.length() >= 2*s)
2067  {
2068    while (noSubset == false)
2069    {
2070      if (T.length() == s)
2071      {
2072        delete [] v;
2073        if (recombination)
2074        {
2075          g= prod (T);
2076          T.removeFirst();
2077          result.append (g/myContent (g));
2078          g /= Lc (g);
2079          appendTestMapDown (result, g, info, source, dest);
2080          return result;
2081        }
2082        else
2083          return CFList (buf/myContent(buf));
2084      }
2085
2086      S= subset (v, s, TT, noSubset);
2087      if (noSubset) break;
2088
2089      g= prod (S);
2090      g /= myContent (g);
2091      if (fdivides (g, buf, quot))
2092      {
2093        buf2= g;
2094        buf2 /= Lc (buf2);
2095        if (!k && beta.level() == 1)
2096        {
2097          if (degree (buf2, alpha) < degMipoBeta)
2098          {
2099            appendTestMapDown (result, buf2, info, source, dest);
2100            buf= quot;
2101            recombination= true;
2102            trueFactor= true;
2103          }
2104        }
2105        else
2106        {
2107          if (!isInExtension (buf2, gamma, k, delta, source, dest))
2108          {
2109            appendTestMapDown (result, buf2, info, source, dest);
2110            buf= quot;
2111            recombination= true;
2112            trueFactor= true;
2113          }
2114        }
2115        if (trueFactor)
2116        {
2117          T= Difference (T, S);
2118
2119          if (T.length() < 2*s || T.length() == s)
2120          {
2121            delete [] v;
2122            buf /= myContent (buf);
2123            buf /= Lc (buf);
2124            appendTestMapDown (result, buf, info, source, dest);
2125            return result;
2126          }
2127          trueFactor= false;
2128          TT= copy (T);
2129          indexUpdate (v, s, T.length(), noSubset);
2130          if (noSubset) break;
2131        }
2132      }
2133    }
2134    s++;
2135    if (T.length() < 2*s || T.length() == s)
2136    {
2137      delete [] v;
2138      appendTestMapDown (result, buf/myContent(buf), info, source, dest);
2139      return result;
2140    }
2141    for (int i= 0; i < T.length(); i++)
2142      v[i]= 0;
2143    noSubset= false;
2144  }
2145  if (T.length() < 2*s)
2146    appendMapDown (result, F/myContent(F), info, source, dest);
2147
2148  delete [] v;
2149  return result;
2150}
2151
2152CFList
2153extFactorize (const CanonicalForm& F, const ExtensionInfo& info);
2154
2155CFList
2156multiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
2157{
2158
2159  if (F.inCoeffDomain())
2160    return CFList (F);
2161
2162  // compress and find main Variable
2163  CFMap N;
2164  CanonicalForm A= myCompress (F, N);
2165
2166  A /= Lc (A); // make monic
2167
2168  Variable alpha= info.getAlpha();
2169  Variable beta= info.getBeta();
2170  CanonicalForm gamma= info.getGamma();
2171  CanonicalForm delta= info.getDelta();
2172  bool extension= info.isInExtension();
2173  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
2174  //univariate case
2175  if (F.isUnivariate())
2176  {
2177    if (extension == false)
2178      return uniFactorizer (F, alpha, GF);
2179    else
2180    {
2181      CFList source, dest;
2182      A= mapDown (F, info, source, dest);
2183      return uniFactorizer (A, beta, GF);
2184    }
2185  }
2186
2187  //bivariate case
2188  if (A.level() == 2)
2189  {
2190    CFList buf= biFactorize (F, info);
2191    return buf;
2192  }
2193
2194  Variable x= Variable (1);
2195  Variable y= Variable (2);
2196
2197  // remove content
2198  CFList contentAi;
2199  CanonicalForm lcmCont= lcmContent (A, contentAi);
2200  A /= lcmCont;
2201
2202  // trivial after content removal
2203  CFList contentAFactors;
2204  if (A.inCoeffDomain())
2205  {
2206    for (CFListIterator i= contentAi; i.hasItem(); i++)
2207    {
2208      if (i.getItem().inCoeffDomain())
2209        continue;
2210      else
2211      {
2212        lcmCont /= i.getItem();
2213        contentAFactors=
2214        Union (multiFactorize (lcmCont, info),
2215               multiFactorize (i.getItem(), info));
2216        break;
2217      }
2218    }
2219    decompress (contentAFactors, N);
2220    normalize (contentAFactors);
2221    return contentAFactors;
2222  }
2223
2224  // factorize content
2225  contentAFactors= multiFactorize (lcmCont, info);
2226
2227  // univariate after content removal
2228  CFList factors;
2229  if (A.isUnivariate ())
2230  {
2231    factors= uniFactorizer (A, alpha, GF);
2232    append (factors, contentAFactors);
2233    decompress (factors, N);
2234    return factors;
2235  }
2236
2237  // check main variable
2238  int swapLevel= 0;
2239  CanonicalForm derivZ;
2240  CanonicalForm gcdDerivZ;
2241  CanonicalForm bufA= A;
2242  Variable z;
2243  for (int i= 1; i <= A.level(); i++)
2244  {
2245    z= Variable (i);
2246    derivZ= deriv (bufA, z);
2247    if (derivZ.isZero())
2248    {
2249      if (i == 1)
2250        swapLevel= 1;
2251      else
2252        continue;
2253    }
2254    else
2255    {
2256      if (swapLevel == 1)
2257      {
2258        swapLevel= i;
2259        bufA= swapvar (A, x, z);
2260      }
2261      gcdDerivZ= gcd (bufA, derivZ);
2262      if (degree (gcdDerivZ) > 0 && !derivZ.isZero())
2263      {
2264        CanonicalForm g= bufA/gcdDerivZ;
2265        CFList factorsG=
2266        Union (multiFactorize (g, info),
2267               multiFactorize (gcdDerivZ, info));
2268        appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
2269        normalize (factorsG);
2270        return factorsG;
2271      }
2272      else
2273      {
2274        A= bufA;
2275        break;
2276      }
2277    }
2278  }
2279
2280
2281  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
2282  bool fail= false;
2283  int swapLevel2= 0;
2284  int level;
2285  int factorNums= 3;
2286  CanonicalForm bivarEval;
2287  CFList biFactors, bufBiFactors;
2288  CanonicalForm evalPoly;
2289  int lift, bufLift;
2290  double logarithm= (double) ilog2 (totaldegree (A));
2291  logarithm /= log2exp;
2292  logarithm= ceil (logarithm);
2293  if (factorNums < (int) logarithm)
2294    factorNums= (int) logarithm;
2295  CFList* bufAeval2= new CFList [A.level() - 2];
2296  CFList* Aeval2= new CFList [A.level() - 2];
2297  int counter;
2298  int differentSecondVar= 0;
2299  // several bivariate factorizations
2300  for (int i= 0; i < factorNums; i++)
2301  {
2302    counter= 0;
2303    bufA= A;
2304    bufAeval= CFList();
2305    bufEvaluation= evalPoints (bufA, bufAeval, alpha, list, GF, fail);
2306    evalPoly= 0;
2307
2308    if (fail && (i == 0))
2309    {
2310      if (!swapLevel)
2311        level= 2;
2312      else
2313        level= swapLevel + 1;
2314
2315      CanonicalForm g;
2316      swapLevel2= newMainVariableSearch (A, Aeval, evaluation, alpha, level, g);
2317
2318      if (!swapLevel2) // need to pass to an extension
2319      {
2320        factors= extFactorize (A, info);
2321        appendSwapDecompress (factors, contentAFactors, N, swapLevel, x);
2322        normalize (factors);
2323        delete [] bufAeval2;
2324        delete [] Aeval2;
2325        return factors;
2326      }
2327      else
2328      {
2329        if (swapLevel2 == -1)
2330        {
2331          CFList factorsG=
2332          Union (multiFactorize (g, info),
2333                 multiFactorize (A/g, info));
2334          appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
2335          normalize (factorsG);
2336          delete [] bufAeval2;
2337          delete [] Aeval2;
2338          return factorsG;
2339        }
2340        fail= false;
2341        bufAeval= Aeval;
2342        bufA= A;
2343        bufEvaluation= evaluation;
2344      }
2345    }
2346    else if (fail && (i > 0))
2347      break;
2348
2349    bivarEval= bufEvaluation.getLast();
2350
2351    evaluationWRTDifferentSecondVars (bufAeval2, bufEvaluation, A);
2352
2353    for (int j= 0; j < A.level() - 2; j++)
2354    {
2355      if (!bufAeval2[j].isEmpty())
2356        counter++;
2357    }
2358
2359    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
2360
2361    TIMING_START (fac_fq_bi_factorizer);
2362    if (!GF && alpha.level() == 1)
2363      bufBiFactors= FpBiSqrfFactorize (bufAeval.getFirst());
2364    else if (GF)
2365      bufBiFactors= GFBiSqrfFactorize (bufAeval.getFirst());
2366    else
2367      bufBiFactors= FqBiSqrfFactorize (bufAeval.getFirst(), alpha);
2368    TIMING_END_AND_PRINT (fac_fq_bi_factorizer,
2369                          "time for bivariate factorization: ");
2370    bufBiFactors.removeFirst();
2371
2372    if (bufBiFactors.length() == 1)
2373    {
2374      if (extension)
2375      {
2376        CFList source, dest;
2377        A= mapDown (A, info, source, dest);
2378      }
2379      factors.append (A);
2380      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
2381                            swapLevel2, x);
2382      normalize (factors);
2383      delete [] bufAeval2;
2384      delete [] Aeval2;
2385      return factors;
2386    }
2387
2388    if (i == 0)
2389    {
2390      Aeval= bufAeval;
2391      evaluation= bufEvaluation;
2392      biFactors= bufBiFactors;
2393      lift= bufLift;
2394      for (int j= 0; j < A.level() - 2; j++)
2395        Aeval2 [j]= bufAeval2 [j];
2396      differentSecondVar= counter;
2397    }
2398    else
2399    {
2400      if (bufBiFactors.length() < biFactors.length() ||
2401          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())) ||
2402          counter > differentSecondVar)
2403      {
2404        Aeval= bufAeval;
2405        evaluation= bufEvaluation;
2406        biFactors= bufBiFactors;
2407        lift= bufLift;
2408        for (int j= 0; j < A.level() - 2; j++)
2409          Aeval2 [j]= bufAeval2 [j];
2410        differentSecondVar= counter;
2411      }
2412    }
2413    int k= 0;
2414    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
2415      evalPoly += j.getItem()*power (x, k);
2416    list.append (evalPoly);
2417  }
2418
2419  delete [] bufAeval2;
2420
2421  sortList (biFactors, x);
2422
2423  int minFactorsLength;
2424  bool irred= false;
2425  factorizationWRTDifferentSecondVars (A, Aeval2, info, minFactorsLength,irred);
2426
2427  if (irred)
2428  {
2429    if (extension)
2430    {
2431      CFList source, dest;
2432      A= mapDown (A, info, source, dest);
2433    }
2434    factors.append (A);
2435    appendSwapDecompress (factors, contentAFactors, N, swapLevel,
2436                          swapLevel2, x);
2437    normalize (factors);
2438    delete [] Aeval2;
2439    return factors;
2440  }
2441
2442  if (minFactorsLength == 0)
2443    minFactorsLength= biFactors.length();
2444  else if (biFactors.length() > minFactorsLength)
2445    refineBiFactors (A, biFactors, Aeval2, evaluation, minFactorsLength);
2446  minFactorsLength= tmin (minFactorsLength, biFactors.length());
2447
2448  if (differentSecondVar == A.level() - 2)
2449  {
2450    bool zeroOccured= false;
2451    for (CFListIterator iter= evaluation; iter.hasItem(); iter++)
2452    {
2453      if (iter.getItem().isZero())
2454      {
2455        zeroOccured= true;
2456        break;
2457      }
2458    }
2459    if (!zeroOccured)
2460    {
2461      factors= sparseHeuristic (A, biFactors, Aeval2, evaluation,
2462                                minFactorsLength);
2463      if (factors.length() == biFactors.length())
2464      {
2465        if (extension)
2466          factors= extNonMonicFactorRecombination (factors, A, info);
2467
2468        appendSwapDecompress (factors, contentAFactors, N, swapLevel,
2469                              swapLevel2, x);
2470        normalize (factors);
2471        delete [] Aeval2;
2472        return factors;
2473      }
2474      else
2475        factors= CFList();
2476      //TODO case where factors.length() > 0
2477    }
2478  }
2479
2480  CFList uniFactors= buildUniFactors (biFactors, evaluation.getLast(), y);
2481
2482  sortByUniFactors (Aeval2, A.level() - 2, uniFactors, evaluation);
2483
2484  CFList * oldAeval= new CFList [A.level() - 2]; //TODO use bufAeval2 for this
2485  for (int i= 0; i < A.level() - 2; i++)
2486    oldAeval[i]= Aeval2[i];
2487
2488  getLeadingCoeffs (A, Aeval2, uniFactors, evaluation);
2489
2490  CFList biFactorsLCs;
2491  for (CFListIterator i= biFactors; i.hasItem(); i++)
2492    biFactorsLCs.append (LC (i.getItem(), 1));
2493
2494  Variable v;
2495  CFList leadingCoeffs= precomputeLeadingCoeff (LC (A, 1), biFactorsLCs, alpha,
2496                                          evaluation, Aeval2, A.level() - 2, v);
2497
2498  if (v.level() != 1)
2499  {
2500    A= swapvar (A, y, v);
2501    int i= A.level();
2502    CanonicalForm evalPoint;
2503    for (CFListIterator iter= evaluation; iter.hasItem(); iter++, i--)
2504    {
2505      if (i == v.level())
2506      {
2507        evalPoint= iter.getItem();
2508        iter.getItem()= evaluation.getLast();
2509        evaluation.removeLast();
2510        evaluation.append (evalPoint);
2511        break;
2512      }
2513    }
2514    for (i= 0; i < A.level() - 2; i++)
2515    {
2516      if (oldAeval[i].isEmpty())
2517        continue;
2518      if (oldAeval[i].getFirst().level() == v.level())
2519      {
2520        CFArray tmp= copy (oldAeval[i]);
2521        oldAeval[i]= biFactors;
2522        for (CFListIterator iter= oldAeval[i]; iter.hasItem(); iter++)
2523          iter.getItem()= swapvar (iter.getItem(), v, y);
2524        for (int ii= 0; ii < tmp.size(); ii++)
2525          tmp[ii]= swapvar (tmp[ii], v, y);
2526        CFArray tmp2= CFArray (tmp.size());
2527        CanonicalForm buf;
2528        for (int ii= 0; ii < tmp.size(); ii++)
2529        {
2530          buf= tmp[ii] (evaluation.getLast(),y);
2531          buf /= Lc (buf);
2532          tmp2[findItem (uniFactors, buf)-1]=tmp[ii];
2533        }
2534        biFactors= CFList();
2535        for (int j= 0; j < tmp2.size(); j++)
2536          biFactors.append (tmp2[j]);
2537      }
2538    }
2539  }
2540
2541  CFListIterator iter;
2542  CanonicalForm oldA= A;
2543  CFList oldBiFactors= biFactors;
2544  if (!leadingCoeffs.getFirst().inCoeffDomain())
2545  {
2546    CanonicalForm tmp= power (leadingCoeffs.getFirst(), biFactors.length() - 1);
2547    A *= tmp;
2548    tmp= leadingCoeffs.getFirst();
2549    iter= evaluation;
2550    for (int i= A.level(); i > 2; i--, iter++)
2551      tmp= tmp (iter.getItem(), i);
2552    if (!tmp.inCoeffDomain())
2553    {
2554      for (CFListIterator i= biFactors; i.hasItem(); i++)
2555      {
2556        i.getItem() *= tmp/LC (i.getItem(), 1);
2557        i.getItem() /= Lc (i.getItem());
2558      }
2559    }
2560  }
2561
2562  leadingCoeffs.removeFirst();
2563
2564  //prepare leading coefficients
2565  CFList* leadingCoeffs2= new CFList [A.level() - 2];
2566  prepareLeadingCoeffs (leadingCoeffs2, A.level(), leadingCoeffs, biFactors,
2567                        evaluation);
2568
2569  Aeval= evaluateAtEval (A, evaluation, 2);
2570  CanonicalForm hh= Lc (Aeval.getFirst());
2571  for (iter= Aeval; iter.hasItem(); iter++)
2572    iter.getItem() /= hh;
2573
2574  A /= hh;
2575
2576  CFList bufFactors= CFList();
2577  if (LucksWangSparseHeuristic (A, biFactors, 2, leadingCoeffs2 [A.level() - 3],
2578                                 factors))
2579  {
2580    int check= biFactors.length();
2581    int * index= new int [factors.length()];
2582    CFList oldFactors= factors;
2583    factors= recoverFactors (A, factors, index);
2584
2585    if (check == factors.length())
2586    {
2587      if (extension)
2588        factors= extNonMonicFactorRecombination (factors, A, info);
2589
2590      if (v.level() != 1)
2591      {
2592        for (iter= factors; iter.hasItem(); iter++)
2593          iter.getItem()= swapvar (iter.getItem(), v, y);
2594      }
2595
2596      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
2597                            swapLevel2, x);
2598      normalize (factors);
2599      delete [] index;
2600      delete [] Aeval2;
2601      return factors;
2602    }
2603    else if (factors.length() > 0)
2604    {
2605      int oneCount= 0;
2606      CFList l;
2607      for (int i= 0; i < check; i++)
2608      {
2609        if (index[i] == 1)
2610        {
2611          iter=biFactors;
2612          for (int j=1; j <= i-oneCount; j++)
2613            iter++;
2614          iter.remove (1);
2615          for (int j= 0; j < A.level() -2; j++)
2616          {
2617            l= leadingCoeffs2[j];
2618            iter= l;
2619            for (int k=1; k <= i-oneCount; k++)
2620              iter++;
2621            iter.remove (1);
2622            leadingCoeffs2[j]=l;
2623          }
2624          oneCount++;
2625        }
2626      }
2627      bufFactors= factors;
2628      factors= CFList();
2629    }
2630    else
2631      factors= CFList();
2632    delete [] index;
2633  }
2634
2635  A= shift2Zero (A, Aeval, evaluation);
2636
2637  for (iter= biFactors; iter.hasItem(); iter++)
2638    iter.getItem()= iter.getItem () (y + evaluation.getLast(), y);
2639
2640  for (int i= 0; i < A.level() - 3; i++)
2641    leadingCoeffs2[i]= CFList();
2642  for (iter= leadingCoeffs2[A.level() - 3]; iter.hasItem(); iter++)
2643  {
2644    iter.getItem()= shift2Zero (iter.getItem(), list, evaluation);
2645    for (int i= A.level() - 4; i > -1; i--)
2646    {
2647      if (i + 1 == A.level() - 3)
2648        leadingCoeffs2[i].append (iter.getItem() (0, i + 4));
2649      else
2650        leadingCoeffs2[i].append (leadingCoeffs2[i+1].getLast() (0, i + 4));
2651    }
2652  }
2653
2654  CFArray Pi;
2655  CFList diophant;
2656  int* liftBounds= new int [A.level() - 1];
2657  int liftBoundsLength= A.level() - 1;
2658  for (int i= 0; i < liftBoundsLength; i++)
2659    liftBounds [i]= degree (A, i + 2) + 1;
2660
2661  Aeval.removeFirst();
2662  bool noOneToOne= false;
2663  factors= nonMonicHenselLift (Aeval, biFactors, leadingCoeffs2, diophant,
2664                               Pi, liftBounds, liftBoundsLength, noOneToOne);
2665
2666  if (!noOneToOne)
2667  {
2668    int check= factors.length();
2669    A= oldA;
2670    factors= recoverFactors (A, factors, evaluation);
2671    if (check != factors.length())
2672      noOneToOne= true;
2673    else
2674      factors= Union (factors, bufFactors);
2675
2676    if (extension && !noOneToOne)
2677      factors= extNonMonicFactorRecombination (factors, A, info);
2678  }
2679  if (noOneToOne)
2680  {
2681    A= shift2Zero (oldA, Aeval, evaluation);
2682    biFactors= oldBiFactors;
2683    for (iter= biFactors; iter.hasItem(); iter++)
2684      iter.getItem()= iter.getItem () (y + evaluation.getLast(), y);
2685    CanonicalForm LCA= LC (Aeval.getFirst(), 1);
2686    CanonicalForm yToLift= power (y, lift);
2687    CFListIterator i= biFactors;
2688    lift= degree (i.getItem(), 2) + degree (LC (i.getItem(), 1)) + 1;
2689    i++;
2690
2691    for (; i.hasItem(); i++)
2692      lift= tmax (lift,
2693                  degree (i.getItem(), 2) + degree (LC (i.getItem(), 1)) + 1);
2694
2695    lift= tmax (degree (Aeval.getFirst() , 2) + 1, lift);
2696
2697    i= biFactors;
2698    yToLift= power (y, lift);
2699    CanonicalForm dummy;
2700    for (; i.hasItem(); i++)
2701    {
2702      LCA= LC (i.getItem(), 1);
2703      extgcd (LCA, yToLift, LCA, dummy);
2704      i.getItem()= mod (i.getItem()*LCA, yToLift);
2705    }
2706
2707    liftBoundsLength= F.level() - 1;
2708    liftBounds= liftingBounds (A, lift);
2709
2710    CFList MOD;
2711    bool earlySuccess;
2712    CFList earlyFactors, liftedFactors;
2713    TIMING_START (fac_fq_hensel_lift);
2714    liftedFactors= henselLiftAndEarly
2715                   (A, MOD, liftBounds, earlySuccess, earlyFactors,
2716                    Aeval, biFactors, evaluation, info);
2717    TIMING_END_AND_PRINT (fac_fq_hensel_lift, "time for hensel lifting: ");
2718
2719    if (!extension)
2720    {
2721      TIMING_START (fac_fq_factor_recombination);
2722      factors= factorRecombination (A, liftedFactors, MOD);
2723      TIMING_END_AND_PRINT (fac_fq_factor_recombination,
2724                            "time for factor recombination: ");
2725    }
2726    else
2727    {
2728      TIMING_START (fac_fq_factor_recombination);
2729      factors= extFactorRecombination (liftedFactors, A, MOD, info, evaluation);
2730      TIMING_END_AND_PRINT (fac_fq_factor_recombination,
2731                            "time for factor recombination: ");
2732    }
2733
2734    if (earlySuccess)
2735      factors= Union (factors, earlyFactors);
2736    if (!extension)
2737    {
2738      for (CFListIterator i= factors; i.hasItem(); i++)
2739      {
2740        int kk= Aeval.getLast().level();
2741        for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
2742        {
2743          if (i.getItem().level() < kk)
2744            continue;
2745          i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
2746        }
2747      }
2748    }
2749  }
2750
2751  if (v.level() != 1)
2752  {
2753    for (CFListIterator iter= factors; iter.hasItem(); iter++)
2754      iter.getItem()= swapvar (iter.getItem(), v, y);
2755  }
2756
2757  swap (factors, swapLevel, swapLevel2, x);
2758  append (factors, contentAFactors);
2759  decompress (factors, N);
2760  normalize (factors);
2761
2762  delete[] liftBounds;
2763
2764  return factors;
2765}
2766
2767/// multivariate factorization over an extension of the initial field
2768CFList
2769extFactorize (const CanonicalForm& F, const ExtensionInfo& info)
2770{
2771  CanonicalForm A= F;
2772
2773  Variable alpha= info.getAlpha();
2774  Variable beta= info.getBeta();
2775  int k= info.getGFDegree();
2776  char cGFName= info.getGFName();
2777  CanonicalForm delta= info.getDelta();
2778  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
2779  Variable w= Variable (1);
2780
2781  CFList factors;
2782  if (!GF && alpha == w)  // we are in F_p
2783  {
2784    CFList factors;
2785    bool extension= true;
2786    int p= getCharacteristic();
2787    if (p*p < (1<<16)) // pass to GF if possible
2788    {
2789      setCharacteristic (getCharacteristic(), 2, 'Z');
2790      ExtensionInfo info= ExtensionInfo (extension);
2791      A= A.mapinto();
2792      factors= multiFactorize (A, info);
2793
2794      CanonicalForm mipo= gf_mipo;
2795      setCharacteristic (getCharacteristic());
2796      Variable vBuf= rootOf (mipo.mapinto());
2797      for (CFListIterator j= factors; j.hasItem(); j++)
2798        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
2799    }
2800    else  // not able to pass to GF, pass to F_p(\alpha)
2801    {
2802      CanonicalForm mipo= randomIrredpoly (2, w);
2803      Variable v= rootOf (mipo);
2804      ExtensionInfo info= ExtensionInfo (v);
2805      factors= multiFactorize (A, info);
2806    }
2807    return factors;
2808  }
2809  else if (!GF && (alpha != w)) // we are in F_p(\alpha)
2810  {
2811    if (k == 1) // need factorization over F_p
2812    {
2813      int extDeg= degree (getMipo (alpha));
2814      extDeg++;
2815      CanonicalForm mipo= randomIrredpoly (extDeg + 1, w);
2816      Variable v= rootOf (mipo);
2817      ExtensionInfo info= ExtensionInfo (v);
2818      factors= multiFactorize (A, info);
2819    }
2820    else
2821    {
2822      if (beta == w)
2823      {
2824        Variable v= chooseExtension (alpha, beta, k);
2825        CanonicalForm primElem, imPrimElem;
2826        bool primFail= false;
2827        Variable vBuf;
2828        primElem= primitiveElement (alpha, vBuf, primFail);
2829        ASSERT (!primFail, "failure in integer factorizer");
2830        if (primFail)
2831          ; //ERROR
2832        else
2833          imPrimElem= mapPrimElem (primElem, vBuf, v);
2834
2835        CFList source, dest;
2836        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
2837                                   source, dest);
2838        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
2839        factors= multiFactorize (bufA, info);
2840      }
2841      else
2842      {
2843        Variable v= chooseExtension (alpha, beta, k);
2844        CanonicalForm primElem, imPrimElem;
2845        bool primFail= false;
2846        Variable vBuf;
2847        ASSERT (!primFail, "failure in integer factorizer");
2848        if (primFail)
2849          ; //ERROR
2850        else
2851          imPrimElem= mapPrimElem (delta, beta, v);
2852
2853        CFList source, dest;
2854        CanonicalForm bufA= mapDown (A, info, source, dest);
2855        source= CFList();
2856        dest= CFList();
2857        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
2858        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
2859        factors= multiFactorize (bufA, info);
2860      }
2861    }
2862    return factors;
2863  }
2864  else // we are in GF (p^k)
2865  {
2866    int p= getCharacteristic();
2867    int extensionDeg= getGFDegree();
2868    bool extension= true;
2869    if (k == 1) // need factorization over F_p
2870    {
2871      extensionDeg++;
2872      if (pow ((double) p, (double) extensionDeg) < (1<<16))
2873      // pass to GF(p^k+1)
2874      {
2875        CanonicalForm mipo= gf_mipo;
2876        setCharacteristic (p);
2877        Variable vBuf= rootOf (mipo.mapinto());
2878        A= GF2FalphaRep (A, vBuf);
2879        setCharacteristic (p, extensionDeg, 'Z');
2880        ExtensionInfo info= ExtensionInfo (extension);
2881        factors= multiFactorize (A.mapinto(), info);
2882      }
2883      else // not able to pass to another GF, pass to F_p(\alpha)
2884      {
2885        CanonicalForm mipo= gf_mipo;
2886        setCharacteristic (p);
2887        Variable vBuf= rootOf (mipo.mapinto());
2888        A= GF2FalphaRep (A, vBuf);
2889        Variable v= chooseExtension (vBuf, beta, k);
2890        ExtensionInfo info= ExtensionInfo (v, extension);
2891        factors= multiFactorize (A, info);
2892      }
2893    }
2894    else // need factorization over GF (p^k)
2895    {
2896      if (pow ((double) p, (double) 2*extensionDeg) < (1<<16))
2897      // pass to GF(p^2k)
2898      {
2899        setCharacteristic (p, 2*extensionDeg, 'Z');
2900        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
2901        factors= multiFactorize (GFMapUp (A, extensionDeg), info);
2902        setCharacteristic (p, extensionDeg, cGFName);
2903      }
2904      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
2905      {
2906        CanonicalForm mipo= gf_mipo;
2907        setCharacteristic (p);
2908        Variable v1= rootOf (mipo.mapinto());
2909        A= GF2FalphaRep (A, v1);
2910        Variable v2= chooseExtension (v1, v1, k);
2911        CanonicalForm primElem, imPrimElem;
2912        bool primFail= false;
2913        Variable vBuf;
2914        primElem= primitiveElement (v1, v1, primFail);
2915        if (primFail)
2916          ; //ERROR
2917        else
2918          imPrimElem= mapPrimElem (primElem, v1, v2);
2919        CFList source, dest;
2920        CanonicalForm bufA= mapUp (A, v1, v2, primElem, imPrimElem,
2921                                     source, dest);
2922        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
2923        factors= multiFactorize (bufA, info);
2924        setCharacteristic (p, k, cGFName);
2925        for (CFListIterator i= factors; i.hasItem(); i++)
2926          i.getItem()= Falpha2GFRep (i.getItem());
2927      }
2928    }
2929    return factors;
2930  }
2931}
2932
2933#endif
2934/* HAVE_NTL */
2935
Note: See TracBrowser for help on using the repository browser.