Context Navigation

source: git/factory/facHensel.cc @ 03c742

spielwiese

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

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

Download in other formats: