Context Navigation

source: git/factory/facHensel.cc @ c770dc

spielwiese

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: