Context Navigation

source: git/factory/facHensel.cc @ c1b52b

spielwiese

Last change on this file since c1b52b was b28747, checked in by Hans Schoenemann <hannes@…>, 4 years ago
fix: factrory with or w/o NTL, FLINT
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	void out_cf(const char s1,const CanonicalForm &f,const char s2);
45
46	TIMING_DEFINE_PRINT (diotime)
47	TIMING_DEFINE_PRINT (product1)
48	TIMING_DEFINE_PRINT (product2)
49	TIMING_DEFINE_PRINT (hensel23)
50	TIMING_DEFINE_PRINT (hensel)
51
52	#if defined (HAVE_NTL) \|\| defined(HAVE_FLINT)
53
54	#if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
55	static
56	CFList productsNTL (const CFList& factors, const CanonicalForm& M)
57	{
58	if (fac_NTL_char != getCharacteristic())
59	{
60	fac_NTL_char= getCharacteristic();
61	zz_p::init (getCharacteristic());
62	}
63	zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
64	zz_pE::init (NTLMipo);
65	zz_pEX prod;
66	vec_zz_pEX v;
67	v.SetLength (factors.length());
68	int j= 0;
69	for (CFListIterator i= factors; i.hasItem(); i++, j++)
70	{
71	if (i.getItem().inCoeffDomain())
72	v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem())));
73	else
74	v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
75	}
76	CFList result;
77	Variable x= Variable (1);
78	for (int j= 0; j < factors.length(); j++)
79	{
80	int k= 0;
81	set(prod);
82	for (int i= 0; i < factors.length(); i++, k++)
83	{
84	if (k == j)
85	continue;
86	prod *= v[i];
87	}
88	result.append (convertNTLzz_pEX2CF (prod, x, M.mvar()));
89	}
90	return result;
91	}
92	#endif
93
94	#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
95	static
96	CFList productsFLINT (const CFList& factors, const CanonicalForm& M)
97	{
98	nmod_poly_t FLINTmipo;
99	fq_nmod_ctx_t fq_con;
100	fq_nmod_poly_t prod;
101	fq_nmod_t buf;
102
103	nmod_poly_init (FLINTmipo, getCharacteristic());
104	convertFacCF2nmod_poly_t (FLINTmipo, M);
105
106	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
107
108	fq_nmod_poly_t * vec=new fq_nmod_poly_t [factors.length()];
109
110	int j= 0;
111
112	for (CFListIterator i= factors; i.hasItem(); i++, j++)
113	{
114	if (i.getItem().inCoeffDomain())
115	{
116	fq_nmod_poly_init (vec[j], fq_con);
117	fq_nmod_init2 (buf, fq_con);
118	convertFacCF2Fq_nmod_t (buf, i.getItem(), fq_con);
119	fq_nmod_poly_set_coeff (vec[j], 0, buf, fq_con);
120	fq_nmod_clear (buf, fq_con);
121	}
122	else
123	convertFacCF2Fq_nmod_poly_t (vec[j], i.getItem(), fq_con);
124	}
125
126	CFList result;
127	Variable x= Variable (1);
128	fq_nmod_poly_init (prod, fq_con);
129	for (j= 0; j < factors.length(); j++)
130	{
131	int k= 0;
132	fq_nmod_poly_one (prod, fq_con);
133	for (int i= 0; i < factors.length(); i++, k++)
134	{
135	if (k == j)
136	continue;
137	fq_nmod_poly_mul (prod, prod, vec[i], fq_con);
138	}
139	result.append (convertFq_nmod_poly_t2FacCF (prod, x, M.mvar(), fq_con));
140	}
141	for (j= 0; j < factors.length(); j++)
142	fq_nmod_poly_clear (vec[j], fq_con);
143
144	nmod_poly_clear (FLINTmipo);
145	fq_nmod_poly_clear (prod, fq_con);
146	fq_nmod_ctx_clear (fq_con);
147	delete [] vec;
148	return result;
149	}
150	#endif
151
152	static
153	void tryDiophantine (CFList& result, const CanonicalForm& F,
154	const CFList& factors, const CanonicalForm& M, bool& fail)
155	{
156	ASSERT (M.isUnivariate(), "expected univariate poly");
157
158	CFList bufFactors= factors;
159	bufFactors.removeFirst();
160	bufFactors.insert (factors.getFirst () (0,2));
161	CanonicalForm inv, leadingCoeff= Lc (F);
162	CFListIterator i= bufFactors;
163	if (bufFactors.getFirst().inCoeffDomain())
164	{
165	if (i.hasItem())
166	i++;
167	}
168	for (; i.hasItem(); i++)
169	{
170	tryInvert (Lc (i.getItem()), M, inv ,fail);
171	if (fail)
172	return;
173	i.getItem()= reduce (i.getItem()*inv, M);
174	}
175	#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
176	bufFactors= productsFLINT (bufFactors, M);
177	#else
178	bufFactors= productsNTL (bufFactors, M);
179	#endif
180
181	CanonicalForm buf1, buf2, buf3, S, T;
182	i= bufFactors;
183	if (i.hasItem())
184	i++;
185	buf1= bufFactors.getFirst();
186	buf2= i.getItem();
187	#ifdef HAVE_NTL
188	Variable x= Variable (1);
189	if (fac_NTL_char != getCharacteristic())
190	{
191	fac_NTL_char= getCharacteristic();
192	zz_p::init (getCharacteristic());
193	}
194	zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
195	zz_pE::init (NTLMipo);
196	zz_pEX NTLbuf1, NTLbuf2, NTLbuf3, NTLS, NTLT;
197	NTLbuf1= convertFacCF2NTLzz_pEX (buf1, NTLMipo);
198	NTLbuf2= convertFacCF2NTLzz_pEX (buf2, NTLMipo);
199	tryNTLXGCD (NTLbuf3, NTLS, NTLT, NTLbuf1, NTLbuf2, fail);
200	if (fail)
201	return;
202	S= convertNTLzz_pEX2CF (NTLS, x, M.mvar());
203	T= convertNTLzz_pEX2CF (NTLT, x, M.mvar());
204	#else
205	tryExtgcd (buf1, buf2, M, buf3, S, T, fail);
206	if (fail)
207	return;
208	#endif
209	result.append (S);
210	result.append (T);
211	if (i.hasItem())
212	i++;
213	for (; i.hasItem(); i++)
214	{
215	#ifdef HAVE_NTL
216	NTLbuf1= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
217	tryNTLXGCD (NTLbuf3, NTLS, NTLT, NTLbuf3, NTLbuf1, fail);
218	if (fail)
219	return;
220	S= convertNTLzz_pEX2CF (NTLS, x, M.mvar());
221	T= convertNTLzz_pEX2CF (NTLT, x, M.mvar());
222	#else
223	buf1= i.getItem();
224	tryExtgcd (buf3, buf1, M, buf3, S, T, fail);
225	if (fail)
226	return;
227	#endif
228	CFListIterator k= factors;
229	for (CFListIterator j= result; j.hasItem(); j++, k++)
230	{
231	j.getItem() *= S;
232	j.getItem()= mod (j.getItem(), k.getItem());
233	j.getItem()= reduce (j.getItem(), M);
234	}
235	result.append (T);
236	}
237	}
238
239	static
240	CFList mapinto (const CFList& L)
241	{
242	CFList result;
243	for (CFListIterator i= L; i.hasItem(); i++)
244	result.append (mapinto (i.getItem()));
245	return result;
246	}
247
248	static
249	int mod (const CFList& L, const CanonicalForm& p)
250	{
251	for (CFListIterator i= L; i.hasItem(); i++)
252	{
253	if (mod (i.getItem(), p) == 0)
254	return 0;
255	}
256	return 1;
257	}
258
259
260	static void
261	chineseRemainder (const CFList & x1, const CanonicalForm & q1,
262	const CFList & x2, const CanonicalForm & q2,
263	CFList & xnew, CanonicalForm & qnew)
264	{
265	ASSERT (x1.length() == x2.length(), "expected lists of equal length");
266	CanonicalForm tmp1, tmp2;
267	CFListIterator j= x2;
268	for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++)
269	{
270	chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2);
271	xnew.append (tmp1);
272	}
273	qnew= tmp2;
274	}
275
276	static
277	CFList Farey (const CFList& L, const CanonicalForm& q)
278	{
279	CFList result;
280	for (CFListIterator i= L; i.hasItem(); i++)
281	result.append (Farey (i.getItem(), q));
282	return result;
283	}
284
285	static
286	CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
287	{
288	CFList result;
289	for (CFListIterator i= L; i.hasItem(); i++)
290	result.append (replacevar (i.getItem(), a, b));
291	return result;
292	}
293
294	CFList
295	modularDiophant (const CanonicalForm& f, const CFList& factors,
296	const CanonicalForm& M)
297	{
298	bool save_rat=!isOn (SW_RATIONAL);
299	On (SW_RATIONAL);
300	CanonicalForm F= f*bCommonDen (f);
301	CFList products= factors;
302	for (CFListIterator i= products; i.hasItem(); i++)
303	{
304	if (products.getFirst().level() == 1)
305	i.getItem() /= Lc (i.getItem());
306	i.getItem() *= bCommonDen (i.getItem());
307	}
308	if (products.getFirst().level() == 1)
309	products.insert (Lc (F));
310	CanonicalForm bound= maxNorm (F);
311	CFList leadingCoeffs;
312	leadingCoeffs.append (lc (F));
313	CanonicalForm dummy;
314	for (CFListIterator i= products; i.hasItem(); i++)
315	{
316	leadingCoeffs.append (lc (i.getItem()));
317	dummy= maxNorm (i.getItem());
318	bound= (dummy > bound) ? dummy : bound;
319	}
320	bound = maxNorm (Lc (F))maxNorm (Lc(F))*bound;
321	bound = boundbound;
322	bound= power (bound, degree (M));
323	bound *= power (CanonicalForm (2),degree (f));
324	CanonicalForm bufBound= bound;
325	int i = cf_getNumBigPrimes() - 1;
326	int p;
327	CFList resultModP, result, newResult;
328	CanonicalForm q (0), newQ;
329	bool fail= false;
330	Variable a= M.mvar();
331	Variable b= Variable (2);
332	setReduce (M.mvar(), false);
333	CanonicalForm mipo= bCommonDen (M)*M;
334	Off (SW_RATIONAL);
335	CanonicalForm modMipo;
336	leadingCoeffs.append (lc (mipo));
337	CFList tmp1, tmp2;
338	bool equal= false;
339	int count= 0;
340	do
341	{
342	p = cf_getBigPrime( i );
343	i--;
344	while ( i >= 0 && mod( leadingCoeffs, p ) == 0)
345	{
346	p = cf_getBigPrime( i );
347	i--;
348	}
349
350	ASSERT (i >= 0, "ran out of primes"); //sic
351
352	setCharacteristic (p);
353	modMipo= mapinto (mipo);
354	modMipo /= lc (modMipo);
355	resultModP= CFList();
356	tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail);
357	setCharacteristic (0);
358	if (fail)
359	{
360	fail= false;
361	continue;
362	}
363
364	if ( q.isZero() )
365	{
366	result= replacevar (mapinto(resultModP), a, b);
367	q= p;
368	}
369	else
370	{
371	result= replacevar (result, a, b);
372	newResult= CFList();
373	chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b),
374	p, newResult, newQ );
375	q= newQ;
376	result= newResult;
377	if (newQ > bound)
378	{
379	count++;
380	tmp1= replacevar (Farey (result, q), b, a);
381	if (tmp2.isEmpty())
382	tmp2= tmp1;
383	else
384	{
385	equal= true;
386	CFListIterator k= tmp1;
387	for (CFListIterator j= tmp2; j.hasItem(); j++, k++)
388	{
389	if (j.getItem() != k.getItem())
390	equal= false;
391	}
392	if (!equal)
393	tmp2= tmp1;
394	}
395	if (count > 2)
396	{
397	bound *= bufBound;
398	equal= false;
399	count= 0;
400	}
401	}
402	if (newQ > bound && equal)
403	{
404	On( SW_RATIONAL );
405	CFList bufResult= result;
406	result= tmp2;
407	setReduce (M.mvar(), true);
408	if (factors.getFirst().level() == 1)
409	{
410	result.removeFirst();
411	CFListIterator j= factors;
412	CanonicalForm denf= bCommonDen (f);
413	for (CFListIterator i= result; i.hasItem(); i++, j++)
414	i.getItem() = Lc (j.getItem())denf;
415	}
416	if (factors.getFirst().level() != 1 &&
417	!bCommonDen (factors.getFirst()).isOne())
418	{
419	CanonicalForm denFirst= bCommonDen (factors.getFirst());
420	for (CFListIterator i= result; i.hasItem(); i++)
421	i.getItem() *= denFirst;
422	}
423
424	CanonicalForm test= 0;
425	CFListIterator jj= factors;
426	for (CFListIterator ii= result; ii.hasItem(); ii++, jj++)
427	test += ii.getItem()*(f/jj.getItem());
428	if (!test.isOne())
429	{
430	bound *= bufBound;
431	equal= false;
432	count= 0;
433	setReduce (M.mvar(), false);
434	result= bufResult;
435	Off (SW_RATIONAL);
436	}
437	else
438	break;
439	}
440	}
441	} while (1);
442	if (save_rat) Off(SW_RATIONAL);
443	return result;
444	}
445
446	void sortList (CFList& list, const Variable& x)
447	{
448	int l= 1;
449	int k= 1;
450	CanonicalForm buf;
451	CFListIterator m;
452	for (CFListIterator i= list; l <= list.length(); i++, l++)
453	{
454	for (CFListIterator j= list; k <= list.length() - l; k++)
455	{
456	m= j;
457	m++;
458	if (degree (j.getItem(), x) > degree (m.getItem(), x))
459	{
460	buf= m.getItem();
461	m.getItem()= j.getItem();
462	j.getItem()= buf;
463	j++;
464	j.getItem()= m.getItem();
465	}
466	else
467	j++;
468	}
469	k= 1;
470	}
471	}
472
473	CFList
474	diophantine (const CanonicalForm& F, const CFList& factors);
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	S= convertNTLZZ_pEX2CF (NTLS, x, gamma);
883	CFListIterator k= bufFactors;
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	/// solve \f$ E=\sum_{i= 1}^r{\sigma_{i}\prod_{j=1, j\neq i}^{r}{f_{j}}} \f$
2120	/// mod M, @a products contains \f$ \prod_{j=1, j\neq i}^{r}{f_{j}} \f$
2121	CFList
2122	diophantine (const CFList& recResult, const CFList& factors,
2123	const CFList& products, const CFList& M, const CanonicalForm& E,
2124	bool& bad)
2125	{
2126	if (M.isEmpty())
2127	{
2128	CFList result;
2129	CFListIterator j= factors;
2130	CanonicalForm buf;
2131	for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2132	{
2133	ASSERT (E.isUnivariate() \|\| E.inCoeffDomain(),
2134	"constant or univariate poly expected");
2135	ASSERT (i.getItem().isUnivariate() \|\| i.getItem().inCoeffDomain(),
2136	"constant or univariate poly expected");
2137	ASSERT (j.getItem().isUnivariate() \|\| j.getItem().inCoeffDomain(),
2138	"constant or univariate poly expected");
2139	buf= mulNTL (E, i.getItem());
2140	result.append (modNTL (buf, j.getItem()));
2141	}
2142	return result;
2143	}
2144	Variable y= M.getLast().mvar();
2145	CFList bufFactors= factors;
2146	for (CFListIterator i= bufFactors; i.hasItem(); i++)
2147	i.getItem()= mod (i.getItem(), y);
2148	CFList bufProducts= products;
2149	for (CFListIterator i= bufProducts; i.hasItem(); i++)
2150	i.getItem()= mod (i.getItem(), y);
2151	CFList buf= M;
2152	buf.removeLast();
2153	CanonicalForm bufE= mod (E, y);
2154	CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2155	bufE, bad);
2156
2157	if (bad)
2158	return CFList();
2159
2160	CanonicalForm e= E;
2161	CFListIterator j= products;
2162	for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2163	e -= j.getItem()*i.getItem();
2164
2165	CFList result= recDiophantine;
2166	int d= degree (M.getLast());
2167	CanonicalForm coeffE;
2168	for (int i= 1; i < d; i++)
2169	{
2170	if (degree (e, y) > 0)
2171	coeffE= e[i];
2172	else
2173	coeffE= 0;
2174	if (!coeffE.isZero())
2175	{
2176	CFListIterator k= result;
2177	recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2178	coeffE, bad);
2179	if (bad)
2180	return CFList();
2181	CFListIterator l= products;
2182	for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2183	{
2184	k.getItem() += j.getItem()*power (y, i);
2185	e -= l.getItem()(j.getItem()power (y, i));
2186	}
2187	}
2188	if (e.isZero())
2189	break;
2190	}
2191	if (!e.isZero())
2192	{
2193	bad= true;
2194	return CFList();
2195	}
2196	return result;
2197	}
2198
2199	#ifdef HAVE_NTL // diophantine
2200	void
2201	nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
2202	CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2203	CFArray& Pi, const CFList& products, int j,
2204	const CFList& MOD, bool& noOneToOne)
2205	{
2206	CanonicalForm E;
2207	CanonicalForm xToJ= power (F.mvar(), j);
2208	Variable x= F.mvar();
2209
2210	// compute the error
2211	#ifdef DEBUGOUTPUT
2212	CanonicalForm test= 1;
2213	for (int i= 0; i < factors.length(); i++)
2214	{
2215	if (i == 0)
2216	test *= mod (bufFactors [i], power (x, j));
2217	else
2218	test *= bufFactors[i];
2219	}
2220	test= mod (test, power (x, j));
2221	test= mod (test, MOD);
2222	CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2223	DEBOUTLN (cerr, "test in nonMonicHenselStep= " << test2);
2224	#endif
2225
2226	if (degree (Pi [factors.length() - 2], x) > 0)
2227	E= F[j] - Pi [factors.length() - 2] [j];
2228	else
2229	E= F[j];
2230
2231	CFArray buf= CFArray (diophant.length());
2232
2233	// actual lifting
2234	TIMING_START (diotime);
2235	CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
2236	noOneToOne);
2237
2238	if (noOneToOne)
2239	return;
2240
2241	int k= 0;
2242	for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2243	{
2244	buf[k]= i.getItem();
2245	bufFactors[k] += xToJ*i.getItem();
2246	}
2247	TIMING_END_AND_PRINT (diotime, "time for dio: ");
2248
2249	// update Pi [0]
2250	TIMING_START (product2);
2251	int degBuf0= degree (bufFactors[0], x);
2252	int degBuf1= degree (bufFactors[1], x);
2253	if (degBuf0 > 0 && degBuf1 > 0)
2254	{
2255	M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2256	if (j + 2 <= M.rows())
2257	M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2258	}
2259	else
2260	M (j + 1, 1)= 0;
2261
2262	CanonicalForm uIZeroJ;
2263	if (degBuf0 > 0 && degBuf1 > 0)
2264	uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2265	mulMod (bufFactors[1] [0], buf[0], MOD);
2266	else if (degBuf0 > 0)
2267	uIZeroJ= mulMod (buf[0], bufFactors[1], MOD) +
2268	mulMod (buf[1], bufFactors[0][0], MOD);
2269	else if (degBuf1 > 0)
2270	uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2271	mulMod (buf[0], bufFactors[1][0], MOD);
2272	else
2273	uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2274	mulMod (buf[0], bufFactors[1], MOD);
2275	Pi [0] += xToJ*uIZeroJ;
2276
2277	CFArray tmp= CFArray (factors.length() - 1);
2278	for (k= 0; k < factors.length() - 1; k++)
2279	tmp[k]= 0;
2280	CFIterator one, two;
2281	one= bufFactors [0];
2282	two= bufFactors [1];
2283	if (degBuf0 > 0 && degBuf1 > 0)
2284	{
2285	while (one.hasTerms() && one.exp() > j) one++;
2286	while (two.hasTerms() && two.exp() > j) two++;
2287	for (k= 1; k <= (j+1)/2; k++)
2288	{
2289	if (k != j - k + 1)
2290	{
2291	if ((one.hasTerms() && one.exp() == j - k + 1) &&
2292	(two.hasTerms() && two.exp() == j - k + 1))
2293	{
2294	tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2295	(bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2296	M (j - k + 2, 1);
2297	one++;
2298	two++;
2299	}
2300	else if (one.hasTerms() && one.exp() == j - k + 1)
2301	{
2302	tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2303	bufFactors[1] [k], MOD) - M (k + 1, 1);
2304	one++;
2305	}
2306	else if (two.hasTerms() && two.exp() == j - k + 1)
2307	{
2308	tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2309	two.coeff()), MOD) - M (k + 1, 1);
2310	two++;
2311	}
2312	}
2313	else
2314	{
2315	tmp[0] += M (k + 1, 1);
2316	}
2317	}
2318	}
2319
2320	if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2321	{
2322	if (j + 2 <= M.rows())
2323	tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2324	(bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
2325	- M(1,1) - M (j + 2,1);
2326	}
2327	else if (degBuf0 >= j + 1)
2328	{
2329	if (degBuf1 > 0)
2330	tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
2331	else
2332	tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
2333	}
2334	else if (degBuf1 >= j + 1)
2335	{
2336	if (degBuf0 > 0)
2337	tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
2338	else
2339	tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
2340	}
2341	Pi [0] += tmp[0]xToJF.mvar();
2342
2343	// update Pi [l]
2344	int degPi, degBuf;
2345	for (int l= 1; l < factors.length() - 1; l++)
2346	{
2347	degPi= degree (Pi [l - 1], x);
2348	degBuf= degree (bufFactors[l + 1], x);
2349	if (degPi > 0 && degBuf > 0)
2350	{
2351	M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2352	if (j + 2 <= M.rows())
2353	M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
2354	MOD);
2355	}
2356	else
2357	M (j + 1, l + 1)= 0;
2358
2359	if (degPi > 0 && degBuf > 0)
2360	uIZeroJ= mulMod (Pi[l - 1] [0], buf[l + 1], MOD) +
2361	mulMod (uIZeroJ, bufFactors[l + 1] [0], MOD);
2362	else if (degPi > 0)
2363	uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD) +
2364	mulMod (Pi[l - 1][0], buf[l + 1], MOD);
2365	else if (degBuf > 0)
2366	uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2367	mulMod (uIZeroJ, bufFactors[l + 1][0], MOD);
2368	else
2369	uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2370	mulMod (uIZeroJ, bufFactors[l + 1], MOD);
2371
2372	Pi [l] += xToJ*uIZeroJ;
2373
2374	one= bufFactors [l + 1];
2375	two= Pi [l - 1];
2376	if (degBuf > 0 && degPi > 0)
2377	{
2378	while (one.hasTerms() && one.exp() > j) one++;
2379	while (two.hasTerms() && two.exp() > j) two++;
2380	for (k= 1; k <= (j+1)/2; k++)
2381	{
2382	if (k != j - k + 1)
2383	{
2384	if ((one.hasTerms() && one.exp() == j - k + 1) &&
2385	(two.hasTerms() && two.exp() == j - k + 1))
2386	{
2387	tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2388	(Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2389	M (j - k + 2, l + 1);
2390	one++;
2391	two++;
2392	}
2393	else if (one.hasTerms() && one.exp() == j - k + 1)
2394	{
2395	tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2396	Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2397	one++;
2398	}
2399	else if (two.hasTerms() && two.exp() == j - k + 1)
2400	{
2401	tmp[l] += mulMod (bufFactors[l + 1] [k],
2402	(Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2403	two++;
2404	}
2405	}
2406	else
2407	tmp[l] += M (k + 1, l + 1);
2408	}
2409	}
2410
2411	if (degPi >= j + 1 && degBuf >= j + 1)
2412	{
2413	if (j + 2 <= M.rows())
2414	tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
2415	(bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
2416	, MOD) - M(1,l+1) - M (j + 2,l+1);
2417	}
2418	else if (degPi >= j + 1)
2419	{
2420	if (degBuf > 0)
2421	tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
2422	else
2423	tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
2424	}
2425	else if (degBuf >= j + 1)
2426	{
2427	if (degPi > 0)
2428	tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
2429	else
2430	tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
2431	}
2432
2433	Pi[l] += tmp[l]xToJF.mvar();
2434	}
2435	TIMING_END_AND_PRINT (product2, "time for product in hensel step: ");
2436	return;
2437	}
2438	#endif
2439
2440	// wrt. Variable (1)
2441	CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
2442	{
2443	if (degree (F, 1) <= 0)
2444	return c;
2445	else
2446	{
2447	CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
2448	result += (swapvar (c, Variable (F.level() + 1), Variable (1))
2449	- LC (result))*power (result.mvar(), degree (result));
2450	return swapvar (result, Variable (F.level() + 1), Variable (1));
2451	}
2452	}
2453
2454	#ifdef HAVE_NTL // nonMonicHenselStep
2455	CFList
2456	nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
2457	diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
2458	const CFList& LCs2, bool& bad)
2459	{
2460	CFList buf= factors;
2461	int k= 0;
2462	int liftBoundBivar= l[k];
2463	CFList bufbuf= factors;
2464	Variable v= Variable (2);
2465
2466	CFList MOD;
2467	MOD.append (power (Variable (2), liftBoundBivar));
2468	CFArray bufFactors= CFArray (factors.length());
2469	k= 0;
2470	CFListIterator j= eval;
2471	j++;
2472	CFListIterator iter1= LCs1;
2473	CFListIterator iter2= LCs2;
2474	iter1++;
2475	iter2++;
2476	bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
2477	bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
2478
2479	CFListIterator i= buf;
2480	i++;
2481	Variable y= j.getItem().mvar();
2482	if (y.level() != 3)
2483	y= Variable (3);
2484
2485	Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
2486	M (1, 1)= Pi[0];
2487	if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2488	Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2489	mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2490	else if (degree (bufFactors[0], y) > 0)
2491	Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2492	else if (degree (bufFactors[1], y) > 0)
2493	Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2494
2495	CFList products;
2496	for (int i= 0; i < bufFactors.size(); i++)
2497	{
2498	if (degree (bufFactors[i], y) > 0)
2499	products.append (eval.getFirst()/bufFactors[i] [0]);
2500	else
2501	products.append (eval.getFirst()/bufFactors[i]);
2502	}
2503
2504	for (int d= 1; d < l[1]; d++)
2505	{
2506	nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
2507	d, MOD, bad);
2508	if (bad)
2509	return CFList();
2510	}
2511	CFList result;
2512	for (k= 0; k < factors.length(); k++)
2513	result.append (bufFactors[k]);
2514	return result;
2515	}
2516	#endif
2517
2518	#ifdef HAVE_NTL // nonMonicHenselStep
2519	CFList
2520	nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
2521	CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2522	int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
2523	)
2524	{
2525	int k= 0;
2526	CFArray bufFactors= CFArray (factors.length());
2527	bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
2528	bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
2529	CFList buf= factors;
2530	Variable y= F.getLast().mvar();
2531	Variable x= F.getFirst().mvar();
2532	CanonicalForm xToLOld= power (x, lOld);
2533	Pi [0]= mod (Pi[0], xToLOld);
2534	M (1, 1)= Pi [0];
2535
2536	if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2537	Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2538	mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2539	else if (degree (bufFactors[0], y) > 0)
2540	Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2541	else if (degree (bufFactors[1], y) > 0)
2542	Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2543
2544	CFList products;
2545	CanonicalForm quot;
2546	for (int i= 0; i < bufFactors.size(); i++)
2547	{
2548	if (degree (bufFactors[i], y) > 0)
2549	{
2550	if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
2551	{
2552	bad= true;
2553	return CFList();
2554	}
2555	products.append (quot);
2556	}
2557	else
2558	{
2559	if (!fdivides (bufFactors[i], F.getFirst(), quot))
2560	{
2561	bad= true;
2562	return CFList();
2563	}
2564	products.append (quot);
2565	}
2566	}
2567
2568	for (int d= 1; d < lNew; d++)
2569	{
2570	nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
2571	d, MOD, bad);
2572	if (bad)
2573	return CFList();
2574	}
2575
2576	CFList result;
2577	for (k= 0; k < factors.length(); k++)
2578	result.append (bufFactors[k]);
2579	return result;
2580	}
2581	#endif
2582
2583	#ifdef HAVE_NTL // nonMonicHenselStep
2584	CFList
2585	nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
2586	lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
2587	const CFArray& Pi, const CFList& diophant, bool& bad)
2588	{
2589	CFList bufDiophant= diophant;
2590	CFList buf= factors;
2591	if (sort)
2592	sortList (buf, Variable (1));
2593	CFArray bufPi= Pi;
2594	CFMatrix M= CFMatrix (l[1], factors.length());
2595	CFList result=
2596	nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
2597	if (bad)
2598	return CFList();
2599
2600	if (eval.length() == 2)
2601	return result;
2602	CFList MOD;
2603	for (int i= 0; i < 2; i++)
2604	MOD.append (power (Variable (i + 2), l[i]));
2605	CFListIterator j= eval;
2606	j++;
2607	CFList bufEval;
2608	bufEval.append (j.getItem());
2609	j++;
2610	CFListIterator jj= LCs1;
2611	CFListIterator jjj= LCs2;
2612	CFList bufLCs1, bufLCs2;
2613	jj++, jjj++;
2614	bufLCs1.append (jj.getItem());
2615	bufLCs2.append (jjj.getItem());
2616	jj++, jjj++;
2617
2618	for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
2619	{
2620	bufEval.append (j.getItem());
2621	bufLCs1.append (jj.getItem());
2622	bufLCs2.append (jjj.getItem());
2623	M= CFMatrix (l[i], factors.length());
2624	result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
2625	l[i - 1], l[i], bufLCs1, bufLCs2, bad);
2626	if (bad)
2627	return CFList();
2628	MOD.append (power (Variable (i + 2), l[i]));
2629	bufEval.removeFirst();
2630	bufLCs1.removeFirst();
2631	bufLCs2.removeFirst();
2632	}
2633	return result;
2634	}
2635	#endif
2636
2637	#ifdef HAVE_NTL // diophantine
2638	CFList
2639	nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
2640	CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
2641	int bivarLiftBound, bool& bad)
2642	{
2643	CFList bufFactors2= factors;
2644
2645	Variable y= Variable (2);
2646	for (CFListIterator i= bufFactors2; i.hasItem(); i++)
2647	i.getItem()= mod (i.getItem(), y);
2648
2649	CanonicalForm bufF= F;
2650	bufF= mod (bufF, y);
2651	bufF= mod (bufF, Variable (3));
2652
2653	TIMING_START (diotime);
2654	diophant= diophantine (bufF, bufFactors2);
2655	TIMING_END_AND_PRINT (diotime, "time for dio in 23: ");
2656
2657	CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
2658
2659	Pi= CFArray (bufFactors2.length() - 1);
2660
2661	CFArray bufFactors= CFArray (bufFactors2.length());
2662	CFListIterator j= LCs;
2663	int i= 0;
2664	for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
2665	bufFactors[i]= replaceLC (k.getItem(), j.getItem());
2666
2667	//initialise Pi
2668	Variable v= Variable (3);
2669	CanonicalForm yToL= power (y, bivarLiftBound);
2670	if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
2671	{
2672	M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
2673	Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
2674	mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
2675	}
2676	else if (degree (bufFactors[0], v) > 0)
2677	{
2678	M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
2679	Pi [0]= M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
2680	}
2681	else if (degree (bufFactors[1], v) > 0)
2682	{
2683	M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
2684	Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
2685	}
2686	else
2687	{
2688	M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
2689	Pi [0]= M (1,1);
2690	}
2691
2692	for (i= 1; i < Pi.size(); i++)
2693	{
2694	if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
2695	{
2696	M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
2697	Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
2698	mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
2699	}
2700	else if (degree (Pi[i-1], v) > 0)
2701	{
2702	M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
2703	Pi [i]= M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
2704	}
2705	else if (degree (bufFactors[i+1], v) > 0)
2706	{
2707	M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
2708	Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
2709	}
2710	else
2711	{
2712	M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
2713	Pi [i]= M (1,i+1);
2714	}
2715	}
2716
2717	CFList products;
2718	bufF= mod (F, Variable (3));
2719	TIMING_START (product1);
2720	for (CFListIterator k= factors; k.hasItem(); k++)
2721	products.append (bufF/k.getItem());
2722	TIMING_END_AND_PRINT (product1, "time for product1 in 23: ");
2723
2724	CFList MOD= CFList (power (v, liftBound));
2725	MOD.insert (yToL);
2726	for (int d= 1; d < liftBound; d++)
2727	{
2728	nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
2729	MOD, bad);
2730	if (bad)
2731	return CFList();
2732	}
2733
2734	CFList result;
2735	for (i= 0; i < factors.length(); i++)
2736	result.append (bufFactors[i]);
2737	return result;
2738	}
2739	#endif
2740
2741	#ifdef HAVE_NTL // nonMonicHenselStep
2742	CFList
2743	nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
2744	CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2745	int& lNew, const CFList& MOD, bool& noOneToOne
2746	)
2747	{
2748
2749	int k= 0;
2750	CFArray bufFactors= CFArray (factors.length());
2751	CFListIterator j= LCs;
2752	for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
2753	bufFactors [k]= replaceLC (i.getItem(), j.getItem());
2754
2755	Variable y= F.getLast().mvar();
2756	Variable x= F.getFirst().mvar();
2757	CanonicalForm xToLOld= power (x, lOld);
2758
2759	Pi [0]= mod (Pi[0], xToLOld);
2760	M (1, 1)= Pi [0];
2761
2762	if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2763	Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2764	mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2765	else if (degree (bufFactors[0], y) > 0)
2766	Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2767	else if (degree (bufFactors[1], y) > 0)
2768	Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2769
2770	for (int i= 1; i < Pi.size(); i++)
2771	{
2772	Pi [i]= mod (Pi [i], xToLOld);
2773	M (1, i + 1)= Pi [i];
2774
2775	if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
2776	Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
2777	mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
2778	else if (degree (Pi[i-1], y) > 0)
2779	Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
2780	else if (degree (bufFactors[i+1], y) > 0)
2781	Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
2782	}
2783
2784	CFList products;
2785	CanonicalForm quot, bufF= F.getFirst();
2786
2787	TIMING_START (product1);
2788	for (int i= 0; i < bufFactors.size(); i++)
2789	{
2790	if (degree (bufFactors[i], y) > 0)
2791	{
2792	if (!fdivides (bufFactors[i] [0], bufF, quot))
2793	{
2794	noOneToOne= true;
2795	return factors;
2796	}
2797	products.append (quot);
2798	}
2799	else
2800	{
2801	if (!fdivides (bufFactors[i], bufF, quot))
2802	{
2803	noOneToOne= true;
2804	return factors;
2805	}
2806	products.append (quot);
2807	}
2808	}
2809	TIMING_END_AND_PRINT (product1, "time for product1 in further hensel:" );
2810
2811	for (int d= 1; d < lNew; d++)
2812	{
2813	nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
2814	products, d, MOD, noOneToOne);
2815	if (noOneToOne)
2816	return CFList();
2817	}
2818
2819	CFList result;
2820	for (k= 0; k < factors.length(); k++)
2821	result.append (bufFactors[k]);
2822	return result;
2823	}
2824	#endif
2825
2826	#ifdef HAVE_NTL // nonMonicHenselLift23
2827	CFList
2828	nonMonicHenselLift (const CFList& eval, const CFList& factors,
2829	CFList* const& LCs, CFList& diophant, CFArray& Pi,
2830	int* liftBound, int length, bool& noOneToOne
2831	)
2832	{
2833	CFList bufDiophant= diophant;
2834	CFList buf= factors;
2835	CFArray bufPi= Pi;
2836	CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
2837	int k= 0;
2838
2839	TIMING_START (hensel23);
2840	CFList result=
2841	nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
2842	liftBound[1], liftBound[0], noOneToOne);
2843	TIMING_END_AND_PRINT (hensel23, "time for 23: ");
2844
2845	if (noOneToOne)
2846	return CFList();
2847
2848	if (eval.length() == 1)
2849	return result;
2850
2851	k++;
2852	CFList MOD;
2853	for (int i= 0; i < 2; i++)
2854	MOD.append (power (Variable (i + 2), liftBound[i]));
2855
2856	CFListIterator j= eval;
2857	CFList bufEval;
2858	bufEval.append (j.getItem());
2859	j++;
2860
2861	for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
2862	{
2863	bufEval.append (j.getItem());
2864	M= CFMatrix (liftBound[i], factors.length() - 1);
2865	TIMING_START (hensel);
2866	result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
2867	liftBound[i-1], liftBound[i], MOD, noOneToOne);
2868	TIMING_END_AND_PRINT (hensel, "time for further hensel: ");
2869	if (noOneToOne)
2870	return result;
2871	MOD.append (power (Variable (i + 2), liftBound[i]));
2872	bufEval.removeFirst();
2873	}
2874
2875	return result;
2876	}
2877	#endif
2878	#endif

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

Download in other formats: