Context Navigation

source: git/factory/facHensel.cc @ f542551

spielwiese

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

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

Download in other formats: