Context Navigation

source: git/factory/facHensel.cc @ f4d7641

fieker-DuValspielwiese

Last change on this file since f4d7641 was b27212, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: wrong index
Property mode set to `100644`
File size: 74.6 KB

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

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

Download in other formats: