source: git/factory/facFqBivar.cc @ 93e7364

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