source: git/factory/facFqBivar.cc @ 0dff6bc

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