source: git/factory/facFqBivar.cc @ 2537fa0

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