source: git/factory/facFqBivar.cc @ 2d35fe

jengelh-datetimespielwiese
Last change on this file since 2d35fe was 2d35fe, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: added more timing info to bivariate factorization
  • Property mode set to 100644
File size: 231.6 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqBivar.cc
5 *
6 * This file provides functions for factorizing a bivariate polynomial over
7 * \f$ F_{p} \f$ , \f$ F_{p}(\alpha ) \f$ or GF, based on "Modern Computer
8 * Algebra, Chapter 15" by J. von zur Gathen & J. Gerhard and "Factoring
9 * multivariate polynomials over a finite field" by L. Bernardin.
10 * Factor Recombination is described in "Factoring polynomials over global
11 * fields" by K. Belabas, M. van Hoeij, J. Klueners, A. Steel
12 *
13 *
14 * @author Martin Lee
15 *
16 **/
17/*****************************************************************************/
18
19#include "config.h"
20
21#include "cf_assert.h"
22#include "debug.h"
23#include "timing.h"
24
25#include "canonicalform.h"
26#include "cf_defs.h"
27#include "cf_map_ext.h"
28#include "cf_random.h"
29#include "facHensel.h"
30#include "facMul.h"
31#include "cf_map.h"
32#include "cf_gcd_smallp.h"
33#include "facFqBivarUtil.h"
34#include "facFqBivar.h"
35#include "cfNewtonPolygon.h"
36#include "algext.h"
37
38#ifdef HAVE_NTL
39#include "NTLconvert.h"
40
41#ifdef HAVE_FLINT
42#include "FLINTconvert.h"
43#endif
44
45TIMING_DEFINE_PRINT(fac_fq_uni_factorizer)
46TIMING_DEFINE_PRINT(fac_fq_bi_hensel_lift)
47TIMING_DEFINE_PRINT(fac_fq_bi_factor_recombination)
48TIMING_DEFINE_PRINT(fac_fq_bi_evaluation)
49TIMING_DEFINE_PRINT(fac_fq_bi_shift_to_zero)
50TIMING_DEFINE_PRINT(fac_fq_logarithmic)
51TIMING_DEFINE_PRINT(fac_fq_compute_lattice_lift)
52TIMING_DEFINE_PRINT(fac_fq_till_reduced)
53TIMING_DEFINE_PRINT(fac_fq_reconstruction)
54TIMING_DEFINE_PRINT(fac_fq_lift)
55TIMING_DEFINE_PRINT(fac_fq_uni_total)
56
57CanonicalForm prodMod0 (const CFList& L, const CanonicalForm& M, const modpk& b)
58{
59  if (L.isEmpty())
60    return 1;
61  else if (L.length() == 1)
62    return mod (L.getFirst()(0, 1) , M);
63  else if (L.length() == 2)
64    return mod (mulNTL (L.getFirst()(0, 1),L.getLast()(0, 1), b), M);
65  else
66  {
67    int l= L.length()/2;
68    CFListIterator i= L;
69    CFList tmp1, tmp2;
70    CanonicalForm buf1, buf2;
71    for (int j= 1; j <= l; j++, i++)
72      tmp1.append (i.getItem());
73    tmp2= Difference (L, tmp1);
74    buf1= prodMod0 (tmp1, M, b);
75    buf2= prodMod0 (tmp2, M, b);
76    return mod (mulNTL (buf1,buf2, b), M);
77  }
78}
79
80CanonicalForm evalPoint (const CanonicalForm& F, CanonicalForm & eval,
81                         const Variable& alpha, CFList& list, const bool& GF,
82                         bool& fail)
83{
84  fail= false;
85  Variable x= Variable(2);
86  Variable y= Variable(1);
87  FFRandom genFF;
88  GFRandom genGF;
89  CanonicalForm random, mipo;
90  double bound;
91  int p= getCharacteristic ();
92  if (alpha.level() != 1)
93  {
94    mipo= getMipo (alpha);
95    int d= degree (mipo);
96    bound= ipower (p, d);
97  }
98  else if (GF)
99  {
100    int d= getGFDegree();
101    bound= ipower (p, d);
102  }
103  else
104    bound= p;
105
106  random= 0;
107  do
108  {
109    if (list.length() >= bound)
110    {
111      fail= true;
112      break;
113    }
114    if (list.isEmpty())
115      random= 0;
116    else if (GF)
117    {
118      if (list.length() == 1)
119        random= getGFGenerator();
120      else
121        random= genGF.generate();
122    }
123    else if (list.length() < p || alpha.level() == 1)
124      random= genFF.generate();
125    else if (alpha != x && list.length() >= p)
126    {
127      if (list.length() == p)
128        random= alpha;
129      else
130      {
131        AlgExtRandomF genAlgExt (alpha);
132        random= genAlgExt.generate();
133      }
134    }
135    if (find (list, random)) continue;
136    eval= F (random, x);
137    if (degree (eval) != degree (F, y))
138    { //leading coeff vanishes
139      if (!find (list, random))
140        list.append (random);
141      continue;
142    }
143    if (degree (gcd (deriv (eval, eval.mvar()), eval), eval.mvar()) > 0)
144    { //evaluated polynomial is not squarefree
145      if (!find (list, random))
146        list.append (random);
147      continue;
148    }
149  } while (find (list, random));
150
151  return random;
152}
153
154CFList
155uniFactorizer (const CanonicalForm& A, const Variable& alpha, const bool& GF)
156{
157  Variable x= A.mvar();
158  if (A.inCoeffDomain())
159    return CFList();
160  ASSERT (A.isUnivariate(),
161          "univariate polynomial expected or constant expected");
162  CFFList factorsA;
163  if (fac_NTL_char != getCharacteristic())
164  {
165    fac_NTL_char= getCharacteristic();
166    zz_p::init (getCharacteristic());
167  }
168  if (GF)
169  {
170    int k= getGFDegree();
171    char cGFName= gf_name;
172    CanonicalForm mipo= gf_mipo;
173    setCharacteristic (getCharacteristic());
174    Variable beta= rootOf (mipo.mapinto());
175    CanonicalForm buf= GF2FalphaRep (A, beta);
176    if (getCharacteristic() > 2)
177    {
178      zz_pX NTLMipo= convertFacCF2NTLzzpX (mipo.mapinto());
179      zz_pE::init (NTLMipo);
180      zz_pEX NTLA= convertFacCF2NTLzz_pEX (buf, NTLMipo);
181      MakeMonic (NTLA);
182      vec_pair_zz_pEX_long NTLFactorsA= CanZass (NTLA);
183      zz_pE multi= to_zz_pE (1);
184      factorsA= convertNTLvec_pair_zzpEX_long2FacCFFList (NTLFactorsA, multi,
185                                                         x, beta);
186    }
187    else
188    {
189      GF2X NTLMipo= convertFacCF2NTLGF2X (mipo.mapinto());
190      GF2E::init (NTLMipo);
191      GF2EX NTLA= convertFacCF2NTLGF2EX (buf, NTLMipo);
192      MakeMonic (NTLA);
193      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
194      GF2E multi= to_GF2E (1);
195      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
196                                                           x, beta);
197    }
198    setCharacteristic (getCharacteristic(), k, cGFName);
199    for (CFFListIterator i= factorsA; i.hasItem(); i++)
200    {
201      buf= i.getItem().factor();
202      buf= Falpha2GFRep (buf);
203      i.getItem()= CFFactor (buf, i.getItem().exp());
204    }
205  }
206  else if (alpha.level() != 1)
207  {
208    if (getCharacteristic() > 2)
209    {
210      zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
211      zz_pE::init (NTLMipo);
212      zz_pEX NTLA= convertFacCF2NTLzz_pEX (A, NTLMipo);
213      MakeMonic (NTLA);
214      vec_pair_zz_pEX_long NTLFactorsA= CanZass (NTLA);
215      zz_pE multi= to_zz_pE (1);
216      factorsA= convertNTLvec_pair_zzpEX_long2FacCFFList (NTLFactorsA, multi,
217                                                           x, alpha);
218    }
219    else
220    {
221      GF2X NTLMipo= convertFacCF2NTLGF2X (getMipo (alpha));
222      GF2E::init (NTLMipo);
223      GF2EX NTLA= convertFacCF2NTLGF2EX (A, NTLMipo);
224      MakeMonic (NTLA);
225      vec_pair_GF2EX_long NTLFactorsA= CanZass (NTLA);
226      GF2E multi= to_GF2E (1);
227      factorsA= convertNTLvec_pair_GF2EX_long2FacCFFList (NTLFactorsA, multi,
228                                                           x, alpha);
229    }
230  }
231  else
232  {
233#ifdef HAVE_FLINT
234    nmod_poly_t FLINTA;
235    convertFacCF2nmod_poly_t (FLINTA, A);
236    nmod_poly_factor_t result;
237    nmod_poly_factor_init (result);
238    mp_limb_t leadingCoeff= nmod_poly_factor (result, FLINTA);
239    factorsA= convertFLINTnmod_poly_factor2FacCFFList (result, leadingCoeff, x);
240    if (factorsA.getFirst().factor().inCoeffDomain())
241      factorsA.removeFirst();
242    nmod_poly_factor_clear (result);
243    nmod_poly_clear (FLINTA);
244#else
245    if (getCharacteristic() > 2)
246    {
247      zz_pX NTLA= convertFacCF2NTLzzpX (A);
248      MakeMonic (NTLA);
249      vec_pair_zz_pX_long NTLFactorsA= CanZass (NTLA);
250      zz_p multi= to_zz_p (1);
251      factorsA= convertNTLvec_pair_zzpX_long2FacCFFList (NTLFactorsA, multi,
252                                                          x);
253    }
254    else
255    {
256      GF2X NTLA= convertFacCF2NTLGF2X (A);
257      vec_pair_GF2X_long NTLFactorsA= CanZass (NTLA);
258      GF2 multi= to_GF2 (1);
259      factorsA= convertNTLvec_pair_GF2X_long2FacCFFList (NTLFactorsA, multi,
260                                                          x);
261    }
262#endif
263  }
264  CFList uniFactors;
265  for (CFFListIterator i= factorsA; i.hasItem(); i++)
266    uniFactors.append (i.getItem().factor());
267  return uniFactors;
268}
269
270/// naive factor recombination as decribed in "Factoring
271/// multivariate polynomials over a finite field" by L Bernardin.
272CFList
273extFactorRecombination (CFList& factors, CanonicalForm& F,
274                        const CanonicalForm& N, const ExtensionInfo& info,
275                        DegreePattern& degs, const CanonicalForm& eval, int s,
276                        int thres)
277{
278  if (factors.length() == 0)
279  {
280    F= 1;
281    return CFList();
282  }
283  if (F.inCoeffDomain())
284    return CFList();
285
286  Variable alpha= info.getAlpha();
287  Variable beta= info.getBeta();
288  CanonicalForm gamma= info.getGamma();
289  CanonicalForm delta= info.getDelta();
290  int k= info.getGFDegree();
291
292  CanonicalForm M= N;
293  int l= degree (N);
294  Variable y= F.mvar();
295  Variable x= Variable (1);
296  CFList source, dest;
297  if (degs.getLength() <= 1 || factors.length() == 1)
298  {
299    CFList result= CFList(mapDown (F(y-eval, y), info, source, dest));
300    F= 1;
301    return result;
302  }
303
304  DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, M) == F " <<
305            (mod (LC (F, 1)*prodMod (factors, M), M)/Lc (mod (LC (F, 1)*prodMod (factors, M), M)) == F/Lc (F)));
306  int degMipoBeta= 1;
307  if (!k && beta.level() != 1)
308    degMipoBeta= degree (getMipo (beta));
309
310  CFList T, S, Diff;
311  T= factors;
312
313  CFList result;
314  CanonicalForm buf, buf2, quot;
315
316  buf= F;
317
318  CanonicalForm g, LCBuf= LC (buf, x);
319  int * v= new int [T.length()];
320  for (int i= 0; i < T.length(); i++)
321    v[i]= 0;
322
323  CFArray TT;
324  DegreePattern bufDegs1, bufDegs2;
325  bufDegs1= degs;
326  int subsetDeg;
327  TT= copy (factors);
328  bool nosubset= false;
329  bool recombination= false;
330  bool trueFactor= false;
331  CanonicalForm test;
332  CanonicalForm buf0= buf (0, x)*LCBuf;
333  while (T.length() >= 2*s && s <= thres)
334  {
335    while (nosubset == false)
336    {
337      if (T.length() == s)
338      {
339        delete [] v;
340        if (recombination)
341        {
342          T.insert (LCBuf);
343          g= prodMod (T, M);
344          T.removeFirst();
345          g /= content(g);
346          g= g (y - eval, y);
347          g /= Lc (g);
348          appendTestMapDown (result, g, info, source, dest);
349          F= 1;
350          return result;
351        }
352        else
353        {
354          appendMapDown (result, F (y - eval, y), info, source, dest);
355          F= 1;
356          return result;
357        }
358      }
359      S= subset (v, s, TT, nosubset);
360      if (nosubset) break;
361      subsetDeg= subsetDegree (S);
362      // skip those combinations that are not possible
363      if (!degs.find (subsetDeg))
364        continue;
365      else
366      {
367        test= prodMod0 (S, M);
368        test *= LCBuf;
369        test = mod (test, M);
370        if (fdivides (test, buf0))
371        {
372          S.insert (LCBuf);
373          g= prodMod (S, M);
374          S.removeFirst();
375          g /= content (g, x);
376          if (fdivides (g, buf, quot))
377          {
378            buf2= g (y - eval, y);
379            buf2 /= Lc (buf2);
380
381            if (!k && beta.level() == 1)
382            {
383              if (degree (buf2, alpha) < degMipoBeta)
384              {
385                buf= quot;
386                LCBuf= LC (buf, x);
387                recombination= true;
388                appendTestMapDown (result, buf2, info, source, dest);
389                trueFactor= true;
390              }
391            }
392            else
393            {
394              if (!isInExtension (buf2, gamma, k, delta, source, dest))
395              {
396                buf= quot;
397                LCBuf= LC (buf, x);
398                recombination= true;
399                appendTestMapDown (result, buf2, info, source, dest);
400                trueFactor= true;
401              }
402            }
403            if (trueFactor)
404            {
405              T= Difference (T, S);
406              l -= degree (g);
407              M= power (y, l);
408              buf0= buf (0, x)*LCBuf;
409
410              // compute new possible degree pattern
411              bufDegs2= DegreePattern (T);
412              bufDegs1.intersect (bufDegs2);
413              bufDegs1.refine ();
414              if (T.length() < 2*s || T.length() == s ||
415                  bufDegs1.getLength() == 1)
416              {
417                delete [] v;
418                if (recombination)
419                {
420                  buf= buf (y-eval,y);
421                  buf /= Lc (buf);
422                  appendTestMapDown (result, buf, info, source,
423                                      dest);
424                  F= 1;
425                  return result;
426                }
427                else
428                {
429                  appendMapDown (result, F (y - eval, y), info, source, dest);
430                  F= 1;
431                  return result;
432                }
433              }
434              trueFactor= false;
435              TT= copy (T);
436              indexUpdate (v, s, T.length(), nosubset);
437              if (nosubset) break;
438            }
439          }
440        }
441      }
442    }
443    s++;
444    if (T.length() < 2*s || T.length() == s)
445    {
446      delete [] v;
447      if (recombination)
448      {
449        buf= buf (y-eval,y);
450        buf /= Lc (buf);
451        appendTestMapDown (result, buf, info, source, dest);
452        F= 1;
453        return result;
454      }
455      else
456      {
457        appendMapDown (result, F (y - eval, y), info, source, dest);
458        F= 1;
459        return result;
460      }
461    }
462    for (int i= 0; i < T.length(); i++)
463      v[i]= 0;
464    nosubset= false;
465  }
466  if (T.length() < 2*s)
467  {
468    appendMapDown (result, F (y - eval, y), info, source, dest);
469    F= 1;
470    delete [] v;
471    return result;
472  }
473
474  if (s > thres)
475  {
476    factors= T;
477    F= buf;
478    degs= bufDegs1;
479  }
480
481  delete [] v;
482  return result;
483}
484
485/// naive factor recombination as decribed in "Factoring
486/// multivariate polynomials over a finite field" by L Bernardin.
487CFList
488factorRecombination (CFList& factors, CanonicalForm& F,
489                     const CanonicalForm& N, DegreePattern& degs, const
490                     CanonicalForm& eval, int s, int thres, const modpk& b,
491                     const CanonicalForm& den
492                    )
493{
494  if (factors.length() == 0)
495  {
496    F= 1;
497    return CFList ();
498  }
499  if (F.inCoeffDomain())
500    return CFList();
501  Variable y= Variable (2);
502  if (degs.getLength() <= 1 || factors.length() == 1)
503  {
504    CFList result= CFList (F(y-eval,y));
505    F= 1;
506    return result;
507  }
508#ifdef DEBUGOUTPUT
509  if (b.getp() == 0)
510    DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, N) == F " <<
511              (mod (LC (F, 1)*prodMod (factors, N),N)/Lc (mod (LC (F, 1)*prodMod (factors, N),N)) == F/Lc(F)));
512  else
513    DEBOUTLN (cerr, "LC (F, 1)*prodMod (factors, N) == F " <<
514              (mod (b(LC (F, 1)*prodMod (factors, N)),N)/Lc (mod (b(LC (F, 1)*prodMod (factors, N)),N)) == F/Lc(F)));
515#endif
516
517  CFList T, S;
518
519  CanonicalForm M= N;
520  int l= degree (N);
521  T= factors;
522  CFList result;
523  Variable x= Variable (1);
524  CanonicalForm denom= den, denQuot;
525  CanonicalForm LCBuf= LC (F, x)*denom;
526  CanonicalForm g, quot, buf= F;
527  int * v= new int [T.length()];
528  for (int i= 0; i < T.length(); i++)
529    v[i]= 0;
530  bool nosubset= false;
531  CFArray TT;
532  DegreePattern bufDegs1, bufDegs2;
533  bufDegs1= degs;
534  int subsetDeg;
535  TT= copy (factors);
536  bool recombination= false;
537  CanonicalForm test;
538  bool isRat= (isOn (SW_RATIONAL) && getCharacteristic() == 0) ||
539               getCharacteristic() > 0;
540  if (!isRat)
541    On (SW_RATIONAL);
542  CanonicalForm buf0= mulNTL (buf (0, x), LCBuf);
543  if (!isRat)
544    Off (SW_RATIONAL);
545  while (T.length() >= 2*s && s <= thres)
546  {
547    while (nosubset == false)
548    {
549      if (T.length() == s)
550      {
551        delete [] v;
552        if (recombination)
553        {
554          T.insert (LCBuf);
555          g= prodMod (T, M);
556          if (b.getp() != 0)
557            g= b(g);
558          T.removeFirst();
559          g /= content (g,x);
560          result.append (g(y-eval,y));
561          F= 1;
562          return result;
563        }
564        else
565        {
566          result= CFList (F(y-eval,y));
567          F= 1;
568          return result;
569        }
570      }
571      S= subset (v, s, TT, nosubset);
572      if (nosubset) break;
573      subsetDeg= subsetDegree (S);
574      // skip those combinations that are not possible
575      if (!degs.find (subsetDeg))
576        continue;
577      else
578      {
579        if (!isRat)
580          On (SW_RATIONAL);
581        test= prodMod0 (S, M);
582        if (!isRat)
583        {
584          test *= bCommonDen (test);
585          Off (SW_RATIONAL);
586        }
587        test= mulNTL (test, LCBuf, b);
588        test= mod (test, M);
589        if (uniFdivides (test, buf0))
590        {
591          if (!isRat)
592            On (SW_RATIONAL);
593          S.insert (LCBuf);
594          g= prodMod (S, M);
595          S.removeFirst();
596          if (!isRat)
597          {
598            g *= bCommonDen(g);
599            Off (SW_RATIONAL);
600          }
601          if (b.getp() != 0)
602            g= b(g);
603          if (!isRat)
604            On (SW_RATIONAL);
605          g /= content (g, x);
606          if (!isRat)
607          {
608            On (SW_RATIONAL);
609            if (!Lc (g).inBaseDomain())
610              g /= Lc (g);
611            g *= bCommonDen (g);
612            Off (SW_RATIONAL);
613            g /= icontent (g);
614            On (SW_RATIONAL);
615          }
616          if (fdivides (g, buf, quot))
617          {
618            denom *= abs (lc (g));
619            recombination= true;
620            result.append (g (y-eval,y));
621            if (b.getp() != 0)
622            {
623              denQuot= bCommonDen (quot);
624              buf= quot*denQuot;
625              Off (SW_RATIONAL);
626              denom /= gcd (denom, denQuot);
627              On (SW_RATIONAL);
628            }
629            else
630              buf= quot;
631            LCBuf= LC (buf, x)*denom;
632            T= Difference (T, S);
633            l -= degree (g);
634            M= power (y, l);
635            buf0= mulNTL (buf (0, x), LCBuf);
636            if (!isRat)
637              Off (SW_RATIONAL);
638            // compute new possible degree pattern
639            bufDegs2= DegreePattern (T);
640            bufDegs1.intersect (bufDegs2);
641            bufDegs1.refine ();
642            if (T.length() < 2*s || T.length() == s ||
643                bufDegs1.getLength() == 1)
644            {
645              delete [] v;
646              if (recombination)
647              {
648                result.append (buf (y-eval,y));
649                F= 1;
650                return result;
651              }
652              else
653              {
654                result= CFList (F (y-eval,y));
655                F= 1;
656                return result;
657              }
658            }
659            TT= copy (T);
660            indexUpdate (v, s, T.length(), nosubset);
661            if (nosubset) break;
662          }
663          if (!isRat)
664            Off (SW_RATIONAL);
665        }
666      }
667    }
668    s++;
669    if (T.length() < 2*s || T.length() == s)
670    {
671      delete [] v;
672      if (recombination)
673      {
674        result.append (buf(y-eval,y));
675        F= 1;
676        return result;
677      }
678      else
679      {
680        result= CFList (F(y-eval,y));
681        F= 1;
682        return result;
683      }
684    }
685    for (int i= 0; i < T.length(); i++)
686      v[i]= 0;
687    nosubset= false;
688  }
689  delete [] v;
690  if (T.length() < 2*s)
691  {
692    result.append (F(y-eval,y));
693    F= 1;
694    return result;
695  }
696
697  if (s > thres)
698  {
699    factors= T;
700    F= buf;
701    degs= bufDegs1;
702  }
703
704  return result;
705}
706
707Variable chooseExtension (const Variable & alpha, const Variable& beta, int k)
708{
709  if (fac_NTL_char != getCharacteristic())
710  {
711    fac_NTL_char= getCharacteristic();
712    zz_p::init (getCharacteristic());
713  }
714  zz_pX NTLIrredpoly;
715  int i=1, m= 2;
716  // extension of F_p needed
717  if (alpha.level() == 1 && beta.level() == 1 && k == 1)
718  {
719    i= 1;
720    m= 2;
721  } //extension of F_p(alpha) needed but want to factorize over F_p
722  else if (alpha.level() != 1 && beta.level() == 1 && k == 1)
723  {
724    i= 1;
725    m= degree (getMipo (alpha)) + 1;
726  } //extension of F_p(alpha) needed for first time
727  else if (alpha.level() != 1 && beta.level() == 1 && k != 1)
728  {
729    i= 2;
730    m= degree (getMipo (alpha));
731  }
732  else if (alpha.level() != 1 && beta.level() != 1 && k != 1)
733  {
734    m= degree (getMipo (beta));
735    i= degree (getMipo (alpha))/m + 1;
736  }
737  BuildIrred (NTLIrredpoly, i*m);
738  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
739  return rootOf (newMipo);
740}
741
742void
743earlyFactorDetection (CFList& reconstructedFactors, CanonicalForm& F, CFList&
744                      factors, int& adaptedLiftBound, int*& factorsFoundIndex,
745                      DegreePattern& degs, bool& success, int deg, const 
746                      CanonicalForm& eval, const modpk& b, CanonicalForm& den)
747{
748  DegreePattern bufDegs1= degs;
749  DegreePattern bufDegs2;
750  CFList T= factors;
751  CanonicalForm buf= F;
752  Variable x= Variable (1);
753  Variable y= Variable (2);
754  CanonicalForm g, quot;
755  CanonicalForm M= power (F.mvar(), deg);
756  adaptedLiftBound= 0;
757  int d= degree (F), l= 0;
758  bool isRat= (isOn (SW_RATIONAL) && getCharacteristic() == 0) || getCharacteristic() > 0;
759  if (!isRat)
760    On (SW_RATIONAL);
761  if (b.getp() != 0)
762    buf *= bCommonDen (buf);
763  CanonicalForm LCBuf= LC (buf, x)*den;
764  CanonicalForm buf0= mulNTL (buf (0,x), LCBuf);
765  CanonicalForm buf1= mulNTL (buf (1,x), LCBuf);
766  if (!isRat)
767    Off (SW_RATIONAL);
768  CanonicalForm test0, test1;
769  CanonicalForm denQuot;
770
771  for (CFListIterator i= factors; i.hasItem(); i++, l++)
772  {
773    if (!bufDegs1.find (degree (i.getItem(), 1)) || factorsFoundIndex[l] == 1)
774      continue;
775    else
776    {
777      test1= mod (mulNTL (i.getItem() (1,x), LCBuf, b), M);
778      if (uniFdivides (test1, buf1))
779      {
780        test0= mod (mulNTL (i.getItem() (0,x), LCBuf, b), M);
781        if (uniFdivides (test0, buf0))
782        {
783          if (!isRat)
784            On (SW_RATIONAL);
785          g= mulMod2 (i.getItem(), LCBuf, M);
786          if (!isRat)
787          {
788            g *= bCommonDen(g);
789            Off (SW_RATIONAL);
790          }
791          if (b.getp() != 0)
792            g= b(g);
793          if (!isRat)
794            On (SW_RATIONAL);
795          g /= content (g, x);
796          if (!isRat)
797          {
798            On (SW_RATIONAL);
799            if (!Lc (g).inBaseDomain())
800              g /= Lc (g);
801            g *= bCommonDen (g);
802            Off (SW_RATIONAL);
803            g /= icontent (g);
804            On (SW_RATIONAL);
805          }
806          if (fdivides (g, buf, quot))
807          {
808            den *= abs (lc (g));
809            reconstructedFactors.append (g (y-eval,y));
810            factorsFoundIndex[l]= 1;
811            if (b.getp() != 0)
812            {
813              denQuot= bCommonDen (quot);
814              buf= quot*denQuot;
815              Off (SW_RATIONAL);
816              den /= gcd (den, denQuot);
817              On (SW_RATIONAL);
818            }
819            else
820              buf= quot;
821            d -= degree (g);
822            LCBuf= LC (buf, x)*den;
823            buf0= mulNTL (buf (0,x), LCBuf);
824            buf1= mulNTL (buf (1,x), LCBuf);
825            if (!isRat)
826              Off (SW_RATIONAL);
827            T= Difference (T, CFList (i.getItem()));
828            F= buf;
829
830            // compute new possible degree pattern
831            bufDegs2= DegreePattern (T);
832            bufDegs1.intersect (bufDegs2);
833            bufDegs1.refine ();
834            if (bufDegs1.getLength() <= 1)
835            {
836              if (!buf.inCoeffDomain())
837              {
838                reconstructedFactors.append (buf (y-eval,y));
839                F= 1;
840              }
841              break;
842            }
843          }
844          if (!isRat)
845            Off (SW_RATIONAL);
846        }
847      }
848    }
849  }
850  adaptedLiftBound= d + 1;
851  if (adaptedLiftBound < deg)
852  {
853    degs= bufDegs1;
854    success= true;
855  }
856  if (bufDegs1.getLength() <= 1)
857    degs= bufDegs1;
858}
859
860void
861earlyFactorDetection (CFList& reconstructedFactors, CanonicalForm& F, CFList&
862                      factors, int& adaptedLiftBound, int*& factorsFoundIndex,
863                      DegreePattern& degs, bool& success, int deg, const
864                      CanonicalForm& eval, const modpk& b)
865{
866  CanonicalForm den= 1;
867  earlyFactorDetection (reconstructedFactors, F, factors, adaptedLiftBound,
868                        factorsFoundIndex, degs, success, deg, eval, b, den);
869}
870
871void
872extEarlyFactorDetection (CFList& reconstructedFactors, CanonicalForm& F, CFList&
873                         factors,int& adaptedLiftBound, int*& factorsFoundIndex,
874                         DegreePattern& degs, bool& success, const
875                         ExtensionInfo& info, const CanonicalForm& eval, int deg
876                        )
877{
878  Variable alpha= info.getAlpha();
879  Variable beta= info.getBeta();
880  CanonicalForm gamma= info.getGamma();
881  CanonicalForm delta= info.getDelta();
882  int k= info.getGFDegree();
883  DegreePattern bufDegs1= degs, bufDegs2;
884  CFList result;
885  CFList T= factors;
886  Variable y= F.mvar();
887  Variable x= Variable (1);
888  CanonicalForm buf= F, LCBuf= LC (buf, x), g, buf2;
889  CanonicalForm M= power (y, deg);
890  adaptedLiftBound= 0;
891  bool trueFactor= false;
892  int d= degree (F), l= 0;
893  CFList source, dest;
894  int degMipoBeta= 1;
895  if (!k && beta.level() != 1)
896    degMipoBeta= degree (getMipo (beta));
897  CanonicalForm quot;
898  for (CFListIterator i= factors; i.hasItem(); i++, l++)
899  {
900    if (!bufDegs1.find (degree (i.getItem(), 1)) || factorsFoundIndex[l] == 1)
901      continue;
902    else
903    {
904      g= mulMod2 (i.getItem(), LCBuf, M);
905      g /= content (g, x);
906      if (fdivides (g, buf, quot))
907      {
908        buf2= g (y - eval, y);
909        buf2 /= Lc (buf2);
910
911        if (!k && beta == x)
912        {
913          if (degree (buf2, alpha) < degMipoBeta)
914          {
915            appendTestMapDown (reconstructedFactors, buf2, info, source, dest);
916            factorsFoundIndex[l]= 1;
917            buf= quot;
918            d -= degree (g);
919            LCBuf= LC (buf, x);
920            trueFactor= true;
921          }
922        }
923        else
924        {
925          if (!isInExtension (buf2, gamma, k, delta, source, dest))
926          {
927            appendTestMapDown (reconstructedFactors, buf2, info, source, dest);
928            factorsFoundIndex[l]= 1;
929            buf= quot;
930            d -= degree (g);
931            LCBuf= LC (buf, x);
932            trueFactor= true;
933          }
934        }
935        if (trueFactor)
936        {
937          T= Difference (T, CFList (i.getItem()));
938          F= buf;
939
940          // compute new possible degree pattern
941          bufDegs2= DegreePattern (T);
942          bufDegs1.intersect (bufDegs2);
943          bufDegs1.refine ();
944          trueFactor= false;
945          if (bufDegs1.getLength() <= 1)
946          {
947            if (!buf.inCoeffDomain())
948            {
949              buf= buf (y - eval, y);
950              buf /= Lc (buf);
951              appendMapDown (reconstructedFactors, buf, info, source, dest);
952              F= 1;
953            }
954            break;
955          }
956        }
957      }
958    }
959  }
960  adaptedLiftBound= d + 1;
961  if (adaptedLiftBound < deg)
962  {
963    degs= bufDegs1;
964    success= true;
965  }
966  if (bufDegs1.getLength() <= 1)
967    degs= bufDegs1;
968}
969
970int*
971getCombinations (int * rightSide, int sizeOfRightSide, int& sizeOfOutput,
972                 int degreeLC)
973{
974  Variable x= Variable (1);
975  int p= getCharacteristic();
976  int d= getGFDegree();
977  char cGFName= gf_name;
978  setCharacteristic(0);
979  CanonicalForm buf= 1;
980  for (int i= 0; i < sizeOfRightSide; i++)
981    buf *= (power (x, rightSide [i]) + 1);
982
983  int j= 0;
984  for (CFIterator i= buf; i.hasTerms(); i++, j++)
985  {
986    if (i.exp() < degreeLC)
987    {
988      j++;
989      break;
990    }
991  }
992
993  ASSERT ( j > 1, "j > 1 expected" );
994
995  int* result = new int  [j - 1];
996  sizeOfOutput= j - 1;
997
998  int i= 0;
999  for (CFIterator m = buf; i < j - 1; i++, m++)
1000    result [i]= m.exp();
1001
1002  if (d > 1)
1003    setCharacteristic (p, d, cGFName);
1004  else
1005    setCharacteristic (p);
1006  return result;
1007}
1008
1009int *
1010getLiftPrecisions (const CanonicalForm& F, int& sizeOfOutput, int degreeLC)
1011{
1012  int sizeOfNewtonPoly;
1013  int ** newtonPolyg= newtonPolygon (F, sizeOfNewtonPoly);
1014  int sizeOfRightSide;
1015  int * rightSide= getRightSide(newtonPolyg, sizeOfNewtonPoly, sizeOfRightSide);
1016  int * result= getCombinations(rightSide, sizeOfRightSide, sizeOfOutput,
1017                                degreeLC);
1018  delete [] rightSide;
1019  for (int i= 0; i < sizeOfNewtonPoly; i++)
1020    delete [] newtonPolyg[i];
1021  delete [] newtonPolyg;
1022  return result;
1023}
1024
1025void
1026deleteFactors (CFList& factors, int* factorsFoundIndex)
1027{
1028  CFList result;
1029  int i= 0;
1030  for (CFListIterator iter= factors; iter.hasItem(); iter++, i++)
1031  {
1032    if (factorsFoundIndex[i] == 1)
1033      continue;
1034    else
1035      result.append (iter.getItem());
1036  }
1037  factors= result;
1038}
1039
1040CFList
1041henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList&
1042                    earlyFactors, DegreePattern& degs, int& liftBound,
1043                    const CFList& uniFactors, const ExtensionInfo& info,
1044                    const CanonicalForm& eval,modpk& b, CanonicalForm& den)
1045{
1046  Variable alpha= info.getAlpha();
1047  Variable beta= info.getBeta();
1048  CanonicalForm gamma= info.getGamma();
1049  CanonicalForm delta= info.getDelta();
1050  bool extension= info.isInExtension();
1051
1052  int sizeOfLiftPre;
1053  int * liftPre= getLiftPrecisions (A, sizeOfLiftPre, degree (LC (A, 1), 2));
1054
1055  Variable x= Variable (1);
1056  Variable y= Variable (2);
1057  CFArray Pi;
1058  CFList diophant;
1059  CFList bufUniFactors= uniFactors;
1060  On (SW_RATIONAL);
1061  CanonicalForm bufA= A;
1062  if (!Lc (A).inBaseDomain())
1063  {
1064    bufA /= Lc (A);
1065    CanonicalForm denBufA= bCommonDen (bufA);
1066    bufA *= denBufA;
1067    Off (SW_RATIONAL);
1068    den /= gcd (den, denBufA);
1069  }
1070  else
1071  {
1072    bufA= A;
1073    Off (SW_RATIONAL);
1074    den /= gcd (den, Lc (A));
1075  }
1076  CanonicalForm lcA0= 0;
1077  bool mipoHasDen= false;
1078  if (getCharacteristic() == 0 && b.getp() != 0)
1079  {
1080    if (alpha.level() == 1)
1081    {
1082      lcA0= lc (A (0, 2));
1083      A *= b.inverse (lcA0);
1084      A= b (A);
1085      for (CFListIterator i= bufUniFactors; i.hasItem(); i++)
1086        i.getItem()= b (i.getItem()*b.inverse (lc (i.getItem())));
1087    }
1088    else
1089    {
1090      lcA0= Lc (A (0,2));
1091      On (SW_RATIONAL);
1092      mipoHasDen= !bCommonDen(getMipo(alpha)).isOne();
1093      Off (SW_RATIONAL);
1094      CanonicalForm lcA0inverse= b.inverse (lcA0);
1095      A *= lcA0inverse;
1096      A= b (A);
1097      // Lc of bufUniFactors is in Z
1098      for (CFListIterator i= bufUniFactors; i.hasItem(); i++)
1099        i.getItem()= b (i.getItem()*b.inverse (lc (i.getItem())));
1100    }
1101  }
1102  bufUniFactors.insert (LC (A, x));
1103  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
1104  earlySuccess= false;
1105  int newLiftBound= 0;
1106
1107  int smallFactorDeg= tmin (11, liftPre [sizeOfLiftPre- 1] + 1);//this is a tunable parameter
1108  int dummy;
1109  int * factorsFoundIndex= new int [uniFactors.length()];
1110  for (int i= 0; i < uniFactors.length(); i++)
1111    factorsFoundIndex [i]= 0;
1112
1113  CFList bufBufUniFactors;
1114  Variable v= alpha;
1115  if (smallFactorDeg >= liftBound || degree (A,y) <= 4)
1116    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M, b, true);
1117  else if (sizeOfLiftPre > 1 && sizeOfLiftPre < 30)
1118  {
1119    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M, b, true);
1120    if (mipoHasDen)
1121    {
1122      for (CFListIterator iter= bufUniFactors; iter.hasItem(); iter++)
1123        if (hasFirstAlgVar (iter.getItem(), v))
1124          break;
1125      if (v != alpha)
1126      {
1127        bufBufUniFactors= bufUniFactors;
1128        for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1129          iter.getItem()= replacevar (iter.getItem(), v, alpha);
1130        A= replacevar (A, alpha, v);
1131      }
1132    }
1133
1134    if (!extension)
1135    {
1136      if (v==alpha)
1137        earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1138                              factorsFoundIndex, degs, earlySuccess,
1139                              smallFactorDeg, eval, b, den);
1140      else
1141        earlyFactorDetection(earlyFactors, bufA, bufBufUniFactors, newLiftBound,
1142                             factorsFoundIndex, degs, earlySuccess,
1143                             smallFactorDeg, eval, b, den);
1144    }
1145    else
1146      extEarlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1147                               factorsFoundIndex, degs, earlySuccess, info,
1148                               eval, smallFactorDeg);
1149    if (degs.getLength() > 1 && !earlySuccess &&
1150        smallFactorDeg != liftPre [sizeOfLiftPre-1] + 1)
1151    {
1152      if (newLiftBound >= liftPre[sizeOfLiftPre-1]+1)
1153      {
1154        bufUniFactors.insert (LC (A, x));
1155        henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
1156                            liftPre[sizeOfLiftPre-1] + 1, Pi, diophant, M, b);
1157        if (v!=alpha)
1158        {
1159          bufBufUniFactors= bufUniFactors;
1160          for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1161            iter.getItem()= replacevar (iter.getItem(), v, alpha);
1162        }
1163        if (!extension)
1164        {
1165          if (v==alpha)
1166          earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1167                                factorsFoundIndex, degs, earlySuccess,
1168                                liftPre[sizeOfLiftPre-1] + 1, eval, b, den);
1169          else
1170          earlyFactorDetection (earlyFactors,bufA,bufBufUniFactors,newLiftBound,
1171                                factorsFoundIndex, degs, earlySuccess,
1172                                liftPre[sizeOfLiftPre-1] + 1, eval, b, den);
1173        }
1174        else
1175          extEarlyFactorDetection (earlyFactors,bufA,bufUniFactors,newLiftBound,
1176                                   factorsFoundIndex, degs, earlySuccess, info,
1177                                   eval, liftPre[sizeOfLiftPre-1] + 1);
1178      }
1179    }
1180    else if (earlySuccess)
1181      liftBound= newLiftBound;
1182
1183    int i= sizeOfLiftPre - 1;
1184    while (degs.getLength() > 1 && !earlySuccess && i - 1 >= 0)
1185    {
1186      if (newLiftBound >= liftPre[i] + 1)
1187      {
1188        bufUniFactors.insert (LC (A, x));
1189        henselLiftResume12 (A, bufUniFactors, liftPre[i] + 1,
1190                            liftPre[i-1] + 1, Pi, diophant, M, b);
1191        if (v!=alpha)
1192        {
1193          bufBufUniFactors= bufUniFactors;
1194          for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1195            iter.getItem()= replacevar (iter.getItem(), v, alpha);
1196        }
1197        if (!extension)
1198        {
1199          if (v==alpha)
1200          earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1201                                factorsFoundIndex, degs, earlySuccess,
1202                                liftPre[i-1] + 1, eval, b, den);
1203          else
1204          earlyFactorDetection (earlyFactors,bufA,bufBufUniFactors,newLiftBound,
1205                                factorsFoundIndex, degs, earlySuccess,
1206                                liftPre[i-1] + 1, eval, b, den);
1207        }
1208        else
1209          extEarlyFactorDetection (earlyFactors,bufA,bufUniFactors,newLiftBound,
1210                                   factorsFoundIndex, degs, earlySuccess, info,
1211                                   eval, liftPre[i-1] + 1);
1212      }
1213      else
1214      {
1215        liftBound= newLiftBound;
1216        break;
1217      }
1218      i--;
1219    }
1220    if (earlySuccess)
1221      liftBound= newLiftBound;
1222    //after here all factors are lifted to liftPre[sizeOfLiftPre-1]
1223  }
1224  else
1225  {
1226    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M, b, true);
1227    if (mipoHasDen)
1228    {
1229      for (CFListIterator iter= bufUniFactors; iter.hasItem(); iter++)
1230        if (hasFirstAlgVar (iter.getItem(), v))
1231          break;
1232      if (v != alpha)
1233      {
1234        bufBufUniFactors= bufUniFactors;
1235        for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1236          iter.getItem()= replacevar (iter.getItem(), v, alpha);
1237        A= replacevar (A, alpha, v);
1238      }
1239    }
1240    if (!extension)
1241    {
1242      if (v==alpha)
1243      earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1244                            factorsFoundIndex, degs, earlySuccess,
1245                            smallFactorDeg, eval, b, den);
1246      else
1247      earlyFactorDetection (earlyFactors, bufA, bufBufUniFactors, newLiftBound,
1248                            factorsFoundIndex, degs, earlySuccess,
1249                            smallFactorDeg, eval, b, den);
1250    }
1251    else
1252      extEarlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1253                               factorsFoundIndex, degs, earlySuccess, info,
1254                               eval, smallFactorDeg);
1255    int i= 1;
1256    while ((degree (A,y)/4)*i + 4 <= smallFactorDeg)
1257      i++;
1258    dummy= tmin (degree (A,y)+1, (degree (A,y)/4)*i+4);
1259    if (degs.getLength() > 1 && !earlySuccess && dummy > smallFactorDeg)
1260    {
1261      bufUniFactors.insert (LC (A, x));
1262      henselLiftResume12 (A, bufUniFactors, smallFactorDeg,
1263                          dummy, Pi, diophant, M, b);
1264      if (v!=alpha)
1265      {
1266        bufBufUniFactors= bufUniFactors;
1267        for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1268          iter.getItem()= replacevar (iter.getItem(), v, alpha);
1269      }
1270      if (!extension)
1271      {
1272        if (v==alpha)
1273        earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1274                              factorsFoundIndex, degs, earlySuccess, dummy,eval,
1275                              b, den);
1276        else
1277        earlyFactorDetection (earlyFactors, bufA,bufBufUniFactors, newLiftBound,
1278                              factorsFoundIndex, degs, earlySuccess, dummy,eval,
1279                              b, den);
1280      }
1281      else
1282        extEarlyFactorDetection (earlyFactors, bufA,bufUniFactors, newLiftBound,
1283                                 factorsFoundIndex, degs, earlySuccess, info,
1284                                 eval, dummy);
1285    }
1286    while (degs.getLength() > 1 && !earlySuccess && i < 4)
1287    {
1288      if (newLiftBound >= dummy)
1289      {
1290        bufUniFactors.insert (LC (A, x));
1291        dummy= tmin (degree (A,y)+1, (degree (A,y)/4)*(i+1)+4);
1292        henselLiftResume12 (A, bufUniFactors, (degree (A,y)/4)*i + 4,
1293                            dummy, Pi, diophant, M, b);
1294        if (v!=alpha)
1295        {
1296          bufBufUniFactors= bufUniFactors;
1297          for (CFListIterator iter= bufBufUniFactors; iter.hasItem(); iter++)
1298            iter.getItem()= replacevar (iter.getItem(), v, alpha);
1299        }
1300        if (!extension)
1301        {
1302          if (v==alpha)
1303          earlyFactorDetection (earlyFactors, bufA, bufUniFactors, newLiftBound,
1304                                factorsFoundIndex, degs, earlySuccess, dummy,
1305                                eval, b, den);
1306          else
1307          earlyFactorDetection (earlyFactors,bufA,bufBufUniFactors,newLiftBound,
1308                                factorsFoundIndex, degs, earlySuccess, dummy,
1309                                eval, b, den);
1310        }
1311        else
1312          extEarlyFactorDetection (earlyFactors,bufA,bufUniFactors,newLiftBound,
1313                                   factorsFoundIndex, degs, earlySuccess, info,
1314                                   eval, dummy);
1315      }
1316      else
1317      {
1318        liftBound= newLiftBound;
1319        break;
1320      }
1321      i++;
1322    }
1323    if (earlySuccess)
1324      liftBound= newLiftBound;
1325  }
1326
1327  A= bufA;
1328  if (earlyFactors.length() > 0 && degs.getLength() > 1)
1329  {
1330    liftBound= degree (A,y) + 1;
1331    earlySuccess= true;
1332    deleteFactors (bufUniFactors, factorsFoundIndex);
1333  }
1334
1335  delete [] factorsFoundIndex;
1336  delete [] liftPre;
1337
1338  return bufUniFactors;
1339}
1340
1341CFList
1342henselLiftAndEarly (CanonicalForm& A, bool& earlySuccess, CFList&
1343                    earlyFactors, DegreePattern& degs, int& liftBound,
1344                    const CFList& uniFactors, const ExtensionInfo& info,
1345                    const CanonicalForm& eval)
1346{
1347  modpk dummy= modpk();
1348  CanonicalForm den= 1;
1349  return henselLiftAndEarly (A, earlySuccess, earlyFactors, degs, liftBound,
1350                             uniFactors, info, eval, dummy, den);
1351}
1352
1353long isReduced (const mat_zz_p& M)
1354{
1355  long i, j, nonZero;
1356  for (i = 1; i <= M.NumRows(); i++)
1357  {
1358    nonZero= 0;
1359    for (j = 1; j <= M.NumCols(); j++)
1360    {
1361      if (!IsZero (M (i,j)))
1362        nonZero++;
1363    }
1364    if (nonZero != 1)
1365      return 0;
1366  }
1367  return 1;
1368}
1369
1370#ifdef HAVE_FLINT
1371long isReduced (const nmod_mat_t M)
1372{
1373  long i, j, nonZero;
1374  for (i = 1; i <= nmod_mat_nrows(M); i++)
1375  {
1376    nonZero= 0;
1377    for (j = 1; j <= nmod_mat_ncols (M); j++)
1378    {
1379      if (!(nmod_mat_entry (M, i-1, j-1)==0))
1380        nonZero++;
1381    }
1382    if (nonZero != 1)
1383      return 0;
1384  }
1385  return 1;
1386}
1387#endif
1388
1389long isReduced (const mat_zz_pE& M)
1390{
1391  long i, j, nonZero;
1392  for (i = 1; i <= M.NumRows(); i++)
1393  {
1394    nonZero= 0;
1395    for (j = 1; j <= M.NumCols(); j++)
1396    {
1397      if (!IsZero (M (i,j)))
1398        nonZero++;
1399    }
1400    if (nonZero != 1)
1401      return 0;
1402  }
1403  return 1;
1404}
1405
1406int * extractZeroOneVecs (const mat_zz_p& M)
1407{
1408  long i, j;
1409  bool nonZeroOne= false;
1410  int * result= new int [M.NumCols()];
1411  for (i = 1; i <= M.NumCols(); i++)
1412  {
1413    for (j = 1; j <= M.NumRows(); j++)
1414    {
1415      if (!(IsOne (M (j,i)) || IsZero (M (j,i))))
1416      {
1417        nonZeroOne= true;
1418        break;
1419      }
1420    }
1421    if (!nonZeroOne)
1422      result [i - 1]= 1;
1423    else
1424      result [i - 1]= 0;
1425    nonZeroOne= false;
1426  }
1427  return result;
1428}
1429
1430#ifdef HAVE_FLINT
1431int * extractZeroOneVecs (const nmod_mat_t M)
1432{
1433  long i, j;
1434  bool nonZeroOne= false;
1435  int * result= new int [nmod_mat_ncols (M)];
1436  for (i = 0; i < nmod_mat_ncols (M); i++)
1437  {
1438    for (j = 0; j < nmod_mat_nrows (M); j++)
1439    {
1440      if (!((nmod_mat_entry (M, j, i) == 1) || (nmod_mat_entry (M, j,i) == 0)))
1441      {
1442        nonZeroOne= true;
1443        break;
1444      }
1445    }
1446    if (!nonZeroOne)
1447      result [i]= 1;
1448    else
1449      result [i]= 0;
1450    nonZeroOne= false;
1451  }
1452  return result;
1453}
1454#endif
1455
1456int * extractZeroOneVecs (const mat_zz_pE& M)
1457{
1458  long i, j;
1459  bool nonZeroOne= false;
1460  int * result= new int [M.NumCols()];
1461  for (i = 1; i <= M.NumCols(); i++)
1462  {
1463    for (j = 1; j <= M.NumRows(); j++)
1464    {
1465      if (!(IsOne (M (j,i)) || IsZero (M (j,i))))
1466      {
1467        nonZeroOne= true;
1468        break;
1469      }
1470    }
1471    if (!nonZeroOne)
1472      result [i - 1]= 1;
1473    else
1474      result [i - 1]= 0;
1475    nonZeroOne= false;
1476  }
1477  return result;
1478}
1479
1480void
1481reconstructionTry (CFList& reconstructedFactors, CanonicalForm& F, const CFList&
1482                   factors, const int liftBound, int& factorsFound, int*&
1483                   factorsFoundIndex, mat_zz_pE& N, const CanonicalForm& eval,
1484                   bool beenInThres
1485                  )
1486{
1487  Variable y= Variable (2);
1488  Variable x= Variable (1);
1489  CanonicalForm yToL= power (y, liftBound);
1490  CanonicalForm bufF= F (y-eval, y);
1491  if (factors.length() == 2)
1492  {
1493    CanonicalForm tmp1, tmp2, tmp3;
1494    tmp1= factors.getFirst();
1495    tmp2= factors.getLast();
1496    tmp1= mulMod2 (tmp1, LC (F,x), yToL);
1497    tmp1 /= content (tmp1, x);
1498    tmp1= tmp1 (y-eval, y);
1499    tmp2= mulMod2 (tmp2, LC (F,x), yToL);
1500    tmp2 /= content (tmp2, x);
1501    tmp2= tmp2 (y-eval, y);
1502    tmp3 = tmp1*tmp2;
1503    if (tmp3/Lc (tmp3) == bufF/Lc (bufF))
1504    {
1505      factorsFound++;
1506      F= 1;
1507      reconstructedFactors.append (tmp1);
1508      reconstructedFactors.append (tmp2);
1509      return;
1510    }
1511  }
1512  CanonicalForm quot, buf;
1513  CFListIterator iter;
1514  for (long i= 1; i <= N.NumCols(); i++)
1515  {
1516    if (factorsFoundIndex [i - 1] == 1)
1517      continue;
1518    iter= factors;
1519    if (beenInThres)
1520    {
1521      int count= 1;
1522      while (count < i)
1523      {
1524        count++;
1525        iter++;
1526      }
1527      buf= iter.getItem();
1528    }
1529    else
1530    {
1531      buf= 1;
1532      for (long j= 1; j <= N.NumRows(); j++, iter++)
1533      {
1534        if (!IsZero (N (j,i)))
1535          buf= mulMod2 (buf, iter.getItem(), yToL);
1536      }
1537    }
1538    buf= mulMod2 (buf, LC (F,x), yToL);
1539    buf /= content (buf, x);
1540    buf= buf (y-eval,y);
1541    if (fdivides (buf, bufF, quot))
1542    {
1543      factorsFoundIndex[i - 1]= 1;
1544      factorsFound++;
1545      bufF= quot;
1546      bufF /= Lc (bufF);
1547      reconstructedFactors.append (buf);
1548    }
1549    if (degree (bufF) <= 0)
1550      return;
1551    if (factorsFound + 1 == N.NumCols())
1552    {
1553      reconstructedFactors.append (bufF);
1554      F= 1;
1555      return;
1556    }
1557  }
1558  if (reconstructedFactors.length() != 0)
1559    F= bufF (y+eval,y);
1560}
1561
1562void
1563reconstructionTry (CFList& reconstructedFactors, CanonicalForm& F, const CFList&
1564                   factors, const int liftBound, int& factorsFound, int*&
1565                   factorsFoundIndex, mat_zz_p& N, const CanonicalForm& eval,
1566                   bool beenInThres
1567                  )
1568{
1569  Variable y= Variable (2);
1570  Variable x= Variable (1);
1571  CanonicalForm yToL= power (y, liftBound);
1572  CanonicalForm bufF= F (y-eval, y);
1573  if (factors.length() == 2)
1574  {
1575    CanonicalForm tmp1, tmp2, tmp3;
1576    tmp1= factors.getFirst();
1577    tmp2= factors.getLast();
1578    tmp1= mulMod2 (tmp1, LC (F,x), yToL);
1579    tmp1 /= content (tmp1, x);
1580    tmp1= tmp1 (y-eval, y);
1581    tmp2= mulMod2 (tmp2, LC (F,x), yToL);
1582    tmp2 /= content (tmp2, x);
1583    tmp2= tmp2 (y-eval,y);
1584    tmp3 = tmp1*tmp2;
1585    if (tmp3/Lc (tmp3) == bufF/Lc (bufF))
1586    {
1587      factorsFound++;
1588      F= 1;
1589      reconstructedFactors.append (tmp1);
1590      reconstructedFactors.append (tmp2);
1591      return;
1592    }
1593  }
1594  CanonicalForm quot, buf;
1595  CFListIterator iter;
1596  for (long i= 1; i <= N.NumCols(); i++)
1597  {
1598    if (factorsFoundIndex [i - 1] == 1)
1599      continue;
1600    iter= factors;
1601    if (beenInThres)
1602    {
1603      int count= 1;
1604      while (count < i)
1605      {
1606        count++;
1607        iter++;
1608      }
1609      buf= iter.getItem();
1610    }
1611    else
1612    {
1613      buf= 1;
1614      for (long j= 1; j <= N.NumRows(); j++, iter++)
1615      {
1616        if (!IsZero (N (j,i)))
1617          buf= mulMod2 (buf, iter.getItem(), yToL);
1618      }
1619    }
1620    buf= mulMod2 (buf, LC (F,x), yToL);
1621    buf /= content (buf, x);
1622    buf= buf (y-eval,y);
1623    if (fdivides (buf, bufF, quot))
1624    {
1625      factorsFoundIndex[i - 1]= 1;
1626      factorsFound++;
1627      bufF= quot;
1628      bufF /= Lc (bufF);
1629      reconstructedFactors.append (buf);
1630    }
1631    if (degree (bufF) <= 0)
1632      return;
1633    if (factorsFound + 1 == N.NumCols())
1634    {
1635      reconstructedFactors.append (bufF);
1636      F=1;
1637      return;
1638    }
1639  }
1640  if (reconstructedFactors.length() != 0)
1641    F= bufF (y+eval,y);
1642}
1643
1644#ifdef HAVE_FLINT
1645void
1646reconstructionTry (CFList& reconstructedFactors, CanonicalForm& F, const CFList&
1647                   factors, const int liftBound, int& factorsFound, int*&
1648                   factorsFoundIndex, nmod_mat_t N, const CanonicalForm& eval,
1649                   bool beenInThres
1650                  )
1651{
1652  Variable y= Variable (2);
1653  Variable x= Variable (1);
1654  CanonicalForm yToL= power (y, liftBound);
1655  CanonicalForm bufF= F (y-eval, y);
1656  if (factors.length() == 2)
1657  {
1658    CanonicalForm tmp1, tmp2, tmp3;
1659    tmp1= factors.getFirst();
1660    tmp2= factors.getLast();
1661    tmp1= mulMod2 (tmp1, LC (F,x), yToL);
1662    tmp1 /= content (tmp1, x);
1663    tmp1= tmp1 (y-eval, y);
1664    tmp2= mulMod2 (tmp2, LC (F,x), yToL);
1665    tmp2 /= content (tmp2, x);
1666    tmp2= tmp2 (y-eval, y);
1667    tmp3 = tmp1*tmp2;
1668    if (tmp3/Lc (tmp3) == bufF/Lc (bufF))
1669    {
1670      factorsFound++;
1671      F= 1;
1672      reconstructedFactors.append (tmp1);
1673      reconstructedFactors.append (tmp2);
1674      return;
1675    }
1676  }
1677  CanonicalForm quot, buf;
1678  CFListIterator iter;
1679  for (long i= 0; i < nmod_mat_ncols (N); i++)
1680  {
1681    if (factorsFoundIndex [i] == 1)
1682      continue;
1683    iter= factors;
1684    if (beenInThres)
1685    {
1686      int count= 0;
1687      while (count < i)
1688      {
1689        count++;
1690        iter++;
1691      }
1692      buf= iter.getItem();
1693    }
1694    else
1695    {
1696      buf= 1;
1697      for (long j= 0; j < nmod_mat_nrows (N); j++, iter++)
1698      {
1699        if (!(nmod_mat_entry (N, j, i) == 0))
1700          buf= mulMod2 (buf, iter.getItem(), yToL);
1701      }
1702    }
1703    buf= mulMod2 (buf, LC (F,x), yToL);
1704    buf /= content (buf, x);
1705    buf= buf (y-eval,y);
1706    if (fdivides (buf, bufF, quot))
1707    {
1708      factorsFoundIndex[i]= 1;
1709      factorsFound++;
1710      bufF= quot;
1711      bufF /= Lc (bufF);
1712      reconstructedFactors.append (buf);
1713    }
1714    if (degree (F) <= 0)
1715      return;
1716    if (factorsFound + 1 == nmod_mat_ncols (N))
1717    {
1718      F= 1;
1719      reconstructedFactors.append (bufF);
1720      return;
1721    }
1722  }
1723  if (reconstructedFactors.length() != 0)
1724    F= bufF (y+eval,y);
1725}
1726#endif
1727
1728CFList
1729reconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs, int
1730                precision, const mat_zz_pE& N, const CanonicalForm& eval
1731               )
1732{
1733  Variable y= Variable (2);
1734  Variable x= Variable (1);
1735  CanonicalForm F= G;
1736  CanonicalForm yToL= power (y, precision);
1737  CanonicalForm quot, buf;
1738  CFList result, factorsConsidered;
1739  CFList bufFactors= factors;
1740  CFListIterator iter;
1741  for (long i= 1; i <= N.NumCols(); i++)
1742  {
1743    if (zeroOneVecs [i - 1] == 0)
1744      continue;
1745    iter= factors;
1746    buf= 1;
1747    factorsConsidered= CFList();
1748    for (long j= 1; j <= N.NumRows(); j++, iter++)
1749    {
1750      if (!IsZero (N (j,i)))
1751      {
1752        factorsConsidered.append (iter.getItem());
1753        buf= mulMod2 (buf, iter.getItem(), yToL);
1754      }
1755    }
1756    buf= mulMod2 (buf, LC (F,x), yToL);
1757    buf /= content (buf, x);
1758    if (fdivides (buf, F, quot))
1759    {
1760      F= quot;
1761      F /= Lc (F);
1762      result.append (buf (y-eval,y));
1763      bufFactors= Difference (bufFactors, factorsConsidered);
1764    }
1765    if (degree (F) <= 0)
1766    {
1767      G= F;
1768      factors= bufFactors;
1769      return result;
1770    }
1771  }
1772  G= F;
1773  factors= bufFactors;
1774  return result;
1775}
1776
1777CFList
1778monicReconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs,
1779                     int precision, const mat_zz_pE& N
1780                    )
1781{
1782  Variable y= Variable (2);
1783  Variable x= Variable (1);
1784  CanonicalForm F= G;
1785  CanonicalForm yToL= power (y, precision);
1786  CanonicalForm quot, buf, buf2;
1787  CFList result;
1788  CFList bufFactors= factors;
1789  CFList factorsConsidered;
1790  CFListIterator iter;
1791  for (long i= 1; i <= N.NumCols(); i++)
1792  {
1793    if (zeroOneVecs [i - 1] == 0)
1794      continue;
1795    iter= factors;
1796    buf= 1;
1797    factorsConsidered= CFList();
1798    for (long j= 1; j <= N.NumRows(); j++, iter++)
1799    {
1800      if (!IsZero (N (j,i)))
1801      {
1802        factorsConsidered.append (iter.getItem());
1803        buf= mulMod2 (buf, iter.getItem(), yToL);
1804      }
1805    }
1806    buf2= buf;
1807    buf= mulMod2 (buf, LC (F,x), yToL);
1808    buf /= content (buf, x);
1809    if (fdivides (buf, F, quot))
1810    {
1811      F= quot;
1812      F /= Lc (F);
1813      result.append (buf2);
1814      bufFactors= Difference (bufFactors, factorsConsidered);
1815    }
1816    if (degree (F) <= 0)
1817    {
1818      G= F;
1819      factors= bufFactors;
1820      return result;
1821    }
1822  }
1823  G= F;
1824  factors= bufFactors;
1825  return result;
1826}
1827
1828CFList
1829extReconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs, int
1830                   precision, const mat_zz_p& N, const ExtensionInfo& info,
1831                   const CanonicalForm& evaluation
1832                  )
1833{
1834  Variable y= Variable (2);
1835  Variable x= Variable (1);
1836  Variable alpha= info.getAlpha();
1837  Variable beta= info.getBeta();
1838  int k= info.getGFDegree();
1839  CanonicalForm gamma= info.getGamma();
1840  CanonicalForm delta= info.getDelta();
1841  CanonicalForm F= G;
1842  CanonicalForm yToL= power (y, precision);
1843  CFList result;
1844  CFList bufFactors= factors;
1845  CFList factorsConsidered;
1846  CanonicalForm buf2, quot, buf;
1847  CFListIterator iter;
1848  for (long i= 1; i <= N.NumCols(); i++)
1849  {
1850    if (zeroOneVecs [i - 1] == 0)
1851      continue;
1852    iter= factors;
1853    buf= 1;
1854    factorsConsidered= CFList();
1855    for (long j= 1; j <= N.NumRows(); j++, iter++)
1856    {
1857      if (!IsZero (N (j,i)))
1858      {
1859        factorsConsidered.append (iter.getItem());
1860        buf= mulMod2 (buf, iter.getItem(), yToL);
1861      }
1862    }
1863    buf= mulMod2 (buf, LC (F,x), yToL);
1864    buf /= content (buf, x);
1865    buf2= buf (y-evaluation, y);
1866    buf2 /= Lc (buf2);
1867    if (!k && beta == x)
1868    {
1869      if (degree (buf2, alpha) < 1)
1870      {
1871        if (fdivides (buf, F, quot))
1872        {
1873          F= quot;
1874          F /= Lc (F);
1875          result.append (buf2);
1876          bufFactors= Difference (bufFactors, factorsConsidered);
1877        }
1878      }
1879    }
1880    else
1881    {
1882      CFList source, dest;
1883
1884      if (!isInExtension (buf2, gamma, k, delta, source, dest))
1885      {
1886        if (fdivides (buf, F, quot))
1887        {
1888          F= quot;
1889          F /= Lc (F);
1890          result.append (buf2);
1891          bufFactors= Difference (bufFactors, factorsConsidered);
1892        }
1893      }
1894    }
1895    if (degree (F) <= 0)
1896    {
1897      G= F;
1898      factors= bufFactors;
1899      return result;
1900    }
1901  }
1902  G= F;
1903  factors= bufFactors;
1904  return result;
1905}
1906
1907#ifdef HAVE_FLINT
1908CFList
1909extReconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs, int
1910                   precision, const nmod_mat_t N, const ExtensionInfo& info,
1911                   const CanonicalForm& evaluation
1912                  )
1913{
1914  Variable y= Variable (2);
1915  Variable x= Variable (1);
1916  Variable alpha= info.getAlpha();
1917  Variable beta= info.getBeta();
1918  int k= info.getGFDegree();
1919  CanonicalForm gamma= info.getGamma();
1920  CanonicalForm delta= info.getDelta();
1921  CanonicalForm F= G;
1922  CanonicalForm yToL= power (y, precision);
1923  CFList result;
1924  CFList bufFactors= factors;
1925  CFList factorsConsidered;
1926  CanonicalForm buf2, quot, buf;
1927  CFListIterator iter;
1928  for (long i= 0; i < nmod_mat_ncols(N); i++)
1929  {
1930    if (zeroOneVecs [i] == 0)
1931      continue;
1932    iter= factors;
1933    buf= 1;
1934    factorsConsidered= CFList();
1935    for (long j= 0; j < nmod_mat_ncols(N); j++, iter++)
1936    {
1937      if (!(nmod_mat_entry (N, j, i) == 0))
1938      {
1939        factorsConsidered.append (iter.getItem());
1940        buf= mulMod2 (buf, iter.getItem(), yToL);
1941      }
1942    }
1943    buf= mulMod2 (buf, LC (F,x), yToL);
1944    buf /= content (buf, x);
1945    buf2= buf (y-evaluation, y);
1946    buf2 /= Lc (buf2);
1947    if (!k && beta == x)
1948    {
1949      if (degree (buf2, alpha) < 1)
1950      {
1951        if (fdivides (buf, F, quot))
1952        {
1953          F= quot;
1954          F /= Lc (F);
1955          result.append (buf2);
1956          bufFactors= Difference (bufFactors, factorsConsidered);
1957        }
1958      }
1959    }
1960    else
1961    {
1962      CFList source, dest;
1963
1964      if (!isInExtension (buf2, gamma, k, delta, source, dest))
1965      {
1966        if (fdivides (buf, F, quot))
1967        {
1968          F= quot;
1969          F /= Lc (F);
1970          result.append (buf2);
1971          bufFactors= Difference (bufFactors, factorsConsidered);
1972        }
1973      }
1974    }
1975    if (degree (F) <= 0)
1976    {
1977      G= F;
1978      factors= bufFactors;
1979      return result;
1980    }
1981  }
1982  G= F;
1983  factors= bufFactors;
1984  return result;
1985}
1986#endif
1987
1988CFList
1989reconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs,
1990                int precision, const mat_zz_p& N, const CanonicalForm& eval)
1991{
1992  Variable y= Variable (2);
1993  Variable x= Variable (1);
1994  CanonicalForm F= G;
1995  CanonicalForm yToL= power (y, precision);
1996  CanonicalForm quot, buf;
1997  CFList result;
1998  CFList bufFactors= factors;
1999  CFList factorsConsidered;
2000  CFListIterator iter;
2001  for (long i= 1; i <= N.NumCols(); i++)
2002  {
2003    if (zeroOneVecs [i - 1] == 0)
2004      continue;
2005    iter= factors;
2006    buf= 1;
2007    factorsConsidered= CFList();
2008    for (long j= 1; j <= N.NumRows(); j++, iter++)
2009    {
2010      if (!IsZero (N (j,i)))
2011      {
2012        factorsConsidered.append (iter.getItem());
2013        buf= mulMod2 (buf, iter.getItem(), yToL);
2014      }
2015    }
2016    buf= mulMod2 (buf, LC (F,x), yToL);
2017    buf /= content (buf, x);
2018    if (fdivides (buf, F, quot))
2019    {
2020      F= quot;
2021      F /= Lc (F);
2022      result.append (buf (y-eval,y));
2023      bufFactors= Difference (bufFactors, factorsConsidered);
2024    }
2025    if (degree (F) <= 0)
2026    {
2027      G= F;
2028      factors= bufFactors;
2029      return result;
2030    }
2031  }
2032  G= F;
2033  factors= bufFactors;
2034  return result;
2035}
2036
2037#ifdef HAVE_FLINT
2038CFList
2039reconstruction (CanonicalForm& G, CFList& factors, int* zeroOneVecs,
2040                int precision, const nmod_mat_t N, const CanonicalForm& eval)
2041{
2042  Variable y= Variable (2);
2043  Variable x= Variable (1);
2044  CanonicalForm F= G;
2045  CanonicalForm yToL= power (y, precision);
2046  CanonicalForm quot, buf;
2047  CFList result;
2048  CFList bufFactors= factors;
2049  CFList factorsConsidered;
2050  CFListIterator iter;
2051  for (long i= 0; i < nmod_mat_ncols (N); i++)
2052  {
2053    if (zeroOneVecs [i] == 0)
2054      continue;
2055    iter= factors;
2056    buf= 1;
2057    factorsConsidered= CFList();
2058    for (long j= 0; j < nmod_mat_nrows (N); j++, iter++)
2059    {
2060      if (!(nmod_mat_entry (N, j, i) == 0))
2061      {
2062        factorsConsidered.append (iter.getItem());
2063        buf= mulMod2 (buf, iter.getItem(), yToL);
2064      }
2065    }
2066    buf= mulMod2 (buf, LC (F,x), yToL);
2067    buf /= content (buf, x);
2068    if (fdivides (buf, F, quot))
2069    {
2070      F= quot;
2071      F /= Lc (F);
2072      result.append (buf (y-eval,y));
2073      bufFactors= Difference (bufFactors, factorsConsidered);
2074    }
2075    if (degree (F) <= 0)
2076    {
2077      G= F;
2078      factors= bufFactors;
2079      return result;
2080    }
2081  }
2082  G= F;
2083  factors= bufFactors;
2084  return result;
2085}
2086#endif
2087
2088void
2089extReconstructionTry (CFList& reconstructedFactors, CanonicalForm& F, const
2090                      CFList& factors, const int liftBound, int& factorsFound,
2091                      int*& factorsFoundIndex, mat_zz_p& N, bool beenInThres,
2092                      const ExtensionInfo& info, const CanonicalForm& evaluation
2093                     )
2094{
2095  Variable y= Variable (2);
2096  Variable x= Variable (1);
2097  Variable alpha= info.getAlpha();
2098  Variable beta= info.getBeta();
2099  int k= info.getGFDegree();
2100  CanonicalForm gamma= info.getGamma();
2101  CanonicalForm delta= info.getDelta();
2102  CanonicalForm yToL= power (y, liftBound);
2103  CFList source, dest;
2104  if (factors.length() == 2)
2105  {
2106    CanonicalForm tmp1, tmp2, tmp3;
2107    tmp1= factors.getFirst();
2108    tmp2= factors.getLast();
2109    tmp1= mulMod2 (tmp1, LC (F,x), yToL);
2110    tmp1 /= content (tmp1, x);
2111    tmp2= mulMod2 (tmp2, LC (F,x), yToL);
2112    tmp2 /= content (tmp2, x);
2113    tmp3 = tmp1*tmp2;
2114    if (tmp3/Lc (tmp3) == F/Lc (F))
2115    {
2116      tmp1= tmp1 (y - evaluation, y);
2117      tmp2= tmp2 (y - evaluation, y);
2118      tmp1 /= Lc (tmp1);
2119      tmp2 /= Lc (tmp2);
2120      if (!k && beta == x && degree (tmp2, alpha) < 1 &&
2121          degree (tmp1, alpha) < 1)
2122      {
2123        factorsFound++;
2124        F= 1;
2125        tmp1= mapDown (tmp1, info, source, dest);
2126        tmp2= mapDown (tmp2, info, source, dest);
2127        reconstructedFactors.append (tmp1);
2128        reconstructedFactors.append (tmp2);
2129        return;
2130      }
2131      else if (!isInExtension (tmp2, gamma, k, delta, source, dest) &&
2132               !isInExtension (tmp1, gamma, k, delta, source, dest))
2133      {
2134        factorsFound++;
2135        F= 1;
2136        tmp1= mapDown (tmp1, info, source, dest);
2137        tmp2= mapDown (tmp2, info, source, dest);
2138        reconstructedFactors.append (tmp1);
2139        reconstructedFactors.append (tmp2);
2140        return;
2141      }
2142    }
2143  }
2144  CanonicalForm quot, buf, buf2;
2145  CFListIterator iter;
2146  for (long i= 1; i <= N.NumCols(); i++)
2147  {
2148    if (factorsFoundIndex [i - 1] == 1)
2149      continue;
2150    iter= factors;
2151    if (beenInThres)
2152    {
2153      int count= 1;
2154      while (count < i)
2155      {
2156        count++;
2157        iter++;
2158      }
2159      buf= iter.getItem();
2160    }
2161    else
2162    {
2163      buf= 1;
2164      for (long j= 1; j <= N.NumRows(); j++, iter++)
2165      {
2166        if (!IsZero (N (j,i)))
2167          buf= mulMod2 (buf, iter.getItem(), yToL);
2168      }
2169    }
2170    buf= mulMod2 (buf, LC (F,x), yToL);
2171    buf /= content (buf, x);
2172    buf2= buf (y - evaluation, y);
2173    buf2 /= Lc (buf2);
2174    if (!k && beta == x)
2175    {
2176      if (degree (buf2, alpha) < 1)
2177      {
2178        if (fdivides (buf, F, quot))
2179        {
2180          factorsFoundIndex[i - 1]= 1;
2181          factorsFound++;
2182          F= quot;
2183          F /= Lc (F);
2184          buf2= mapDown (buf2, info, source, dest);
2185          reconstructedFactors.append (buf2);
2186        }
2187      }
2188    }
2189    else
2190    {
2191      if (!isInExtension (buf2, gamma, k, delta, source, dest))
2192      {
2193        if (fdivides (buf, F, quot))
2194        {
2195          factorsFoundIndex[i - 1]= 1;
2196          factorsFound++;
2197          F= quot;
2198          F /= Lc (F);
2199          buf2= mapDown (buf2, info, source, dest);
2200          reconstructedFactors.append (buf2);
2201        }
2202      }
2203    }
2204    if (degree (F) <= 0)
2205      return;
2206    if (factorsFound + 1 == N.NumCols())
2207    {
2208      CanonicalForm tmp= F (y - evaluation, y);
2209      tmp= mapDown (tmp, info, source, dest);
2210      reconstructedFactors.append (tmp);
2211      return;
2212    }
2213  }
2214}
2215
2216#ifdef HAVE_FLINT
2217void
2218extReconstructionTry (CFList& reconstructedFactors, CanonicalForm& F, const
2219                      CFList& factors, const int liftBound, int& factorsFound,
2220                      int*& factorsFoundIndex, nmod_mat_t N, bool beenInThres,
2221                      const ExtensionInfo& info, const CanonicalForm& evaluation
2222                     )
2223{
2224  Variable y= Variable (2);
2225  Variable x= Variable (1);
2226  Variable alpha= info.getAlpha();
2227  Variable beta= info.getBeta();
2228  int k= info.getGFDegree();
2229  CanonicalForm gamma= info.getGamma();
2230  CanonicalForm delta= info.getDelta();
2231  CanonicalForm yToL= power (y, liftBound);
2232  CFList source, dest;
2233  if (factors.length() == 2)
2234  {
2235    CanonicalForm tmp1, tmp2, tmp3;
2236    tmp1= factors.getFirst();
2237    tmp2= factors.getLast();
2238    tmp1= mulMod2 (tmp1, LC (F,x), yToL);
2239    tmp1 /= content (tmp1, x);
2240    tmp2= mulMod2 (tmp2, LC (F,x), yToL);
2241    tmp2 /= content (tmp2, x);
2242    tmp3 = tmp1*tmp2;
2243    if (tmp3/Lc (tmp3) == F/Lc (F))
2244    {
2245      tmp1= tmp1 (y - evaluation, y);
2246      tmp2= tmp2 (y - evaluation, y);
2247      tmp1 /= Lc (tmp1);
2248      tmp2 /= Lc (tmp2);
2249      if (!k && beta == x && degree (tmp2, alpha) < 1 &&
2250          degree (tmp1, alpha) < 1)
2251      {
2252        factorsFound++;
2253        F= 1;
2254        tmp1= mapDown (tmp1, info, source, dest);
2255        tmp2= mapDown (tmp2, info, source, dest);
2256        reconstructedFactors.append (tmp1);
2257        reconstructedFactors.append (tmp2);
2258        return;
2259      }
2260      else if (!isInExtension (tmp2, gamma, k, delta, source, dest) &&
2261               !isInExtension (tmp1, gamma, k, delta, source, dest))
2262      {
2263        factorsFound++;
2264        F= 1;
2265        tmp1= mapDown (tmp1, info, source, dest);
2266        tmp2= mapDown (tmp2, info, source, dest);
2267        reconstructedFactors.append (tmp1);
2268        reconstructedFactors.append (tmp2);
2269        return;
2270      }
2271    }
2272  }
2273  CanonicalForm quot, buf, buf2;
2274  CFListIterator iter;
2275  for (long i= 0; i < nmod_mat_ncols (N); i++)
2276  {
2277    if (factorsFoundIndex [i] == 1)
2278      continue;
2279    iter= factors;
2280    if (beenInThres)
2281    {
2282      int count= 0;
2283      while (count < i)
2284      {
2285        count++;
2286        iter++;
2287      }
2288      buf= iter.getItem();
2289    }
2290    else
2291    {
2292      buf= 1;
2293      for (long j= 0; j < nmod_mat_nrows (N); j++, iter++)
2294      {
2295        if (!(nmod_mat_entry (N, j, i) == 0))
2296          buf= mulMod2 (buf, iter.getItem(), yToL);
2297      }
2298    }
2299    buf= mulMod2 (buf, LC (F,x), yToL);
2300    buf /= content (buf, x);
2301    buf2= buf (y - evaluation, y);
2302    buf2 /= Lc (buf2);
2303    if (!k && beta == x)
2304    {
2305      if (degree (buf2, alpha) < 1)
2306      {
2307        if (fdivides (buf, F, quot))
2308        {
2309          factorsFoundIndex[i]= 1;
2310          factorsFound++;
2311          F= quot;
2312          F /= Lc (F);
2313          buf2= mapDown (buf2, info, source, dest);
2314          reconstructedFactors.append (buf2);
2315        }
2316      }
2317    }
2318    else
2319    {
2320      if (!isInExtension (buf2, gamma, k, delta, source, dest))
2321      {
2322        if (fdivides (buf, F, quot))
2323        {
2324          factorsFoundIndex[i]= 1;
2325          factorsFound++;
2326          F= quot;
2327          F /= Lc (F);
2328          buf2= mapDown (buf2, info, source, dest);
2329          reconstructedFactors.append (buf2);
2330        }
2331      }
2332    }
2333    if (degree (F) <= 0)
2334      return;
2335    if (factorsFound + 1 == nmod_mat_nrows (N))
2336    {
2337      CanonicalForm tmp= F (y - evaluation, y);
2338      tmp= mapDown (tmp, info, source, dest);
2339      reconstructedFactors.append (tmp);
2340      return;
2341    }
2342  }
2343}
2344#endif
2345
2346//over Fp
2347int
2348liftAndComputeLattice (const CanonicalForm& F, int* bounds, int sizeBounds, int
2349                       start, int liftBound, int minBound, CFList& factors,
2350                       mat_zz_p& NTLN, CFList& diophant, CFMatrix& M, CFArray&
2351                       Pi, CFArray& bufQ, bool& irreducible
2352                      )
2353{
2354  CanonicalForm LCF= LC (F, 1);
2355  CFArray *A= new CFArray [factors.length() - 1];
2356  bool wasInBounds= false;
2357  bool hitBound= false;
2358  int l= (minBound+1)*2;
2359  int stepSize= 2;
2360  int oldL= l/2;
2361  bool reduced= false;
2362  mat_zz_p NTLK, *NTLC;
2363  CFMatrix C;
2364  CFArray buf;
2365  CFListIterator j;
2366  CanonicalForm truncF;
2367  Variable y= F.mvar();
2368  while (l <= liftBound)
2369  {
2370    TIMING_START (fac_fq_compute_lattice_lift);
2371    if (start)
2372    {
2373      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
2374      start= 0;
2375    }
2376    else
2377    {
2378      if (wasInBounds)
2379        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
2380      else
2381        henselLift12 (F, factors, l, Pi, diophant, M);
2382    }
2383    TIMING_END_AND_PRINT (fac_fq_compute_lattice_lift,
2384                          "time to lift in compute lattice: ");
2385
2386    factors.insert (LCF);
2387    j= factors;
2388    j++;
2389
2390    truncF= mod (F, power (y, l));
2391    TIMING_START (fac_fq_logarithmic);
2392    for (int i= 0; i < factors.length() - 1; i++, j++)
2393    {
2394      if (!wasInBounds)
2395        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
2396      else
2397        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
2398                                     bufQ[i]);
2399    }
2400    TIMING_END_AND_PRINT (fac_fq_logarithmic,
2401                          "time to compute logarithmic derivative: ");
2402
2403    for (int i= 0; i < sizeBounds; i++)
2404    {
2405      if (bounds [i] + 1 <= l/2)
2406      {
2407        wasInBounds= true;
2408        int k= tmin (bounds [i] + 1, l/2);
2409        C= CFMatrix (l - k, factors.length() - 1);
2410        for (int ii= 0; ii < factors.length() - 1; ii++)
2411        {
2412          if (A[ii].size() - 1 >= i)
2413          {
2414            buf= getCoeffs (A[ii] [i], k);
2415            writeInMatrix (C, buf, ii + 1, 0);
2416          }
2417        }
2418        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
2419        NTLK= (*NTLC)*NTLN;
2420        transpose (NTLK, NTLK);
2421        kernel (NTLK, NTLK);
2422        transpose (NTLK, NTLK);
2423        NTLN *= NTLK;
2424
2425        if (NTLN.NumCols() == 1)
2426        {
2427          irreducible= true;
2428          break;
2429        }
2430        if (isReduced (NTLN) && l > (minBound+1)*2)
2431        {
2432          reduced= true;
2433          break;
2434        }
2435      }
2436    }
2437
2438    if (irreducible)
2439      break;
2440    if (reduced)
2441      break;
2442    oldL= l;
2443    l += stepSize;
2444    stepSize *= 2;
2445    if (l > liftBound)
2446    {
2447      if (!hitBound)
2448      {
2449        l= liftBound;
2450        hitBound= true;
2451      }
2452      else
2453        break;
2454    }
2455  }
2456  delete [] A;
2457  if (!wasInBounds)
2458  {
2459    if (start)
2460      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
2461    else
2462      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
2463    factors.insert (LCF);
2464  }
2465  return l;
2466}
2467
2468#ifdef HAVE_FLINT
2469int
2470liftAndComputeLattice (const CanonicalForm& F, int* bounds, int sizeBounds, int
2471                       start, int liftBound, int minBound, CFList& factors,
2472                       nmod_mat_t FLINTN, CFList& diophant, CFMatrix& M,CFArray&
2473                       Pi, CFArray& bufQ, bool& irreducible
2474                      )
2475{
2476  CanonicalForm LCF= LC (F, 1);
2477  CFArray *A= new CFArray [factors.length() - 1];
2478  bool wasInBounds= false;
2479  bool hitBound= false;
2480  int l= (minBound+1)*2;
2481  int stepSize= 2;
2482  int oldL= l/2;
2483  bool reduced= false;
2484  long rank;
2485  nmod_mat_t FLINTK, FLINTC, null;
2486  CFMatrix C;
2487  CFArray buf;
2488  CFListIterator j;
2489  CanonicalForm truncF;
2490  Variable y= F.mvar();
2491  while (l <= liftBound)
2492  {
2493    TIMING_START (fac_fq_compute_lattice_lift);
2494    if (start)
2495    {
2496      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
2497      start= 0;
2498    }
2499    else
2500    {
2501      if (wasInBounds)
2502        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
2503      else
2504        henselLift12 (F, factors, l, Pi, diophant, M);
2505    }
2506    TIMING_END_AND_PRINT (fac_fq_compute_lattice_lift,
2507                          "time to lift in compute lattice: ");
2508
2509    factors.insert (LCF);
2510    j= factors;
2511    j++;
2512
2513    truncF= mod (F, power (y, l));
2514    TIMING_START (fac_fq_logarithmic);
2515    for (int i= 0; i < factors.length() - 1; i++, j++)
2516    {
2517      if (!wasInBounds)
2518        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
2519      else
2520        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
2521                                     bufQ[i]);
2522    }
2523    TIMING_END_AND_PRINT (fac_fq_logarithmic,
2524                          "time to compute logarithmic derivative: ");
2525
2526    for (int i= 0; i < sizeBounds; i++)
2527    {
2528      if (bounds [i] + 1 <= l/2)
2529      {
2530        wasInBounds= true;
2531        int k= tmin (bounds [i] + 1, l/2);
2532        C= CFMatrix (l - k, factors.length() - 1);
2533        for (int ii= 0; ii < factors.length() - 1; ii++)
2534        {
2535          if (A[ii].size() - 1 >= i)
2536          {
2537            buf= getCoeffs (A[ii] [i], k);
2538            writeInMatrix (C, buf, ii + 1, 0);
2539          }
2540        }
2541
2542        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
2543        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
2544                       getCharacteristic());
2545        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
2546        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
2547                       getCharacteristic());
2548        rank= nmod_mat_nullspace (null, FLINTK);
2549        nmod_mat_clear (FLINTK);
2550        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
2551        nmod_mat_clear (FLINTC);
2552        nmod_mat_init_set (FLINTC, FLINTN);
2553        nmod_mat_clear (FLINTN);
2554        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
2555                       getCharacteristic());
2556        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
2557
2558        nmod_mat_clear (FLINTC);
2559        nmod_mat_window_clear (FLINTK);
2560        nmod_mat_clear (null);
2561        if (nmod_mat_ncols (FLINTN) == 1)
2562        {
2563          irreducible= true;
2564          break;
2565        }
2566        if (isReduced (FLINTN) && l > (minBound+1)*2)
2567        {
2568          reduced= true;
2569          break;
2570        }
2571      }
2572    }
2573
2574    if (irreducible)
2575      break;
2576    if (reduced)
2577      break;
2578    oldL= l;
2579    l += stepSize;
2580    stepSize *= 2;
2581    if (l > liftBound)
2582    {
2583      if (!hitBound)
2584      {
2585        l= liftBound;
2586        hitBound= true;
2587      }
2588      else
2589        break;
2590    }
2591  }
2592  delete [] A;
2593  if (!wasInBounds)
2594  {
2595    if (start)
2596      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
2597    else
2598      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
2599    factors.insert (LCF);
2600  }
2601  return l;
2602}
2603#endif
2604
2605//over field extension
2606int
2607extLiftAndComputeLattice (const CanonicalForm& F, int* bounds, int sizeBounds,
2608                          int liftBound, int minBound, int start, CFList&
2609                          factors, mat_zz_p& NTLN, CFList& diophant,
2610                          CFMatrix& M, CFArray& Pi, CFArray& bufQ, bool&
2611                          irreducible, const CanonicalForm& evaluation, const
2612                          ExtensionInfo& info, CFList& source, CFList& dest
2613                         )
2614{
2615  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
2616  CanonicalForm LCF= LC (F, 1);
2617  CFArray *A= new CFArray [factors.length() - 1];
2618  bool wasInBounds= false;
2619  bool hitBound= false;
2620  int degMipo;
2621  Variable alpha;
2622  alpha= info.getAlpha();
2623  degMipo= degree (getMipo (alpha));
2624
2625  Variable gamma= info.getBeta();
2626  CanonicalForm primElemAlpha= info.getGamma();
2627  CanonicalForm imPrimElemAlpha= info.getDelta();
2628
2629  int stepSize= 2;
2630  int l= ((minBound+1)/degMipo+1)*2;
2631  l= tmax (l, 2);
2632  if (start > l)
2633    l= start;
2634  int oldL= l/2;
2635  bool reduced= false;
2636  Variable y= F.mvar();
2637  Variable x= Variable (1);
2638  CanonicalForm powX, imBasis, truncF;
2639  CFMatrix Mat, C;
2640  CFArray buf;
2641  CFIterator iter;
2642  mat_zz_p* NTLMat, *NTLC, NTLK;
2643  CFListIterator j;
2644  while (l <= liftBound)
2645  {
2646    TIMING_START (fac_fq_compute_lattice_lift);
2647    if (start)
2648    {
2649      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
2650      start= 0;
2651    }
2652    else
2653    {
2654      if (wasInBounds)
2655        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
2656      else
2657        henselLift12 (F, factors, l, Pi, diophant, M);
2658    }
2659    TIMING_END_AND_PRINT (fac_fq_compute_lattice_lift,
2660                          "time to lift in compute lattice: ");
2661
2662    factors.insert (LCF);
2663
2664    if (GF)
2665      setCharacteristic (getCharacteristic());
2666
2667    powX= power (y-gamma, l);
2668    Mat= CFMatrix (l*degMipo, l*degMipo);
2669    for (int i= 0; i < l*degMipo; i++)
2670    {
2671      imBasis= mod (power (y, i), powX);
2672      imBasis= imBasis (power (y, degMipo), y);
2673      imBasis= imBasis (y, gamma);
2674      iter= imBasis;
2675      for (; iter.hasTerms(); iter++)
2676        Mat (iter.exp()+ 1, i+1)= iter.coeff();
2677    }
2678
2679    NTLMat= convertFacCFMatrix2NTLmat_zz_p (Mat);
2680    *NTLMat= inv (*NTLMat);
2681
2682    if (GF)
2683      setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2684
2685    j= factors;
2686    j++;
2687
2688    truncF= mod (F, power (y, l));
2689    TIMING_START (fac_fq_logarithmic);
2690    for (int i= 0; i < factors.length() - 1; i++, j++)
2691    {
2692      if (!wasInBounds)
2693        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
2694      else
2695        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
2696                                     bufQ[i]);
2697    }
2698    TIMING_END_AND_PRINT (fac_fq_logarithmic,
2699                          "time to compute logarithmic derivative: ");
2700
2701    for (int i= 0; i < sizeBounds; i++)
2702    {
2703      if (bounds [i] + 1 <= (l/2)*degMipo)
2704      {
2705        wasInBounds= true;
2706        int k= tmin (bounds [i] + 1, (l/2)*degMipo);
2707        C= CFMatrix (l*degMipo - k, factors.length() - 1);
2708
2709        for (int ii= 0; ii < factors.length() - 1; ii++)
2710        {
2711          if (A[ii].size() - 1 >= i)
2712          {
2713            if (GF)
2714            {
2715              A [ii] [i]= A [ii] [i] (y-evaluation, y);
2716              setCharacteristic (getCharacteristic());
2717              A[ii] [i]= GF2FalphaRep (A[ii] [i], alpha);
2718              if (alpha != gamma)
2719                A [ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
2720                                     gamma, source, dest
2721                                    );
2722              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
2723            }
2724            else
2725            {
2726              A [ii] [i]= A [ii] [i] (y-evaluation, y);
2727              if (alpha != gamma)
2728                A[ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
2729                                    gamma, source, dest
2730                                   );
2731              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
2732            }
2733            writeInMatrix (C, buf, ii + 1, 0);
2734          }
2735          if (GF)
2736            setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2737        }
2738
2739        if (GF)
2740          setCharacteristic(getCharacteristic());
2741
2742        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
2743        NTLK= (*NTLC)*NTLN;
2744        transpose (NTLK, NTLK);
2745        kernel (NTLK, NTLK);
2746        transpose (NTLK, NTLK);
2747        NTLN *= NTLK;
2748
2749        if (GF)
2750          setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2751
2752        if (NTLN.NumCols() == 1)
2753        {
2754          irreducible= true;
2755          break;
2756        }
2757        if (isReduced (NTLN))
2758        {
2759          reduced= true;
2760          break;
2761        }
2762      }
2763    }
2764
2765    if (NTLN.NumCols() == 1)
2766    {
2767      irreducible= true;
2768      break;
2769    }
2770    if (reduced)
2771      break;
2772    oldL= l;
2773    l += stepSize;
2774    stepSize *= 2;
2775    if (l > liftBound)
2776    {
2777      if (!hitBound)
2778      {
2779        l= liftBound;
2780        hitBound= true;
2781      }
2782      else
2783        break;
2784    }
2785  }
2786  delete [] A;
2787  if (!wasInBounds)
2788  {
2789    if (start)
2790      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
2791    else
2792      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
2793    factors.insert (LCF);
2794  }
2795  return l;
2796}
2797
2798/*#ifdef HAVE_FLINT
2799//over field extension
2800int
2801extLiftAndComputeLattice (const CanonicalForm& F, int* bounds, int sizeBounds,
2802                          int liftBound, int minBound, int start, CFList&
2803                          factors, nmod_mat_t FLINTN, CFList& diophant,
2804                          CFMatrix& M, CFArray& Pi, CFArray& bufQ, bool&
2805                          irreducible, const CanonicalForm& evaluation, const
2806                          ExtensionInfo& info, CFList& source, CFList& dest
2807                         )
2808{
2809  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
2810  CanonicalForm LCF= LC (F, 1);
2811  CFArray *A= new CFArray [factors.length() - 1];
2812  bool wasInBounds= false;
2813  bool hitBound= false;
2814  int degMipo;
2815  Variable alpha;
2816  alpha= info.getAlpha();
2817  degMipo= degree (getMipo (alpha));
2818
2819  Variable gamma= info.getBeta();
2820  CanonicalForm primElemAlpha= info.getGamma();
2821  CanonicalForm imPrimElemAlpha= info.getDelta();
2822
2823  int stepSize= 2;
2824  int l= ((minBound+1)/degMipo+1)*2;
2825  l= tmax (l, 2);
2826  if (start > l)
2827    l= start;
2828  int oldL= l/2;
2829  bool reduced= false;
2830  Variable y= F.mvar();
2831  Variable x= Variable (1);
2832  CanonicalForm powX, imBasis, truncF;
2833  CFMatrix Mat, C;
2834  CFArray buf;
2835  CFIterator iter;
2836  long rank;
2837  nmod_mat_t FLINTMat, FLINTMatInv, FLINTC, FLINTK, null;
2838  CFListIterator j;
2839  while (l <= liftBound)
2840  {
2841    if (start)
2842    {
2843      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
2844      start= 0;
2845    }
2846    else
2847    {
2848      if (wasInBounds)
2849        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
2850      else
2851        henselLift12 (F, factors, l, Pi, diophant, M);
2852    }
2853
2854    factors.insert (LCF);
2855
2856    if (GF)
2857      setCharacteristic (getCharacteristic());
2858
2859    powX= power (y-gamma, l);
2860    Mat= CFMatrix (l*degMipo, l*degMipo);
2861    for (int i= 0; i < l*degMipo; i++)
2862    {
2863      imBasis= mod (power (y, i), powX);
2864      imBasis= imBasis (power (y, degMipo), y);
2865      imBasis= imBasis (y, gamma);
2866      iter= imBasis;
2867      for (; iter.hasTerms(); iter++)
2868        Mat (iter.exp()+ 1, i+1)= iter.coeff();
2869    }
2870
2871    convertFacCFMatrix2nmod_mat_t (FLINTMat, Mat);
2872    nmod_mat_init (FLINTMatInv, nmod_mat_nrows (FLINTMat),
2873                   nmod_mat_nrows (FLINTMat), getCharacteristic());
2874    nmod_mat_inv (FLINTMatInv, FLINTMat);
2875
2876    if (GF)
2877      setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2878
2879    j= factors;
2880    j++;
2881
2882    truncF= mod (F, power (y, l));
2883    for (int i= 0; i < factors.length() - 1; i++, j++)
2884    {
2885      if (!wasInBounds)
2886        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
2887      else
2888        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
2889                                     bufQ[i]);
2890    }
2891
2892    for (int i= 0; i < sizeBounds; i++)
2893    {
2894      if (bounds [i] + 1 <= (l/2)*degMipo)
2895      {
2896        wasInBounds= true;
2897        int k= tmin (bounds [i] + 1, (l/2)*degMipo);
2898        C= CFMatrix (l*degMipo - k, factors.length() - 1);
2899
2900        for (int ii= 0; ii < factors.length() - 1; ii++)
2901        {
2902          if (A[ii].size() - 1 >= i)
2903          {
2904            if (GF)
2905            {
2906              A [ii] [i]= A [ii] [i] (y-evaluation, y);
2907              setCharacteristic (getCharacteristic());
2908              A[ii] [i]= GF2FalphaRep (A[ii] [i], alpha);
2909              if (alpha != gamma)
2910                A [ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
2911                                     gamma, source, dest
2912                                    );
2913              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, FLINTMatInv); //TODO
2914            }
2915            else
2916            {
2917              A [ii] [i]= A [ii] [i] (y-evaluation, y);
2918              if (alpha != gamma)
2919                A[ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
2920                                    gamma, source, dest
2921                                   );
2922              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, FLINTMatInv); //TODO
2923            }
2924            writeInMatrix (C, buf, ii + 1, 0);
2925          }
2926          if (GF)
2927            setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2928        }
2929
2930        if (GF)
2931          setCharacteristic(getCharacteristic());
2932
2933        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
2934        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
2935                       getCharacteristic());
2936        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
2937        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
2938                       getCharacteristic());
2939        rank= nmod_mat_nullspace (null, FLINTK);
2940        nmod_mat_clear (FLINTK);
2941        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
2942        nmod_mat_clear (FLINTC);
2943        nmod_mat_init_set (FLINTC, FLINTN);
2944        nmod_mat_clear (FLINTN);
2945        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
2946                       getCharacteristic());
2947        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
2948
2949        nmod_mat_clear (FLINTC);
2950        nmod_mat_window_clear (FLINTK);
2951        nmod_mat_clear (null);
2952
2953        if (GF)
2954          setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
2955
2956        if (nmod_mat_ncols (FLINTN) == 1)
2957        {
2958          irreducible= true;
2959          break;
2960        }
2961        if (isReduced (FLINTN))
2962        {
2963          reduced= true;
2964          break;
2965        }
2966      }
2967    }
2968
2969    nmod_mat_clear (FLINTMat);
2970    nmod_mat_clear (FLINTMatInv);
2971
2972    if (nmod_mat_ncols (FLINTN) == 1)
2973    {
2974      irreducible= true;
2975      break;
2976    }
2977    if (reduced)
2978      break;
2979    oldL= l;
2980    l += stepSize;
2981    stepSize *= 2;
2982    if (l > liftBound)
2983    {
2984      if (!hitBound)
2985      {
2986        l= liftBound;
2987        hitBound= true;
2988      }
2989      else
2990        break;
2991    }
2992  }
2993  delete [] A;
2994  if (!wasInBounds)
2995  {
2996    if (start)
2997      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
2998    else
2999      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
3000    factors.insert (LCF);
3001  }
3002  return l;
3003}
3004#endif*/
3005
3006// over Fq
3007int
3008liftAndComputeLattice (const CanonicalForm& F, int* bounds, int sizeBounds,
3009                       int start, int liftBound, int minBound, CFList& factors,
3010                       mat_zz_pE& NTLN, CFList& diophant, CFMatrix& M, CFArray&
3011                       Pi, CFArray& bufQ, bool& irreducible
3012                      )
3013{
3014  CanonicalForm LCF= LC (F, 1);
3015  CFArray *A= new CFArray [factors.length() - 1];
3016  bool wasInBounds= false;
3017  bool hitBound= false;
3018  int l= (minBound+1)*2;
3019  int stepSize= 2;
3020  int oldL= l/2;
3021  bool reduced= false;
3022  CFListIterator j;
3023  mat_zz_pE* NTLC, NTLK;
3024  CFArray buf;
3025  CFMatrix C;
3026  Variable y= F.mvar();
3027  CanonicalForm truncF;
3028  while (l <= liftBound)
3029  {
3030    TIMING_START (fac_fq_compute_lattice_lift);
3031    if (start)
3032    {
3033      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
3034      start= 0;
3035    }
3036    else
3037    {
3038      if (wasInBounds)
3039        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
3040      else
3041        henselLift12 (F, factors, l, Pi, diophant, M);
3042    }
3043    TIMING_END_AND_PRINT (fac_fq_compute_lattice_lift,
3044                          "time to lift in compute lattice: ");
3045
3046    factors.insert (LCF);
3047    j= factors;
3048    j++;
3049
3050    truncF= mod (F, power (y,l));
3051    TIMING_START (fac_fq_logarithmic);
3052    for (int i= 0; i < factors.length() - 1; i++, j++)
3053    {
3054      if (l == (minBound+1)*2)
3055      {
3056        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
3057      }
3058      else
3059      {
3060        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
3061                                     bufQ[i]
3062                                    );
3063      }
3064    }
3065    TIMING_END_AND_PRINT (fac_fq_logarithmic,
3066                          "time to compute logarithmic derivative: ");
3067
3068    for (int i= 0; i < sizeBounds; i++)
3069    {
3070      if (bounds [i] + 1 <= l/2)
3071      {
3072        wasInBounds= true;
3073        int k= tmin (bounds [i] + 1, l/2);
3074        C= CFMatrix (l - k, factors.length() - 1);
3075        for (int ii= 0; ii < factors.length() - 1; ii++)
3076        {
3077
3078          if (A[ii].size() - 1 >= i)
3079          {
3080            buf= getCoeffs (A[ii] [i], k);
3081            writeInMatrix (C, buf, ii + 1, 0);
3082          }
3083        }
3084
3085        NTLC= convertFacCFMatrix2NTLmat_zz_pE(C);
3086        NTLK= (*NTLC)*NTLN;
3087        transpose (NTLK, NTLK);
3088        kernel (NTLK, NTLK);
3089        transpose (NTLK, NTLK);
3090        NTLN *= NTLK;
3091
3092        if (NTLN.NumCols() == 1)
3093        {
3094          irreducible= true;
3095          break;
3096        }
3097        if (isReduced (NTLN) && l > (minBound+1)*2)
3098        {
3099          reduced= true;
3100          break;
3101        }
3102      }
3103    }
3104
3105    if (NTLN.NumCols() == 1)
3106    {
3107      irreducible= true;
3108      break;
3109    }
3110    if (reduced)
3111      break;
3112    oldL= l;
3113    l += stepSize;
3114    stepSize *= 2;
3115    if (l > liftBound)
3116    {
3117      if (!hitBound)
3118      {
3119        l= liftBound;
3120        hitBound= true;
3121      }
3122      else
3123        break;
3124    }
3125  }
3126  delete [] A;
3127  if (!wasInBounds)
3128  {
3129    if (start)
3130      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
3131    else
3132      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
3133    factors.insert (LCF);
3134  }
3135  return l;
3136}
3137
3138#ifdef HAVE_FLINT
3139int
3140liftAndComputeLatticeFq2Fp (const CanonicalForm& F, int* bounds, int sizeBounds,
3141                            int start, int liftBound, int minBound, CFList&
3142                            factors, nmod_mat_t FLINTN, CFList& diophant,
3143                            CFMatrix& M, CFArray& Pi, CFArray& bufQ, bool&
3144                            irreducible, const Variable& alpha
3145                           )
3146#else
3147int
3148liftAndComputeLatticeFq2Fp (const CanonicalForm& F, int* bounds, int sizeBounds,
3149                            int start, int liftBound, int minBound, CFList&
3150                            factors, mat_zz_p& NTLN, CFList& diophant, CFMatrix&
3151                            M, CFArray& Pi, CFArray& bufQ, bool& irreducible,
3152                            const Variable& alpha
3153                           )
3154#endif
3155{
3156  CanonicalForm LCF= LC (F, 1);
3157  CFArray *A= new CFArray [factors.length() - 1];
3158  bool wasInBounds= false;
3159  int l= (minBound+1)*2;
3160  int oldL= l/2;
3161  int stepSize= 2;
3162  bool hitBound= false;
3163  int extensionDeg= degree (getMipo (alpha));
3164  bool reduced= false;
3165  CFListIterator j;
3166  CFMatrix C;
3167  CFArray buf;
3168#ifdef HAVE_FLINT
3169  long rank;
3170  nmod_mat_t FLINTC, FLINTK, null;
3171#else
3172  mat_zz_p* NTLC, NTLK;
3173#endif
3174  Variable y= F.mvar();
3175  CanonicalForm truncF;
3176  while (l <= liftBound)
3177  {
3178    TIMING_START (fac_fq_compute_lattice_lift);
3179    if (start)
3180    {
3181      henselLiftResume12 (F, factors, start, l, Pi, diophant, M);
3182      start= 0;
3183    }
3184    else
3185    {
3186      if (wasInBounds)
3187        henselLiftResume12 (F, factors, oldL, l, Pi, diophant, M);
3188      else
3189        henselLift12 (F, factors, l, Pi, diophant, M);
3190    }
3191    TIMING_END_AND_PRINT (fac_fq_compute_lattice_lift,
3192                          "time to lift in compute lattice: ");
3193
3194    factors.insert (LCF);
3195    j= factors;
3196    j++;
3197
3198    truncF= mod (F, power (y,l));
3199    TIMING_START (fac_fq_logarithmic);
3200    for (int i= 0; i < factors.length() - 1; i++, j++)
3201    {
3202      if (l == (minBound+1)*2)
3203      {
3204        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ[i]);
3205      }
3206      else
3207      {
3208        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
3209                                     bufQ[i]
3210                                    );
3211      }
3212    }
3213    TIMING_END_AND_PRINT (fac_fq_logarithmic,
3214                          "time to compute logarithmic derivative: ");
3215
3216    for (int i= 0; i < sizeBounds; i++)
3217    {
3218      if (bounds [i] + 1 <= l/2)
3219      {
3220        wasInBounds= true;
3221        int k= tmin (bounds [i] + 1, l/2);
3222        C= CFMatrix ((l - k)*extensionDeg, factors.length() - 1);
3223        for (int ii= 0; ii < factors.length() - 1; ii++)
3224        {
3225          if (A[ii].size() - 1 >= i)
3226          {
3227            buf= getCoeffs (A[ii] [i], k, alpha);
3228            writeInMatrix (C, buf, ii + 1, 0);
3229          }
3230        }
3231
3232#ifdef HAVE_FLINT
3233        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
3234        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
3235                       getCharacteristic());
3236        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
3237        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
3238                       getCharacteristic());
3239        rank= nmod_mat_nullspace (null, FLINTK);
3240        nmod_mat_clear (FLINTK);
3241        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
3242        nmod_mat_clear (FLINTC);
3243        nmod_mat_init_set (FLINTC, FLINTN);
3244        nmod_mat_clear (FLINTN);
3245        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
3246                       getCharacteristic());
3247        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
3248
3249        nmod_mat_clear (FLINTC);
3250        nmod_mat_window_clear (FLINTK);
3251        nmod_mat_clear (null);
3252#else
3253        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
3254        NTLK= (*NTLC)*NTLN;
3255        transpose (NTLK, NTLK);
3256        kernel (NTLK, NTLK);
3257        transpose (NTLK, NTLK);
3258        NTLN *= NTLK;
3259#endif
3260
3261#ifdef HAVE_FLINT
3262        if (nmod_mat_nrows (FLINTN) == 1)
3263#else
3264        if (NTLN.NumCols() == 1)
3265#endif
3266        {
3267          irreducible= true;
3268          break;
3269        }
3270#ifdef HAVE_FLINT
3271        if (isReduced (FLINTN) && l > (minBound+1)*2)
3272#else
3273        if (isReduced (NTLN) && l > (minBound+1)*2)
3274#endif
3275        {
3276          reduced= true;
3277          break;
3278        }
3279      }
3280    }
3281
3282#ifdef HAVE_FLINT
3283    if (nmod_mat_ncols (FLINTN) == 1)
3284#else
3285    if (NTLN.NumCols() == 1)
3286#endif
3287    {
3288      irreducible= true;
3289      break;
3290    }
3291    if (reduced)
3292      break;
3293    oldL= l;
3294    l += stepSize;
3295    stepSize *= 2;
3296    if (l > liftBound)
3297    {
3298      if (!hitBound)
3299      {
3300        l= liftBound;
3301        hitBound= true;
3302      }
3303      else
3304        break;
3305    }
3306  }
3307  delete [] A;
3308  if (!wasInBounds)
3309  {
3310    if (start)
3311      henselLiftResume12 (F, factors, start, degree (F) + 1, Pi, diophant, M);
3312    else
3313      henselLift12 (F, factors, degree (F) + 1, Pi, diophant, M);
3314    factors.insert (LCF);
3315  }
3316  return l;
3317}
3318
3319CFList
3320increasePrecision (CanonicalForm& F, CFList& factors, int factorsFound,
3321                   int oldNumCols, int oldL, int precision,
3322                   const CanonicalForm& eval
3323                  )
3324{
3325  int d;
3326  bool isIrreducible= false;
3327  int* bounds= computeBounds (F, d, isIrreducible);
3328  Variable y= F.mvar();
3329  if (isIrreducible)
3330  {
3331    delete [] bounds;
3332    CanonicalForm G= F;
3333    F= 1;
3334    return CFList (G (y-eval, y));
3335  }
3336  CFArray * A= new CFArray [factors.length()];
3337  CFArray bufQ= CFArray (factors.length());
3338#ifdef HAVE_FLINT
3339  nmod_mat_t FLINTN;
3340  nmod_mat_init (FLINTN,factors.length(),factors.length(), getCharacteristic());
3341  for (long i=factors.length()-1; i >= 0; i--)
3342    nmod_mat_entry (FLINTN, i, i)= 1;
3343#else
3344  mat_zz_p NTLN;
3345  ident (NTLN, factors.length());
3346#endif
3347  int minBound= bounds[0];
3348  for (int i= 1; i < d; i++)
3349  {
3350    if (bounds[i] != 0)
3351      minBound= tmin (minBound, bounds[i]);
3352  }
3353  int l= tmax (2*(minBound + 1), oldL);
3354  int oldL2= l/2;
3355  int stepSize= 2;
3356  bool useOldQs= false;
3357  bool hitBound= false;
3358  CFListIterator j;
3359  CFMatrix C;
3360  CFArray buf;
3361#ifdef HAVE_FLINT
3362  long rank;
3363  nmod_mat_t FLINTC, FLINTK, null;
3364#else
3365  mat_zz_p* NTLC, NTLK;
3366#endif
3367  CanonicalForm truncF;
3368  while (l <= precision)
3369  {
3370    j= factors;
3371    truncF= mod (F, power (y,l));
3372    if (useOldQs)
3373    {
3374      for (int i= 0; i < factors.length(); i++, j++)
3375        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL2, bufQ[i],
3376                                     bufQ[i]
3377                                    );
3378    }
3379    else
3380    {
3381      for (int i= 0; i < factors.length(); i++, j++)
3382        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
3383    }
3384    useOldQs= true;
3385    for (int i= 0; i < d; i++)
3386    {
3387      if (bounds [i] + 1 <= l/2)
3388      {
3389        int k= tmin (bounds [i] + 1, l/2);
3390        C= CFMatrix (l - k, factors.length());
3391        for (int ii= 0; ii < factors.length(); ii++)
3392        {
3393          if (A[ii].size() - 1 >= i)
3394          {
3395            buf= getCoeffs (A[ii] [i], k);
3396            writeInMatrix (C, buf, ii + 1, 0);
3397          }
3398        }
3399#ifdef HAVE_FLINT
3400        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
3401        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
3402                       getCharacteristic());
3403        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
3404        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
3405                       getCharacteristic());
3406        rank= nmod_mat_nullspace (null, FLINTK);
3407        nmod_mat_clear (FLINTK);
3408        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
3409        nmod_mat_clear (FLINTC);
3410        nmod_mat_init_set (FLINTC, FLINTN);
3411        nmod_mat_clear (FLINTN);
3412        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
3413                       getCharacteristic());
3414        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
3415
3416        nmod_mat_clear (FLINTC);
3417        nmod_mat_window_clear (FLINTK);
3418        nmod_mat_clear (null);
3419#else
3420        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
3421        NTLK= (*NTLC)*NTLN;
3422        transpose (NTLK, NTLK);
3423        kernel (NTLK, NTLK);
3424        transpose (NTLK, NTLK);
3425        NTLN *= NTLK;
3426#endif
3427#ifdef HAVE_FLINT
3428        if (nmod_mat_ncols (FLINTN) == 1)
3429        {
3430          nmod_mat_clear (FLINTN);
3431#else
3432        if (NTLN.NumCols() == 1)
3433        {
3434#endif
3435          delete [] A;
3436          delete [] bounds;
3437          CanonicalForm G= F;
3438          F= 1;
3439          return CFList (G (y-eval,y));
3440        }
3441      }
3442    }
3443
3444#ifdef HAVE_FLINT
3445    if (nmod_mat_ncols (FLINTN) < oldNumCols - factorsFound)
3446    {
3447      if (isReduced (FLINTN))
3448      {
3449        int * factorsFoundIndex= new int [nmod_mat_ncols (FLINTN)];
3450        for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
3451#else
3452    if (NTLN.NumCols() < oldNumCols - factorsFound)
3453    {
3454      if (isReduced (NTLN))
3455      {
3456        int * factorsFoundIndex= new int [NTLN.NumCols()];
3457        for (long i= 0; i < NTLN.NumCols(); i++)
3458#endif
3459          factorsFoundIndex[i]= 0;
3460        int factorsFound2= 0;
3461        CFList result;
3462        CanonicalForm bufF= F;
3463#ifdef HAVE_FLINT
3464        reconstructionTry (result, bufF, factors, degree (F) + 1, factorsFound2,
3465                           factorsFoundIndex, FLINTN, eval, false
3466                          );
3467        if (result.length() == nmod_mat_ncols (FLINTN))
3468        {
3469          nmod_mat_clear (FLINTN);
3470#else
3471        reconstructionTry (result, bufF, factors, degree (F) + 1, factorsFound2,
3472                           factorsFoundIndex, NTLN, eval, false
3473                          );
3474        if (result.length() == NTLN.NumCols())
3475        {
3476#endif
3477          delete [] factorsFoundIndex;
3478          delete [] A;
3479          delete [] bounds;
3480          F= 1;
3481          return result;
3482        }
3483        delete [] factorsFoundIndex;
3484      }
3485      else if (l == precision)
3486      {
3487        CanonicalForm bufF= F;
3488#ifdef HAVE_FLINT
3489        int * zeroOne= extractZeroOneVecs (FLINTN);
3490        CFList result= reconstruction (bufF,factors,zeroOne,precision,FLINTN, eval);
3491        nmod_mat_clear (FLINTN);
3492#else
3493        int * zeroOne= extractZeroOneVecs (NTLN);
3494        CFList result= reconstruction (bufF, factors, zeroOne, precision, NTLN, eval);
3495#endif
3496        F= bufF;
3497        delete [] zeroOne;
3498        delete [] A;
3499        delete [] bounds;
3500        return result;
3501      }
3502    }
3503    oldL2= l;
3504    l += stepSize;
3505    stepSize *= 2;
3506    if (l > precision)
3507    {
3508      if (!hitBound)
3509      {
3510        l= precision;
3511        hitBound= true;
3512      }
3513      else
3514        break;
3515    }
3516  }
3517#ifdef HAVE_FLINT
3518  nmod_mat_clear (FLINTN);
3519#endif
3520  delete [] bounds;
3521  delete [] A;
3522  return CFList();
3523}
3524
3525CFList
3526increasePrecision (CanonicalForm& F, CFList& factors, int factorsFound,
3527                   int oldNumCols, int oldL, const Variable&,
3528                   int precision, const CanonicalForm& eval
3529                  )
3530{
3531  int d;
3532  bool isIrreducible= false;
3533  Variable y= F.mvar();
3534  int* bounds= computeBounds (F, d, isIrreducible);
3535  if (isIrreducible)
3536  {
3537    delete [] bounds;
3538    CanonicalForm G= F;
3539    F= 1;
3540    return CFList (G (y-eval,y));
3541  }
3542  CFArray * A= new CFArray [factors.length()];
3543  CFArray bufQ= CFArray (factors.length());
3544  mat_zz_pE NTLN;
3545  ident (NTLN, factors.length());
3546  int minBound= bounds[0];
3547  for (int i= 1; i < d; i++)
3548  {
3549    if (bounds[i] != 0)
3550      minBound= tmin (minBound, bounds[i]);
3551  }
3552  int l= tmax (2*(minBound + 1), oldL);
3553  int oldL2= l/2;
3554  int stepSize= 2;
3555  bool useOldQs= false;
3556  bool hitBound= false;
3557  CFListIterator j;
3558  CFMatrix C;
3559  mat_zz_pE* NTLC, NTLK;
3560  CFArray buf;
3561  CanonicalForm truncF;
3562  while (l <= precision)
3563  {
3564    j= factors;
3565    truncF= mod (F, power (y,l));
3566    if (useOldQs)
3567    {
3568      for (int i= 0; i < factors.length(); i++, j++)
3569        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL2, bufQ[i],
3570                                     bufQ[i]
3571                                    );
3572    }
3573    else
3574    {
3575      for (int i= 0; i < factors.length(); i++, j++)
3576        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
3577    }
3578    useOldQs= true;
3579    for (int i= 0; i < d; i++)
3580    {
3581      if (bounds [i] + 1 <= l/2)
3582      {
3583        int k= tmin (bounds [i] + 1, l/2);
3584        C= CFMatrix (l - k, factors.length());
3585        for (int ii= 0; ii < factors.length(); ii++)
3586        {
3587          if (A[ii].size() - 1 >= i)
3588          {
3589            buf= getCoeffs (A[ii] [i], k);
3590            writeInMatrix (C, buf, ii + 1, 0);
3591          }
3592        }
3593        NTLC= convertFacCFMatrix2NTLmat_zz_pE(C);
3594        NTLK= (*NTLC)*NTLN;
3595        transpose (NTLK, NTLK);
3596        kernel (NTLK, NTLK);
3597        transpose (NTLK, NTLK);
3598        NTLN *= NTLK;
3599        if (NTLN.NumCols() == 1)
3600        {
3601          delete [] A;
3602          delete [] bounds;
3603          CanonicalForm G= F;
3604          F= 1;
3605          return CFList (G (y-eval,y));
3606        }
3607      }
3608    }
3609
3610    if (NTLN.NumCols() < oldNumCols - factorsFound)
3611    {
3612      if (isReduced (NTLN))
3613      {
3614        int * factorsFoundIndex= new int [NTLN.NumCols()];
3615        for (long i= 0; i < NTLN.NumCols(); i++)
3616          factorsFoundIndex[i]= 0;
3617        int factorsFound2= 0;
3618        CFList result;
3619        CanonicalForm bufF= F;
3620        reconstructionTry (result, bufF, factors, degree (F) + 1, factorsFound2,
3621                           factorsFoundIndex, NTLN, eval, false);
3622        if (result.length() == NTLN.NumCols())
3623        {
3624          delete [] factorsFoundIndex;
3625          delete [] A;
3626          delete [] bounds;
3627          F= 1;
3628          return result;
3629        }
3630        delete [] factorsFoundIndex;
3631      }
3632      else if (l == precision)
3633      {
3634        CanonicalForm bufF= F;
3635        int * zeroOne= extractZeroOneVecs (NTLN);
3636        CFList result= reconstruction (bufF, factors, zeroOne, precision, NTLN, eval);
3637        F= bufF;
3638        delete [] zeroOne;
3639        delete [] A;
3640        delete [] bounds;
3641        return result;
3642      }
3643    }
3644    oldL2= l;
3645    l += stepSize;
3646    stepSize *= 2;
3647    if (l > precision)
3648    {
3649      if (!hitBound)
3650      {
3651        l= precision;
3652        hitBound= true;
3653      }
3654      else
3655        break;
3656    }
3657  }
3658  delete [] bounds;
3659  delete [] A;
3660  return CFList();
3661}
3662
3663//over field extension
3664CFList
3665extIncreasePrecision (CanonicalForm& F, CFList& factors, int factorsFound,
3666                      int oldNumCols, int oldL, const CanonicalForm& evaluation,
3667                      const ExtensionInfo& info, CFList& source, CFList& dest,
3668                      int precision
3669                     )
3670{
3671  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
3672  int degMipo= degree (getMipo (info.getAlpha()));
3673  Variable alpha= info.getAlpha();
3674  int d;
3675  bool isIrreducible= false;
3676  int* bounds= computeBounds (F, d, isIrreducible);
3677  if (isIrreducible)
3678  {
3679    delete [] bounds;
3680    Variable y= Variable (2);
3681    CanonicalForm tmp= F (y - evaluation, y);
3682    CFList source, dest;
3683    tmp= mapDown (tmp, info, source, dest);
3684    F= 1;
3685    return CFList (tmp);
3686  }
3687
3688  CFArray * A= new CFArray [factors.length()];
3689  CFArray bufQ= CFArray (factors.length());
3690  if (fac_NTL_char != getCharacteristic())
3691  {
3692    fac_NTL_char= getCharacteristic();
3693    zz_p::init (getCharacteristic());
3694  }
3695  mat_zz_p NTLN;
3696  ident (NTLN, factors.length());
3697  int minBound= bounds[0];
3698  for (int i= 1; i < d; i++)
3699  {
3700    if (bounds[i] != 0)
3701      minBound= tmin (minBound, bounds[i]);
3702  }
3703  int l= tmax (oldL, 2*((minBound+1)/degMipo+1));
3704  int oldL2= l/2;
3705  int stepSize= 2;
3706  bool useOldQs= false;
3707  bool hitBound= false;
3708  Variable gamma= info.getBeta();
3709  CanonicalForm primElemAlpha= info.getGamma();
3710  CanonicalForm imPrimElemAlpha= info.getDelta();
3711  CFListIterator j;
3712  Variable y= F.mvar();
3713  CanonicalForm powX, imBasis, truncF;
3714  CFMatrix Mat, C;
3715  CFIterator iter;
3716  mat_zz_p* NTLMat,*NTLC, NTLK;
3717  CFArray buf;
3718  while (l <= precision)
3719  {
3720    j= factors;
3721    if (GF)
3722      setCharacteristic (getCharacteristic());
3723    powX= power (y-gamma, l);
3724    Mat= CFMatrix (l*degMipo, l*degMipo);
3725    for (int i= 0; i < l*degMipo; i++)
3726    {
3727      imBasis= mod (power (y, i), powX);
3728      imBasis= imBasis (power (y, degMipo), y);
3729      imBasis= imBasis (y, gamma);
3730      iter= imBasis;
3731      for (; iter.hasTerms(); iter++)
3732          Mat (iter.exp()+ 1, i+1)= iter.coeff();
3733    }
3734
3735    NTLMat= convertFacCFMatrix2NTLmat_zz_p (Mat);
3736    *NTLMat= inv (*NTLMat);
3737    if (GF)
3738      setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
3739
3740    truncF= mod (F, power (y, l));
3741    if (useOldQs)
3742    {
3743      for (int i= 0; i < factors.length(); i++, j++)
3744        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL2, bufQ[i],
3745                                     bufQ[i]
3746                                    );
3747    }
3748    else
3749    {
3750      for (int i= 0; i < factors.length(); i++, j++)
3751        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
3752    }
3753    useOldQs= true;
3754    for (int i= 0; i < d; i++)
3755    {
3756      if (bounds [i] + 1 <= (l/2)*degMipo)
3757      {
3758        int k= tmin (bounds [i] + 1, (l/2)*degMipo);
3759        C= CFMatrix (l*degMipo - k, factors.length());
3760        for (int ii= 0; ii < factors.length(); ii++)
3761        {
3762          if (A[ii].size() - 1 >= i)
3763          {
3764            if (GF)
3765            {
3766              A[ii] [i]= A [ii] [i] (y-evaluation, y);
3767              setCharacteristic (getCharacteristic());
3768              A[ii] [i]= GF2FalphaRep (A[ii] [i], alpha);
3769              if (alpha != gamma)
3770                A [ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
3771                                     gamma, source, dest
3772                                    );
3773              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
3774            }
3775            else
3776            {
3777              A [ii] [i]= A [ii] [i] (y-evaluation, y);
3778              if (alpha != gamma)
3779                A[ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
3780                                    gamma, source, dest
3781                                   );
3782              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
3783            }
3784            writeInMatrix (C, buf, ii + 1, 0);
3785          }
3786          if (GF)
3787            setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
3788        }
3789
3790        if (GF)
3791          setCharacteristic(getCharacteristic());
3792
3793        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
3794        NTLK= (*NTLC)*NTLN;
3795        transpose (NTLK, NTLK);
3796        kernel (NTLK, NTLK);
3797        transpose (NTLK, NTLK);
3798        NTLN *= NTLK;
3799
3800        if (GF)
3801          setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
3802
3803        if (NTLN.NumCols() == 1)
3804        {
3805          Variable y= Variable (2);
3806          CanonicalForm tmp= F (y - evaluation, y);
3807          CFList source, dest;
3808          tmp= mapDown (tmp, info, source, dest);
3809          delete [] A;
3810          delete [] bounds;
3811          F= 1;
3812          return CFList (tmp);
3813        }
3814      }
3815    }
3816
3817    if (NTLN.NumCols() < oldNumCols - factorsFound)
3818    {
3819      if (isReduced (NTLN))
3820      {
3821        int * factorsFoundIndex= new int [NTLN.NumCols()];
3822        for (long i= 0; i < NTLN.NumCols(); i++)
3823          factorsFoundIndex[i]= 0;
3824        int factorsFound2= 0;
3825        CFList result;
3826        CanonicalForm bufF= F;
3827        extReconstructionTry (result, bufF, factors,degree (F)+1, factorsFound2,
3828                              factorsFoundIndex, NTLN, false, info, evaluation
3829                             );
3830        if (result.length() == NTLN.NumCols())
3831        {
3832          delete [] factorsFoundIndex;
3833          delete [] A;
3834          delete [] bounds;
3835          F= 1;
3836          return result;
3837        }
3838        delete [] factorsFoundIndex;
3839      }
3840      else if (l == precision)
3841      {
3842        CanonicalForm bufF= F;
3843        int * zeroOne= extractZeroOneVecs (NTLN);
3844        CFList result= extReconstruction (bufF, factors, zeroOne, precision,
3845                                          NTLN, info, evaluation
3846                                         );
3847        F= bufF;
3848        delete [] zeroOne;
3849        delete [] A;
3850        delete [] bounds;
3851        return result;
3852      }
3853    }
3854    oldL2= l;
3855    l += stepSize;
3856    stepSize *= 2;
3857    if (l > precision)
3858    {
3859      if (!hitBound)
3860      {
3861        hitBound= true;
3862        l= precision;
3863      }
3864      else
3865        break;
3866    }
3867  }
3868  delete [] bounds;
3869  delete [] A;
3870  return CFList();
3871}
3872
3873CFList
3874increasePrecision2 (const CanonicalForm& F, CFList& factors,
3875                    const Variable& alpha, int precision)
3876{
3877  int d;
3878  bool isIrreducible= false;
3879  int* bounds= computeBounds (F, d, isIrreducible);
3880  if (isIrreducible)
3881  {
3882    delete [] bounds;
3883    return CFList (F);
3884  }
3885  CFArray * A= new CFArray [factors.length()];
3886  CFArray bufQ= CFArray (factors.length());
3887  if (fac_NTL_char != getCharacteristic())
3888  {
3889    fac_NTL_char= getCharacteristic();
3890    zz_p::init (getCharacteristic());
3891  }
3892  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
3893  zz_pE::init (NTLMipo);
3894  mat_zz_pE NTLN;
3895  ident (NTLN, factors.length());
3896  int minBound= bounds[0];
3897  for (int i= 1; i < d; i++)
3898  {
3899    if (bounds[i] != 0)
3900      minBound= tmin (minBound, bounds[i]);
3901  }
3902  int l= tmin (2*(minBound + 1), precision);
3903  int oldL= l/2;
3904  int stepSize= 2;
3905  bool useOldQs= false;
3906  bool hitBound= false;
3907  CFListIterator j;
3908  CFMatrix C;
3909  CFArray buf;
3910  mat_zz_pE* NTLC, NTLK;
3911  Variable y= F.mvar();
3912  CanonicalForm truncF;
3913  while (l <= precision)
3914  {
3915    j= factors;
3916    truncF= mod (F, power (y, l));
3917    if (useOldQs)
3918    {
3919      for (int i= 0; i < factors.length(); i++, j++)
3920        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i], bufQ[i]);
3921    }
3922    else
3923    {
3924      for (int i= 0; i < factors.length(); i++, j++)
3925        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
3926    }
3927    useOldQs= true;
3928    for (int i= 0; i < d; i++)
3929    {
3930      if (bounds [i] + 1 <= l/2)
3931      {
3932        int k= tmin (bounds [i] + 1, l/2);
3933        C= CFMatrix (l - k, factors.length());
3934        for (int ii= 0; ii < factors.length(); ii++)
3935        {
3936          if (A[ii].size() - 1 >= i)
3937          {
3938            buf= getCoeffs (A[ii] [i], k);
3939            writeInMatrix (C, buf, ii + 1, 0);
3940          }
3941        }
3942        NTLC= convertFacCFMatrix2NTLmat_zz_pE(C);
3943        NTLK= (*NTLC)*NTLN;
3944        transpose (NTLK, NTLK);
3945        kernel (NTLK, NTLK);
3946        transpose (NTLK, NTLK);
3947        NTLN *= NTLK;
3948        if (NTLN.NumCols() == 1)
3949        {
3950          delete [] A;
3951          delete [] bounds;
3952          return CFList (F);
3953        }
3954      }
3955    }
3956
3957    if (isReduced (NTLN) || l == precision)
3958    {
3959      CanonicalForm bufF= F;
3960      int * zeroOne= extractZeroOneVecs (NTLN);
3961      CFList bufFactors= factors;
3962      CFList result= monicReconstruction (bufF, factors, zeroOne, precision,
3963                                          NTLN
3964                                         );
3965      if (result.length() != NTLN.NumCols() && l != precision)
3966        factors= bufFactors;
3967      if (result.length() == NTLN.NumCols())
3968      {
3969        delete [] zeroOne;
3970        delete [] A;
3971        delete [] bounds;
3972        return result;
3973      }
3974      if (l == precision)
3975      {
3976        delete [] zeroOne;
3977        delete [] A;
3978        delete [] bounds;
3979        return Union (result, factors);
3980      }
3981    }
3982    oldL= l;
3983    l += stepSize;
3984    stepSize *= 2;
3985    if (l > precision)
3986    {
3987      if (!hitBound)
3988      {
3989        l= precision;
3990        hitBound= true;
3991      }
3992      else
3993        break;
3994    }
3995  }
3996  delete [] bounds;
3997  delete [] A;
3998  return CFList();
3999}
4000
4001CFList
4002increasePrecisionFq2Fp (CanonicalForm& F, CFList& factors, int factorsFound,
4003                        int oldNumCols, int oldL, const Variable& alpha,
4004                        int precision, const CanonicalForm& eval
4005                       )
4006{
4007  int d;
4008  bool isIrreducible= false;
4009  Variable y= F.mvar();
4010  int* bounds= computeBounds (F, d, isIrreducible);
4011  if (isIrreducible)
4012  {
4013    delete [] bounds;
4014    CanonicalForm G= F;
4015    F= 1;
4016    return CFList (G (y-eval,y));
4017  }
4018  int extensionDeg= degree (getMipo (alpha));
4019  CFArray * A= new CFArray [factors.length()];
4020  CFArray bufQ= CFArray (factors.length());
4021#ifdef HAVE_FLINT
4022  nmod_mat_t FLINTN;
4023  nmod_mat_init (FLINTN,factors.length(),factors.length(), getCharacteristic());
4024  for (long i=factors.length()-1; i >= 0; i--)
4025    nmod_mat_entry (FLINTN, i, i)= 1;
4026#else
4027  mat_zz_p NTLN;
4028  ident (NTLN, factors.length());
4029#endif
4030  int minBound= bounds[0];
4031  for (int i= 1; i < d; i++)
4032  {
4033    if (bounds[i] != 0)
4034      minBound= tmin (minBound, bounds[i]);
4035  }
4036  int l= tmax (2*(minBound + 1), oldL);
4037  int oldL2= l/2;
4038  int stepSize= 2;
4039  bool useOldQs= false;
4040  bool hitBound= false;
4041  CFListIterator j;
4042  CFMatrix C;
4043#ifdef HAVE_FLINT
4044  long rank;
4045  nmod_mat_t FLINTC, FLINTK, null;
4046#else
4047  mat_zz_p* NTLC, NTLK;
4048#endif
4049  CFArray buf;
4050  CanonicalForm truncF;
4051  while (l <= precision)
4052  {
4053    j= factors;
4054    truncF= mod (F, power (y, l));
4055    if (useOldQs)
4056    {
4057      for (int i= 0; i < factors.length(); i++, j++)
4058        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL2, bufQ[i],
4059                                     bufQ[i]
4060                                    );
4061    }
4062    else
4063    {
4064      for (int i= 0; i < factors.length(); i++, j++)
4065        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
4066    }
4067    useOldQs= true;
4068    for (int i= 0; i < d; i++)
4069    {
4070      if (bounds [i] + 1 <= l/2)
4071      {
4072        int k= tmin (bounds [i] + 1, l/2);
4073        C= CFMatrix ((l - k)*extensionDeg, factors.length());
4074        for (int ii= 0; ii < factors.length(); ii++)
4075        {
4076          if (A[ii].size() - 1 >= i)
4077          {
4078            buf= getCoeffs (A[ii] [i], k, alpha);
4079            writeInMatrix (C, buf, ii + 1, 0);
4080          }
4081        }
4082#ifdef HAVE_FLINT
4083        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
4084        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
4085                       getCharacteristic());
4086        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
4087        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
4088                       getCharacteristic());
4089        rank= nmod_mat_nullspace (null, FLINTK);
4090        nmod_mat_clear (FLINTK);
4091        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
4092        nmod_mat_clear (FLINTC);
4093        nmod_mat_init_set (FLINTC, FLINTN);
4094        nmod_mat_clear (FLINTN);
4095        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
4096                       getCharacteristic());
4097        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
4098
4099        nmod_mat_clear (FLINTC);
4100        nmod_mat_window_clear (FLINTK);
4101        nmod_mat_clear (null);
4102#else
4103        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
4104        NTLK= (*NTLC)*NTLN;
4105        transpose (NTLK, NTLK);
4106        kernel (NTLK, NTLK);
4107        transpose (NTLK, NTLK);
4108        NTLN *= NTLK;
4109#endif
4110#ifdef HAVE_FLINT
4111        if (nmod_mat_ncols (FLINTN) == 1)
4112        {
4113          nmod_mat_clear (FLINTN);
4114#else
4115        if (NTLN.NumCols() == 1)
4116        {
4117#endif
4118          delete [] A;
4119          delete [] bounds;
4120          CanonicalForm G= F;
4121          F= 1;
4122          return CFList (G (y-eval,y));
4123        }
4124      }
4125    }
4126
4127#ifdef HAVE_FLINT
4128    if (nmod_mat_ncols (FLINTN) < oldNumCols - factorsFound)
4129    {
4130      if (isReduced (FLINTN))
4131      {
4132        int * factorsFoundIndex= new int [nmod_mat_ncols (FLINTN)];
4133        for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
4134#else
4135    if (NTLN.NumCols() < oldNumCols - factorsFound)
4136    {
4137      if (isReduced (NTLN))
4138      {
4139        int * factorsFoundIndex= new int [NTLN.NumCols()];
4140        for (long i= 0; i < NTLN.NumCols(); i++)
4141#endif
4142          factorsFoundIndex[i]= 0;
4143        int factorsFound2= 0;
4144        CFList result;
4145        CanonicalForm bufF= F;
4146#ifdef HAVE_FLINT
4147        reconstructionTry (result, bufF, factors, degree (F) + 1, factorsFound2,
4148                           factorsFoundIndex, FLINTN, eval, false
4149                          );
4150        if (result.length() == nmod_mat_ncols (FLINTN))
4151        {
4152          nmod_mat_clear (FLINTN);
4153#else
4154        reconstructionTry (result, bufF, factors, degree (F) + 1, factorsFound2,
4155                           factorsFoundIndex, NTLN, eval, false
4156                          );
4157        if (result.length() == NTLN.NumCols())
4158        {
4159#endif
4160          delete [] factorsFoundIndex;
4161          delete [] A;
4162          delete [] bounds;
4163          F= 1;
4164          return result;
4165        }
4166        delete [] factorsFoundIndex;
4167      }
4168      else if (l == precision)
4169      {
4170        CanonicalForm bufF= F;
4171#ifdef HAVE_FLINT
4172        int * zeroOne= extractZeroOneVecs (FLINTN);
4173        CFList result= reconstruction (bufF,factors,zeroOne,precision,FLINTN, eval);
4174        nmod_mat_clear (FLINTN);
4175#else
4176        int * zeroOne= extractZeroOneVecs (NTLN);
4177        CFList result= reconstruction (bufF, factors, zeroOne, precision, NTLN, eval);
4178#endif
4179        F= bufF;
4180        delete [] zeroOne;
4181        delete [] A;
4182        delete [] bounds;
4183        return result;
4184      }
4185    }
4186    oldL2= l;
4187    l += stepSize;
4188    stepSize *= 2;
4189    if (l > precision)
4190    {
4191      if (!hitBound)
4192      {
4193        hitBound= true;
4194        l= precision;
4195      }
4196      else
4197        break;
4198    }
4199  }
4200#ifdef HAVE_FLINT
4201  nmod_mat_clear (FLINTN);
4202#endif
4203  delete [] bounds;
4204  delete [] A;
4205  return CFList();
4206}
4207
4208#ifdef HAVE_FLINT
4209CFList
4210increasePrecision (CanonicalForm& F, CFList& factors, int oldL, int
4211                   l, int d, int* bounds, CFArray& bufQ, nmod_mat_t FLINTN,
4212                   const CanonicalForm& eval
4213                  )
4214#else
4215CFList
4216increasePrecision (CanonicalForm& F, CFList& factors, int oldL, int
4217                   l, int d, int* bounds, CFArray& bufQ, mat_zz_p& NTLN,
4218                   const CanonicalForm& eval
4219                  )
4220#endif
4221{
4222  CFList result= CFList();
4223  CFArray * A= new CFArray [factors.length()];
4224  int oldL2= oldL/2;
4225  bool hitBound= false;
4226#ifdef HAVE_FLINT
4227  if (nmod_mat_nrows (FLINTN) != factors.length()) //refined factors
4228  {
4229    nmod_mat_clear (FLINTN);
4230    nmod_mat_init(FLINTN,factors.length(),factors.length(),getCharacteristic());
4231    for (long i=factors.length()-1; i >= 0; i--)
4232      nmod_mat_entry (FLINTN, i, i)= 1;
4233    bufQ= CFArray (factors.length());
4234  }
4235#else
4236  if (NTLN.NumRows() != factors.length()) //refined factors
4237  {
4238    ident (NTLN, factors.length());
4239    bufQ= CFArray (factors.length());
4240  }
4241#endif
4242  bool useOldQs= false;
4243  CFListIterator j;
4244  CFMatrix C;
4245  CFArray buf;
4246#ifdef HAVE_FLINT
4247  long rank;
4248  nmod_mat_t FLINTC, FLINTK, null;
4249#else
4250  mat_zz_p* NTLC, NTLK;
4251#endif
4252  CanonicalForm bufF, truncF;
4253  CFList bufUniFactors;
4254  Variable y= F.mvar();
4255  while (oldL <= l)
4256  {
4257    j= factors;
4258    truncF= mod (F, power (y, oldL));
4259    if (useOldQs)
4260    {
4261      for (int i= 0; i < factors.length(); i++, j++)
4262        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, oldL2, bufQ[i],
4263                                     bufQ[i]
4264                                    );
4265    }
4266    else
4267    {
4268      for (int i= 0; i < factors.length(); i++, j++)
4269        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, bufQ [i]);
4270    }
4271    useOldQs= true;
4272
4273    for (int i= 0; i < d; i++)
4274    {
4275      if (bounds [i] + 1 <= oldL/2)
4276      {
4277        int k= tmin (bounds [i] + 1, oldL/2);
4278        C= CFMatrix (oldL - k, factors.length());
4279        for (int ii= 0; ii < factors.length(); ii++)
4280        {
4281          if (A[ii].size() - 1 >= i)
4282          {
4283            buf= getCoeffs (A[ii] [i], k);
4284            writeInMatrix (C, buf, ii + 1, 0);
4285          }
4286        }
4287#ifdef HAVE_FLINT
4288        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
4289        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
4290                       getCharacteristic());
4291        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
4292        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
4293                       getCharacteristic());
4294        rank= nmod_mat_nullspace (null, FLINTK);
4295        nmod_mat_clear (FLINTK);
4296        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
4297        nmod_mat_clear (FLINTC);
4298        nmod_mat_init_set (FLINTC, FLINTN);
4299        nmod_mat_clear (FLINTN);
4300        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
4301                       getCharacteristic());
4302        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
4303
4304        nmod_mat_clear (FLINTC);
4305        nmod_mat_window_clear (FLINTK);
4306        nmod_mat_clear (null);
4307#else
4308        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
4309        NTLK= (*NTLC)*NTLN;
4310        transpose (NTLK, NTLK);
4311        kernel (NTLK, NTLK);
4312        transpose (NTLK, NTLK);
4313        NTLN *= NTLK;
4314#endif
4315#ifdef HAVE_FLINT
4316        if (nmod_mat_ncols (FLINTN) == 1)
4317#else
4318        if (NTLN.NumCols() == 1)
4319#endif
4320        {
4321          delete [] A;
4322          return CFList (F (y-eval,y));
4323        }
4324      }
4325    }
4326#ifdef HAVE_FLINT
4327    if (nmod_mat_ncols (FLINTN) == 1)
4328#else
4329    if (NTLN.NumCols() == 1)
4330#endif
4331    {
4332      delete [] A;
4333      return CFList (F (y-eval,y));
4334    }
4335    int * zeroOneVecs;
4336#ifdef HAVE_FLINT
4337    zeroOneVecs= extractZeroOneVecs (FLINTN);
4338#else
4339    zeroOneVecs= extractZeroOneVecs (NTLN);
4340#endif
4341    bufF= F;
4342    bufUniFactors= factors;
4343#ifdef HAVE_FLINT
4344    result= reconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, FLINTN, eval);
4345#else
4346    result= reconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, NTLN, eval);
4347#endif
4348    delete [] zeroOneVecs;
4349    if (degree (bufF) + 1 + degree (LC (bufF, 1)) < oldL && result.length() > 0)
4350    {
4351      F= bufF;
4352      factors= bufUniFactors;
4353      delete [] A;
4354      return result;
4355    }
4356
4357    result= CFList();
4358    oldL2= oldL;
4359    oldL *= 2;
4360    if (oldL > l)
4361    {
4362      if (!hitBound)
4363      {
4364        oldL= l;
4365        hitBound= true;
4366      }
4367      else
4368        break;
4369    }
4370  }
4371  delete [] A;
4372  return result;
4373}
4374
4375CFList
4376increasePrecision (CanonicalForm& F, CFList& factors, int oldL, int
4377                   l, int d, int* bounds, CFArray& bufQ, mat_zz_pE& NTLN,
4378                   const CanonicalForm& eval
4379                  )
4380{
4381  CFList result= CFList();
4382  CFArray * A= new CFArray [factors.length()];
4383  int oldL2= oldL/2;
4384  bool hitBound= false;
4385  bool useOldQs= false;
4386  if (NTLN.NumRows() != factors.length()) //refined factors
4387    ident (NTLN, factors.length());
4388  CFListIterator j;
4389  CFMatrix C;
4390  CFArray buf;
4391  mat_zz_pE* NTLC, NTLK;
4392  CanonicalForm bufF, truncF;
4393  CFList bufUniFactors;
4394  Variable y= F.mvar();
4395  while (oldL <= l)
4396  {
4397    j= factors;
4398    truncF= mod (F, power (y, oldL));
4399    if (useOldQs)
4400    {
4401      for (int i= 0; i < factors.length(); i++, j++)
4402        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, oldL2, bufQ[i],
4403                                     bufQ[i]
4404                                    );
4405    }
4406    else
4407    {
4408      for (int i= 0; i < factors.length(); i++, j++)
4409        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, bufQ [i]);
4410    }
4411    useOldQs= true;
4412
4413    for (int i= 0; i < d; i++)
4414    {
4415      if (bounds [i] + 1 <= oldL/2)
4416      {
4417        int k= tmin (bounds [i] + 1, oldL/2);
4418        C= CFMatrix (oldL - k, factors.length());
4419        for (int ii= 0; ii < factors.length(); ii++)
4420        {
4421          if (A[ii].size() - 1 >= i)
4422          {
4423            buf= getCoeffs (A[ii] [i], k);
4424            writeInMatrix (C, buf, ii + 1, 0);
4425          }
4426        }
4427        NTLC= convertFacCFMatrix2NTLmat_zz_pE(C);
4428        NTLK= (*NTLC)*NTLN;
4429        transpose (NTLK, NTLK);
4430        kernel (NTLK, NTLK);
4431        transpose (NTLK, NTLK);
4432        NTLN *= NTLK;
4433        if (NTLN.NumCols() == 1)
4434        {
4435          delete [] A;
4436          return CFList (F (y-eval,y));
4437        }
4438      }
4439    }
4440    if (NTLN.NumCols() == 1)
4441    {
4442      delete [] A;
4443      return CFList (F (y-eval,y));
4444    }
4445
4446    int * zeroOneVecs;
4447    zeroOneVecs= extractZeroOneVecs (NTLN);
4448    bufF= F;
4449    bufUniFactors= factors;
4450    result= reconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, NTLN, eval);
4451    delete [] zeroOneVecs;
4452    if (degree (bufF) + 1 + degree (LC (bufF, 1)) < l && result.length() > 0)
4453    {
4454      F= bufF;
4455      factors= bufUniFactors;
4456      delete [] A;
4457      return result;
4458    }
4459
4460    result= CFList();
4461    oldL2= oldL;
4462    oldL *= 2;
4463    if (oldL > l)
4464    {
4465      if (!hitBound)
4466      {
4467        oldL= l;
4468        hitBound= true;
4469      }
4470      else
4471        break;
4472    }
4473  }
4474  delete [] A;
4475  return result;
4476}
4477
4478//over field extension
4479CFList
4480extIncreasePrecision (CanonicalForm& F, CFList& factors, int oldL, int l, int d,
4481                      int* bounds, CFArray& bufQ, mat_zz_p& NTLN, const
4482                      CanonicalForm& evaluation, const ExtensionInfo& info,
4483                      CFList& source, CFList& dest
4484                     )
4485{
4486  CFList result= CFList();
4487  CFArray * A= new CFArray [factors.length()];
4488  int oldL2= oldL/2; //be careful
4489  bool hitBound= false;
4490  bool useOldQs= false;
4491  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
4492  int degMipo= degree (getMipo (info.getAlpha()));
4493  Variable alpha= info.getAlpha();
4494
4495  Variable gamma= info.getBeta();
4496  CanonicalForm primElemAlpha= info.getGamma();
4497  CanonicalForm imPrimElemAlpha= info.getDelta();
4498  if (NTLN.NumRows() != factors.length()) //refined factors
4499    ident (NTLN, factors.length());
4500  Variable y= F.mvar();
4501  CFListIterator j;
4502  CanonicalForm powX, imBasis, bufF, truncF;
4503  CFMatrix Mat, C;
4504  CFIterator iter;
4505  mat_zz_p* NTLMat;
4506  CFArray buf;
4507  mat_zz_p* NTLC, NTLK;
4508  CFList bufUniFactors;
4509  while (oldL <= l)
4510  {
4511    j= factors;
4512    if (GF)
4513      setCharacteristic (getCharacteristic());
4514
4515    powX= power (y-gamma, oldL);
4516    Mat= CFMatrix (oldL*degMipo, oldL*degMipo);
4517    for (int i= 0; i < oldL*degMipo; i++)
4518    {
4519      imBasis= mod (power (y, i), powX);
4520      imBasis= imBasis (power (y, degMipo), y);
4521      imBasis= imBasis (y, gamma);
4522      iter= imBasis;
4523      for (; iter.hasTerms(); iter++)
4524        Mat (iter.exp()+ 1, i+1)= iter.coeff();
4525    }
4526
4527    NTLMat= convertFacCFMatrix2NTLmat_zz_p (Mat);
4528    *NTLMat= inv (*NTLMat);
4529    if (GF)
4530      setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
4531
4532    truncF= mod (F, power (y, oldL));
4533    if (useOldQs)
4534    {
4535      for (int i= 0; i < factors.length(); i++, j++)
4536        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, oldL2, bufQ[i],
4537                                     bufQ[i]);
4538    }
4539    else
4540    {
4541      for (int i= 0; i < factors.length(); i++, j++)
4542        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, bufQ [i]);
4543    }
4544    useOldQs= true;
4545
4546    for (int i= 0; i < d; i++)
4547    {
4548      if (bounds [i] + 1 <= oldL/2)
4549      {
4550        int k= tmin (bounds [i] + 1, oldL/2);
4551        C= CFMatrix (oldL*degMipo - k, factors.length());
4552        for (int ii= 0; ii < factors.length(); ii++)
4553        {
4554          if (A[ii].size() - 1 >= i)
4555          {
4556            if (GF)
4557            {
4558              A [ii] [i]= A [ii] [i] (y-evaluation, y);
4559              setCharacteristic (getCharacteristic());
4560              A[ii] [i]= GF2FalphaRep (A[ii] [i], alpha);
4561              if (alpha != gamma)
4562                A [ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
4563                                     gamma, source, dest
4564                                    );
4565              buf= getCoeffs (A[ii] [i], k, oldL, degMipo, gamma, 0, *NTLMat);
4566            }
4567            else
4568            {
4569              A [ii] [i]= A [ii] [i] (y-evaluation, y);
4570              if (alpha != gamma)
4571                A[ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
4572                                    gamma, source, dest
4573                                   );
4574              buf= getCoeffs (A[ii] [i], k, oldL, degMipo, gamma, 0, *NTLMat);
4575            }
4576            writeInMatrix (C, buf, ii + 1, 0);
4577          }
4578          if (GF)
4579            setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
4580        }
4581
4582        if (GF)
4583          setCharacteristic(getCharacteristic());
4584
4585        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
4586        NTLK= (*NTLC)*NTLN;
4587        transpose (NTLK, NTLK);
4588        kernel (NTLK, NTLK);
4589        transpose (NTLK, NTLK);
4590        NTLN *= NTLK;
4591
4592        if (GF)
4593          setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
4594
4595        if (NTLN.NumCols() == 1)
4596        {
4597          Variable y= Variable (2);
4598          CanonicalForm tmp= F (y - evaluation, y);
4599          CFList source, dest;
4600          tmp= mapDown (tmp, info, source, dest);
4601          delete [] A;
4602          return CFList (tmp);
4603        }
4604      }
4605    }
4606    if (NTLN.NumCols() == 1)
4607    {
4608      Variable y= Variable (2);
4609      CanonicalForm tmp= F (y - evaluation, y);
4610      CFList source, dest;
4611      tmp= mapDown (tmp, info, source, dest);
4612      delete [] A;
4613      return CFList (tmp);
4614    }
4615
4616    int * zeroOneVecs;
4617    zeroOneVecs= extractZeroOneVecs (NTLN);
4618    bufF= F;
4619    bufUniFactors= factors;
4620    result= extReconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, NTLN,
4621                               info, evaluation
4622                              );
4623    delete [] zeroOneVecs;
4624    if (degree (bufF) + 1 + degree (LC (bufF, 1)) < l && result.length() > 0)
4625    {
4626      F= bufF;
4627      factors= bufUniFactors;
4628      return result;
4629    }
4630
4631    result= CFList();
4632    oldL2= oldL;
4633    oldL *= 2;
4634    if (oldL > l)
4635    {
4636      if (!hitBound)
4637      {
4638        oldL= l;
4639        hitBound= true;
4640      }
4641      else
4642        break;
4643    }
4644  }
4645  delete [] A;
4646  return result;
4647}
4648
4649#ifdef HAVE_FLINT
4650CFList
4651increasePrecisionFq2Fp (CanonicalForm& F, CFList& factors, int oldL, int l,
4652                        int d, int* bounds, CFArray& bufQ, nmod_mat_t FLINTN,
4653                        const Variable& alpha, const CanonicalForm& eval
4654                       )
4655#else
4656CFList
4657increasePrecisionFq2Fp (CanonicalForm& F, CFList& factors, int oldL, int l,
4658                        int d, int* bounds, CFArray& bufQ, mat_zz_p& NTLN,
4659                        const Variable& alpha, const CanonicalForm& eval
4660                       )
4661#endif
4662{
4663  CFList result= CFList();
4664  CFArray * A= new CFArray [factors.length()];
4665  int extensionDeg= degree (getMipo (alpha));
4666  int oldL2= oldL/2;
4667  bool hitBound= false;
4668  bool useOldQs= false;
4669#ifdef HAVE_FLINT
4670  if (nmod_mat_nrows (FLINTN) != factors.length()) //refined factors
4671  {
4672    nmod_mat_clear (FLINTN);
4673    nmod_mat_init(FLINTN,factors.length(),factors.length(),getCharacteristic());
4674    for (long i=factors.length()-1; i >= 0; i--)
4675      nmod_mat_entry (FLINTN, i, i)= 1;
4676  }
4677#else
4678  if (NTLN.NumRows() != factors.length()) //refined factors
4679    ident (NTLN, factors.length());
4680#endif
4681  CFListIterator j;
4682  CFMatrix C;
4683  CFArray buf;
4684#ifdef HAVE_FLINT
4685  long rank;
4686  nmod_mat_t FLINTC, FLINTK, null;
4687#else
4688  mat_zz_p* NTLC, NTLK;
4689#endif
4690  CanonicalForm bufF, truncF;
4691  CFList bufUniFactors;
4692  Variable y= F.mvar();
4693  while (oldL <= l)
4694  {
4695    j= factors;
4696    truncF= mod (F, power (y, oldL));
4697    if (useOldQs)
4698    {
4699      for (int i= 0; i < factors.length(); i++, j++)
4700        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, oldL2, bufQ[i],
4701                                     bufQ[i]
4702                                    );
4703    }
4704    else
4705    {
4706      for (int i= 0; i < factors.length(); i++, j++)
4707        A[i]= logarithmicDerivative (truncF, j.getItem(), oldL, bufQ [i]);
4708    }
4709    useOldQs= true;
4710
4711    for (int i= 0; i < d; i++)
4712    {
4713      if (bounds [i] + 1 <= oldL/2)
4714      {
4715        int k= tmin (bounds [i] + 1, oldL/2);
4716        C= CFMatrix ((oldL - k)*extensionDeg, factors.length());
4717        for (int ii= 0; ii < factors.length(); ii++)
4718        {
4719          if (A[ii].size() - 1 >= i)
4720          {
4721            buf= getCoeffs (A[ii] [i], k, alpha);
4722            writeInMatrix (C, buf, ii + 1, 0);
4723          }
4724        }
4725#ifdef HAVE_FLINT
4726        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
4727        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
4728                       getCharacteristic());
4729        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
4730        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
4731                       getCharacteristic());
4732        rank= nmod_mat_nullspace (null, FLINTK);
4733        nmod_mat_clear (FLINTK);
4734        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
4735        nmod_mat_clear (FLINTC);
4736        nmod_mat_init_set (FLINTC, FLINTN);
4737        nmod_mat_clear (FLINTN);
4738        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
4739                       getCharacteristic());
4740        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
4741
4742        nmod_mat_clear (FLINTC);
4743        nmod_mat_window_clear (FLINTK);
4744        nmod_mat_clear (null);
4745#else
4746        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
4747        NTLK= (*NTLC)*NTLN;
4748        transpose (NTLK, NTLK);
4749        kernel (NTLK, NTLK);
4750        transpose (NTLK, NTLK);
4751        NTLN *= NTLK;
4752#endif
4753#ifdef HAVE_FLINT
4754        if (nmod_mat_ncols (FLINTN) == 1)
4755#else
4756        if (NTLN.NumCols() == 1)
4757#endif
4758        {
4759          delete [] A;
4760          return CFList (F(y-eval,y));
4761        }
4762      }
4763    }
4764
4765    int * zeroOneVecs;
4766#ifdef HAVE_FLINT
4767    zeroOneVecs= extractZeroOneVecs (FLINTN);
4768#else
4769    zeroOneVecs= extractZeroOneVecs (NTLN);
4770#endif
4771
4772    bufF= F;
4773    bufUniFactors= factors;
4774#ifdef HAVE_FLINT
4775    result= reconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, FLINTN, eval);
4776#else
4777    result= reconstruction (bufF, bufUniFactors, zeroOneVecs, oldL, NTLN, eval);
4778#endif
4779    delete [] zeroOneVecs;
4780    if (degree (bufF) + 1 + degree (LC (bufF, 1)) < l && result.length() > 0)
4781    {
4782      F= bufF;
4783      factors= bufUniFactors;
4784      delete [] A;
4785      return result;
4786    }
4787
4788    result= CFList();
4789    oldL2= oldL;
4790    oldL *= 2;
4791    if (oldL > l)
4792    {
4793      if (!hitBound)
4794      {
4795        oldL= l;
4796        hitBound= true;
4797      }
4798      else
4799        break;
4800    }
4801  }
4802  delete [] A;
4803  return result;
4804}
4805
4806#ifdef HAVE_FLINT
4807CFList
4808furtherLiftingAndIncreasePrecision (CanonicalForm& F, CFList&
4809                                    factors, int l, int liftBound, int d, int*
4810                                    bounds, nmod_mat_t FLINTN, CFList& diophant,
4811                                    CFMatrix& M, CFArray& Pi, CFArray& bufQ,
4812                                    const CanonicalForm& eval
4813                                   )
4814#else
4815CFList
4816furtherLiftingAndIncreasePrecision (CanonicalForm& F, CFList&
4817                                    factors, int l, int liftBound, int d, int*
4818                                    bounds, mat_zz_p& NTLN, CFList& diophant,
4819                                    CFMatrix& M, CFArray& Pi, CFArray& bufQ,
4820                                    const CanonicalForm& eval
4821                                   )
4822#endif
4823{
4824  CanonicalForm LCF= LC (F, 1);
4825  CFList result;
4826  bool irreducible= false;
4827  CFList bufFactors= factors;
4828  CFArray *A = new CFArray [bufFactors.length()];
4829  bool useOldQs= false;
4830  bool hitBound= false;
4831  int oldL= l;
4832  int stepSize= 8; //TODO choose better step size?
4833  l += tmax (tmin (8, degree (F) + 1 + degree (LC (F, 1))-l), 2);
4834#ifdef HAVE_FLINT
4835  if (nmod_mat_nrows (FLINTN) != factors.length()) //refined factors
4836  {
4837    nmod_mat_clear (FLINTN);
4838    nmod_mat_init(FLINTN,factors.length(),factors.length(),getCharacteristic());
4839    for (long i=factors.length()-1; i >= 0; i--)
4840      nmod_mat_entry (FLINTN, i, i)= 1;
4841  }
4842#else
4843  if (NTLN.NumRows() != factors.length()) //refined factors
4844    ident (NTLN, factors.length());
4845#endif
4846  CFListIterator j;
4847  CFMatrix C;
4848  CFArray buf;
4849#ifdef HAVE_FLINT
4850  long rank;
4851  nmod_mat_t FLINTC, FLINTK, null;
4852#else
4853  mat_zz_p* NTLC, NTLK;
4854#endif
4855  CanonicalForm bufF, truncF;
4856  Variable y= F.mvar();
4857  while (l <= liftBound)
4858  {
4859    bufFactors.insert (LCF);
4860    henselLiftResume12 (F, bufFactors, oldL, l, Pi, diophant, M);
4861    bufFactors.insert (LCF);
4862    bufFactors.removeFirst();
4863    j= bufFactors;
4864    truncF= mod (F, power (y, l));
4865    if (useOldQs)
4866    {
4867      for (int i= 0; i < bufFactors.length(); i++, j++)
4868        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
4869                                     bufQ[i]);
4870    }
4871    else
4872    {
4873      for (int i= 0; i < bufFactors.length(); i++, j++)
4874        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
4875    }
4876    for (int i= 0; i < d; i++)
4877    {
4878      if (bounds [i] + 1 <= l/2)
4879      {
4880        int k= tmin (bounds [i] + 1, l/2);
4881        C= CFMatrix (l - k, bufFactors.length());
4882        for (int ii= 0; ii < bufFactors.length(); ii++)
4883        {
4884          if (A[ii].size() - 1 >= i)
4885          {
4886            buf= getCoeffs (A[ii] [i], k);
4887            writeInMatrix (C, buf, ii + 1, 0);
4888          }
4889        }
4890#ifdef HAVE_FLINT
4891        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
4892        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
4893                       getCharacteristic());
4894        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
4895        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
4896                       getCharacteristic());
4897        rank= nmod_mat_nullspace (null, FLINTK);
4898        nmod_mat_clear (FLINTK);
4899        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
4900        nmod_mat_clear (FLINTC);
4901        nmod_mat_init_set (FLINTC, FLINTN);
4902        nmod_mat_clear (FLINTN);
4903        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
4904                       getCharacteristic());
4905        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
4906
4907        nmod_mat_clear (FLINTC);
4908        nmod_mat_window_clear (FLINTK);
4909        nmod_mat_clear (null);
4910#else
4911        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
4912        NTLK= (*NTLC)*NTLN;
4913        transpose (NTLK, NTLK);
4914        kernel (NTLK, NTLK);
4915        transpose (NTLK, NTLK);
4916        NTLN *= NTLK;
4917#endif
4918#ifdef HAVE_FLINT
4919        if (nmod_mat_ncols (FLINTN) == 1)
4920#else
4921        if (NTLN.NumCols() == 1)
4922#endif
4923        {
4924          irreducible= true;
4925          break;
4926        }
4927      }
4928    }
4929
4930#ifdef HAVE_FLINT
4931    if (nmod_mat_ncols (FLINTN) == 1)
4932#else
4933    if (NTLN.NumCols() == 1)
4934#endif
4935    {
4936      irreducible= true;
4937      break;
4938    }
4939
4940#ifdef HAVE_FLINT
4941    int * zeroOneVecs= extractZeroOneVecs (FLINTN);
4942#else
4943    int * zeroOneVecs= extractZeroOneVecs (NTLN);
4944#endif
4945    bufF= F;
4946#ifdef HAVE_FLINT
4947    result= reconstruction (bufF, bufFactors, zeroOneVecs, l, FLINTN, eval);
4948#else
4949    result= reconstruction (bufF, bufFactors, zeroOneVecs, l, NTLN, eval);
4950#endif
4951    delete [] zeroOneVecs;
4952    if (result.length() > 0 && degree (bufF) + 1 + degree (LC (bufF, 1)) <= l)
4953    {
4954      F= bufF;
4955      factors= bufFactors;
4956      delete [] A;
4957      return result;
4958    }
4959
4960#ifdef HAVE_FLINT
4961    if (isReduced (FLINTN))
4962#else
4963    if (isReduced (NTLN))
4964#endif
4965    {
4966      int factorsFound= 0;
4967      bufF= F;
4968#ifdef HAVE_FLINT
4969      int* factorsFoundIndex= new int [nmod_mat_ncols (FLINTN)];
4970      for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
4971#else
4972      int* factorsFoundIndex= new int [NTLN.NumCols()];
4973      for (long i= 0; i < NTLN.NumCols(); i++)
4974#endif
4975        factorsFoundIndex[i]= 0;
4976#ifdef HAVE_FLINT
4977      if (l < liftBound)
4978        reconstructionTry (result, bufF, bufFactors, l, factorsFound,
4979                           factorsFoundIndex, FLINTN, eval, false
4980                          );
4981      else
4982        reconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
4983                           degree (LCF), factorsFound, factorsFoundIndex,
4984                           FLINTN, eval, false
4985                          );
4986
4987      if (nmod_mat_ncols (FLINTN) == result.length())
4988#else
4989      if (l < liftBound)
4990        reconstructionTry (result, bufF, bufFactors, l, factorsFound,
4991                           factorsFoundIndex, NTLN, eval, false
4992                          );
4993      else
4994        reconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
4995                           degree (LCF), factorsFound, factorsFoundIndex,
4996                           NTLN, eval, false
4997                          );
4998
4999      if (NTLN.NumCols() == result.length())
5000#endif
5001      {
5002        delete [] A;
5003        delete [] factorsFoundIndex;
5004        return result;
5005      }
5006      delete [] factorsFoundIndex;
5007    }
5008    result= CFList();
5009    oldL= l;
5010    stepSize *= 2;
5011    l += stepSize;
5012    if (l > liftBound)
5013    {
5014      if (!hitBound)
5015      {
5016        l= liftBound;
5017        hitBound= true;
5018      }
5019      else
5020        break;
5021    }
5022  }
5023  if (irreducible)
5024  {
5025    delete [] A;
5026    return CFList (F (y-eval,y));
5027  }
5028  delete [] A;
5029  factors= bufFactors;
5030  return CFList();
5031}
5032
5033//Fq
5034CFList
5035furtherLiftingAndIncreasePrecision (CanonicalForm& F, CFList&
5036                                    factors, int l, int liftBound, int d, int*
5037                                    bounds, mat_zz_pE& NTLN, CFList& diophant,
5038                                    CFMatrix& M, CFArray& Pi, CFArray& bufQ,
5039                                    const CanonicalForm& eval
5040                                   )
5041{
5042  CanonicalForm LCF= LC (F, 1);
5043  CFList result;
5044  bool irreducible= false;
5045  CFList bufFactors= factors;
5046  CFArray *A = new CFArray [bufFactors.length()];
5047  bool useOldQs= false;
5048  bool hitBound= false;
5049  int oldL= l;
5050  int stepSize= 8; //TODO choose better step size?
5051  l += tmax (tmin (8, degree (F) + 1 + degree (LC (F, 1))-l), 2);
5052  if (NTLN.NumRows() != factors.length()) //refined factors
5053    ident (NTLN, factors.length());
5054  CFListIterator j;
5055  CFArray buf;
5056  mat_zz_pE* NTLC, NTLK;
5057  CanonicalForm bufF, truncF;
5058  Variable y= F.mvar();
5059  while (l <= liftBound)
5060  {
5061    bufFactors.insert (LCF);
5062    henselLiftResume12 (F, bufFactors, oldL, l, Pi, diophant, M);
5063    j= bufFactors;
5064    truncF= mod (F, power (y, l));
5065    if (useOldQs)
5066    {
5067      for (int i= 0; i < bufFactors.length(); i++, j++)
5068        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
5069                                     bufQ[i]);
5070    }
5071    else
5072    {
5073      for (int i= 0; i < bufFactors.length(); i++, j++)
5074        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
5075    }
5076    for (int i= 0; i < d; i++)
5077    {
5078      if (bounds [i] + 1 <= l/2)
5079      {
5080        int k= tmin (bounds [i] + 1, l/2);
5081        CFMatrix C= CFMatrix (l - k, bufFactors.length());
5082        for (int ii= 0; ii < bufFactors.length(); ii++)
5083        {
5084          if (A[ii].size() - 1 >= i)
5085          {
5086            buf= getCoeffs (A[ii] [i], k);
5087            writeInMatrix (C, buf, ii + 1, 0);
5088          }
5089        }
5090        NTLC= convertFacCFMatrix2NTLmat_zz_pE(C);
5091        NTLK= (*NTLC)*NTLN;
5092        transpose (NTLK, NTLK);
5093        kernel (NTLK, NTLK);
5094        transpose (NTLK, NTLK);
5095        NTLN *= NTLK;
5096        if (NTLN.NumCols() == 1)
5097        {
5098          irreducible= true;
5099          break;
5100        }
5101      }
5102    }
5103    if (NTLN.NumCols() == 1)
5104    {
5105      irreducible= true;
5106      break;
5107    }
5108
5109    int * zeroOneVecs= extractZeroOneVecs (NTLN);
5110    bufF= F;
5111    result= reconstruction (bufF, bufFactors, zeroOneVecs, l, NTLN, eval);
5112    delete [] zeroOneVecs;
5113    if (result.length() > 0 && degree (bufF) + 1 + degree (LC (bufF, 1)) <= l)
5114    {
5115      F= bufF;
5116      factors= bufFactors;
5117      delete [] A;
5118      return result;
5119    }
5120
5121    if (isReduced (NTLN))
5122    {
5123      int factorsFound= 0;
5124      bufF= F;
5125      int* factorsFoundIndex= new int [NTLN.NumCols()];
5126      for (long i= 0; i < NTLN.NumCols(); i++)
5127        factorsFoundIndex[i]= 0;
5128      if (l < liftBound)
5129        reconstructionTry (result, bufF, bufFactors, l, factorsFound,
5130                           factorsFoundIndex, NTLN, eval, false
5131                          );
5132      else
5133        reconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
5134                           degree (LCF), factorsFound, factorsFoundIndex,
5135                           NTLN, eval, false
5136                          );
5137      if (NTLN.NumCols() == result.length())
5138      {
5139        delete [] A;
5140        delete [] factorsFoundIndex;
5141        return result;
5142      }
5143      delete [] factorsFoundIndex;
5144    }
5145    result= CFList();
5146    oldL= l;
5147    stepSize *= 2;
5148    l += stepSize;
5149    if (l > liftBound)
5150    {
5151      if (!hitBound)
5152      {
5153        l= liftBound;
5154        hitBound= true;
5155      }
5156      else
5157        break;
5158    }
5159  }
5160  if (irreducible)
5161  {
5162    delete [] A;
5163    return CFList (F (y-eval,y));
5164  }
5165  delete [] A;
5166  factors= bufFactors;
5167  return CFList();
5168}
5169
5170//over field extension
5171CFList
5172extFurtherLiftingAndIncreasePrecision (CanonicalForm& F, CFList& factors, int l,
5173                                       int liftBound, int d, int* bounds,
5174                                       mat_zz_p& NTLN, CFList& diophant,
5175                                       CFMatrix& M, CFArray& Pi, CFArray& bufQ,
5176                                       const CanonicalForm& evaluation, const
5177                                       ExtensionInfo& info, CFList& source,
5178                                       CFList& dest
5179                                      )
5180{
5181  CanonicalForm LCF= LC (F, 1);
5182  CFList result;
5183  bool irreducible= false;
5184  CFList bufFactors= factors;
5185  CFArray *A = new CFArray [bufFactors.length()];
5186  bool useOldQs= false;
5187  bool hitBound= false;
5188  bool GF= (CFFactory::gettype()==GaloisFieldDomain);
5189  int degMipo= degree (getMipo (info.getAlpha()));
5190  Variable alpha= info.getAlpha();
5191  int oldL= l; //be careful
5192  int stepSize= 8;
5193  l += tmax (tmin (8, degree (F) + 1 + degree (LC (F, 1))-l),2);
5194  Variable gamma= info.getBeta();
5195  CanonicalForm primElemAlpha= info.getGamma();
5196  CanonicalForm imPrimElemAlpha= info.getDelta();
5197  if (NTLN.NumRows() != factors.length()) //refined factors
5198    ident (NTLN, factors.length());
5199  Variable y= F.mvar();
5200  CanonicalForm powX, imBasis, bufF, truncF;
5201  CFMatrix Mat, C;
5202  CFIterator iter;
5203  mat_zz_p* NTLMat,*NTLC, NTLK;
5204  CFListIterator j;
5205  CFArray buf;
5206  while (l <= liftBound)
5207  {
5208    bufFactors.insert (LCF);
5209    henselLiftResume12 (F, bufFactors, oldL, l, Pi, diophant, M);
5210
5211    if (GF)
5212      setCharacteristic (getCharacteristic());
5213
5214    powX= power (y-gamma, l);
5215    Mat= CFMatrix (l*degMipo, l*degMipo);
5216    for (int i= 0; i < l*degMipo; i++)
5217    {
5218
5219      imBasis= mod (power (y, i), powX);
5220      imBasis= imBasis (power (y, degMipo), y);
5221      imBasis= imBasis (y, gamma);
5222      iter= imBasis;
5223      for (; iter.hasTerms(); iter++)
5224        Mat (iter.exp()+ 1, i+1)= iter.coeff();
5225    }
5226
5227    NTLMat= convertFacCFMatrix2NTLmat_zz_p (Mat);
5228    *NTLMat= inv (*NTLMat);
5229
5230    if (GF)
5231      setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
5232
5233    j= bufFactors;
5234    truncF= mod (F, power (y, l));
5235    if (useOldQs)
5236    {
5237      for (int i= 0; i < bufFactors.length(); i++, j++)
5238        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
5239                                     bufQ[i]);
5240    }
5241    else
5242    {
5243      for (int i= 0; i < bufFactors.length(); i++, j++)
5244        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
5245    }
5246    for (int i= 0; i < d; i++)
5247    {
5248      if (bounds [i] + 1 <= l/2)
5249      {
5250        int k= tmin (bounds [i] + 1, l/2);
5251        C= CFMatrix (l*degMipo - k, bufFactors.length());
5252        for (int ii= 0; ii < bufFactors.length(); ii++)
5253        {
5254          if (A[ii].size() - 1 >= i)
5255          {
5256            if (GF)
5257            {
5258              A [ii] [i]= A [ii] [i] (y-evaluation, y);
5259              setCharacteristic (getCharacteristic());
5260              A[ii] [i]= GF2FalphaRep (A[ii] [i], alpha);
5261              if (alpha != gamma)
5262                A [ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
5263                                     gamma, source, dest
5264                                    );
5265              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
5266            }
5267            else
5268            {
5269              A [ii] [i]= A [ii] [i] (y-evaluation, y);
5270              if (alpha != gamma)
5271                A[ii] [i]= mapDown (A[ii] [i], imPrimElemAlpha, primElemAlpha,
5272                                    gamma, source, dest
5273                                   );
5274              buf= getCoeffs (A[ii] [i], k, l, degMipo, gamma, 0, *NTLMat);
5275            }
5276            writeInMatrix (C, buf, ii + 1, 0);
5277          }
5278          if (GF)
5279            setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
5280        }
5281
5282        if (GF)
5283          setCharacteristic(getCharacteristic());
5284        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
5285        NTLK= (*NTLC)*NTLN;
5286        transpose (NTLK, NTLK);
5287        kernel (NTLK, NTLK);
5288        transpose (NTLK, NTLK);
5289        NTLN *= NTLK;
5290        if (GF)
5291          setCharacteristic (getCharacteristic(), degMipo, info.getGFName());
5292
5293        if (NTLN.NumCols() == 1)
5294        {
5295          irreducible= true;
5296          break;
5297        }
5298      }
5299    }
5300    if (NTLN.NumCols() == 1)
5301    {
5302      irreducible= true;
5303      break;
5304    }
5305
5306    int * zeroOneVecs= extractZeroOneVecs (NTLN);
5307    bufF= F;
5308    result= extReconstruction (bufF, bufFactors, zeroOneVecs, l, NTLN, info,
5309                               evaluation
5310                              );
5311    delete [] zeroOneVecs;
5312    if (result.length() > 0 && degree (bufF) + 1 + degree (LC (bufF, 1)) <= l)
5313    {
5314      F= bufF;
5315      factors= bufFactors;
5316      delete [] A;
5317      return result;
5318    }
5319
5320    if (isReduced (NTLN))
5321    {
5322      int factorsFound= 0;
5323      bufF= F;
5324      int* factorsFoundIndex= new int [NTLN.NumCols()];
5325      for (long i= 0; i < NTLN.NumCols(); i++)
5326        factorsFoundIndex[i]= 0;
5327      if (l < degree (bufF) + 1 + degree (LCF))
5328        extReconstructionTry (result, bufF, bufFactors, l, factorsFound,
5329                              factorsFoundIndex, NTLN, false, info, evaluation
5330                             );
5331      else
5332        extReconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
5333                              degree (LCF), factorsFound, factorsFoundIndex,
5334                              NTLN, false, info, evaluation
5335                             );
5336      if (NTLN.NumCols() == result.length())
5337      {
5338        delete [] A;
5339        delete [] factorsFoundIndex;
5340        return result;
5341      }
5342      delete [] factorsFoundIndex;
5343    }
5344    result= CFList();
5345    oldL= l;
5346    stepSize *= 2;
5347    l += stepSize;
5348    if (l > liftBound)
5349    {
5350      if (!hitBound)
5351      {
5352        l= liftBound;
5353        hitBound= true;
5354      }
5355      else
5356        break;
5357    }
5358  }
5359  if (irreducible)
5360  {
5361    delete [] A;
5362    Variable y= Variable (2);
5363    CanonicalForm tmp= F (y - evaluation, y);
5364    CFList source, dest;
5365    tmp= mapDown (tmp, info, source, dest);
5366    return CFList (tmp);
5367  }
5368  delete [] A;
5369  factors= bufFactors;
5370  return CFList();
5371}
5372
5373#ifdef HAVE_FLINT
5374CFList
5375furtherLiftingAndIncreasePrecisionFq2Fp (CanonicalForm& F, CFList& factors, int
5376                                         l, int liftBound, int d, int* bounds,
5377                                         nmod_mat_t FLINTN, CFList& diophant,
5378                                         CFMatrix& M, CFArray& Pi, CFArray& bufQ,
5379                                         const Variable& alpha,
5380                                         const CanonicalForm& eval
5381                                        )
5382#else
5383CFList
5384furtherLiftingAndIncreasePrecisionFq2Fp (CanonicalForm& F, CFList& factors, int
5385                                         l, int liftBound, int d, int* bounds,
5386                                         mat_zz_p& NTLN, CFList& diophant,
5387                                         CFMatrix& M, CFArray& Pi, CFArray& bufQ,
5388                                         const Variable& alpha,
5389                                         const CanonicalForm& eval
5390                                        )
5391#endif
5392{
5393  CanonicalForm LCF= LC (F, 1);
5394  CFList result;
5395  bool irreducible= false;
5396  CFList bufFactors= factors;
5397  CFArray *A = new CFArray [bufFactors.length()];
5398  bool useOldQs= false;
5399  int extensionDeg= degree (getMipo (alpha));
5400  bool hitBound= false;
5401  int oldL= l;
5402  int stepSize= 8; //TODO choose better step size?
5403  l += tmax (tmin (8, degree (F) + 1 + degree (LC (F, 1))-l), 2);
5404#ifdef HAVE_FLINT
5405  if (nmod_mat_nrows (FLINTN) != factors.length()) //refined factors
5406  {
5407    nmod_mat_clear (FLINTN);
5408    nmod_mat_init(FLINTN,factors.length(),factors.length(),getCharacteristic());
5409    for (long i=factors.length()-1; i >= 0; i--)
5410      nmod_mat_entry (FLINTN, i, i)= 1;
5411  }
5412#else
5413  if (NTLN.NumRows() != factors.length()) //refined factors
5414    ident (NTLN, factors.length());
5415#endif
5416  CFListIterator j;
5417  CFMatrix C;
5418#ifdef HAVE_FLINT
5419  long rank;
5420  nmod_mat_t FLINTC, FLINTK, null;
5421#else
5422  mat_zz_p* NTLC, NTLK;
5423#endif
5424  CanonicalForm bufF, truncF;
5425  Variable y= F.mvar();
5426  while (l <= liftBound)
5427  {
5428    bufFactors.insert (LCF);
5429    henselLiftResume12 (F, bufFactors, oldL, l, Pi, diophant, M);
5430    j= bufFactors;
5431    truncF= mod (F, power (y, l));
5432    if (useOldQs)
5433    {
5434      for (int i= 0; i < bufFactors.length(); i++, j++)
5435        A[i]= logarithmicDerivative (truncF, j.getItem(), l, oldL, bufQ[i],
5436                                     bufQ[i]);
5437    }
5438    else
5439    {
5440      for (int i= 0; i < bufFactors.length(); i++, j++)
5441        A[i]= logarithmicDerivative (truncF, j.getItem(), l, bufQ [i]);
5442    }
5443    for (int i= 0; i < d; i++)
5444    {
5445      if (bounds [i] + 1 <= l/2)
5446      {
5447        int k= tmin (bounds [i] + 1, l/2);
5448        C= CFMatrix ((l - k)*extensionDeg, bufFactors.length());
5449        for (int ii= 0; ii < bufFactors.length(); ii++)
5450        {
5451          CFArray buf;
5452          if (A[ii].size() - 1 >= i)
5453          {
5454            buf= getCoeffs (A[ii] [i], k, alpha);
5455            writeInMatrix (C, buf, ii + 1, 0);
5456          }
5457        }
5458#ifdef HAVE_FLINT
5459        convertFacCFMatrix2nmod_mat_t (FLINTC, C);
5460        nmod_mat_init (FLINTK, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTN),
5461                       getCharacteristic());
5462        nmod_mat_mul (FLINTK, FLINTC, FLINTN);
5463        nmod_mat_init (null, nmod_mat_ncols (FLINTK), nmod_mat_ncols (FLINTK),
5464                       getCharacteristic());
5465        rank= nmod_mat_nullspace (null, FLINTK);
5466        nmod_mat_clear (FLINTK);
5467        nmod_mat_window_init (FLINTK, null, 0, 0, nmod_mat_nrows(null), rank);
5468        nmod_mat_clear (FLINTC);
5469        nmod_mat_init_set (FLINTC, FLINTN);
5470        nmod_mat_clear (FLINTN);
5471        nmod_mat_init (FLINTN, nmod_mat_nrows (FLINTC), nmod_mat_ncols (FLINTK),
5472                       getCharacteristic());
5473        nmod_mat_mul (FLINTN, FLINTC, FLINTK); //no aliasing allowed!!
5474
5475        nmod_mat_clear (FLINTC);
5476        nmod_mat_window_clear (FLINTK);
5477        nmod_mat_clear (null);
5478#else
5479        NTLC= convertFacCFMatrix2NTLmat_zz_p(C);
5480        NTLK= (*NTLC)*NTLN;
5481        transpose (NTLK, NTLK);
5482        kernel (NTLK, NTLK);
5483        transpose (NTLK, NTLK);
5484        NTLN *= NTLK;
5485#endif
5486#ifdef HAVE_FLINT
5487        if (nmod_mat_ncols (FLINTN) == 1)
5488#else
5489        if (NTLN.NumCols() == 1)
5490#endif
5491        {
5492          irreducible= true;
5493          break;
5494        }
5495      }
5496    }
5497#ifdef HAVE_FLINT
5498    if (nmod_mat_ncols (FLINTN) == 1)
5499#else
5500    if (NTLN.NumCols() == 1)
5501#endif
5502    {
5503      irreducible= true;
5504      break;
5505    }
5506
5507#ifdef HAVE_FLINT
5508    int * zeroOneVecs= extractZeroOneVecs (FLINTN);
5509#else
5510    int * zeroOneVecs= extractZeroOneVecs (NTLN);
5511#endif
5512    CanonicalForm bufF= F;
5513#ifdef HAVE_FLINT
5514    result= reconstruction (bufF, bufFactors, zeroOneVecs, l, FLINTN, eval);
5515#else
5516    result= reconstruction (bufF, bufFactors, zeroOneVecs, l, NTLN,  eval);
5517#endif
5518    delete [] zeroOneVecs;
5519    if (result.length() > 0 && degree (bufF) + 1 + degree (LC (bufF, 1)) <= l)
5520    {
5521      F= bufF;
5522      factors= bufFactors;
5523      delete [] A;
5524      return result;
5525    }
5526
5527#ifdef HAVE_FLINT
5528    if (isReduced (FLINTN))
5529#else
5530    if (isReduced (NTLN))
5531#endif
5532    {
5533      int factorsFound= 0;
5534      bufF= F;
5535#ifdef HAVE_FLINT
5536      int* factorsFoundIndex= new int [nmod_mat_ncols (FLINTN)];
5537      for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
5538#else
5539      int* factorsFoundIndex= new int [NTLN.NumCols()];
5540      for (long i= 0; i < NTLN.NumCols(); i++)
5541#endif
5542        factorsFoundIndex[i]= 0;
5543#ifdef HAVE_FLINT
5544      if (l < degree (bufF) + 1 + degree (LCF))
5545        reconstructionTry (result, bufF, bufFactors, l, factorsFound,
5546                           factorsFoundIndex, FLINTN, eval, false
5547                          );
5548      else
5549        reconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
5550                           degree (LCF), factorsFound, factorsFoundIndex,
5551                           FLINTN, eval, false
5552                          );
5553      if (nmod_mat_ncols (FLINTN) == result.length())
5554#else
5555      if (l < degree (bufF) + 1 + degree (LCF))
5556        reconstructionTry (result, bufF, bufFactors, l, factorsFound,
5557                           factorsFoundIndex, NTLN, eval, false
5558                          );
5559      else
5560        reconstructionTry (result, bufF, bufFactors, degree (bufF) + 1 +
5561                           degree (LCF), factorsFound, factorsFoundIndex,
5562                           NTLN, eval, false
5563                          );
5564      if (NTLN.NumCols() == result.length())
5565#endif
5566      {
5567        delete [] A;
5568        delete [] factorsFoundIndex;
5569        return result;
5570      }
5571      delete [] factorsFoundIndex;
5572    }
5573    result= CFList();
5574    oldL= l;
5575    stepSize *= 2;
5576    l += stepSize;
5577    if (l > liftBound)
5578    {
5579      if (!hitBound)
5580      {
5581        l= liftBound;
5582        hitBound= true;
5583      }
5584      else
5585        break;
5586    }
5587  }
5588  if (irreducible)
5589  {
5590    delete [] A;
5591    return CFList (F (y-eval,y));
5592  }
5593  delete [] A;
5594  factors= bufFactors;
5595  return CFList();
5596}
5597
5598void
5599refineAndRestartLift (const CanonicalForm& F, const mat_zz_p& NTLN, int
5600                      liftBound, int l, CFList& factors, CFMatrix& M, CFArray&
5601                      Pi, CFList& diophant
5602                     )
5603{
5604  CFList bufFactors;
5605  Variable y= Variable (2);
5606  CanonicalForm LCF= LC (F, 1);
5607  CFListIterator iter;
5608  CanonicalForm buf;
5609  for (long i= 1; i <= NTLN.NumCols(); i++)
5610  {
5611    iter= factors;
5612    buf= 1;
5613    for (long j= 1; j <= NTLN.NumRows(); j++, iter++)
5614    {
5615      if (!IsZero (NTLN (j,i)))
5616        buf= mulNTL (buf, mod (iter.getItem(), y));
5617    }
5618    bufFactors.append (buf);
5619  }
5620  factors= bufFactors;
5621  M= CFMatrix (liftBound, factors.length());
5622  Pi= CFArray();
5623  diophant= CFList();
5624  factors.insert (LCF);
5625  henselLift12 (F, factors, l, Pi, diophant, M);
5626}
5627
5628#ifdef HAVE_FLINT
5629void
5630refineAndRestartLift (const CanonicalForm& F, const nmod_mat_t FLINTN, int
5631                      liftBound, int l, CFList& factors, CFMatrix& M, CFArray&
5632                      Pi, CFList& diophant
5633                     )
5634{
5635  CFList bufFactors;
5636  Variable y= Variable (2);
5637  CanonicalForm LCF= LC (F, 1);
5638  CFListIterator iter;
5639  CanonicalForm buf;
5640  for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
5641  {
5642    iter= factors;
5643    buf= 1;
5644    for (long j= 0; j < nmod_mat_nrows (FLINTN); j++, iter++)
5645    {
5646      if (!(nmod_mat_entry (FLINTN,j,i) == 0))
5647        buf= mulNTL (buf, mod (iter.getItem(), y));
5648    }
5649    bufFactors.append (buf);
5650  }
5651  factors= bufFactors;
5652  M= CFMatrix (liftBound, factors.length());
5653  Pi= CFArray();
5654  diophant= CFList();
5655  factors.insert (LCF);
5656  henselLift12 (F, factors, l, Pi, diophant, M);
5657}
5658#endif
5659
5660void
5661refineAndRestartLift (const CanonicalForm& F, const mat_zz_pE& NTLN, int
5662                      liftBound, int l, CFList& factors, CFMatrix& M, CFArray&
5663                      Pi, CFList& diophant
5664                     )
5665{
5666  CFList bufFactors;
5667  Variable y= Variable (2);
5668  CanonicalForm LCF= LC (F, 1);
5669  CFListIterator iter;
5670  CanonicalForm buf;
5671  for (long i= 1; i <= NTLN.NumCols(); i++)
5672  {
5673    iter= factors;
5674    buf= 1;
5675    for (long j= 1; j <= NTLN.NumRows(); j++, iter++)
5676    {
5677      if (!IsZero (NTLN (j,i)))
5678        buf= mulNTL (buf, mod (iter.getItem(), y));
5679    }
5680    bufFactors.append (buf);
5681  }
5682  factors= bufFactors;
5683  M= CFMatrix (liftBound, factors.length());
5684  Pi= CFArray();
5685  diophant= CFList();
5686  factors.insert (LCF);
5687  henselLift12 (F, factors, l, Pi, diophant, M);
5688}
5689
5690#ifdef HAVE_FLINT
5691CFList
5692earlyReconstructionAndLifting (const CanonicalForm& F, const nmod_mat_t N,
5693                               CanonicalForm& bufF, CFList& factors, int& l,
5694                               int& factorsFound, bool beenInThres, CFMatrix& M,
5695                               CFArray& Pi, CFList& diophant, bool symmetric,
5696                               const CanonicalForm& evaluation
5697                              )
5698#else
5699CFList
5700earlyReconstructionAndLifting (const CanonicalForm& F, const mat_zz_p& N,
5701                               CanonicalForm& bufF, CFList& factors, int& l,
5702                               int& factorsFound, bool beenInThres, CFMatrix& M,
5703                               CFArray& Pi, CFList& diophant, bool symmetric,
5704                               const CanonicalForm& evaluation
5705                              )
5706#endif
5707{
5708  int sizeOfLiftPre;
5709  int * liftPre= getLiftPrecisions (F, sizeOfLiftPre, degree (LC (F, 1), 2));
5710
5711  Variable y= F.mvar();
5712  factorsFound= 0;
5713  CanonicalForm LCF= LC (F, 1);
5714  CFList result;
5715  int smallFactorDeg= tmin (11, liftPre [sizeOfLiftPre- 1] + 1);
5716#ifdef HAVE_FLINT
5717  nmod_mat_t FLINTN;
5718  nmod_mat_init_set (FLINTN, N);
5719  int * factorsFoundIndex= new int [nmod_mat_ncols (FLINTN)];
5720  for (long i= 0; i < nmod_mat_ncols (FLINTN); i++)
5721#else
5722  mat_zz_p NTLN= N;
5723  int * factorsFoundIndex= new int [NTLN.NumCols()];
5724  for (long i= 0; i < NTLN.NumCols(); i++)
5725#endif
5726    factorsFoundIndex [i]= 0;
5727
5728  if (degree (F) + 1 > smallFactorDeg)
5729  {
5730    if (l < smallFactorDeg)
5731    {
5732      TIMING_START (fac_fq_lift);
5733      factors.insert (LCF);
5734      henselLiftResume12 (F, factors, l, smallFactorDeg, Pi, diophant, M);
5735      TIMING_END_AND_PRINT (fac_fq_lift, "time to lift in reconstruction0: ");
5736      l= smallFactorDeg;
5737    }
5738#ifdef HAVE_FLINT
5739    TIMING_START (fac_fq_reconstruction);
5740    reconstructionTry (result, bufF, factors, smallFactorDeg, factorsFound,
5741                       factorsFoundIndex, FLINTN, evaluation, beenInThres
5742                      );
5743    TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct0: ");
5744    if (result.length() == nmod_mat_ncols (FLINTN))
5745    {
5746      nmod_mat_clear (FLINTN);
5747#else
5748    TIMING_START (fac_fq_reconstruction);
5749    reconstructionTry (result, bufF, factors, smallFactorDeg, factorsFound,
5750                       factorsFoundIndex, NTLN, evaluation, beenInThres
5751                      );
5752    TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct0: ");
5753    if (result.length() == NTLN.NumCols())
5754    {
5755#endif
5756      delete [] liftPre;
5757      delete [] factorsFoundIndex;
5758      return result;
5759    }
5760  }
5761
5762  int i= sizeOfLiftPre - 1;
5763  int dummy= 1;
5764  if (sizeOfLiftPre > 1 && sizeOfLiftPre < 30)
5765  {
5766    while (i > 0)
5767    {
5768      if (l < liftPre[i-1] + 1)
5769      {
5770        factors.insert (LCF);
5771        TIMING_START (fac_fq_lift);
5772        henselLiftResume12 (F, factors, l, liftPre[i-1] + 1, Pi, diophant, M);
5773        TIMING_END_AND_PRINT (fac_fq_lift, "time to lift in reconstruction1: ");
5774        l= liftPre[i-1] + 1;
5775      }
5776      else
5777      {
5778        i--;
5779        if (i != 0)
5780          continue;
5781      }
5782#ifdef HAVE_FLINT
5783      TIMING_START (fac_fq_reconstruction);
5784      reconstructionTry (result, bufF, factors, l, factorsFound,
5785                         factorsFoundIndex, FLINTN, evaluation, beenInThres
5786                        );
5787      TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct1: ");
5788      if (result.length() == nmod_mat_ncols (FLINTN))
5789      {
5790        nmod_mat_clear (FLINTN);
5791#else
5792      TIMING_START (fac_fq_reconstruction);
5793      reconstructionTry (result, bufF, factors, l, factorsFound,
5794                         factorsFoundIndex, NTLN, evaluation, beenInThres
5795                        );
5796      TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct1: ");
5797      if (result.length() == NTLN.NumCols())
5798      {
5799#endif
5800        delete [] liftPre;
5801        delete [] factorsFoundIndex;
5802        return result;
5803      }
5804      i--;
5805    }
5806  }
5807  else
5808  {
5809    i= 1;
5810    while (((degree (F,y)/4)*i+1) + 4 <= smallFactorDeg)
5811      i++;
5812    while (i < 5)
5813    {
5814      dummy= tmin (degree (F,y)+1, ((degree (F,y)/4)+1)*i+4);
5815      if (l < dummy)
5816      {
5817        factors.insert (LCF);
5818        TIMING_START (fac_fq_lift);
5819        henselLiftResume12 (F, factors, l, dummy, Pi, diophant, M);
5820        TIMING_END_AND_PRINT (fac_fq_lift, "time to lift in reconstruction2: ");
5821        l= dummy;
5822        if (i == 1 && degree (F)%4==0 && symmetric && factors.length() == 2 &&
5823            LC (F,1).inCoeffDomain() &&
5824           (degree (factors.getFirst(), 1) == degree (factors.getLast(),1)))
5825        {
5826          Variable x= Variable (1);
5827          CanonicalForm g, h, gg, hh, multiplier1, multiplier2, check1, check2;
5828          int m= degree (F)/4+1;
5829          g= factors.getFirst();
5830          h= factors.getLast();
5831          g= mod (g, power (y,m));
5832          h= mod (h, power (y,m));
5833          g= g (y-evaluation, y);
5834          h= h (y-evaluation, y);
5835          gg= mod (swapvar (g,x,y),power (x,m));
5836          gg= gg (y + evaluation, y);
5837          multiplier1= factors.getLast()[m-1][0]/gg[m-1][0];
5838          gg= div (gg, power (y,m));
5839          gg= gg*power (y,m);
5840          hh= mod (swapvar (h,x,y),power (x,m));
5841          hh= hh (y + evaluation, y);
5842          multiplier2= factors.getFirst()[m-1][0]/hh[m-1][0];
5843          hh= div (hh, power (y,m));
5844          hh= hh*power (y,m);
5845          gg= multiplier1*gg+mod (factors.getLast(), power (y,m));
5846          hh= multiplier2*hh+mod (factors.getFirst(), power (y,m));
5847          check1= gg (y-evaluation,y);
5848          check2= hh (y-evaluation,y);
5849          CanonicalForm oldcheck1= check1;
5850          check1= swapvar (check1, x, y);
5851          if (check1/Lc (check1) == check2/Lc (check2))
5852          {
5853#ifdef HAVE_FLINT
5854            nmod_mat_clear (FLINTN);
5855#endif
5856            result.append (oldcheck1);
5857            result.append (check2);
5858            delete [] liftPre;
5859            delete [] factorsFoundIndex;
5860            return result;
5861          }
5862        }
5863      }
5864      else
5865      {
5866        i++;
5867        if (i < 5)
5868          continue;
5869      }
5870#ifdef HAVE_FLINT
5871      TIMING_START (fac_fq_reconstruction);
5872      reconstructionTry (result, bufF, factors, l, factorsFound,
5873                         factorsFoundIndex, FLINTN, evaluation, beenInThres
5874                        );
5875      TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct2: ");
5876      if (result.length() == nmod_mat_ncols (FLINTN))
5877      {
5878        nmod_mat_clear (FLINTN);
5879#else
5880      TIMING_START (fac_fq_reconstruction);
5881      reconstructionTry (result, bufF, factors, l, factorsFound,
5882                         factorsFoundIndex, NTLN, evaluation, beenInThres
5883                        );
5884      TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct2: ");
5885      if (result.length() == NTLN.NumCols())
5886      {
5887#endif
5888        delete [] liftPre;
5889        delete [] factorsFoundIndex;
5890        return result;
5891      }
5892      i++;
5893    }
5894  }
5895
5896#ifdef HAVE_FLINT
5897  nmod_mat_clear (FLINTN);
5898#endif
5899  delete [] liftPre;
5900  delete [] factorsFoundIndex;
5901  return result;
5902}
5903
5904CFList
5905earlyReconstructionAndLifting (const CanonicalForm& F, const mat_zz_pE& N,
5906                               CanonicalForm& bufF, CFList& factors, int& l,
5907                               int& factorsFound, bool beenInThres, CFMatrix& M,
5908                               CFArray& Pi, CFList& diophant, bool symmetric,
5909                               const CanonicalForm& evaluation
5910                              )
5911{
5912  int sizeOfLiftPre;
5913  int * liftPre= getLiftPrecisions (F, sizeOfLiftPre, degree (LC (F, 1), 2));
5914  Variable y= F.mvar();
5915  factorsFound= 0;
5916  CanonicalForm LCF= LC (F, 1);
5917  CFList result;
5918  int smallFactorDeg= 11;
5919  mat_zz_pE NTLN= N;
5920  int * factorsFoundIndex= new int [NTLN.NumCols()];
5921  for (long i= 0; i < NTLN.NumCols(); i++)
5922    factorsFoundIndex [i]= 0;
5923
5924  if (degree (F) + 1 > smallFactorDeg)
5925  {
5926    if (l < smallFactorDeg)
5927    {
5928      TIMING_START (fac_fq_lift);
5929      factors.insert (LCF);
5930      henselLiftResume12 (F, factors, l, smallFactorDeg, Pi, diophant, M);
5931      TIMING_END_AND_PRINT (fac_fq_lift, "time to lift in reconstruction0: ");
5932      l= smallFactorDeg;
5933    }
5934    TIMING_START (fac_fq_reconstruction);
5935    reconstructionTry (result, bufF, factors, smallFactorDeg, factorsFound,
5936                       factorsFoundIndex, NTLN, evaluation, beenInThres
5937                      );
5938    TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct0: ");
5939    if (result.length() == NTLN.NumCols())
5940    {
5941      delete [] liftPre;
5942      delete [] factorsFoundIndex;
5943      return result;
5944    }
5945  }
5946
5947  int i= sizeOfLiftPre - 1;
5948  int dummy= 1;
5949  if (sizeOfLiftPre > 1 && sizeOfLiftPre < 30)
5950  {
5951    while (i > 0)
5952    {
5953      if (l < liftPre[i-1] + 1)
5954      {
5955        factors.insert (LCF);
5956        TIMING_START (fac_fq_lift);
5957        henselLiftResume12 (F, factors, l, liftPre[i-1] + 1, Pi, diophant, M);
5958        TIMING_END_AND_PRINT (fac_fq_lift, "time to lift in reconstruction1: ");
5959        l= liftPre[i-1] + 1;
5960      }
5961      else
5962      {
5963        i--;
5964        if (i != 0)
5965          continue;
5966      }
5967      TIMING_START (fac_fq_reconstruction);
5968      reconstructionTry (result, bufF, factors, l, factorsFound,
5969                         factorsFoundIndex, NTLN, evaluation, beenInThres
5970                        );
5971      TIMING_END_AND_PRINT (fac_fq_reconstruction, "time to reconstruct1: ");
5972      if (result.length() == NTLN.NumCols())