Context Navigation

source: git/factory/facHensel.cc @ 650f2d8

spielwiese

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

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

Download in other formats: