source: git/factory/facFqBivar.cc @ 514ade

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