source: git/factory/facHensel.cc @ d80a1a

spielwiese
Last change on this file since d80a1a was d80a1a, checked in by Martin Lee <martinlee84@…>, 13 years ago
bug fix for trac ticket #341 git-svn-id: file:///usr/local/Singular/svn/trunk@14255 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 78.2 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facHensel.cc
5 *
6 * This file implements functions to lift factors via Hensel lifting and
7 * functions for modular multiplication and division with remainder.
8 *
9 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
10 * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
11 * remainder is described in "Fast Recursive Division" by C. Burnikel and
12 * J. Ziegler. Karatsuba multiplication is described in "Modern Computer
13 * Algebra" by J. von zur Gathen and J. Gerhard.
14 *
15 * @author Martin Lee
16 *
17 * @internal @version \$Id$
18 *
19 **/
20/*****************************************************************************/
21
22#include "assert.h"
23#include "debug.h"
24#include "timing.h"
25
26#include "facHensel.h"
27#include "cf_util.h"
28#include "fac_util.h"
29#include "cf_algorithm.h"
30
31#ifdef HAVE_NTL
32#include <NTL/lzz_pEX.h>
33#include "NTLconvert.h"
34
35CanonicalForm
36mulNTL (const CanonicalForm& F, const CanonicalForm& G)
37{
38  if (F.inCoeffDomain() || G.inCoeffDomain())
39    return F*G;
40  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
41  ASSERT (F.level() == G.level(), "expected polys of same level");
42  if (CFFactory::gettype() == GaloisFieldDomain)
43    return F*G;
44  zz_p::init (getCharacteristic());
45  Variable alpha;
46  CanonicalForm result;
47  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
48  {
49    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
50    zz_pE::init (NTLMipo);
51    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
52    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
53    mul (NTLF, NTLF, NTLG);
54    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
55  }
56  else
57  {
58    zz_pX NTLF= convertFacCF2NTLzzpX (F);
59    zz_pX NTLG= convertFacCF2NTLzzpX (G);
60    mul (NTLF, NTLF, NTLG);
61    result= convertNTLzzpX2CF(NTLF, F.mvar());
62  }
63  return result;
64}
65
66CanonicalForm
67modNTL (const CanonicalForm& F, const CanonicalForm& G)
68{
69  if (F.inCoeffDomain() && G.isUnivariate())
70    return F;
71  else if (F.inCoeffDomain() && G.inCoeffDomain())
72    return mod (F, G);
73  else if (F.isUnivariate() && G.inCoeffDomain())
74    return mod (F,G);
75
76  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
77  ASSERT (F.level() == G.level(), "expected polys of same level");
78  if (CFFactory::gettype() == GaloisFieldDomain)
79    return mod (F, G);
80  zz_p::init (getCharacteristic());
81  Variable alpha;
82  CanonicalForm result;
83  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
84  {
85    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
86    zz_pE::init (NTLMipo);
87    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
88    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
89    rem (NTLF, NTLF, NTLG);
90    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
91  }
92  else
93  {
94    zz_pX NTLF= convertFacCF2NTLzzpX (F);
95    zz_pX NTLG= convertFacCF2NTLzzpX (G);
96    rem (NTLF, NTLF, NTLG);
97    result= convertNTLzzpX2CF(NTLF, F.mvar());
98  }
99  return result;
100}
101
102CanonicalForm
103divNTL (const CanonicalForm& F, const CanonicalForm& G)
104{
105  if (F.inCoeffDomain() && G.isUnivariate())
106    return F;
107  else if (F.inCoeffDomain() && G.inCoeffDomain())
108    return div (F, G);
109  else if (F.isUnivariate() && G.inCoeffDomain())
110    return div (F,G);
111
112  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
113  ASSERT (F.level() == G.level(), "expected polys of same level");
114  if (CFFactory::gettype() == GaloisFieldDomain)
115    return div (F, G);
116  zz_p::init (getCharacteristic());
117  Variable alpha;
118  CanonicalForm result;
119  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
120  {
121    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
122    zz_pE::init (NTLMipo);
123    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
124    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
125    div (NTLF, NTLF, NTLG);
126    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
127  }
128  else
129  {
130    zz_pX NTLF= convertFacCF2NTLzzpX (F);
131    zz_pX NTLG= convertFacCF2NTLzzpX (G);
132    div (NTLF, NTLF, NTLG);
133    result= convertNTLzzpX2CF(NTLF, F.mvar());
134  }
135  return result;
136}
137
138/*
139void
140divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
141           CanonicalForm& R)
142{
143  if (F.inCoeffDomain() && G.isUnivariate())
144  {
145    R= F;
146    Q= 0;
147  }
148  else if (F.inCoeffDomain() && G.inCoeffDomain())
149  {
150    divrem (F, G, Q, R);
151    return;
152  }
153  else if (F.isUnivariate() && G.inCoeffDomain())
154  {
155    divrem (F, G, Q, R);
156    return;
157  }
158
159  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
160  ASSERT (F.level() == G.level(), "expected polys of same level");
161  if (CFFactory::gettype() == GaloisFieldDomain)
162  {
163    divrem (F, G, Q, R);
164    return;
165  }
166  zz_p::init (getCharacteristic());
167  Variable alpha;
168  CanonicalForm result;
169  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
170  {
171    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
172    zz_pE::init (NTLMipo);
173    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
174    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
175    zz_pEX NTLQ;
176    zz_pEX NTLR;
177    DivRem (NTLQ, NTLR, NTLF, NTLG);
178    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
179    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
180    return;
181  }
182  else
183  {
184    zz_pX NTLF= convertFacCF2NTLzzpX (F);
185    zz_pX NTLG= convertFacCF2NTLzzpX (G);
186    zz_pX NTLQ;
187    zz_pX NTLR;
188    DivRem (NTLQ, NTLR, NTLF, NTLG);
189    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
190    R= convertNTLzzpX2CF(NTLR, F.mvar());
191    return;
192  }
193}*/
194
195CanonicalForm mod (const CanonicalForm& F, const CFList& M)
196{
197  CanonicalForm A= F;
198  for (CFListIterator i= M; i.hasItem(); i++)
199    A= mod (A, i.getItem());
200  return A;
201}
202
203zz_pX kronSubFp (const CanonicalForm& A, int d)
204{
205  int degAy= degree (A);
206  zz_pX result;
207  result.rep.SetLength (d*(degAy + 1));
208
209  zz_p *resultp;
210  resultp= result.rep.elts();
211  zz_pX buf;
212  zz_p *bufp;
213
214  for (CFIterator i= A; i.hasTerms(); i++)
215  {
216    if (i.coeff().inCoeffDomain())
217      buf= convertFacCF2NTLzzpX (i.coeff());
218    else
219      buf= convertFacCF2NTLzzpX (i.coeff());
220
221    int k= i.exp()*d;
222    bufp= buf.rep.elts();
223    int bufRepLength= (int) buf.rep.length();
224    for (int j= 0; j < bufRepLength; j++)
225      resultp [j + k]= bufp [j];
226  }
227  result.normalize();
228
229  return result;
230}
231
232zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
233{
234  int degAy= degree (A);
235  zz_pEX result;
236  result.rep.SetLength (d*(degAy + 1));
237
238  Variable v;
239  zz_pE *resultp;
240  resultp= result.rep.elts();
241  zz_pEX buf1;
242  zz_pE *buf1p;
243  zz_pX buf2;
244  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
245
246  for (CFIterator i= A; i.hasTerms(); i++)
247  {
248    if (i.coeff().inCoeffDomain())
249    {
250      buf2= convertFacCF2NTLzzpX (i.coeff());
251      buf1= to_zz_pEX (to_zz_pE (buf2));
252    }
253    else
254      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
255
256    int k= i.exp()*d;
257    buf1p= buf1.rep.elts();
258    int buf1RepLength= (int) buf1.rep.length();
259    for (int j= 0; j < buf1RepLength; j++)
260      resultp [j + k]= buf1p [j];
261  }
262  result.normalize();
263
264  return result;
265}
266
267void
268kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
269                const Variable& alpha)
270{
271  int degAy= degree (A);
272  subA1.rep.SetLength ((long) d*(degAy + 1));
273  subA2.rep.SetLength ((long) d*(degAy + 1));
274
275  Variable v;
276  zz_pE *subA1p;
277  zz_pE *subA2p;
278  subA1p= subA1.rep.elts();
279  subA2p= subA2.rep.elts();
280  zz_pEX buf;
281  zz_pE *bufp;
282  zz_pX buf2;
283  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
284
285  for (CFIterator i= A; i.hasTerms(); i++)
286  {
287    if (i.coeff().inCoeffDomain())
288    {
289      buf2= convertFacCF2NTLzzpX (i.coeff());
290      buf= to_zz_pEX (to_zz_pE (buf2));
291    }
292    else
293      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
294
295    int k= i.exp()*d;
296    int kk= (degAy - i.exp())*d;
297    bufp= buf.rep.elts();
298    int bufRepLength= (int) buf.rep.length();
299    for (int j= 0; j < bufRepLength; j++)
300    {
301      subA1p [j + k]= bufp [j];
302      subA2p [j + kk]= bufp [j];
303    }
304  }
305  subA1.normalize();
306  subA2.normalize();
307}
308
309void
310kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
311{
312  int degAy= degree (A);
313  subA1.rep.SetLength ((long) d*(degAy + 1));
314  subA2.rep.SetLength ((long) d*(degAy + 1));
315
316  Variable v;
317  zz_p *subA1p;
318  zz_p *subA2p;
319  subA1p= subA1.rep.elts();
320  subA2p= subA2.rep.elts();
321  zz_pX buf;
322  zz_p *bufp;
323
324  for (CFIterator i= A; i.hasTerms(); i++)
325  {
326    buf= convertFacCF2NTLzzpX (i.coeff());
327
328    int k= i.exp()*d;
329    int kk= (degAy - i.exp())*d;
330    bufp= buf.rep.elts();
331    int bufRepLength= (int) buf.rep.length();
332    for (int j= 0; j < bufRepLength; j++)
333    {
334      subA1p [j + k]= bufp [j];
335      subA2p [j + kk]= bufp [j];
336    }
337  }
338  subA1.normalize();
339  subA2.normalize();
340}
341
342CanonicalForm
343reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
344              const Variable& alpha)
345{
346  Variable y= Variable (2);
347  Variable x= Variable (1);
348
349  zz_pEX f= F;
350  zz_pEX g= G;
351  int degf= deg(f);
352  int degg= deg(g);
353
354  zz_pEX buf1;
355  zz_pEX buf2;
356  zz_pEX buf3;
357
358  zz_pE *buf1p;
359  zz_pE *buf2p;
360  zz_pE *buf3p;
361  if (f.rep.length() < (long) d*(k+1)) //zero padding
362    f.rep.SetLength ((long)d*(k+1));
363
364  zz_pE *gp= g.rep.elts();
365  zz_pE *fp= f.rep.elts();
366  CanonicalForm result= 0;
367  int i= 0;
368  int lf= 0;
369  int lg= d*k;
370  int degfSubLf= degf;
371  int deggSubLg= degg-lg;
372  int repLengthBuf2;
373  int repLengthBuf1;
374  int ii;
375  zz_pE zzpEZero= zz_pE();
376
377  while (degf >= lf || lg >= 0)
378  {
379    if (degfSubLf >= d)
380      repLengthBuf1= d;
381    else if (degfSubLf < 0)
382      repLengthBuf1= 0;
383    else
384      repLengthBuf1= degfSubLf + 1;
385    buf1.rep.SetLength((long) repLengthBuf1);
386
387    buf1p= buf1.rep.elts();
388    for (int ind= 0; ind < repLengthBuf1; ind++)
389      buf1p [ind]= fp [ind + lf];
390    buf1.normalize();
391
392    repLengthBuf1= buf1.rep.length();
393
394    if (deggSubLg >= d - 1)
395      repLengthBuf2= d - 1;
396    else if (deggSubLg < 0)
397      repLengthBuf2= 0;
398    else
399      repLengthBuf2= deggSubLg + 1;
400
401    buf2.rep.SetLength ((long) repLengthBuf2);
402    buf2p= buf2.rep.elts();
403    for (int ind= 0; ind < repLengthBuf2; ind++)
404    {
405      buf2p [ind]= gp [ind + lg];
406    }
407    buf2.normalize();
408
409    repLengthBuf2= buf2.rep.length();
410
411    buf3.rep.SetLength((long) repLengthBuf2 + d);
412    buf3p= buf3.rep.elts();
413    buf2p= buf2.rep.elts();
414    buf1p= buf1.rep.elts();
415    for (int ind= 0; ind < repLengthBuf1; ind++)
416      buf3p [ind]= buf1p [ind];
417    for (int ind= repLengthBuf1; ind < d; ind++)
418      buf3p [ind]= zzpEZero;
419    for (int ind= 0; ind < repLengthBuf2; ind++)
420      buf3p [ind + d]= buf2p [ind];
421    buf3.normalize();
422
423    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
424    i++;
425
426
427    lf= i*d;
428    degfSubLf= degf - lf;
429
430    lg= d*(k-i);
431    deggSubLg= degg - lg;
432
433    buf1p= buf1.rep.elts();
434
435    if (lg >= 0 && deggSubLg > 0)
436    {
437      if (repLengthBuf2 > degfSubLf + 1)
438        degfSubLf= repLengthBuf2 - 1;
439      int tmp= tmin (repLengthBuf1, deggSubLg);
440      for (int ind= 0; ind < tmp; ind++)
441        gp [ind + lg] -= buf1p [ind];
442    }
443
444    if (lg < 0)
445      break;
446
447    buf2p= buf2.rep.elts();
448    if (degfSubLf >= 0)
449    {
450      for (int ind= 0; ind < repLengthBuf2; ind++)
451        fp [ind + lf] -= buf2p [ind];
452    }
453  }
454
455  return result;
456}
457
458CanonicalForm
459reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
460{
461  Variable y= Variable (2);
462  Variable x= Variable (1);
463
464  zz_pX f= F;
465  zz_pX g= G;
466  int degf= deg(f);
467  int degg= deg(g);
468
469  zz_pX buf1;
470  zz_pX buf2;
471  zz_pX buf3;
472
473  zz_p *buf1p;
474  zz_p *buf2p;
475  zz_p *buf3p;
476
477  if (f.rep.length() < (long) d*(k+1)) //zero padding
478    f.rep.SetLength ((long)d*(k+1));
479
480  zz_p *gp= g.rep.elts();
481  zz_p *fp= f.rep.elts();
482  CanonicalForm result= 0;
483  int i= 0;
484  int lf= 0;
485  int lg= d*k;
486  int degfSubLf= degf;
487  int deggSubLg= degg-lg;
488  int repLengthBuf2;
489  int repLengthBuf1;
490  int ii;
491  zz_p zzpZero= zz_p();
492  while (degf >= lf || lg >= 0)
493  {
494    if (degfSubLf >= d)
495      repLengthBuf1= d;
496    else if (degfSubLf < 0)
497      repLengthBuf1= 0;
498    else
499      repLengthBuf1= degfSubLf + 1;
500    buf1.rep.SetLength((long) repLengthBuf1);
501
502    buf1p= buf1.rep.elts();
503    for (int ind= 0; ind < repLengthBuf1; ind++)
504      buf1p [ind]= fp [ind + lf];
505    buf1.normalize();
506
507    repLengthBuf1= buf1.rep.length();
508   
509    if (deggSubLg >= d - 1)
510      repLengthBuf2= d - 1;
511    else if (deggSubLg < 0)
512      repLengthBuf2= 0;
513    else
514      repLengthBuf2= deggSubLg + 1;
515
516    buf2.rep.SetLength ((long) repLengthBuf2);
517    buf2p= buf2.rep.elts();
518    for (int ind= 0; ind < repLengthBuf2; ind++)
519      buf2p [ind]= gp [ind + lg];
520
521    buf2.normalize();
522
523    repLengthBuf2= buf2.rep.length();
524
525
526    buf3.rep.SetLength((long) repLengthBuf2 + d);
527    buf3p= buf3.rep.elts();
528    buf2p= buf2.rep.elts();
529    buf1p= buf1.rep.elts();
530    for (int ind= 0; ind < repLengthBuf1; ind++)
531      buf3p [ind]= buf1p [ind];
532    for (int ind= repLengthBuf1; ind < d; ind++)
533      buf3p [ind]= zzpZero;
534    for (int ind= 0; ind < repLengthBuf2; ind++)
535      buf3p [ind + d]= buf2p [ind];
536    buf3.normalize();
537
538    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
539    i++;
540
541
542    lf= i*d;
543    degfSubLf= degf - lf;
544
545    lg= d*(k-i);
546    deggSubLg= degg - lg;
547
548    buf1p= buf1.rep.elts();
549
550    if (lg >= 0 && deggSubLg > 0)
551    {
552      if (repLengthBuf2 > degfSubLf + 1)
553        degfSubLf= repLengthBuf2 - 1;
554      int tmp= tmin (repLengthBuf1, deggSubLg);
555      for (int ind= 0; ind < tmp; ind++)
556        gp [ind + lg] -= buf1p [ind];
557    }
558    if (lg < 0)
559      break;
560
561    buf2p= buf2.rep.elts();
562    if (degfSubLf >= 0)
563    {
564      for (int ind= 0; ind < repLengthBuf2; ind++)
565        fp [ind + lf] -= buf2p [ind];
566    }
567  }
568
569  return result;
570}
571
572CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
573{
574  Variable y= Variable (2);
575  Variable x= Variable (1);
576
577  zz_pEX f= F;
578  zz_pE *fp= f.rep.elts();
579
580  zz_pEX buf;
581  zz_pE *bufp;
582  CanonicalForm result= 0;
583  int i= 0;
584  int degf= deg(f);
585  int k= 0;
586  int degfSubK;
587  int repLength;
588  while (degf >= k)
589  {
590    degfSubK= degf - k;
591    if (degfSubK >= d)
592      repLength= d;
593    else
594      repLength= degfSubK + 1;
595
596    buf.rep.SetLength ((long) repLength);
597    bufp= buf.rep.elts();
598    for (int j= 0; j < repLength; j++)
599      bufp [j]= fp [j + k];
600    buf.normalize();
601
602    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
603    i++;
604    k= d*i;
605  }
606
607  return result;
608}
609
610CanonicalForm reverseSubstFp (const zz_pX& F, int d)
611{
612  Variable y= Variable (2);
613  Variable x= Variable (1);
614
615  zz_pX f= F;
616  zz_p *fp= f.rep.elts();
617
618  zz_pX buf;
619  zz_p *bufp;
620  CanonicalForm result= 0;
621  int i= 0;
622  int degf= deg(f);
623  int k= 0;
624  int degfSubK;
625  int repLength;
626  while (degf >= k)
627  {
628    degfSubK= degf - k;
629    if (degfSubK >= d)
630      repLength= d;
631    else
632      repLength= degfSubK + 1;
633
634    buf.rep.SetLength ((long) repLength);
635    bufp= buf.rep.elts();
636    for (int j= 0; j < repLength; j++)
637      bufp [j]= fp [j + k];
638    buf.normalize();
639
640    result += convertNTLzzpX2CF (buf, x)*power (y, i);
641    i++;
642    k= d*i;
643  }
644
645  return result;
646}
647
648// assumes input to be reduced mod M and to be an element of Fq not Fp
649CanonicalForm
650mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const 
651                  CanonicalForm& M)
652{
653  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
654  int d2= tmax (degree (F, 2), degree (G, 2));
655
656  zz_pX F1, F2;
657  kronSubRecipro (F1, F2, F, d1);
658  zz_pX G1, G2;
659  kronSubRecipro (G1, G2, G, d1);
660
661  int k= d1*degree (M);
662  MulTrunc (F1, F1, G1, (long) k);
663
664  mul (F2, F2, G2);
665  F2 >>= k;
666
667  return reverseSubst (F1, F2, d1, d2);
668}
669
670//Kronecker substitution
671CanonicalForm
672mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
673              CanonicalForm& M)
674{
675  CanonicalForm A= F;
676  CanonicalForm B= G;
677
678  int p= getCharacteristic();
679
680  int degAx= degree (A, 1);
681  int degAy= degree (A, 2);
682  int degBx= degree (B, 1);
683  int degBy= degree (B, 2);
684  int d1= degAx + 1 + degBx;
685  int d2= tmax (degree (A, 2), degree (B, 2));
686
687  if (d1 > 128 && d2 > 160 && (degAy == degBy))
688    return mulMod2NTLFpReci (A, B, M);
689
690  zz_pX NTLA= kronSubFp (A, d1);
691  zz_pX NTLB= kronSubFp (B, d1);
692
693  int k= d1*degree (M);
694  MulTrunc (NTLA, NTLA, NTLB, (long) k);
695
696  A= reverseSubstFp (NTLA, d1);
697
698  return A;
699}
700
701// assumes input to be reduced mod M and to be an element of Fq not Fp
702CanonicalForm
703mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const 
704                  CanonicalForm& M, const Variable& alpha)
705{
706  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
707  int d2= tmax (degree (F, 2), degree (G, 2));
708
709  zz_pEX F1, F2;
710  kronSubRecipro (F1, F2, F, d1, alpha);
711  zz_pEX G1, G2;
712  kronSubRecipro (G1, G2, G, d1, alpha);
713
714  int k1= d1*degree (M);
715  MulTrunc (F1, F1, G1, (long) k1);
716
717  mul (F2, F2, G2);
718
719  F2 >>= k1;
720
721  CanonicalForm result= reverseSubst (F1, F2, d1, d2, alpha);
722
723  return result;
724}
725
726CanonicalForm
727mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
728              CanonicalForm& M)
729{
730  Variable alpha;
731  CanonicalForm A= F;
732  CanonicalForm B= G;
733
734  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
735  {
736    int degAx= degree (A, 1);
737    int degAy= degree (A, 2);
738    int degBx= degree (B, 1);
739    int degBy= degree (B, 2);
740    int d1= degAx + degBx + 1;
741    int d2= tmax (degree (A, 2), degree (B, 2));
742    zz_p::init (getCharacteristic());
743    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
744    zz_pE::init (NTLMipo);
745
746    int degMipo= degree (getMipo (alpha));
747    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy))
748      return mulMod2NTLFqReci (A, B, M, alpha);
749
750    zz_pEX NTLA= kronSub (A, d1, alpha);
751    zz_pEX NTLB= kronSub (B, d1, alpha);
752
753    int k= d1*degree (M);
754
755    MulTrunc (NTLA, NTLA, NTLB, (long) k);
756
757    A= reverseSubst (NTLA, d1, alpha);
758
759    return A;
760  }
761  else
762    return mulMod2NTLFp (A, B, M);
763}
764
765CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
766                       const CanonicalForm& M)
767{
768  if (A.isZero() || B.isZero())
769    return 0;
770
771  ASSERT (M.isUnivariate(), "M must be univariate");
772
773  CanonicalForm F= mod (A, M);
774  CanonicalForm G= mod (B, M);
775  if (F.inCoeffDomain() || G.inCoeffDomain())
776    return F*G;
777  Variable y= M.mvar();
778  int degF= degree (F, y);
779  int degG= degree (G, y);
780
781  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) && 
782      (F.level() == G.level()))
783  { 
784    CanonicalForm result= mulNTL (F, G);
785    return mod (result, M);
786  }
787  else if (degF <= 1 && degG <= 1)
788  { 
789    CanonicalForm result= F*G;
790    return mod (result, M);
791  }
792
793  int sizeF= size (F);
794  int sizeG= size (G);
795 
796  int fallBackToNaive= 50;
797  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
798    return mod (F*G, M);
799
800  if (CFFactory::gettype() != GaloisFieldDomain &&
801      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
802    return mulMod2NTLFq (F, G, M);
803
804  int m= (int) ceil (degree (M)/2.0);
805  if (degF >= m || degG >= m)
806  {
807    CanonicalForm MLo= power (y, m);
808    CanonicalForm MHi= power (y, degree (M) - m);
809    CanonicalForm F0= mod (F, MLo);
810    CanonicalForm F1= div (F, MLo);
811    CanonicalForm G0= mod (G, MLo);
812    CanonicalForm G1= div (G, MLo);
813    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
814    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
815    CanonicalForm F0G0= mulMod2 (F0, G0, M);
816    return F0G0 + MLo*(F0G1 + F1G0);
817  }
818  else
819  {
820    m= (int) ceil (tmax (degF, degG)/2.0);
821    CanonicalForm yToM= power (y, m);
822    CanonicalForm F0= mod (F, yToM);
823    CanonicalForm F1= div (F, yToM);
824    CanonicalForm G0= mod (G, yToM);
825    CanonicalForm G1= div (G, yToM);
826    CanonicalForm H00= mulMod2 (F0, G0, M);
827    CanonicalForm H11= mulMod2 (F1, G1, M);
828    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
829    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
830  }
831  DEBOUTLN (cerr, "fatal end in mulMod2");
832}
833
834CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
835                      const CFList& MOD)
836{
837  if (A.isZero() || B.isZero())
838    return 0;
839
840  if (MOD.length() == 1)
841    return mulMod2 (A, B, MOD.getLast());
842
843  CanonicalForm M= MOD.getLast();
844  CanonicalForm F= mod (A, M);
845  CanonicalForm G= mod (B, M);
846  if (F.inCoeffDomain() || G.inCoeffDomain())
847    return F*G;
848  Variable y= M.mvar();
849  int degF= degree (F, y);
850  int degG= degree (G, y);
851
852  if ((degF <= 1 && F.level() <= M.level()) &&
853      (degG <= 1 && G.level() <= M.level()))
854  {
855    CFList buf= MOD;
856    buf.removeLast();
857    if (degF == 1 && degG == 1)
858    {
859      CanonicalForm F0= mod (F, y);
860      CanonicalForm F1= div (F, y);
861      CanonicalForm G0= mod (G, y);
862      CanonicalForm G1= div (G, y);
863      if (degree (M) > 2)
864      {
865        CanonicalForm H00= mulMod (F0, G0, buf);
866        CanonicalForm H11= mulMod (F1, G1, buf);
867        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
868        return H11*y*y + (H01 - H00 - H11)*y + H00;
869      }
870      else //here degree (M) == 2
871      {
872        buf.append (y);
873        CanonicalForm F0G1= mulMod (F0, G1, buf);
874        CanonicalForm F1G0= mulMod (F1, G0, buf);
875        CanonicalForm F0G0= mulMod (F0, G0, MOD);
876        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
877        return result;
878      }
879    }
880    else if (degF == 1 && degG == 0)
881      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
882    else if (degF == 0 && degG == 1)
883      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
884    else
885      return mulMod (F, G, buf);
886  }
887  int m= (int) ceil (degree (M)/2.0);
888  if (degF >= m || degG >= m)
889  {
890    CanonicalForm MLo= power (y, m);
891    CanonicalForm MHi= power (y, degree (M) - m);
892    CanonicalForm F0= mod (F, MLo);
893    CanonicalForm F1= div (F, MLo);
894    CanonicalForm G0= mod (G, MLo);
895    CanonicalForm G1= div (G, MLo);
896    CFList buf= MOD;
897    buf.removeLast();
898    buf.append (MHi);
899    CanonicalForm F0G1= mulMod (F0, G1, buf);
900    CanonicalForm F1G0= mulMod (F1, G0, buf);
901    CanonicalForm F0G0= mulMod (F0, G0, MOD);
902    return F0G0 + MLo*(F0G1 + F1G0);
903  }
904  else
905  {
906    m= (int) ceil (tmax (degF, degG)/2.0);
907    CanonicalForm yToM= power (y, m);
908    CanonicalForm F0= mod (F, yToM);
909    CanonicalForm F1= div (F, yToM);
910    CanonicalForm G0= mod (G, yToM);
911    CanonicalForm G1= div (G, yToM);
912    CanonicalForm H00= mulMod (F0, G0, MOD);
913    CanonicalForm H11= mulMod (F1, G1, MOD);
914    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
915    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
916  }
917  DEBOUTLN (cerr, "fatal end in mulMod");
918}
919
920CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
921{
922  if (L.isEmpty())
923    return 1;
924  int l= L.length();
925  if (l == 1)
926    return mod (L.getFirst(), M);
927  else if (l == 2) {
928    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
929    return result;
930  }
931  else
932  {
933    l /= 2;
934    CFList tmp1, tmp2;
935    CFListIterator i= L;
936    CanonicalForm buf1, buf2;
937    for (int j= 1; j <= l; j++, i++)
938      tmp1.append (i.getItem());
939    tmp2= Difference (L, tmp1);
940    buf1= prodMod (tmp1, M);
941    buf2= prodMod (tmp2, M);
942    CanonicalForm result= mulMod2 (buf1, buf2, M);
943    return result;
944  }
945}
946
947CanonicalForm prodMod (const CFList& L, const CFList& M)
948{
949  if (L.isEmpty())
950    return 1;
951  else if (L.length() == 1)
952    return L.getFirst();
953  else if (L.length() == 2)
954    return mulMod (L.getFirst(), L.getLast(), M);
955  else
956  {
957    int l= L.length()/2;
958    CFListIterator i= L;
959    CFList tmp1, tmp2;
960    CanonicalForm buf1, buf2;
961    for (int j= 1; j <= l; j++, i++)
962      tmp1.append (i.getItem());
963    tmp2= Difference (L, tmp1);
964    buf1= prodMod (tmp1, M);
965    buf2= prodMod (tmp2, M);
966    return mulMod (buf1, buf2, M);
967  }
968}
969
970
971CanonicalForm reverse (const CanonicalForm& F, int d)
972{
973  if (d == 0)
974    return F;
975  CanonicalForm A= F;
976  Variable y= Variable (2);
977  Variable x= Variable (1);
978  if (degree (A, x) > 0)
979  {
980    A= swapvar (A, x, y);
981    CanonicalForm result= 0;
982    CFIterator i= A;
983    while (d - i.exp() < 0)
984      i++;
985
986    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
987      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
988    return result;
989  }
990  else
991    return A*power (x, d);
992}
993
994CanonicalForm
995newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
996{
997  int l= ilog2(n);
998
999  CanonicalForm g= mod (F, M)[0] [0];
1000
1001  ASSERT (!g.isZero(), "expected a unit");
1002
1003  Variable alpha;
1004
1005  if (!g.isOne())
1006    g = 1/g;
1007  Variable x= Variable (1);
1008  CanonicalForm result;
1009  int exp= 0;
1010  if (n & 1)
1011  {
1012    result= g;
1013    exp= 1;
1014  }
1015  CanonicalForm h;
1016
1017  for (int i= 1; i <= l; i++)
1018  {
1019    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
1020    h= mod (h, power (x, (1 << i)) - 1);
1021    h= div (h, power (x, (1 << (i - 1))));
1022    h= mod (h, M);
1023    g -= power (x, (1 << (i - 1)))*
1024         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
1025
1026    if (n & (1 << i))
1027    {
1028      if (exp)
1029      {
1030        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
1031        h= mod (h, power (x, exp + (1 << i)) - 1);
1032        h= div (h, power (x, exp));
1033        h= mod (h, M);
1034        result -= power(x, exp)*mod (mulMod2 (g, h, M),
1035                                       power (x, (1 << i)));
1036        exp += (1 << i);
1037      }
1038      else
1039      {
1040        exp= (1 << i);
1041        result= g;
1042      }
1043    }
1044  }
1045
1046  return result;
1047}
1048
1049CanonicalForm
1050newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
1051           M)
1052{
1053  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
1054  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
1055
1056  CanonicalForm A= mod (F, M);
1057  CanonicalForm B= mod (G, M);
1058
1059  Variable x= Variable (1);
1060  int degA= degree (A, x);
1061  int degB= degree (B, x);
1062  int m= degA - degB;
1063  if (m < 0)
1064    return 0;
1065
1066  CanonicalForm Q;
1067  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1068  {
1069    CanonicalForm R;
1070    divrem2 (A, B, Q, R, M);
1071  }
1072  else
1073  {
1074    CanonicalForm R= reverse (A, degA);
1075    CanonicalForm revB= reverse (B, degB);
1076    revB= newtonInverse (revB, m + 1, M);
1077    Q= mulMod2 (R, revB, M);
1078    Q= mod (Q, power (x, m + 1));
1079    Q= reverse (Q, m);
1080  }
1081
1082  return Q;
1083}
1084
1085void
1086newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1087              CanonicalForm& R, const CanonicalForm& M)
1088{
1089  CanonicalForm A= mod (F, M);
1090  CanonicalForm B= mod (G, M);
1091  Variable x= Variable (1);
1092  int degA= degree (A, x);
1093  int degB= degree (B, x);
1094  int m= degA - degB;
1095
1096  if (m < 0)
1097  {
1098    R= A;
1099    Q= 0;
1100    return;
1101  }
1102
1103  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1104  {
1105     divrem2 (A, B, Q, R, M);
1106  }
1107  else
1108  {
1109    R= reverse (A, degA);
1110
1111    CanonicalForm revB= reverse (B, degB);
1112    revB= newtonInverse (revB, m + 1, M);
1113    Q= mulMod2 (R, revB, M);
1114
1115    Q= mod (Q, power (x, m + 1));
1116    Q= reverse (Q, m);
1117
1118    R= A - mulMod2 (Q, B, M);
1119  }
1120}
1121
1122static inline
1123CFList split (const CanonicalForm& F, const int m, const Variable& x)
1124{
1125  CanonicalForm A= F;
1126  CanonicalForm buf= 0;
1127  bool swap= false;
1128  if (degree (A, x) <= 0)
1129    return CFList(A);
1130  else if (x.level() != A.level())
1131  {
1132    swap= true;
1133    A= swapvar (A, x, A.mvar());
1134  }
1135
1136  int j= (int) floor ((double) degree (A)/ m);
1137  CFList result;
1138  CFIterator i= A;
1139  for (; j >= 0; j--)
1140  {
1141    while (i.hasTerms() && i.exp() - j*m >= 0)
1142    {
1143      if (swap)
1144        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
1145      else
1146        buf += i.coeff()*power (x, i.exp() - j*m);
1147      i++;
1148    }
1149    if (swap)
1150      result.append (swapvar (buf, x, F.mvar()));
1151    else
1152      result.append (buf);
1153    buf= 0;
1154  }
1155  return result;
1156}
1157
1158static inline
1159void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1160               CanonicalForm& R, const CFList& M);
1161
1162static inline
1163void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1164               CanonicalForm& R, const CFList& M)
1165{
1166  CanonicalForm A= mod (F, M);
1167  CanonicalForm B= mod (G, M);
1168  Variable x= Variable (1);
1169  int degB= degree (B, x);
1170  int degA= degree (A, x);
1171  if (degA < degB)
1172  {
1173    Q= 0;
1174    R= A;
1175    return;
1176  }
1177  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
1178  if (degB < 1)
1179  {
1180    divrem (A, B, Q, R);
1181    Q= mod (Q, M);
1182    R= mod (R, M);
1183    return;
1184  }
1185
1186  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
1187  CFList splitA= split (A, m, x);
1188  CFList splitB= split (B, m, x);
1189  if (splitA.length() == 3)
1190    splitA.insert (0);
1191  if (splitA.length() == 2)
1192  {
1193    splitA.insert (0);
1194    splitA.insert (0);
1195  }
1196  if (splitA.length() == 1)
1197  {
1198    splitA.insert (0);
1199    splitA.insert (0);
1200    splitA.insert (0);
1201  }
1202
1203  CanonicalForm xToM= power (x, m);
1204
1205  CFListIterator i= splitA;
1206  CanonicalForm H= i.getItem();
1207  i++;
1208  H *= xToM;
1209  H += i.getItem();
1210  i++;
1211  H *= xToM;
1212  H += i.getItem();
1213  i++;
1214
1215  divrem32 (H, B, Q, R, M);
1216
1217  CFList splitR= split (R, m, x);
1218  if (splitR.length() == 1)
1219    splitR.insert (0);
1220
1221  H= splitR.getFirst();
1222  H *= xToM;
1223  H += splitR.getLast();
1224  H *= xToM;
1225  H += i.getItem();
1226
1227  CanonicalForm bufQ;
1228  divrem32 (H, B, bufQ, R, M);
1229
1230  Q *= xToM;
1231  Q += bufQ;
1232  return;
1233}
1234
1235static inline
1236void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1237               CanonicalForm& R, const CFList& M)
1238{
1239  CanonicalForm A= mod (F, M);
1240  CanonicalForm B= mod (G, M);
1241  Variable x= Variable (1);
1242  int degB= degree (B, x);
1243  int degA= degree (A, x);
1244  if (degA < degB)
1245  {
1246    Q= 0;
1247    R= A;
1248    return;
1249  }
1250  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
1251  if (degB < 1)
1252  {
1253    divrem (A, B, Q, R);
1254    Q= mod (Q, M);
1255    R= mod (R, M);
1256    return;
1257  }
1258  int m= (int) ceil ((double) (degB + 1)/ 2.0);
1259
1260  CFList splitA= split (A, m, x);
1261  CFList splitB= split (B, m, x);
1262
1263  if (splitA.length() == 2)
1264  {
1265    splitA.insert (0);
1266  }
1267  if (splitA.length() == 1)
1268  {
1269    splitA.insert (0);
1270    splitA.insert (0);
1271  }
1272  CanonicalForm xToM= power (x, m);
1273
1274  CanonicalForm H;
1275  CFListIterator i= splitA;
1276  i++;
1277
1278  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
1279  {
1280    H= splitA.getFirst()*xToM + i.getItem();
1281    divrem21 (H, splitB.getFirst(), Q, R, M);
1282  }
1283  else
1284  {
1285    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
1286       splitB.getFirst()*xToM;
1287    Q= xToM - 1;
1288  }
1289
1290  H= mulMod (Q, splitB.getLast(), M);
1291
1292  R= R*xToM + splitA.getLast() - H;
1293
1294  while (degree (R, x) >= degB)
1295  {
1296    xToM= power (x, degree (R, x) - degB);
1297    Q += LC (R, x)*xToM;
1298    R -= mulMod (LC (R, x), B, M)*xToM;
1299    Q= mod (Q, M);
1300    R= mod (R, M);
1301  }
1302
1303  return;
1304}
1305
1306void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1307              CanonicalForm& R, const CanonicalForm& M)
1308{
1309  CanonicalForm A= mod (F, M);
1310  CanonicalForm B= mod (G, M);
1311  Variable x= Variable (1);
1312  int degB= degree (B, x);
1313  if (degB > degree (A, x))
1314  {
1315    Q= 0;
1316    R= A;
1317    return;
1318  }
1319
1320  CFList splitA= split (A, degB, x);
1321
1322  CanonicalForm xToDegB= power (x, degB);
1323  CanonicalForm H, bufQ;
1324  Q= 0;
1325  CFListIterator i= splitA;
1326  H= i.getItem()*xToDegB;
1327  i++;
1328  H += i.getItem();
1329  CFList buf;
1330  while (i.hasItem())
1331  {
1332    buf= CFList (M);
1333    divrem21 (H, B, bufQ, R, buf);
1334    i++;
1335    if (i.hasItem())
1336      H= R*xToDegB + i.getItem();
1337    Q *= xToDegB;
1338    Q += bufQ;
1339  }
1340  return;
1341}
1342
1343void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1344             CanonicalForm& R, const CFList& MOD)
1345{
1346  CanonicalForm A= mod (F, MOD);
1347  CanonicalForm B= mod (G, MOD);
1348  Variable x= Variable (1);
1349  int degB= degree (B, x);
1350  if (degB > degree (A, x))
1351  {
1352    Q= 0;
1353    R= A;
1354    return;
1355  }
1356
1357  if (degB <= 0)
1358  {
1359    divrem (A, B, Q, R);
1360    Q= mod (Q, MOD);
1361    R= mod (R, MOD);
1362    return;
1363  }
1364  CFList splitA= split (A, degB, x);
1365
1366  CanonicalForm xToDegB= power (x, degB);
1367  CanonicalForm H, bufQ;
1368  Q= 0;
1369  CFListIterator i= splitA;
1370  H= i.getItem()*xToDegB;
1371  i++;
1372  H += i.getItem();
1373  while (i.hasItem())
1374  {
1375    divrem21 (H, B, bufQ, R, MOD);
1376    i++;
1377    if (i.hasItem())
1378      H= R*xToDegB + i.getItem();
1379    Q *= xToDegB;
1380    Q += bufQ;
1381  }
1382  return;
1383}
1384
1385void sortList (CFList& list, const Variable& x)
1386{
1387  int l= 1;
1388  int k= 1;
1389  CanonicalForm buf;
1390  CFListIterator m;
1391  for (CFListIterator i= list; l <= list.length(); i++, l++)
1392  {
1393    for (CFListIterator j= list; k <= list.length() - l; k++)
1394    {
1395      m= j;
1396      m++;
1397      if (degree (j.getItem(), x) > degree (m.getItem(), x))
1398      {
1399        buf= m.getItem();
1400        m.getItem()= j.getItem();
1401        j.getItem()= buf;
1402        j++;
1403        j.getItem()= m.getItem();
1404      }
1405      else
1406        j++;
1407    }
1408    k= 1;
1409  }
1410}
1411
1412static inline
1413CFList diophantine (const CanonicalForm& F, const CFList& factors)
1414{
1415  CanonicalForm buf1, buf2, buf3, S, T;
1416  CFListIterator i= factors;
1417  CFList result;
1418  if (i.hasItem())
1419    i++;
1420  buf1= F/factors.getFirst();
1421  buf2= divNTL (F, i.getItem());
1422  buf3= extgcd (buf1, buf2, S, T);
1423  result.append (S);
1424  result.append (T);
1425  if (i.hasItem())
1426    i++;
1427  for (; i.hasItem(); i++)
1428  {
1429    buf1= divNTL (F, i.getItem());
1430    buf3= extgcd (buf3, buf1, S, T);
1431    CFListIterator k= factors;
1432    for (CFListIterator j= result; j.hasItem(); j++, k++)
1433    {
1434      j.getItem()= mulNTL (j.getItem(), S);
1435      j.getItem()= modNTL (j.getItem(), k.getItem());
1436    }
1437    result.append (T);
1438  }
1439  return result;
1440}
1441
1442void
1443henselStep12 (const CanonicalForm& F, const CFList& factors,
1444              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1445              CFArray& Pi, int j)
1446{
1447  CanonicalForm E;
1448  CanonicalForm xToJ= power (F.mvar(), j);
1449  Variable x= F.mvar();
1450  // compute the error
1451  if (j == 1)
1452    E= F[j];
1453  else
1454  {
1455    if (degree (Pi [factors.length() - 2], x) > 0)
1456      E= F[j] - Pi [factors.length() - 2] [j];
1457    else
1458      E= F[j];
1459  }
1460
1461  CFArray buf= CFArray (diophant.length());
1462  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1463  int k= 0;
1464  CanonicalForm remainder;
1465  // actual lifting
1466  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1467  {
1468    if (degree (bufFactors[k], x) > 0)
1469    {
1470      if (k > 0)
1471        remainder= modNTL (E, bufFactors[k] [0]);
1472      else
1473        remainder= E;
1474    }
1475    else
1476      remainder= modNTL (E, bufFactors[k]);
1477
1478    buf[k]= mulNTL (i.getItem(), remainder);
1479    if (degree (bufFactors[k], x) > 0)
1480      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1481    else
1482      buf[k]= modNTL (buf[k], bufFactors[k]);
1483  }
1484  for (k= 1; k < factors.length(); k++)
1485    bufFactors[k] += xToJ*buf[k];
1486
1487  // update Pi [0]
1488  int degBuf0= degree (bufFactors[0], x);
1489  int degBuf1= degree (bufFactors[1], x);
1490  if (degBuf0 > 0 && degBuf1 > 0)
1491    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1492  CanonicalForm uIZeroJ;
1493  if (j == 1)
1494  {
1495    if (degBuf0 > 0 && degBuf1 > 0)
1496      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1497                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1498    else if (degBuf0 > 0)
1499      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1500    else if (degBuf1 > 0)
1501      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1502    else
1503      uIZeroJ= 0;
1504    Pi [0] += xToJ*uIZeroJ;
1505  }
1506  else
1507  {
1508    if (degBuf0 > 0 && degBuf1 > 0)
1509      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1510                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1511    else if (degBuf0 > 0)
1512      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1513    else if (degBuf1 > 0)
1514      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1515    else
1516      uIZeroJ= 0;
1517    Pi [0] += xToJ*uIZeroJ;
1518  }
1519  CFArray tmp= CFArray (factors.length() - 1);
1520  for (k= 0; k < factors.length() - 1; k++)
1521    tmp[k]= 0;
1522  CFIterator one, two;
1523  one= bufFactors [0];
1524  two= bufFactors [1];
1525  if (degBuf0 > 0 && degBuf1 > 0)
1526  {
1527    for (k= 1; k <= (int) ceil (j/2.0); k++)
1528    {
1529      if (k != j - k + 1)
1530      {
1531        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
1532        {
1533          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
1534                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1535          one++;
1536          two++;
1537        }
1538        else if (one.exp() == j - k + 1)
1539        {
1540          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
1541                     M (k + 1, 1);
1542          one++;
1543        }
1544        else if (two.exp() == j - k + 1)
1545        {
1546          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
1547                    M (k + 1, 1);
1548          two++;
1549        }
1550      }
1551      else
1552      {
1553        tmp[0] += M (k + 1, 1);
1554      }
1555    }
1556  }
1557  Pi [0] += tmp[0]*xToJ*F.mvar();
1558
1559  // update Pi [l]
1560  int degPi, degBuf;
1561  for (int l= 1; l < factors.length() - 1; l++)
1562  {
1563    degPi= degree (Pi [l - 1], x);
1564    degBuf= degree (bufFactors[l + 1], x);
1565    if (degPi > 0 && degBuf > 0)
1566      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1567    if (j == 1)
1568    {
1569      if (degPi > 0 && degBuf > 0)
1570        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1571                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
1572                  M (1, l + 1));
1573      else if (degPi > 0)
1574        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
1575      else if (degBuf > 0)
1576        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
1577    }
1578    else
1579    {
1580      if (degPi > 0 && degBuf > 0)
1581      {
1582        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1583        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
1584      }
1585      else if (degPi > 0)
1586        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
1587      else if (degBuf > 0)
1588      {
1589        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1590        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
1591      }
1592      Pi[l] += xToJ*uIZeroJ;
1593    }
1594    one= bufFactors [l + 1];
1595    two= Pi [l - 1];
1596    if (two.exp() == j + 1)
1597    {
1598      if (degBuf > 0 && degPi > 0)
1599      {
1600          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
1601          two++;
1602      }
1603      else if (degPi > 0)
1604      {
1605          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
1606          two++;
1607      }
1608    }
1609    if (degBuf > 0 && degPi > 0)
1610    {
1611      for (k= 1; k <= (int) ceil (j/2.0); k++)
1612      {
1613        if (k != j - k + 1)
1614        {
1615          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
1616          {
1617            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
1618                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1619            one++;
1620            two++;
1621          }
1622          else if (one.exp() == j - k + 1)
1623          {
1624            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
1625                       M (k + 1, l + 1);
1626            one++;
1627          }
1628          else if (two.exp() == j - k + 1)
1629          {
1630            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
1631                      M (k + 1, l + 1);
1632            two++;
1633          }
1634        }
1635        else
1636          tmp[l] += M (k + 1, l + 1);
1637      }
1638    }
1639    Pi[l] += tmp[l]*xToJ*F.mvar();
1640  }
1641  return;
1642}
1643
1644void
1645henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1646              CFList& diophant, CFMatrix& M)
1647{
1648  sortList (factors, Variable (1));
1649  Pi= CFArray (factors.length() - 1);
1650  CFListIterator j= factors;
1651  diophant= diophantine (F[0], factors);
1652  DEBOUTLN (cerr, "diophant= " << diophant);
1653  j++;
1654  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
1655  M (1, 1)= Pi [0];
1656  int i= 1;
1657  if (j.hasItem())
1658    j++;
1659  for (j; j.hasItem(); j++, i++)
1660  {
1661    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
1662    M (1, i + 1)= Pi [i];
1663  }
1664  CFArray bufFactors= CFArray (factors.length());
1665  i= 0;
1666  for (CFListIterator k= factors; k.hasItem(); i++, k++)
1667  {
1668    if (i == 0)
1669      bufFactors[i]= mod (k.getItem(), F.mvar());
1670    else
1671      bufFactors[i]= k.getItem();
1672  }
1673  for (i= 1; i < l; i++)
1674    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
1675
1676  CFListIterator k= factors;
1677  for (i= 0; i < factors.length (); i++, k++)
1678    k.getItem()= bufFactors[i];
1679  factors.removeFirst();
1680  return;
1681}
1682
1683void
1684henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
1685                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
1686{
1687  CFArray bufFactors= CFArray (factors.length());
1688  int i= 0;
1689  CanonicalForm xToStart= power (F.mvar(), start);
1690  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1691  {
1692    if (i == 0)
1693      bufFactors[i]= mod (k.getItem(), xToStart);
1694    else
1695      bufFactors[i]= k.getItem();
1696  }
1697  for (i= start; i < end; i++)
1698    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
1699
1700  CFListIterator k= factors;
1701  for (i= 0; i < factors.length(); k++, i++)
1702    k.getItem()= bufFactors [i];
1703  factors.removeFirst();
1704  return;
1705}
1706
1707static inline
1708CFList
1709biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
1710{
1711  Variable y= F.mvar();
1712  CFList result;
1713  if (y.level() == 1)
1714  {
1715    result= diophantine (F, factors);
1716    return result;
1717  }
1718  else
1719  {
1720    CFList buf= factors;
1721    for (CFListIterator i= buf; i.hasItem(); i++)
1722      i.getItem()= mod (i.getItem(), y);
1723    CanonicalForm A= mod (F, y);
1724    int bufD= 1;
1725    CFList recResult= biDiophantine (A, buf, bufD);
1726    CanonicalForm e= 1;
1727    CFList p;
1728    CFArray bufFactors= CFArray (factors.length());
1729    CanonicalForm yToD= power (y, d);
1730    int k= 0;
1731    for (CFListIterator i= factors; i.hasItem(); i++, k++)
1732    {
1733      bufFactors [k]= i.getItem();
1734    }
1735    CanonicalForm b;
1736    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1737    {
1738      b= 1;
1739      if (fdivides (bufFactors[k], F))
1740        b= F/bufFactors[k];
1741      else
1742      {
1743        for (int l= 0; l < factors.length(); l++)
1744        {
1745          if (l == k)
1746            continue;
1747          else
1748          {
1749            b= mulMod2 (b, bufFactors[l], yToD);
1750          }
1751        }
1752      }
1753      p.append (b);
1754    }
1755
1756    CFListIterator j= p;
1757    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1758      e -= i.getItem()*j.getItem();
1759
1760    if (e.isZero())
1761      return recResult;
1762    CanonicalForm coeffE;
1763    CFList s;
1764    result= recResult;
1765    CanonicalForm g;
1766    for (int i= 1; i < d; i++)
1767    {
1768      if (degree (e, y) > 0)
1769        coeffE= e[i];
1770      else
1771        coeffE= 0;
1772      if (!coeffE.isZero())
1773      {
1774        CFListIterator k= result;
1775        CFListIterator l= p;
1776        int ii= 0;
1777        j= recResult;
1778        for (; j.hasItem(); j++, k++, l++, ii++)
1779        {
1780          g= coeffE*j.getItem();
1781          if (degree (bufFactors[ii], y) <= 0)
1782            g= mod (g, bufFactors[ii]);
1783          else
1784            g= mod (g, bufFactors[ii][0]);
1785          k.getItem() += g*power (y, i);
1786          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
1787          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
1788                    mod (e, power (y, i + 1)));
1789        }
1790      }
1791      if (e.isZero())
1792        break;
1793    }
1794
1795    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
1796
1797#ifdef DEBUGOUTPUT
1798    CanonicalForm test= 0;
1799    j= p;
1800    for (CFListIterator i= result; i.hasItem(); i++, j++)
1801      test += mod (i.getItem()*j.getItem(), power (y, d));
1802    DEBOUTLN (cerr, "test= " << test);
1803#endif
1804    return result;
1805  }
1806}
1807
1808static inline
1809CFList
1810multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
1811                     const CFList& recResult, const CFList& M, const int d)
1812{
1813  Variable y= F.mvar();
1814  CFList result;
1815  CFListIterator i;
1816  CanonicalForm e= 1;
1817  CFListIterator j= factors;
1818  CFList p;
1819  CFArray bufFactors= CFArray (factors.length());
1820  CanonicalForm yToD= power (y, d);
1821  int k= 0;
1822  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1823    bufFactors [k]= i.getItem();
1824  CanonicalForm b;
1825  CFList buf= M;
1826  buf.removeLast();
1827  buf.append (yToD);
1828  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1829  {
1830    b= 1;
1831    if (fdivides (bufFactors[k], F))
1832      b= F/bufFactors[k];
1833    else
1834    {
1835      for (int l= 0; l < factors.length(); l++)
1836      {
1837        if (l == k)
1838          continue;
1839        else
1840        {
1841          b= mulMod (b, bufFactors[l], buf);
1842        }
1843      }
1844    }
1845    p.append (b);
1846  }
1847  j= p;
1848  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1849    e -= mulMod (i.getItem(), j.getItem(), M);
1850
1851  if (e.isZero())
1852    return recResult;
1853  CanonicalForm coeffE;
1854  CFList s;
1855  result= recResult;
1856  CanonicalForm g;
1857  for (int i= 1; i < d; i++)
1858  {
1859    if (degree (e, y) > 0)
1860      coeffE= e[i];
1861    else
1862      coeffE= 0;
1863    if (!coeffE.isZero())
1864    {
1865      CFListIterator k= result;
1866      CFListIterator l= p;
1867      j= recResult;
1868      int ii= 0;
1869      CanonicalForm dummy;
1870      for (; j.hasItem(); j++, k++, l++, ii++)
1871      {
1872        g= mulMod (coeffE, j.getItem(), M);
1873        if (degree (bufFactors[ii], y) <= 0)
1874          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
1875                  g, M);
1876        else
1877          divrem (g, bufFactors[ii][0], dummy, g, M);
1878        k.getItem() += g*power (y, i);
1879        e -= mulMod (g*power (y, i), l.getItem(), M);
1880      }
1881    }
1882
1883    if (e.isZero())
1884      break;
1885  }
1886
1887#ifdef DEBUGOUTPUT
1888    CanonicalForm test= 0;
1889    j= p;
1890    for (CFListIterator i= result; i.hasItem(); i++, j++)
1891      test += mod (i.getItem()*j.getItem(), power (y, d));
1892    DEBOUTLN (cerr, "test= " << test);
1893#endif
1894  return result;
1895}
1896
1897static inline
1898void
1899henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
1900            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
1901            const CFList& MOD)
1902{
1903  CanonicalForm E;
1904  CanonicalForm xToJ= power (F.mvar(), j);
1905  Variable x= F.mvar();
1906  // compute the error
1907  if (j == 1)
1908  {
1909    E= F[j];
1910#ifdef DEBUGOUTPUT
1911    CanonicalForm test= 1;
1912    for (int i= 0; i < factors.length(); i++)
1913    {
1914      if (i == 0)
1915        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
1916      else
1917        test= mulMod (test, bufFactors[i], MOD);
1918    }
1919    CanonicalForm test2= mod (F-test, xToJ);
1920
1921    test2= mod (test2, MOD);
1922    DEBOUTLN (cerr, "test= " << test2);
1923#endif
1924  }
1925  else
1926  {
1927#ifdef DEBUGOUTPUT
1928    CanonicalForm test= 1;
1929    for (int i= 0; i < factors.length(); i++)
1930    {
1931      if (i == 0)
1932        test *= mod (bufFactors [i], power (x, j));
1933      else
1934        test *= bufFactors[i];
1935    }
1936    test= mod (test, power (x, j));
1937    test= mod (test, MOD);
1938    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1939    DEBOUTLN (cerr, "test= " << test2);
1940#endif
1941
1942    if (degree (Pi [factors.length() - 2], x) > 0)
1943      E= F[j] - Pi [factors.length() - 2] [j];
1944    else
1945      E= F[j];
1946  }
1947
1948  CFArray buf= CFArray (diophant.length());
1949  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1950  int k= 0;
1951  // actual lifting
1952  CanonicalForm dummy, rest1;
1953  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1954  {
1955    if (degree (bufFactors[k], x) > 0)
1956    {
1957      if (k > 0)
1958        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
1959      else
1960        rest1= E;
1961    }
1962    else
1963      divrem (E, bufFactors[k], dummy, rest1, MOD);
1964
1965    buf[k]= mulMod (i.getItem(), rest1, MOD);
1966
1967    if (degree (bufFactors[k], x) > 0)
1968      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
1969    else
1970      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
1971  }
1972  for (k= 1; k < factors.length(); k++)
1973    bufFactors[k] += xToJ*buf[k];
1974
1975  // update Pi [0]
1976  int degBuf0= degree (bufFactors[0], x);
1977  int degBuf1= degree (bufFactors[1], x);
1978  if (degBuf0 > 0 && degBuf1 > 0)
1979    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
1980  CanonicalForm uIZeroJ;
1981  if (j == 1)
1982  {
1983    if (degBuf0 > 0 && degBuf1 > 0)
1984      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1985                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1986    else if (degBuf0 > 0)
1987      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
1988    else if (degBuf1 > 0)
1989      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1990    else
1991      uIZeroJ= 0;
1992    Pi [0] += xToJ*uIZeroJ;
1993  }
1994  else
1995  {
1996    if (degBuf0 > 0 && degBuf1 > 0)
1997      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1998                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1999    else if (degBuf0 > 0)
2000      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
2001    else if (degBuf1 > 0)
2002      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2003    else
2004      uIZeroJ= 0;
2005    Pi [0] += xToJ*uIZeroJ;
2006  }
2007
2008  CFArray tmp= CFArray (factors.length() - 1);
2009  for (k= 0; k < factors.length() - 1; k++)
2010    tmp[k]= 0;
2011  CFIterator one, two;
2012  one= bufFactors [0];
2013  two= bufFactors [1];
2014  if (degBuf0 > 0 && degBuf1 > 0)
2015  {
2016    for (k= 1; k <= (int) ceil (j/2.0); k++)
2017    {
2018      if (k != j - k + 1)
2019      {
2020        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2021        {
2022          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2023                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2024                    M (j - k + 2, 1);
2025          one++;
2026          two++;
2027        }
2028        else if (one.exp() == j - k + 1)
2029        {
2030          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2031                    bufFactors[1] [k], MOD) - M (k + 1, 1);
2032          one++;
2033        }
2034        else if (two.exp() == j - k + 1)
2035        {
2036          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2037                    two.coeff()), MOD) - M (k + 1, 1);
2038          two++;
2039        }
2040      }
2041      else
2042      {
2043        tmp[0] += M (k + 1, 1);
2044      }
2045    }
2046  }
2047  Pi [0] += tmp[0]*xToJ*F.mvar();
2048
2049  // update Pi [l]
2050  int degPi, degBuf;
2051  for (int l= 1; l < factors.length() - 1; l++)
2052  {
2053    degPi= degree (Pi [l - 1], x);
2054    degBuf= degree (bufFactors[l + 1], x);
2055    if (degPi > 0 && degBuf > 0)
2056      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2057    if (j == 1)
2058    {
2059      if (degPi > 0 && degBuf > 0)
2060        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
2061                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
2062                  M (1, l + 1));
2063      else if (degPi > 0)
2064        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
2065      else if (degBuf > 0)
2066        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
2067    }
2068    else
2069    {
2070      if (degPi > 0 && degBuf > 0)
2071      {
2072        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2073        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
2074      }
2075      else if (degPi > 0)
2076        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
2077      else if (degBuf > 0)
2078      {
2079        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2080        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
2081      }
2082      Pi[l] += xToJ*uIZeroJ;
2083    }
2084    one= bufFactors [l + 1];
2085    two= Pi [l - 1];
2086    if (two.exp() == j + 1)
2087    {
2088      if (degBuf > 0 && degPi > 0)
2089      {
2090          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
2091          two++;
2092      }
2093      else if (degPi > 0)
2094      {
2095          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
2096          two++;
2097      }
2098    }
2099    if (degBuf > 0 && degPi > 0)
2100    {
2101      for (k= 1; k <= (int) ceil (j/2.0); k++)
2102      {
2103        if (k != j - k + 1)
2104        {
2105          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2106          {
2107            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2108                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2109                      M (j - k + 2, l + 1);
2110            one++;
2111            two++;
2112          }
2113          else if (one.exp() == j - k + 1)
2114          {
2115            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2116                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2117            one++;
2118          }
2119          else if (two.exp() == j - k + 1)
2120          {
2121            tmp[l] += mulMod (bufFactors[l + 1] [k],
2122                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2123            two++;
2124          }
2125        }
2126        else
2127          tmp[l] += M (k + 1, l + 1);
2128      }
2129    }
2130    Pi[l] += tmp[l]*xToJ*F.mvar();
2131  }
2132
2133  return;
2134}
2135
2136CFList
2137henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
2138              diophant, CFArray& Pi, CFMatrix& M)
2139{
2140  CFList buf= factors;
2141  int k= 0;
2142  int liftBound;
2143  int liftBoundBivar= l[k];
2144  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
2145  CFList MOD;
2146  MOD.append (power (Variable (2), liftBoundBivar));
2147  CFArray bufFactors= CFArray (factors.length());
2148  k= 0;
2149  CFListIterator j= eval;
2150  j++;
2151  buf.removeFirst();
2152  buf.insert (LC (j.getItem(), 1));
2153  for (CFListIterator i= buf; i.hasItem(); i++, k++)
2154    bufFactors[k]= i.getItem();
2155  Pi= CFArray (factors.length() - 1);
2156  CFListIterator i= buf;
2157  i++;
2158  Variable y= j.getItem().mvar();
2159  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
2160  M (1, 1)= Pi [0];
2161  k= 1;
2162  if (i.hasItem())
2163    i++;
2164  for (i; i.hasItem(); i++, k++)
2165  {
2166    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
2167    M (1, k + 1)= Pi [k];
2168  }
2169
2170  for (int d= 1; d < l[1]; d++)
2171    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
2172  CFList result;
2173  for (k= 1; k < factors.length(); k++)
2174    result.append (bufFactors[k]);
2175  return result;
2176}
2177
2178void
2179henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
2180                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
2181                  const CFList& MOD)
2182{
2183  CFArray bufFactors= CFArray (factors.length());
2184  int i= 0;
2185  CanonicalForm xToStart= power (F.mvar(), start);
2186  for (CFListIterator k= factors; k.hasItem(); k++, i++)
2187  {
2188    if (i == 0)
2189      bufFactors[i]= mod (k.getItem(), xToStart);
2190    else
2191      bufFactors[i]= k.getItem();
2192  }
2193  for (i= start; i < end; i++)
2194    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
2195
2196  CFListIterator k= factors;
2197  for (i= 0; i < factors.length(); k++, i++)
2198    k.getItem()= bufFactors [i];
2199  factors.removeFirst();
2200  return;
2201}
2202
2203CFList
2204henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
2205            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
2206            lNew)
2207{
2208  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
2209  int k= 0;
2210  CFArray bufFactors= CFArray (factors.length());
2211  for (CFListIterator i= factors; i.hasItem(); i++, k++)
2212  {
2213    if (k == 0)
2214      bufFactors[k]= LC (F.getLast(), 1);
2215    else
2216      bufFactors[k]= i.getItem();
2217  }
2218  CFList buf= factors;
2219  buf.removeFirst();
2220  buf.insert (LC (F.getLast(), 1));
2221  CFListIterator i= buf;
2222  i++;
2223  Variable y= F.getLast().mvar();
2224  Variable x= F.getFirst().mvar();
2225  CanonicalForm xToLOld= power (x, lOld);
2226  Pi [0]= mod (Pi[0], xToLOld);
2227  M (1, 1)= Pi [0];
2228  k= 1;
2229  if (i.hasItem())
2230    i++;
2231  for (i; i.hasItem(); i++, k++)
2232  {
2233    Pi [k]= mod (Pi [k], xToLOld);
2234    M (1, k + 1)= Pi [k];
2235  }
2236
2237  for (int d= 1; d < lNew; d++)
2238    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
2239  CFList result;
2240  for (k= 1; k < factors.length(); k++)
2241    result.append (bufFactors[k]);
2242  return result;
2243}
2244
2245CFList
2246henselLift (const CFList& eval, const CFList& factors, const int* l, const int
2247            lLength)
2248{
2249  CFList diophant;
2250  CFList buf= factors;
2251  buf.insert (LC (eval.getFirst(), 1));
2252  sortList (buf, Variable (1));
2253  CFArray Pi;
2254  CFMatrix M= CFMatrix (l[1], factors.length());
2255  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
2256  if (eval.length() == 2)
2257    return result;
2258  CFList MOD;
2259  for (int i= 0; i < 2; i++)
2260    MOD.append (power (Variable (i + 2), l[i]));
2261  CFListIterator j= eval;
2262  j++;
2263  CFList bufEval;
2264  bufEval.append (j.getItem());
2265  j++;
2266
2267  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
2268  {
2269    result.insert (LC (bufEval.getFirst(), 1));
2270    bufEval.append (j.getItem());
2271    M= CFMatrix (l[i], factors.length());
2272    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
2273    MOD.append (power (Variable (i + 2), l[i]));
2274    bufEval.removeFirst();
2275  }
2276  return result;
2277}
2278
2279void
2280henselStep122 (const CanonicalForm& F, const CFList& factors,
2281              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2282              CFArray& Pi, int j, const CFArray& LCs)
2283{
2284  Variable x= F.mvar();
2285  CanonicalForm xToJ= power (x, j);
2286
2287  CanonicalForm E;
2288  // compute the error
2289  if (degree (Pi [factors.length() - 2], x) > 0)
2290    E= F[j] - Pi [factors.length() - 2] [j];
2291  else
2292    E= F[j];
2293
2294  CFArray buf= CFArray (diophant.length());
2295
2296  int k= 0;
2297  CanonicalForm remainder;
2298  // actual lifting
2299  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
2300  {
2301    if (degree (bufFactors[k], x) > 0)
2302      remainder= modNTL (E, bufFactors[k] [0]);
2303    else
2304      remainder= modNTL (E, bufFactors[k]);
2305    buf[k]= mulNTL (i.getItem(), remainder);
2306    if (degree (bufFactors[k], x) > 0)
2307      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
2308    else
2309      buf[k]= modNTL (buf[k], bufFactors[k]); 
2310  }
2311
2312  for (k= 0; k < factors.length(); k++)
2313    bufFactors[k] += xToJ*buf[k];
2314
2315  // update Pi [0]
2316  int degBuf0= degree (bufFactors[0], x);
2317  int degBuf1= degree (bufFactors[1], x);
2318  if (degBuf0 > 0 && degBuf1 > 0)
2319  {
2320    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
2321    if (j + 2 <= M.rows())
2322      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
2323  }
2324  CanonicalForm uIZeroJ;
2325  if (degBuf0 > 0 && degBuf1 > 0)
2326    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
2327  else if (degBuf0 > 0)
2328    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
2329  else if (degBuf1 > 0)
2330    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
2331  else
2332    uIZeroJ= 0;
2333  Pi [0] += xToJ*uIZeroJ;
2334
2335  CFArray tmp= CFArray (factors.length() - 1);
2336  for (k= 0; k < factors.length() - 1; k++)
2337    tmp[k]= 0;
2338  CFIterator one, two;
2339  one= bufFactors [0];
2340  two= bufFactors [1];
2341  if (degBuf0 > 0 && degBuf1 > 0)
2342  {
2343    while (one.hasTerms() && one.exp() > j) one++;
2344    while (two.hasTerms() && two.exp() > j) two++;
2345    for (k= 1; k <= (int) ceil (j/2.0); k++)
2346    {
2347      if (k != j - k + 1)
2348      {
2349        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2350        {
2351          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
2352                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
2353          one++;
2354          two++;
2355        }
2356        else if (one.exp() == j - k + 1)
2357        {
2358          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
2359                      M (k + 1, 1);
2360          one++;
2361        }
2362        else if (two.exp() == j - k + 1)
2363        {
2364          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
2365                    M (k + 1, 1);
2366          two++;
2367        }
2368      }
2369      else
2370        tmp[0] += M (k + 1, 1);
2371    }
2372  }
2373
2374  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2375  {
2376    if (j + 2 <= M.rows())
2377      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2378                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
2379                         - M(1,1) - M (j + 2,1);
2380  }
2381  else if (degBuf0 >= j + 1)
2382  {
2383    if (degBuf1 > 0)
2384      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
2385    else
2386      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
2387  }
2388  else if (degBuf1 >= j + 1)
2389  {
2390    if (degBuf0 > 0)
2391      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
2392    else
2393      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
2394  }
2395
2396  Pi [0] += tmp[0]*xToJ*F.mvar();
2397
2398  /*// update Pi [l]
2399  int degPi, degBuf;
2400  for (int l= 1; l < factors.length() - 1; l++)
2401  {
2402    degPi= degree (Pi [l - 1], x);
2403    degBuf= degree (bufFactors[l + 1], x);
2404    if (degPi > 0 && degBuf > 0)
2405      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
2406    if (j == 1)
2407    {
2408      if (degPi > 0 && degBuf > 0)
2409        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
2410                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
2411                  M (1, l + 1));
2412      else if (degPi > 0)
2413        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
2414      else if (degBuf > 0)
2415        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
2416    }
2417    else
2418    {
2419      if (degPi > 0 && degBuf > 0)
2420      {
2421        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
2422        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
2423      }
2424      else if (degPi > 0)
2425        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
2426      else if (degBuf > 0)
2427      {
2428        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
2429        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
2430      }
2431      Pi[l] += xToJ*uIZeroJ;
2432    }
2433    one= bufFactors [l + 1];
2434    two= Pi [l - 1];
2435    if (two.exp() == j + 1)
2436    {
2437      if (degBuf > 0 && degPi > 0)
2438      {
2439          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
2440          two++;
2441      }
2442      else if (degPi > 0)
2443      {
2444          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
2445          two++;
2446      }
2447    }
2448    if (degBuf > 0 && degPi > 0)
2449    {
2450      for (k= 1; k <= (int) ceil (j/2.0); k++)
2451      {
2452        if (k != j - k + 1)
2453        {
2454          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2455          {
2456            tmp[l] += mulNTL ((bufFactors[l + 1][k]+one.coeff()),(Pi[l-1] [k] +
2457                        two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
2458            one++;
2459            two++;
2460          }
2461          else if (one.exp() == j - k + 1)
2462          {
2463            tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()),Pi[l - 1] [k]) -
2464                        M (k + 1, l + 1);
2465            one++;
2466          }
2467          else if (two.exp() == j - k + 1)
2468          {
2469            tmp[l] += mulNTL (bufFactors[l+1][k],(Pi[l-1] [k] + two.coeff())) -
2470                      M (k + 1, l + 1);
2471            two++;
2472          }
2473        }
2474        else
2475          tmp[l] += M (k + 1, l + 1);
2476      }
2477    }
2478    Pi[l] += tmp[l]*xToJ*F.mvar();
2479  }*/
2480  return;
2481}
2482
2483void
2484henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
2485              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
2486{
2487  if (sort)
2488    sortList (factors, Variable (1));
2489  Pi= CFArray (factors.length() - 2);
2490  CFList bufFactors2= factors;
2491  bufFactors2.removeFirst();
2492  diophant.removeFirst();
2493  CFListIterator iter= diophant;
2494  CanonicalForm s,t;
2495  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
2496  diophant= CFList();
2497  diophant.append (t);
2498  diophant.append (s);
2499  DEBOUTLN (cerr, "diophant= " << diophant);
2500
2501  CFArray bufFactors= CFArray (bufFactors2.length());
2502  int i= 0;
2503  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
2504    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
2505
2506  Variable x= F.mvar();
2507  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
2508  {
2509    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
2510    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
2511                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
2512  }
2513  else if (degree (bufFactors[0], x) > 0)
2514  {
2515    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
2516    Pi [0]= M (1, 1) +
2517            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
2518  }
2519  else if (degree (bufFactors[1], x) > 0)
2520  {
2521    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
2522    Pi [0]= M (1, 1) +
2523            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
2524  }
2525  else
2526  {
2527    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
2528    Pi [0]= M (1, 1);
2529  }
2530
2531  for (i= 1; i < l; i++)
2532    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
2533
2534  factors= CFList();
2535  for (i= 0; i < bufFactors.size(); i++)
2536    factors.append (bufFactors[i]);
2537  return;
2538}
2539
2540/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
2541/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
2542static inline
2543CFList
2544diophantine (const CFList& recResult, const CFList& factors, 
2545             const CFList& products, const CFList& M, const CanonicalForm& E)
2546{
2547  if (M.isEmpty())
2548  {
2549    CFList result;
2550    CFListIterator j= factors;
2551    CanonicalForm buf;
2552    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2553    {
2554      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
2555              "constant or univariate poly expected");
2556      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
2557              "constant or univariate poly expected");
2558      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
2559              "constant or univariate poly expected");
2560      buf= mulNTL (E, i.getItem());
2561      result.append (modNTL (buf, j.getItem()));
2562    }
2563    return result;
2564  }
2565  Variable y= M.getLast().mvar();
2566  CFList bufFactors= factors;
2567  for (CFListIterator i= bufFactors; i.hasItem(); i++)
2568    i.getItem()= mod (i.getItem(), y);
2569  CFList bufProducts= products;
2570  for (CFListIterator i= bufProducts; i.hasItem(); i++)
2571    i.getItem()= mod (i.getItem(), y);
2572  CFList buf= M;
2573  buf.removeLast();
2574  CanonicalForm bufE= mod (E, y);
2575  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2576                                      bufE);
2577
2578  CanonicalForm e= E;
2579  CFListIterator j= products;
2580  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2581    e -= mulMod (i.getItem(), j.getItem(), M);
2582
2583  CFList result= recDiophantine;
2584  int d= degree (M.getLast());
2585  CanonicalForm coeffE;
2586  for (int i= 1; i < d; i++)
2587  {
2588    if (degree (e, y) > 0)
2589      coeffE= e[i];
2590    else
2591      coeffE= 0;
2592    if (!coeffE.isZero())
2593    {
2594      CFListIterator k= result;
2595      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2596                                   coeffE);
2597      CFListIterator l= products;
2598      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2599      {
2600        k.getItem() += j.getItem()*power (y, i);
2601        e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
2602      }
2603    }
2604    if (e.isZero())
2605      break;
2606  }
2607  return result;
2608}
2609
2610// non monic case
2611static inline
2612CFList
2613multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
2614                     const CFList& recResult, const CFList& M, const int d)
2615{
2616  Variable y= F.mvar();
2617  CFList result;
2618  CFListIterator i;
2619  CanonicalForm e= 1;
2620  CFListIterator j= factors;
2621  CFList p;
2622  CFArray bufFactors= CFArray (factors.length());
2623  CanonicalForm yToD= power (y, d);
2624  int k= 0;
2625  for (CFListIterator i= factors; i.hasItem(); i++, k++)
2626    bufFactors [k]= i.getItem();
2627  CanonicalForm b;
2628  CFList buf= M;
2629  buf.removeLast();
2630  buf.append (yToD);
2631  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
2632  {
2633    b= 1;
2634    if (fdivides (bufFactors[k], F))
2635      b= F/bufFactors[k];
2636    else 
2637    {
2638      for (int l= 0; l < factors.length(); l++)
2639      {
2640        if (l == k)
2641          continue;
2642        else
2643        {
2644          b= mulMod (b, bufFactors[l], buf);
2645        }
2646      }
2647    }
2648    p.append (b);
2649  }
2650  j= p;
2651  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2652    e -= mulMod (i.getItem(), j.getItem(), M);
2653  if (e.isZero())
2654    return recResult;
2655  CanonicalForm coeffE;
2656  result= recResult;
2657  CanonicalForm g;
2658  buf.removeLast();
2659  for (int i= 1; i < d; i++)
2660  {
2661    if (degree (e, y) > 0)
2662      coeffE= e[i];
2663    else
2664      coeffE= 0;
2665    if (!coeffE.isZero())
2666    {
2667      CFListIterator k= result;
2668      CFListIterator l= p;
2669      j= recResult;
2670      int ii= 0;
2671      CanonicalForm dummy;
2672      CFList recDiophantine;
2673      CFList buf2, buf3;
2674      buf2= factors;
2675      buf3= p;
2676      for (CFListIterator iii= buf2; iii.hasItem(); iii++)
2677        iii.getItem()= mod (iii.getItem(), y);
2678      for (CFListIterator iii= buf3; iii.hasItem(); iii++)
2679        iii.getItem()= mod (iii.getItem(), y);
2680      recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
2681      CFListIterator iter= recDiophantine;
2682      for (; j.hasItem(); j++, k++, l++, ii++, iter++)
2683      {
2684        k.getItem() += iter.getItem()*power (y, i);
2685        e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
2686      }
2687    }
2688    if (e.isZero())
2689      break;
2690  }
2691
2692#ifdef DEBUGOUTPUT
2693    CanonicalForm test= 0;
2694    j= p; 
2695    for (CFListIterator i= result; i.hasItem(); i++, j++)
2696      test += mod (i.getItem()*j.getItem(), power (y, d));
2697    DEBOUTLN (cerr, "test= " << test);
2698#endif
2699  return result;
2700}
2701
2702// so far only tested for two! factor Hensel lifting
2703static inline
2704void
2705henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
2706            const CFList& diophant, CFMatrix& M, CFArray& Pi,
2707            const CFList& products, int j, const CFList& MOD)
2708{
2709  CanonicalForm E;
2710  CanonicalForm xToJ= power (F.mvar(), j);
2711  Variable x= F.mvar();
2712
2713  // compute the error
2714#ifdef DEBUGOUTPUT
2715    CanonicalForm test= 1;
2716    for (int i= 0; i < factors.length(); i++)
2717    {
2718      if (i == 0)
2719        test *= mod (bufFactors [i], power (x, j));
2720      else
2721        test *= bufFactors[i];
2722    }
2723    test= mod (test, power (x, j));
2724    test= mod (test, MOD);
2725    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2726    DEBOUTLN (cerr, "test= " << test2);
2727#endif
2728
2729  if (degree (Pi [factors.length() - 2], x) > 0)
2730    E= F[j] - Pi [factors.length() - 2] [j];
2731  else
2732    E= F[j];
2733
2734  CFArray buf= CFArray (diophant.length());
2735
2736  // actual lifting
2737  CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
2738
2739  int k= 0;
2740  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2741  {
2742    buf[k]= i.getItem();
2743    bufFactors[k] += xToJ*i.getItem();
2744  }
2745
2746  // update Pi [0]
2747  int degBuf0= degree (bufFactors[0], x);
2748  int degBuf1= degree (bufFactors[1], x);
2749  if (degBuf0 > 0 && degBuf1 > 0)
2750  {
2751    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2752    if (j + 2 <= M.rows())
2753      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2754  }
2755  CanonicalForm uIZeroJ;
2756
2757  if (degBuf0 > 0 && degBuf1 > 0)
2758    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2759             mulMod (bufFactors[1] [0], buf[0], MOD);
2760  else if (degBuf0 > 0)
2761    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
2762  else if (degBuf1 > 0)
2763    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2764  else
2765    uIZeroJ= 0;
2766  Pi [0] += xToJ*uIZeroJ; 
2767
2768  CFArray tmp= CFArray (factors.length() - 1);
2769  for (k= 0; k < factors.length() - 1; k++)
2770    tmp[k]= 0;
2771  CFIterator one, two;
2772  one= bufFactors [0];
2773  two= bufFactors [1];
2774  if (degBuf0 > 0 && degBuf1 > 0)
2775  {
2776    while (one.hasTerms() && one.exp() > j) one++;
2777    while (two.hasTerms() && two.exp() > j) two++;
2778    for (k= 1; k <= (int) ceil (j/2.0); k++)
2779    {
2780      if (k != j - k + 1)
2781      {
2782        if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2783        {
2784          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2785                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2786                    M (j - k + 2, 1);
2787          one++;
2788          two++;
2789        }
2790        else if (one.exp() == j - k + 1)
2791        {
2792          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2793                    bufFactors[1] [k], MOD) - M (k + 1, 1);
2794          one++;
2795        }
2796        else if (two.exp() == j - k + 1)
2797        {
2798          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2799                    two.coeff()), MOD) - M (k + 1, 1);
2800          two++;
2801        }
2802      }
2803      else
2804      {
2805        tmp[0] += M (k + 1, 1);
2806      }
2807    }
2808  }
2809
2810  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2811  {
2812    if (j + 2 <= M.rows())
2813      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2814                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
2815                         - M(1,1) - M (j + 2,1);
2816  }
2817  else if (degBuf0 >= j + 1)
2818  {
2819    if (degBuf1 > 0)
2820      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
2821    else
2822      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
2823  }
2824  else if (degBuf1 >= j + 1)
2825  {
2826    if (degBuf0 > 0)
2827      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
2828    else
2829      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
2830  }
2831  Pi [0] += tmp[0]*xToJ*F.mvar();
2832
2833  // update Pi [l]
2834  int degPi, degBuf;
2835  for (int l= 1; l < factors.length() - 1; l++)
2836  {
2837    degPi= degree (Pi [l - 1], x);
2838    degBuf= degree (bufFactors[l + 1], x);
2839    if (degPi > 0 && degBuf > 0)
2840      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2841    if (j == 1)
2842    {
2843      if (degPi > 0 && degBuf > 0)
2844        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
2845                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
2846                  M (1, l + 1));
2847      else if (degPi > 0)
2848        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
2849      else if (degBuf > 0)
2850        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
2851    }
2852    else
2853    {
2854      if (degPi > 0 && degBuf > 0)
2855      {
2856        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2857        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
2858      }
2859      else if (degPi > 0)
2860        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
2861      else if (degBuf > 0)
2862      {
2863        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2864        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
2865      }
2866      Pi[l] += xToJ*uIZeroJ;
2867    }
2868    one= bufFactors [l + 1];
2869    two= Pi [l - 1];
2870    if (two.exp() == j + 1)
2871    {
2872      if (degBuf > 0 && degPi > 0)
2873      {
2874          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
2875          two++;
2876      }
2877      else if (degPi > 0)
2878      {
2879          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
2880          two++;
2881      }
2882    }
2883    if (degBuf > 0 && degPi > 0)
2884    {
2885      for (k= 1; k <= (int) ceil (j/2.0); k++)
2886      {
2887        if (k != j - k + 1)
2888        {
2889          if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
2890          {
2891            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2892                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2893                      M (j - k + 2, l + 1);
2894            one++;
2895            two++;
2896          }
2897          else if (one.exp() == j - k + 1)
2898          {
2899            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2900                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2901            one++;
2902          }
2903          else if (two.exp() == j - k + 1)
2904          {
2905            tmp[l] += mulMod (bufFactors[l + 1] [k],
2906                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2907            two++;
2908          }
2909        }
2910        else
2911          tmp[l] += M (k + 1, l + 1);
2912      }
2913    }
2914    Pi[l] += tmp[l]*xToJ*F.mvar();
2915  }
2916
2917  return;
2918}
2919
2920// wrt. Variable (1)
2921CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
2922{
2923  if (degree (F, 1) <= 0)
2924   return c;
2925  else
2926  {
2927    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
2928    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
2929              - LC (result))*power (result.mvar(), degree (result));
2930    return swapvar (result, Variable (F.level() + 1), Variable (1));
2931  }
2932}
2933
2934// so far only for two factor Hensel lifting
2935CFList
2936henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
2937              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
2938              const CFList& LCs2)
2939{
2940  CFList buf= factors;
2941  int k= 0;
2942  int liftBoundBivar= l[k];
2943  CFList bufbuf= factors;
2944  Variable v= Variable (2);
2945
2946  CFList MOD;
2947  MOD.append (power (Variable (2), liftBoundBivar));
2948  CFArray bufFactors= CFArray (factors.length());
2949  k= 0;
2950  CFListIterator j= eval;
2951  j++;
2952  CFListIterator iter1= LCs1;
2953  CFListIterator iter2= LCs2;
2954  iter1++;
2955  iter2++;
2956  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
2957  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
2958
2959  CFListIterator i= buf;
2960  i++;
2961  Variable y= j.getItem().mvar();
2962  if (y.level() != 3)
2963    y= Variable (3);
2964
2965  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
2966  M (1, 1)= Pi[0];
2967  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2968    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2969                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2970  else if (degree (bufFactors[0], y) > 0)
2971    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2972  else if (degree (bufFactors[1], y) > 0)
2973    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2974
2975  CFList products;
2976  for (int i= 0; i < bufFactors.size(); i++)
2977  {
2978    if (degree (bufFactors[i], y) > 0)
2979      products.append (M (1, 1)/bufFactors[i] [0]);
2980    else
2981      products.append (M (1, 1)/bufFactors[i]);
2982  }
2983
2984  for (int d= 1; d < l[1]; d++)
2985    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
2986  CFList result;
2987  for (k= 0; k < factors.length(); k++)
2988    result.append (bufFactors[k]);
2989  return result;
2990}
2991
2992
2993CFList
2994henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
2995            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
2996            lNew, const CFList& LCs1, const CFList& LCs2)
2997{
2998  int k= 0;
2999  CFArray bufFactors= CFArray (factors.length());
3000  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
3001  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
3002  CFList buf= factors;
3003  Variable y= F.getLast().mvar();
3004  Variable x= F.getFirst().mvar();
3005  CanonicalForm xToLOld= power (x, lOld);
3006  Pi [0]= mod (Pi[0], xToLOld);
3007  M (1, 1)= Pi [0];
3008
3009  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
3010    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) + 
3011                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
3012  else if (degree (bufFactors[0], y) > 0)
3013    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
3014  else if (degree (bufFactors[1], y) > 0)
3015    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
3016
3017  CFList products;
3018  for (int i= 0; i < bufFactors.size(); i++)
3019  {
3020    if (degree (bufFactors[i], y) > 0)
3021    {
3022      ASSERT (fdivides (bufFactors[i][0], M (1, 1)), "expected exact division");
3023      products.append (M (1, 1)/bufFactors[i] [0]);
3024    }
3025    else
3026    {
3027      ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
3028      products.append (M (1, 1)/bufFactors[i]);
3029    }
3030  }
3031
3032  for (int d= 1; d < lNew; d++)
3033    henselStep2 (F.getLast(),buf,bufFactors, diophant, M, Pi, products, d, MOD);
3034
3035  CFList result;
3036  for (k= 0; k < factors.length(); k++)
3037    result.append (bufFactors[k]);
3038  return result;
3039}
3040
3041// so far only for two! factor Hensel lifting
3042CFList
3043henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
3044             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
3045             const CFArray& Pi, const CFList& diophant)
3046{
3047  CFList bufDiophant= diophant;
3048  CFList buf= factors;
3049  if (sort)
3050    sortList (buf, Variable (1));
3051  CFArray bufPi= Pi;
3052  CFMatrix M= CFMatrix (l[1], factors.length());
3053  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
3054  if (eval.length() == 2)
3055    return result;
3056  CFList MOD;
3057  for (int i= 0; i < 2; i++)
3058    MOD.append (power (Variable (i + 2), l[i]));
3059  CFListIterator j= eval;
3060  j++;
3061  CFList bufEval;
3062  bufEval.append (j.getItem());
3063  j++;
3064  CFListIterator jj= LCs1;
3065  CFListIterator jjj= LCs2;
3066  CFList bufLCs1, bufLCs2;
3067  jj++, jjj++;
3068  bufLCs1.append (jj.getItem());
3069  bufLCs2.append (jjj.getItem());
3070  jj++, jjj++;
3071
3072  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
3073  {
3074    bufEval.append (j.getItem());
3075    bufLCs1.append (jj.getItem());
3076    bufLCs2.append (jjj.getItem());
3077    M= CFMatrix (l[i], factors.length());
3078    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
3079                         l[i], bufLCs1, bufLCs2);
3080    MOD.append (power (Variable (i + 2), l[i]));
3081    bufEval.removeFirst();
3082    bufLCs1.removeFirst();
3083    bufLCs2.removeFirst();
3084  }
3085  return result;
3086}
3087
3088#endif
3089/* HAVE_NTL */
3090
Note: See TracBrowser for help on using the repository browser.