Context Navigation

source: git/factory/cfModGcd.cc @ 767da0

spielwiese

Last change on this file since 767da0 was 4b24c19, checked in by Hans Schoenemann <hannes@…>, 4 years ago
Revert "use: chineseRemainderCached" This reverts commit 0aeeeeaed39a3ac804f59d6a4f58df7e1572b687.
Property mode set to `100644`
File size: 116.8 KB

Line
1	// -- c++ --
2	//*****************************************************************************
3	/** @file cfModGcd.cc
4	*
5	* This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
6	* \f$ F_{p}(\alpha ) \f$, GF or Z based on Alg. 7.1. and 7.2. as described in
7	* "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn via modular
8	* computations. And sparse modular variants as described in Zippel
9	* "Effective Polynomial Computation", deKleine, Monagan, Wittkopf
10	* "Algorithms for the non-monic case of the sparse modular GCD algorithm" and
11	* Javadi "A new solution to the normalization problem"
12	*
13	* @author Martin Lee
14	* @date 22.10.2009
15	*
16	* @par Copyright:
17	* (c) by The SINGULAR Team, see LICENSE file
18	*
19	**/
20	//*****************************************************************************
21
22
23	#include "config.h"
24
25
26	#include "cf_assert.h"
27	#include "debug.h"
28	#include "timing.h"
29
30	#include "canonicalform.h"
31	#include "cfGcdUtil.h"
32	#include "cf_map.h"
33	#include "cf_util.h"
34	#include "cf_irred.h"
35	#include "templates/ftmpl_functions.h"
36	#include "cf_random.h"
37	#include "cf_reval.h"
38	#include "facHensel.h"
39	#include "cf_iter.h"
40	#include "cfNewtonPolygon.h"
41	#include "cf_algorithm.h"
42	#include "cf_primes.h"
43
44	// inline helper functions:
45	#include "cf_map_ext.h"
46
47	#ifdef HAVE_NTL
48	#include "NTLconvert.h"
49	#endif
50
51	#ifdef HAVE_FLINT
52	#include "FLINTconvert.h"
53	#endif
54
55	#include "cfModGcd.h"
56
57	TIMING_DEFINE_PRINT(gcd_recursion)
58	TIMING_DEFINE_PRINT(newton_interpolation)
59	TIMING_DEFINE_PRINT(termination_test)
60	TIMING_DEFINE_PRINT(ez_p_compress)
61	TIMING_DEFINE_PRINT(ez_p_hensel_lift)
62	TIMING_DEFINE_PRINT(ez_p_eval)
63	TIMING_DEFINE_PRINT(ez_p_content)
64
65	bool
66	terminationTest (const CanonicalForm& F, const CanonicalForm& G,
67	const CanonicalForm& coF, const CanonicalForm& coG,
68	const CanonicalForm& cand)
69	{
70	CanonicalForm LCCand= abs (LC (cand));
71	if (LCCand*abs (LC (coF)) == abs (LC (F)))
72	{
73	if (LCCand*abs (LC (coG)) == abs (LC (G)))
74	{
75	if (abs (cand)*abs (coF) == abs (F))
76	{
77	if (abs (cand)*abs (coG) == abs (G))
78	return true;
79	}
80	return false;
81	}
82	return false;
83	}
84	return false;
85	}
86
87	#if defined(HAVE_NTL)\|\| defined(HAVE_FLINT)
88
89	static const double log2exp= 1.442695041;
90
91	/// compressing two polynomials F and G, M is used for compressing,
92	/// N to reverse the compression
93	int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
94	CFMap & N, bool topLevel)
95	{
96	int n= tmax (F.level(), G.level());
97	int * degsf= NEW_ARRAY(int,n + 1);
98	int * degsg= NEW_ARRAY(int,n + 1);
99
100	for (int i = n; i >= 0; i--)
101	degsf[i]= degsg[i]= 0;
102
103	degsf= degrees (F, degsf);
104	degsg= degrees (G, degsg);
105
106	int both_non_zero= 0;
107	int f_zero= 0;
108	int g_zero= 0;
109	int both_zero= 0;
110
111	if (topLevel)
112	{
113	for (int i= 1; i <= n; i++)
114	{
115	if (degsf[i] != 0 && degsg[i] != 0)
116	{
117	both_non_zero++;
118	continue;
119	}
120	if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
121	{
122	f_zero++;
123	continue;
124	}
125	if (degsg[i] == 0 && degsf[i] && i <= F.level())
126	{
127	g_zero++;
128	continue;
129	}
130	}
131
132	if (both_non_zero == 0)
133	{
134	DELETE_ARRAY(degsf);
135	DELETE_ARRAY(degsg);
136	return 0;
137	}
138
139	// map Variables which do not occur in both polynomials to higher levels
140	int k= 1;
141	int l= 1;
142	for (int i= 1; i <= n; i++)
143	{
144	if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
145	{
146	if (k + both_non_zero != i)
147	{
148	M.newpair (Variable (i), Variable (k + both_non_zero));
149	N.newpair (Variable (k + both_non_zero), Variable (i));
150	}
151	k++;
152	}
153	if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
154	{
155	if (l + g_zero + both_non_zero != i)
156	{
157	M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
158	N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
159	}
160	l++;
161	}
162	}
163
164	// sort Variables x_{i} in increasing order of
165	// min(deg_{x_{i}}(f),deg_{x_{i}}(g))
166	int m= tmax (F.level(), G.level());
167	int min_max_deg;
168	k= both_non_zero;
169	l= 0;
170	int i= 1;
171	while (k > 0)
172	{
173	if (degsf [i] != 0 && degsg [i] != 0)
174	min_max_deg= tmax (degsf[i], degsg[i]);
175	else
176	min_max_deg= 0;
177	while (min_max_deg == 0)
178	{
179	i++;
180	if (degsf [i] != 0 && degsg [i] != 0)
181	min_max_deg= tmax (degsf[i], degsg[i]);
182	else
183	min_max_deg= 0;
184	}
185	for (int j= i + 1; j <= m; j++)
186	{
187	if (degsf[j] != 0 && degsg [j] != 0 &&
188	tmax (degsf[j],degsg[j]) <= min_max_deg)
189	{
190	min_max_deg= tmax (degsf[j], degsg[j]);
191	l= j;
192	}
193	}
194	if (l != 0)
195	{
196	if (l != k)
197	{
198	M.newpair (Variable (l), Variable(k));
199	N.newpair (Variable (k), Variable(l));
200	degsf[l]= 0;
201	degsg[l]= 0;
202	l= 0;
203	}
204	else
205	{
206	degsf[l]= 0;
207	degsg[l]= 0;
208	l= 0;
209	}
210	}
211	else if (l == 0)
212	{
213	if (i != k)
214	{
215	M.newpair (Variable (i), Variable (k));
216	N.newpair (Variable (k), Variable (i));
217	degsf[i]= 0;
218	degsg[i]= 0;
219	}
220	else
221	{
222	degsf[i]= 0;
223	degsg[i]= 0;
224	}
225	i++;
226	}
227	k--;
228	}
229	}
230	else
231	{
232	//arrange Variables such that no gaps occur
233	for (int i= 1; i <= n; i++)
234	{
235	if (degsf[i] == 0 && degsg[i] == 0)
236	{
237	both_zero++;
238	continue;
239	}
240	else
241	{
242	if (both_zero != 0)
243	{
244	M.newpair (Variable (i), Variable (i - both_zero));
245	N.newpair (Variable (i - both_zero), Variable (i));
246	}
247	}
248	}
249	}
250
251	DELETE_ARRAY(degsf);
252	DELETE_ARRAY(degsg);
253
254	return 1;
255	}
256
257	static inline CanonicalForm
258	uni_content (const CanonicalForm & F);
259
260	CanonicalForm
261	uni_content (const CanonicalForm& F, const Variable& x)
262	{
263	if (F.inCoeffDomain())
264	return F.genOne();
265	if (F.level() == x.level() && F.isUnivariate())
266	return F;
267	if (F.level() != x.level() && F.isUnivariate())
268	return F.genOne();
269
270	if (x.level() != 1)
271	{
272	CanonicalForm f= swapvar (F, x, Variable (1));
273	CanonicalForm result= uni_content (f);
274	return swapvar (result, x, Variable (1));
275	}
276	else
277	return uni_content (F);
278	}
279
280	/// compute the content of F, where F is considered as an element of
281	/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
282	static inline CanonicalForm
283	uni_content (const CanonicalForm & F)
284	{
285	if (F.inBaseDomain())
286	return F.genOne();
287	if (F.level() == 1 && F.isUnivariate())
288	return F;
289	if (F.level() != 1 && F.isUnivariate())
290	return F.genOne();
291	if (degree (F,1) == 0) return F.genOne();
292
293	int l= F.level();
294	if (l == 2)
295	return content(F);
296	else
297	{
298	CanonicalForm pol, c = 0;
299	CFIterator i = F;
300	for (; i.hasTerms(); i++)
301	{
302	pol= i.coeff();
303	pol= uni_content (pol);
304	c= gcd (c, pol);
305	if (c.isOne())
306	return c;
307	}
308	return c;
309	}
310	}
311
312	CanonicalForm
313	extractContents (const CanonicalForm& F, const CanonicalForm& G,
314	CanonicalForm& contentF, CanonicalForm& contentG,
315	CanonicalForm& ppF, CanonicalForm& ppG, const int d)
316	{
317	CanonicalForm uniContentF, uniContentG, gcdcFcG;
318	contentF= 1;
319	contentG= 1;
320	ppF= F;
321	ppG= G;
322	CanonicalForm result= 1;
323	for (int i= 1; i <= d; i++)
324	{
325	uniContentF= uni_content (F, Variable (i));
326	uniContentG= uni_content (G, Variable (i));
327	gcdcFcG= gcd (uniContentF, uniContentG);
328	contentF *= uniContentF;
329	contentG *= uniContentG;
330	ppF /= uniContentF;
331	ppG /= uniContentG;
332	result *= gcdcFcG;
333	}
334	return result;
335	}
336
337	/// compute the leading coefficient of F, where F is considered as an element
338	/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on
339	/// \f$ Mon (x_{2},\ldots ,x_{n}) \f$ is dp.
340	static inline
341	CanonicalForm uni_lcoeff (const CanonicalForm& F)
342	{
343	if (F.level() > 1)
344	{
345	Variable x= Variable (2);
346	int deg= totaldegree (F, x, F.mvar());
347	for (CFIterator i= F; i.hasTerms(); i++)
348	{
349	if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
350	return uni_lcoeff (i.coeff());
351	}
352	}
353	return F;
354	}
355
356	/// Newton interpolation - Incremental algorithm.
357	/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
358	/// computes the interpolation polynomial assuming that
359	/// the polynomials in u are the results of evaluating the variabe x
360	/// of the unknown polynomial at the alpha values.
361	/// This incremental version receives only the values of alpha_n and u_n and
362	/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
363	/// the polynomial interpolating in all the points.
364	/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
365	static inline CanonicalForm
366	newtonInterp(const CanonicalForm & alpha, const CanonicalForm & u,
367	const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
368	const Variable & x)
369	{
370	CanonicalForm interPoly;
371
372	interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
373	*newtonPoly;
374	return interPoly;
375	}
376
377	/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
378	/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
379	/// fail if there are no field elements left which have not been used before
380	static inline CanonicalForm
381	randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
382	bool & fail)
383	{
384	fail= false;
385	Variable x= F.mvar();
386	AlgExtRandomF genAlgExt (alpha);
387	FFRandom genFF;
388	CanonicalForm random, mipo;
389	mipo= getMipo (alpha);
390	int p= getCharacteristic ();
391	int d= degree (mipo);
392	double bound= pow ((double) p, (double) d);
393	do
394	{
395	if (list.length() == bound)
396	{
397	fail= true;
398	break;
399	}
400	if (list.length() < p)
401	{
402	random= genFF.generate();
403	while (find (list, random))
404	random= genFF.generate();
405	}
406	else
407	{
408	random= genAlgExt.generate();
409	while (find (list, random))
410	random= genAlgExt.generate();
411	}
412	if (F (random, x) == 0)
413	{
414	list.append (random);
415	continue;
416	}
417	} while (find (list, random));
418	return random;
419	}
420
421	static inline
422	Variable chooseExtension (const Variable & alpha)
423	{
424	int i, m;
425	// extension of F_p needed
426	if (alpha.level() == 1)
427	{
428	i= 1;
429	m= 2;
430	} //extension of F_p(alpha)
431	if (alpha.level() != 1)
432	{
433	i= 4;
434	m= degree (getMipo (alpha));
435	}
436	#ifdef HAVE_FLINT
437	nmod_poly_t Irredpoly;
438	nmod_poly_init(Irredpoly,getCharacteristic());
439	nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
440	CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
441	nmod_poly_clear(Irredpoly);
442	#else
443	if (fac_NTL_char != getCharacteristic())
444	{
445	fac_NTL_char= getCharacteristic();
446	zz_p::init (getCharacteristic());
447	}
448	zz_pX NTLIrredpoly;
449	BuildIrred (NTLIrredpoly, i*m);
450	CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
451	#endif
452	return rootOf (newMipo);
453	}
454
455	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
456	CanonicalForm
457	modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
458	CanonicalForm& coF, CanonicalForm& coG,
459	Variable & alpha, CFList& l, bool& topLevel);
460	#endif
461
462	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
463	CanonicalForm
464	modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
465	Variable & alpha, CFList& l, bool& topLevel)
466	{
467	CanonicalForm dummy1, dummy2;
468	CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
469	topLevel);
470	return result;
471	}
472	#endif
473
474	/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
475	/// l and topLevel are only used internally, output is monic
476	/// based on Alg. 7.2. as described in "Algorithms for
477	/// Computer Algebra" by Geddes, Czapor, Labahn
478	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
479	CanonicalForm
480	modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
481	CanonicalForm& coF, CanonicalForm& coG,
482	Variable & alpha, CFList& l, bool& topLevel)
483	{
484	CanonicalForm A= F;
485	CanonicalForm B= G;
486	if (F.isZero() && degree(G) > 0)
487	{
488	coF= 0;
489	coG= Lc (G);
490	return G/Lc(G);
491	}
492	else if (G.isZero() && degree (F) > 0)
493	{
494	coF= Lc (F);
495	coG= 0;
496	return F/Lc(F);
497	}
498	else if (F.isZero() && G.isZero())
499	{
500	coF= coG= 0;
501	return F.genOne();
502	}
503	if (F.inBaseDomain() \|\| G.inBaseDomain())
504	{
505	coF= F;
506	coG= G;
507	return F.genOne();
508	}
509	if (F.isUnivariate() && fdivides(F, G, coG))
510	{
511	coF= Lc (F);
512	return F/Lc(F);
513	}
514	if (G.isUnivariate() && fdivides(G, F, coF))
515	{
516	coG= Lc (G);
517	return G/Lc(G);
518	}
519	if (F == G)
520	{
521	coF= coG= Lc (F);
522	return F/Lc(F);
523	}
524
525	CFMap M,N;
526	int best_level= myCompress (A, B, M, N, topLevel);
527
528	if (best_level == 0)
529	{
530	coF= F;
531	coG= G;
532	return B.genOne();
533	}
534
535	A= M(A);
536	B= M(B);
537
538	Variable x= Variable(1);
539
540	//univariate case
541	if (A.isUnivariate() && B.isUnivariate())
542	{
543	CanonicalForm result= gcd (A, B);
544	coF= N (A/result);
545	coG= N (B/result);
546	return N (result);
547	}
548
549	CanonicalForm cA, cB; // content of A and B
550	CanonicalForm ppA, ppB; // primitive part of A and B
551	CanonicalForm gcdcAcB;
552
553	cA = uni_content (A);
554	cB = uni_content (B);
555	gcdcAcB= gcd (cA, cB);
556	ppA= A/cA;
557	ppB= B/cB;
558
559	CanonicalForm lcA, lcB; // leading coefficients of A and B
560	CanonicalForm gcdlcAlcB;
561
562	lcA= uni_lcoeff (ppA);
563	lcB= uni_lcoeff (ppB);
564
565	gcdlcAlcB= gcd (lcA, lcB);
566
567	int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
568
569	if (d == 0)
570	{
571	coF= N (ppA*(cA/gcdcAcB));
572	coG= N (ppB*(cB/gcdcAcB));
573	return N(gcdcAcB);
574	}
575
576	int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
577	if (d0 < d)
578	d= d0;
579	if (d == 0)
580	{
581	coF= N (ppA*(cA/gcdcAcB));
582	coG= N (ppB*(cB/gcdcAcB));
583	return N(gcdcAcB);
584	}
585
586	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
587	CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
588	coG_m, ppCoF, ppCoG;
589
590	newtonPoly= 1;
591	m= gcdlcAlcB;
592	G_m= 0;
593	coF= 0;
594	coG= 0;
595	H= 0;
596	bool fail= false;
597	topLevel= false;
598	bool inextension= false;
599	Variable V_buf= alpha, V_buf4= alpha;
600	CanonicalForm prim_elem, im_prim_elem;
601	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
602	CFList source, dest;
603	int bound1= degree (ppA, 1);
604	int bound2= degree (ppB, 1);
605	do
606	{
607	random_element= randomElement (mlcAlcB, V_buf, l, fail);
608	if (fail)
609	{
610	source= CFList();
611	dest= CFList();
612
613	Variable V_buf3= V_buf;
614	V_buf= chooseExtension (V_buf);
615	bool prim_fail= false;
616	Variable V_buf2;
617	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
618	if (V_buf4 == alpha)
619	prim_elem_alpha= prim_elem;
620
621	if (V_buf3 != V_buf4)
622	{
623	m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
624	G_m= mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
625	coF_m= mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
626	coG_m= mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
627	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
628	source, dest);
629	ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
630	ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
631	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
632	source, dest);
633	lcA= mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
634	lcB= mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
635	for (CFListIterator i= l; i.hasItem(); i++)
636	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
637	source, dest);
638	}
639
640	ASSERT (!prim_fail, "failure in integer factorizer");
641	if (prim_fail)
642	; //ERROR
643	else
644	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
645
646	if (V_buf4 == alpha)
647	im_prim_elem_alpha= im_prim_elem;
648	else
649	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
650	im_prim_elem, source, dest);
651	DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
652	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
653	inextension= true;
654	for (CFListIterator i= l; i.hasItem(); i++)
655	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
656	im_prim_elem, source, dest);
657	m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
658	G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
659	coF_m= mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
660	coG_m= mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
661	newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
662	source, dest);
663	ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
664	ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
665	gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
666	source, dest);
667	lcA= mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
668	lcB= mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
669
670	fail= false;
671	random_element= randomElement (mlcAlcB, V_buf, l, fail );
672	DEBOUTLN (cerr, "fail= " << fail);
673	CFList list;
674	TIMING_START (gcd_recursion);
675	G_random_element=
676	modGCDFq (ppA (random_element, x), ppB (random_element, x),
677	coF_random_element, coG_random_element, V_buf,
678	list, topLevel);
679	TIMING_END_AND_PRINT (gcd_recursion,
680	"time for recursive call: ");
681	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
682	V_buf4= V_buf;
683	}
684	else
685	{
686	CFList list;
687	TIMING_START (gcd_recursion);
688	G_random_element=
689	modGCDFq (ppA(random_element, x), ppB(random_element, x),
690	coF_random_element, coG_random_element, V_buf,
691	list, topLevel);
692	TIMING_END_AND_PRINT (gcd_recursion,
693	"time for recursive call: ");
694	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
695	}
696
697	if (!G_random_element.inCoeffDomain())
698	d0= totaldegree (G_random_element, Variable(2),
699	Variable (G_random_element.level()));
700	else
701	d0= 0;
702
703	if (d0 == 0)
704	{
705	if (inextension)
706	{
707	CFList u, v;
708	ppA= mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
709	ppB= mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
710	prune1 (alpha);
711	}
712	coF= N (ppA*(cA/gcdcAcB));
713	coG= N (ppB*(cB/gcdcAcB));
714	return N(gcdcAcB);
715	}
716	if (d0 > d)
717	{
718	if (!find (l, random_element))
719	l.append (random_element);
720	continue;
721	}
722
723	G_random_element=
724	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
725	* G_random_element;
726	coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
727	*coF_random_element;
728	coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
729	*coG_random_element;
730
731	if (!G_random_element.inCoeffDomain())
732	d0= totaldegree (G_random_element, Variable(2),
733	Variable (G_random_element.level()));
734	else
735	d0= 0;
736
737	if (d0 < d)
738	{
739	m= gcdlcAlcB;
740	newtonPoly= 1;
741	G_m= 0;
742	d= d0;
743	coF_m= 0;
744	coG_m= 0;
745	}
746
747	TIMING_START (newton_interpolation);
748	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
749	coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
750	coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
751	TIMING_END_AND_PRINT (newton_interpolation,
752	"time for newton interpolation: ");
753
754	//termination test
755	if ((uni_lcoeff (H) == gcdlcAlcB) \|\| (G_m == H))
756	{
757	TIMING_START (termination_test);
758	if (gcdlcAlcB.isOne())
759	cH= 1;
760	else
761	cH= uni_content (H);
762	ppH= H/cH;
763	ppH /= Lc (ppH);
764	CanonicalForm lcppH= gcdlcAlcB/cH;
765	CanonicalForm ccoF= lcppH/Lc (lcppH);
766	CanonicalForm ccoG= lcppH/Lc (lcppH);
767	ppCoF= coF/ccoF;
768	ppCoG= coG/ccoG;
769	if (inextension)
770	{
771	if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
772	(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
773	terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) \|\|
774	(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
775	{
776	CFList u, v;
777	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
778	ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
779	ppCoF= mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
780	ppCoG= mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
781	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
782	coF= N ((cA/gcdcAcB)*ppCoF);
783	coG= N ((cB/gcdcAcB)*ppCoG);
784	TIMING_END_AND_PRINT (termination_test,
785	"time for successful termination test Fq: ");
786	prune1 (alpha);
787	return N(gcdcAcB*ppH);
788	}
789	}
790	else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
791	(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
792	terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) \|\|
793	(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
794	{
795	coF= N ((cA/gcdcAcB)*ppCoF);
796	coG= N ((cB/gcdcAcB)*ppCoG);
797	TIMING_END_AND_PRINT (termination_test,
798	"time for successful termination test Fq: ");
799	return N(gcdcAcB*ppH);
800	}
801	TIMING_END_AND_PRINT (termination_test,
802	"time for unsuccessful termination test Fq: ");
803	}
804
805	G_m= H;
806	coF_m= coF;
807	coG_m= coG;
808	newtonPoly= newtonPoly*(x - random_element);
809	m= m*(x - random_element);
810	if (!find (l, random_element))
811	l.append (random_element);
812	} while (1);
813	}
814	#endif
815
816	/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
817	/// univariate polynomial, returns fail if there are no field elements left
818	/// which have not been used before
819	static inline
820	CanonicalForm
821	GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
822	{
823	fail= false;
824	Variable x= F.mvar();
825	GFRandom genGF;
826	CanonicalForm random;
827	int p= getCharacteristic();
828	int d= getGFDegree();
829	int bound= ipower (p, d);
830	do
831	{
832	if (list.length() == bound)
833	{
834	fail= true;
835	break;
836	}
837	if (list.length() < 1)
838	random= 0;
839	else
840	{
841	random= genGF.generate();
842	while (find (list, random))
843	random= genGF.generate();
844	}
845	if (F (random, x) == 0)
846	{
847	list.append (random);
848	continue;
849	}
850	} while (find (list, random));
851	return random;
852	}
853
854	CanonicalForm
855	modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
856	CanonicalForm& coF, CanonicalForm& coG,
857	CFList& l, bool& topLevel);
858
859	CanonicalForm
860	modGCDGF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
861	bool& topLevel)
862	{
863	CanonicalForm dummy1, dummy2;
864	CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
865	return result;
866	}
867
868	/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
869	/// Computer Algebra" by Geddes, Czapor, Labahn
870	/// Usually this algorithm will be faster than modGCDFq since GF has
871	/// faster field arithmetics, however it might fail if the input is large since
872	/// the size of the base field is bounded by 2^16, output is monic
873	CanonicalForm
874	modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
875	CanonicalForm& coF, CanonicalForm& coG,
876	CFList& l, bool& topLevel)
877	{
878	CanonicalForm A= F;
879	CanonicalForm B= G;
880	if (F.isZero() && degree(G) > 0)
881	{
882	coF= 0;
883	coG= Lc (G);
884	return G/Lc(G);
885	}
886	else if (G.isZero() && degree (F) > 0)
887	{
888	coF= Lc (F);
889	coG= 0;
890	return F/Lc(F);
891	}
892	else if (F.isZero() && G.isZero())
893	{
894	coF= coG= 0;
895	return F.genOne();
896	}
897	if (F.inBaseDomain() \|\| G.inBaseDomain())
898	{
899	coF= F;
900	coG= G;
901	return F.genOne();
902	}
903	if (F.isUnivariate() && fdivides(F, G, coG))
904	{
905	coF= Lc (F);
906	return F/Lc(F);
907	}
908	if (G.isUnivariate() && fdivides(G, F, coF))
909	{
910	coG= Lc (G);
911	return G/Lc(G);
912	}
913	if (F == G)
914	{
915	coF= coG= Lc (F);
916	return F/Lc(F);
917	}
918
919	CFMap M,N;
920	int best_level= myCompress (A, B, M, N, topLevel);
921
922	if (best_level == 0)
923	{
924	coF= F;
925	coG= G;
926	return B.genOne();
927	}
928
929	A= M(A);
930	B= M(B);
931
932	Variable x= Variable(1);
933
934	//univariate case
935	if (A.isUnivariate() && B.isUnivariate())
936	{
937	CanonicalForm result= gcd (A, B);
938	coF= N (A/result);
939	coG= N (B/result);
940	return N (result);
941	}
942
943	CanonicalForm cA, cB; // content of A and B
944	CanonicalForm ppA, ppB; // primitive part of A and B
945	CanonicalForm gcdcAcB;
946
947	cA = uni_content (A);
948	cB = uni_content (B);
949	gcdcAcB= gcd (cA, cB);
950	ppA= A/cA;
951	ppB= B/cB;
952
953	CanonicalForm lcA, lcB; // leading coefficients of A and B
954	CanonicalForm gcdlcAlcB;
955
956	lcA= uni_lcoeff (ppA);
957	lcB= uni_lcoeff (ppB);
958
959	gcdlcAlcB= gcd (lcA, lcB);
960
961	int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
962	if (d == 0)
963	{
964	coF= N (ppA*(cA/gcdcAcB));
965	coG= N (ppB*(cB/gcdcAcB));
966	return N(gcdcAcB);
967	}
968	int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
969	if (d0 < d)
970	d= d0;
971	if (d == 0)
972	{
973	coF= N (ppA*(cA/gcdcAcB));
974	coG= N (ppB*(cB/gcdcAcB));
975	return N(gcdcAcB);
976	}
977
978	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
979	CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
980	coG_m, ppCoF, ppCoG;
981
982	newtonPoly= 1;
983	m= gcdlcAlcB;
984	G_m= 0;
985	coF= 0;
986	coG= 0;
987	H= 0;
988	bool fail= false;
989	topLevel= false;
990	bool inextension= false;
991	int p=-1;
992	int k= getGFDegree();
993	int kk;
994	int expon;
995	char gf_name_buf= gf_name;
996	int bound1= degree (ppA, 1);
997	int bound2= degree (ppB, 1);
998	do
999	{
1000	random_element= GFRandomElement (mlcAlcB, l, fail);
1001	if (fail)
1002	{
1003	p= getCharacteristic();
1004	expon= 2;
1005	kk= getGFDegree();
1006	if (ipower (p, kk*expon) < (1 << 16))
1007	setCharacteristic (p, kk*(int)expon, 'b');
1008	else
1009	{
1010	expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
1011	ASSERT (expon >= 2, "not enough points in modGCDGF");
1012	setCharacteristic (p, (int)(kk*expon), 'b');
1013	}
1014	inextension= true;
1015	fail= false;
1016	for (CFListIterator i= l; i.hasItem(); i++)
1017	i.getItem()= GFMapUp (i.getItem(), kk);
1018	m= GFMapUp (m, kk);
1019	G_m= GFMapUp (G_m, kk);
1020	newtonPoly= GFMapUp (newtonPoly, kk);
1021	coF_m= GFMapUp (coF_m, kk);
1022	coG_m= GFMapUp (coG_m, kk);
1023	ppA= GFMapUp (ppA, kk);
1024	ppB= GFMapUp (ppB, kk);
1025	gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
1026	lcA= GFMapUp (lcA, kk);
1027	lcB= GFMapUp (lcB, kk);
1028	random_element= GFRandomElement (mlcAlcB, l, fail);
1029	DEBOUTLN (cerr, "fail= " << fail);
1030	CFList list;
1031	TIMING_START (gcd_recursion);
1032	G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1033	coF_random_element, coG_random_element,
1034	list, topLevel);
1035	TIMING_END_AND_PRINT (gcd_recursion,
1036	"time for recursive call: ");
1037	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1038	}
1039	else
1040	{
1041	CFList list;
1042	TIMING_START (gcd_recursion);
1043	G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1044	coF_random_element, coG_random_element,
1045	list, topLevel);
1046	TIMING_END_AND_PRINT (gcd_recursion,
1047	"time for recursive call: ");
1048	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1049	}
1050
1051	if (!G_random_element.inCoeffDomain())
1052	d0= totaldegree (G_random_element, Variable(2),
1053	Variable (G_random_element.level()));
1054	else
1055	d0= 0;
1056
1057	if (d0 == 0)
1058	{
1059	if (inextension)
1060	{
1061	ppA= GFMapDown (ppA, k);
1062	ppB= GFMapDown (ppB, k);
1063	setCharacteristic (p, k, gf_name_buf);
1064	}
1065	coF= N (ppA*(cA/gcdcAcB));
1066	coG= N (ppB*(cB/gcdcAcB));
1067	return N(gcdcAcB);
1068	}
1069	if (d0 > d)
1070	{
1071	if (!find (l, random_element))
1072	l.append (random_element);
1073	continue;
1074	}
1075
1076	G_random_element=
1077	(gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
1078	G_random_element;
1079
1080	coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1081	*coF_random_element;
1082	coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1083	*coG_random_element;
1084
1085	if (!G_random_element.inCoeffDomain())
1086	d0= totaldegree (G_random_element, Variable(2),
1087	Variable (G_random_element.level()));
1088	else
1089	d0= 0;
1090
1091	if (d0 < d)
1092	{
1093	m= gcdlcAlcB;
1094	newtonPoly= 1;
1095	G_m= 0;
1096	d= d0;
1097	coF_m= 0;
1098	coG_m= 0;
1099	}
1100
1101	TIMING_START (newton_interpolation);
1102	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1103	coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1104	coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1105	TIMING_END_AND_PRINT (newton_interpolation,
1106	"time for newton interpolation: ");
1107
1108	//termination test
1109	if ((uni_lcoeff (H) == gcdlcAlcB) \|\| (G_m == H))
1110	{
1111	TIMING_START (termination_test);
1112	if (gcdlcAlcB.isOne())
1113	cH= 1;
1114	else
1115	cH= uni_content (H);
1116	ppH= H/cH;
1117	ppH /= Lc (ppH);
1118	CanonicalForm lcppH= gcdlcAlcB/cH;
1119	CanonicalForm ccoF= lcppH/Lc (lcppH);
1120	CanonicalForm ccoG= lcppH/Lc (lcppH);
1121	ppCoF= coF/ccoF;
1122	ppCoG= coG/ccoG;
1123	if (inextension)
1124	{
1125	if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1126	(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1127	terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) \|\|
1128	(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1129	{
1130	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1131	ppH= GFMapDown (ppH, k);
1132	ppCoF= GFMapDown (ppCoF, k);
1133	ppCoG= GFMapDown (ppCoG, k);
1134	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1135	coF= N ((cA/gcdcAcB)*ppCoF);
1136	coG= N ((cB/gcdcAcB)*ppCoG);
1137	setCharacteristic (p, k, gf_name_buf);
1138	TIMING_END_AND_PRINT (termination_test,
1139	"time for successful termination GF: ");
1140	return N(gcdcAcB*ppH);
1141	}
1142	}
1143	else
1144	{
1145	if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1146	(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1147	terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) \|\|
1148	(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1149	{
1150	coF= N ((cA/gcdcAcB)*ppCoF);
1151	coG= N ((cB/gcdcAcB)*ppCoG);
1152	TIMING_END_AND_PRINT (termination_test,
1153	"time for successful termination GF: ");
1154	return N(gcdcAcB*ppH);
1155	}
1156	}
1157	TIMING_END_AND_PRINT (termination_test,
1158	"time for unsuccessful termination GF: ");
1159	}
1160
1161	G_m= H;
1162	coF_m= coF;
1163	coG_m= coG;
1164	newtonPoly= newtonPoly*(x - random_element);
1165	m= m*(x - random_element);
1166	if (!find (l, random_element))
1167	l.append (random_element);
1168	} while (1);
1169	}
1170
1171	static inline
1172	CanonicalForm
1173	FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1174	{
1175	fail= false;
1176	Variable x= F.mvar();
1177	FFRandom genFF;
1178	CanonicalForm random;
1179	int p= getCharacteristic();
1180	int bound= p;
1181	do
1182	{
1183	if (list.length() == bound)
1184	{
1185	fail= true;
1186	break;
1187	}
1188	if (list.length() < 1)
1189	random= 0;
1190	else
1191	{
1192	random= genFF.generate();
1193	while (find (list, random))
1194	random= genFF.generate();
1195	}
1196	if (F (random, x) == 0)
1197	{
1198	list.append (random);
1199	continue;
1200	}
1201	} while (find (list, random));
1202	return random;
1203	}
1204
1205	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
1206	CanonicalForm
1207	modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1208	CanonicalForm& coF, CanonicalForm& coG,
1209	bool& topLevel, CFList& l);
1210	#endif
1211
1212	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
1213	CanonicalForm
1214	modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1215	bool& topLevel, CFList& l)
1216	{
1217	CanonicalForm dummy1, dummy2;
1218	CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1219	return result;
1220	}
1221	#endif
1222
1223	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
1224	CanonicalForm
1225	modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1226	CanonicalForm& coF, CanonicalForm& coG,
1227	bool& topLevel, CFList& l)
1228	{
1229	CanonicalForm A= F;
1230	CanonicalForm B= G;
1231	if (F.isZero() && degree(G) > 0)
1232	{
1233	coF= 0;
1234	coG= Lc (G);
1235	return G/Lc(G);
1236	}
1237	else if (G.isZero() && degree (F) > 0)
1238	{
1239	coF= Lc (F);
1240	coG= 0;
1241	return F/Lc(F);
1242	}
1243	else if (F.isZero() && G.isZero())
1244	{
1245	coF= coG= 0;
1246	return F.genOne();
1247	}
1248	if (F.inBaseDomain() \|\| G.inBaseDomain())
1249	{
1250	coF= F;
1251	coG= G;
1252	return F.genOne();
1253	}
1254	if (F.isUnivariate() && fdivides(F, G, coG))
1255	{
1256	coF= Lc (F);
1257	return F/Lc(F);
1258	}
1259	if (G.isUnivariate() && fdivides(G, F, coF))
1260	{
1261	coG= Lc (G);
1262	return G/Lc(G);
1263	}
1264	if (F == G)
1265	{
1266	coF= coG= Lc (F);
1267	return F/Lc(F);
1268	}
1269	CFMap M,N;
1270	int best_level= myCompress (A, B, M, N, topLevel);
1271
1272	if (best_level == 0)
1273	{
1274	coF= F;
1275	coG= G;
1276	return B.genOne();
1277	}
1278
1279	A= M(A);
1280	B= M(B);
1281
1282	Variable x= Variable (1);
1283
1284	//univariate case
1285	if (A.isUnivariate() && B.isUnivariate())
1286	{
1287	CanonicalForm result= gcd (A, B);
1288	coF= N (A/result);
1289	coG= N (B/result);
1290	return N (result);
1291	}
1292
1293	CanonicalForm cA, cB; // content of A and B
1294	CanonicalForm ppA, ppB; // primitive part of A and B
1295	CanonicalForm gcdcAcB;
1296
1297	cA = uni_content (A);
1298	cB = uni_content (B);
1299	gcdcAcB= gcd (cA, cB);
1300	ppA= A/cA;
1301	ppB= B/cB;
1302
1303	CanonicalForm lcA, lcB; // leading coefficients of A and B
1304	CanonicalForm gcdlcAlcB;
1305	lcA= uni_lcoeff (ppA);
1306	lcB= uni_lcoeff (ppB);
1307
1308	gcdlcAlcB= gcd (lcA, lcB);
1309
1310	int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1311	int d0;
1312
1313	if (d == 0)
1314	{
1315	coF= N (ppA*(cA/gcdcAcB));
1316	coG= N (ppB*(cB/gcdcAcB));
1317	return N(gcdcAcB);
1318	}
1319
1320	d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1321
1322	if (d0 < d)
1323	d= d0;
1324
1325	if (d == 0)
1326	{
1327	coF= N (ppA*(cA/gcdcAcB));
1328	coG= N (ppB*(cB/gcdcAcB));
1329	return N(gcdcAcB);
1330	}
1331
1332	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1333	CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1334	coF_m, coG_m, ppCoF, ppCoG;
1335
1336	newtonPoly= 1;
1337	m= gcdlcAlcB;
1338	H= 0;
1339	coF= 0;
1340	coG= 0;
1341	G_m= 0;
1342	Variable alpha, V_buf, cleanUp;
1343	bool fail= false;
1344	bool inextension= false;
1345	topLevel= false;
1346	CFList source, dest;
1347	int bound1= degree (ppA, 1);
1348	int bound2= degree (ppB, 1);
1349	do
1350	{
1351	if (inextension)
1352	random_element= randomElement (mlcAlcB, V_buf, l, fail);
1353	else
1354	random_element= FpRandomElement (mlcAlcB, l, fail);
1355
1356	if (!fail && !inextension)
1357	{
1358	CFList list;
1359	TIMING_START (gcd_recursion);
1360	G_random_element=
1361	modGCDFp (ppA (random_element,x), ppB (random_element,x),
1362	coF_random_element, coG_random_element, topLevel,
1363	list);
1364	TIMING_END_AND_PRINT (gcd_recursion,
1365	"time for recursive call: ");
1366	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1367	}
1368	else if (!fail && inextension)
1369	{
1370	CFList list;
1371	TIMING_START (gcd_recursion);
1372	G_random_element=
1373	modGCDFq (ppA (random_element, x), ppB (random_element, x),
1374	coF_random_element, coG_random_element, V_buf,
1375	list, topLevel);
1376	TIMING_END_AND_PRINT (gcd_recursion,
1377	"time for recursive call: ");
1378	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1379	}
1380	else if (fail && !inextension)
1381	{
1382	source= CFList();
1383	dest= CFList();
1384	CFList list;
1385	CanonicalForm mipo;
1386	int deg= 2;
1387	bool initialized= false;
1388	do
1389	{
1390	mipo= randomIrredpoly (deg, x);
1391	if (initialized)
1392	setMipo (alpha, mipo);
1393	else
1394	alpha= rootOf (mipo);
1395	inextension= true;
1396	initialized= true;
1397	fail= false;
1398	random_element= randomElement (mlcAlcB, alpha, l, fail);
1399	deg++;
1400	} while (fail);
1401	list= CFList();
1402	V_buf= alpha;
1403	cleanUp= alpha;
1404	TIMING_START (gcd_recursion);
1405	G_random_element=
1406	modGCDFq (ppA (random_element, x), ppB (random_element, x),
1407	coF_random_element, coG_random_element, alpha,
1408	list, topLevel);
1409	TIMING_END_AND_PRINT (gcd_recursion,
1410	"time for recursive call: ");
1411	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1412	}
1413	else if (fail && inextension)
1414	{
1415	source= CFList();
1416	dest= CFList();
1417
1418	Variable V_buf3= V_buf;
1419	V_buf= chooseExtension (V_buf);
1420	bool prim_fail= false;
1421	Variable V_buf2;
1422	CanonicalForm prim_elem, im_prim_elem;
1423	prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1424
1425	if (V_buf3 != alpha)
1426	{
1427	m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1428	G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1429	coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1430	coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1431	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1432	source, dest);
1433	ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1434	ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1435	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1436	dest);
1437	lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1438	lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1439	for (CFListIterator i= l; i.hasItem(); i++)
1440	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1441	source, dest);
1442	}
1443
1444	ASSERT (!prim_fail, "failure in integer factorizer");
1445	if (prim_fail)
1446	; //ERROR
1447	else
1448	im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1449
1450	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1451	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1452
1453	for (CFListIterator i= l; i.hasItem(); i++)
1454	i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1455	im_prim_elem, source, dest);
1456	m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1457	G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1458	coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1459	coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1460	newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1461	source, dest);
1462	ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1463	ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1464	gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1465	source, dest);
1466	lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1467	lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1468	fail= false;
1469	random_element= randomElement (mlcAlcB, V_buf, l, fail );
1470	DEBOUTLN (cerr, "fail= " << fail);
1471	CFList list;
1472	TIMING_START (gcd_recursion);
1473	G_random_element=
1474	modGCDFq (ppA (random_element, x), ppB (random_element, x),
1475	coF_random_element, coG_random_element, V_buf,
1476	list, topLevel);
1477	TIMING_END_AND_PRINT (gcd_recursion,
1478	"time for recursive call: ");
1479	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1480	alpha= V_buf;
1481	}
1482
1483	if (!G_random_element.inCoeffDomain())
1484	d0= totaldegree (G_random_element, Variable(2),
1485	Variable (G_random_element.level()));
1486	else
1487	d0= 0;
1488
1489	if (d0 == 0)
1490	{
1491	if (inextension)
1492	prune (cleanUp);
1493	coF= N (ppA*(cA/gcdcAcB));
1494	coG= N (ppB*(cB/gcdcAcB));
1495	return N(gcdcAcB);
1496	}
1497
1498	if (d0 > d)
1499	{
1500	if (!find (l, random_element))
1501	l.append (random_element);
1502	continue;
1503	}
1504
1505	G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1506	*G_random_element;
1507
1508	coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1509	*coF_random_element;
1510	coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1511	*coG_random_element;
1512
1513	if (!G_random_element.inCoeffDomain())
1514	d0= totaldegree (G_random_element, Variable(2),
1515	Variable (G_random_element.level()));
1516	else
1517	d0= 0;
1518
1519	if (d0 < d)
1520	{
1521	m= gcdlcAlcB;
1522	newtonPoly= 1;
1523	G_m= 0;
1524	d= d0;
1525	coF_m= 0;
1526	coG_m= 0;
1527	}
1528
1529	TIMING_START (newton_interpolation);
1530	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1531	coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1532	coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1533	TIMING_END_AND_PRINT (newton_interpolation,
1534	"time for newton_interpolation: ");
1535
1536	//termination test
1537	if ((uni_lcoeff (H) == gcdlcAlcB) \|\| (G_m == H))
1538	{
1539	TIMING_START (termination_test);
1540	if (gcdlcAlcB.isOne())
1541	cH= 1;
1542	else
1543	cH= uni_content (H);
1544	ppH= H/cH;
1545	ppH /= Lc (ppH);
1546	CanonicalForm lcppH= gcdlcAlcB/cH;
1547	CanonicalForm ccoF= lcppH/Lc (lcppH);
1548	CanonicalForm ccoG= lcppH/Lc (lcppH);
1549	ppCoF= coF/ccoF;
1550	ppCoG= coG/ccoG;
1551	DEBOUTLN (cerr, "ppH= " << ppH);
1552	if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1553	(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1554	terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) \|\|
1555	(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1556	{
1557	if (inextension)
1558	prune (cleanUp);
1559	coF= N ((cA/gcdcAcB)*ppCoF);
1560	coG= N ((cB/gcdcAcB)*ppCoG);
1561	TIMING_END_AND_PRINT (termination_test,
1562	"time for successful termination Fp: ");
1563	return N(gcdcAcB*ppH);
1564	}
1565	TIMING_END_AND_PRINT (termination_test,
1566	"time for unsuccessful termination Fp: ");
1567	}
1568
1569	G_m= H;
1570	coF_m= coF;
1571	coG_m= coG;
1572	newtonPoly= newtonPoly*(x - random_element);
1573	m= m*(x - random_element);
1574	if (!find (l, random_element))
1575	l.append (random_element);
1576	} while (1);
1577	}
1578	#endif
1579
1580	CFArray
1581	solveVandermonde (const CFArray& M, const CFArray& A)
1582	{
1583	int r= M.size();
1584	ASSERT (A.size() == r, "vector does not have right size");
1585
1586	if (r == 1)
1587	{
1588	CFArray result= CFArray (1);
1589	result [0]= A [0] / M [0];
1590	return result;
1591	}
1592	// check solvability
1593	bool notDistinct= false;
1594	for (int i= 0; i < r - 1; i++)
1595	{
1596	for (int j= i + 1; j < r; j++)
1597	{
1598	if (M [i] == M [j])
1599	{
1600	notDistinct= true;
1601	break;
1602	}
1603	}
1604	}
1605	if (notDistinct)
1606	return CFArray();
1607
1608	CanonicalForm master= 1;
1609	Variable x= Variable (1);
1610	for (int i= 0; i < r; i++)
1611	master *= x - M [i];
1612	CFList Pj;
1613	CanonicalForm tmp;
1614	for (int i= 0; i < r; i++)
1615	{
1616	tmp= master/(x - M [i]);
1617	tmp /= tmp (M [i], 1);
1618	Pj.append (tmp);
1619	}
1620	CFArray result= CFArray (r);
1621
1622	CFListIterator j= Pj;
1623	for (int i= 1; i <= r; i++, j++)
1624	{
1625	tmp= 0;
1626	for (int l= 0; l < A.size(); l++)
1627	tmp += A[l]*j.getItem()[l];
1628	result[i - 1]= tmp;
1629	}
1630	return result;
1631	}
1632
1633	CFArray
1634	solveGeneralVandermonde (const CFArray& M, const CFArray& A)
1635	{
1636	int r= M.size();
1637	ASSERT (A.size() == r, "vector does not have right size");
1638	if (r == 1)
1639	{
1640	CFArray result= CFArray (1);
1641	result [0]= A[0] / M [0];
1642	return result;
1643	}
1644	// check solvability
1645	bool notDistinct= false;
1646	for (int i= 0; i < r - 1; i++)
1647	{
1648	for (int j= i + 1; j < r; j++)
1649	{
1650	if (M [i] == M [j])
1651	{
1652	notDistinct= true;
1653	break;
1654	}
1655	}
1656	}
1657	if (notDistinct)
1658	return CFArray();
1659
1660	CanonicalForm master= 1;
1661	Variable x= Variable (1);
1662	for (int i= 0; i < r; i++)
1663	master *= x - M [i];
1664	master *= x;
1665	CFList Pj;
1666	CanonicalForm tmp;
1667	for (int i= 0; i < r; i++)
1668	{
1669	tmp= master/(x - M [i]);
1670	tmp /= tmp (M [i], 1);
1671	Pj.append (tmp);
1672	}
1673
1674	CFArray result= CFArray (r);
1675
1676	CFListIterator j= Pj;
1677	for (int i= 1; i <= r; i++, j++)
1678	{
1679	tmp= 0;
1680
1681	for (int l= 1; l <= A.size(); l++)
1682	tmp += A[l - 1]*j.getItem()[l];
1683	result[i - 1]= tmp;
1684	}
1685	return result;
1686	}
1687
1688	/// M in row echolon form, rk rank of M
1689	CFArray
1690	readOffSolution (const CFMatrix& M, const long rk)
1691	{
1692	CFArray result= CFArray (rk);
1693	CanonicalForm tmp1, tmp2, tmp3;
1694	for (int i= rk; i >= 1; i--)
1695	{
1696	tmp3= 0;
1697	tmp1= M (i, M.columns());
1698	for (int j= M.columns() - 1; j >= 1; j--)
1699	{
1700	tmp2= M (i, j);
1701	if (j == i)
1702	break;
1703	else
1704	tmp3 += tmp2*result[j - 1];
1705	}
1706	result[i - 1]= (tmp1 - tmp3)/tmp2;
1707	}
1708	return result;
1709	}
1710
1711	CFArray
1712	readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
1713	{
1714	CFArray result= CFArray (M.rows());
1715	CanonicalForm tmp1, tmp2, tmp3;
1716	int k;
1717	for (int i= M.rows(); i >= 1; i--)
1718	{
1719	tmp3= 0;
1720	tmp1= L[i - 1];
1721	k= 0;
1722	for (int j= M.columns(); j >= 1; j--, k++)
1723	{
1724	tmp2= M (i, j);
1725	if (j == i)
1726	break;
1727	else
1728	{
1729	if (k > partialSol.size() - 1)
1730	tmp3 += tmp2*result[j - 1];
1731	else
1732	tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1733	}
1734	}
1735	result [i - 1]= (tmp1 - tmp3)/tmp2;
1736	}
1737	return result;
1738	}
1739
1740	long
1741	gaussianElimFp (CFMatrix& M, CFArray& L)
1742	{
1743	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1744	CFMatrix *N;
1745	N= new CFMatrix (M.rows(), M.columns() + 1);
1746
1747	for (int i= 1; i <= M.rows(); i++)
1748	for (int j= 1; j <= M.columns(); j++)
1749	(*N) (i, j)= M (i, j);
1750
1751	int j= 1;
1752	for (int i= 0; i < L.size(); i++, j++)
1753	(*N) (j, M.columns() + 1)= L[i];
1754	#ifdef HAVE_FLINT
1755	nmod_mat_t FLINTN;
1756	convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1757	long rk= nmod_mat_rref (FLINTN);
1758
1759	delete N;
1760	N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1761	nmod_mat_clear (FLINTN);
1762	#else
1763	int p= getCharacteristic ();
1764	if (fac_NTL_char != p)
1765	{
1766	fac_NTL_char= p;
1767	zz_p::init (p);
1768	}
1769	mat_zz_p NTLN= convertFacCFMatrix2NTLmat_zz_p(N);
1770	delete N;
1771	long rk= gauss (*NTLN);
1772
1773	N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1774	delete NTLN;
1775	#endif
1776
1777	L= CFArray (M.rows());
1778	for (int i= 0; i < M.rows(); i++)
1779	L[i]= (*N) (i + 1, M.columns() + 1);
1780	M= (*N) (1, M.rows(), 1, M.columns());
1781	delete N;
1782	return rk;
1783	}
1784
1785	long
1786	gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
1787	{
1788	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1789	CFMatrix *N;
1790	N= new CFMatrix (M.rows(), M.columns() + 1);
1791
1792	for (int i= 1; i <= M.rows(); i++)
1793	for (int j= 1; j <= M.columns(); j++)
1794	(*N) (i, j)= M (i, j);
1795
1796	int j= 1;
1797	for (int i= 0; i < L.size(); i++, j++)
1798	(*N) (j, M.columns() + 1)= L[i];
1799	#ifdef HAVE_FLINT
1800	// convert mipo
1801	nmod_poly_t mipo1;
1802	convertFacCF2nmod_poly_t(mipo1,getMipo(alpha));
1803	fq_nmod_ctx_t ctx;
1804	fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1805	nmod_poly_clear(mipo1);
1806	// convert matrix
1807	fq_nmod_mat_t FLINTN;
1808	convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1809	// rank
1810	long rk= fq_nmod_mat_rref (FLINTN,ctx);
1811	// clean up
1812	fq_nmod_mat_clear (FLINTN,ctx);
1813	fq_nmod_ctx_clear(ctx);
1814	#elif defined(HAVE_NTL)
1815	int p= getCharacteristic ();
1816	if (fac_NTL_char != p)
1817	{
1818	fac_NTL_char= p;
1819	zz_p::init (p);
1820	}
1821	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1822	zz_pE::init (NTLMipo);
1823	mat_zz_pE NTLN= convertFacCFMatrix2NTLmat_zz_pE(N);
1824	long rk= gauss (*NTLN);
1825	N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1826	delete NTLN;
1827	#else
1828	factoryError("NTL/FLINT missing: gaussianElimFq");
1829	#endif
1830	delete N;
1831
1832	M= (*N) (1, M.rows(), 1, M.columns());
1833	L= CFArray (M.rows());
1834	for (int i= 0; i < M.rows(); i++)
1835	L[i]= (*N) (i + 1, M.columns() + 1);
1836
1837	delete N;
1838	return rk;
1839	}
1840
1841	CFArray
1842	solveSystemFp (const CFMatrix& M, const CFArray& L)
1843	{
1844	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1845	CFMatrix *N;
1846	N= new CFMatrix (M.rows(), M.columns() + 1);
1847
1848	for (int i= 1; i <= M.rows(); i++)
1849	for (int j= 1; j <= M.columns(); j++)
1850	(*N) (i, j)= M (i, j);
1851
1852	int j= 1;
1853	for (int i= 0; i < L.size(); i++, j++)
1854	(*N) (j, M.columns() + 1)= L[i];
1855
1856	#ifdef HAVE_FLINT
1857	nmod_mat_t FLINTN;
1858	convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1859	long rk= nmod_mat_rref (FLINTN);
1860	#else
1861	int p= getCharacteristic ();
1862	if (fac_NTL_char != p)
1863	{
1864	fac_NTL_char= p;
1865	zz_p::init (p);
1866	}
1867	mat_zz_p NTLN= convertFacCFMatrix2NTLmat_zz_p(N);
1868	long rk= gauss (*NTLN);
1869	#endif
1870	delete N;
1871	if (rk != M.columns())
1872	{
1873	#ifdef HAVE_FLINT
1874	nmod_mat_clear (FLINTN);
1875	#else
1876	delete NTLN;
1877	#endif
1878	return CFArray();
1879	}
1880	#ifdef HAVE_FLINT
1881	N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1882	nmod_mat_clear (FLINTN);
1883	#else
1884	N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1885	delete NTLN;
1886	#endif
1887	CFArray A= readOffSolution (*N, rk);
1888
1889	delete N;
1890	return A;
1891	}
1892
1893	CFArray
1894	solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1895	{
1896	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1897	CFMatrix *N;
1898	N= new CFMatrix (M.rows(), M.columns() + 1);
1899
1900	for (int i= 1; i <= M.rows(); i++)
1901	for (int j= 1; j <= M.columns(); j++)
1902	(*N) (i, j)= M (i, j);
1903	int j= 1;
1904	for (int i= 0; i < L.size(); i++, j++)
1905	(*N) (j, M.columns() + 1)= L[i];
1906	#ifdef HAVE_FLINT
1907	// convert mipo
1908	nmod_poly_t mipo1;
1909	convertFacCF2nmod_poly_t(mipo1,getMipo(alpha));
1910	fq_nmod_ctx_t ctx;
1911	fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1912	nmod_poly_clear(mipo1);
1913	// convert matrix
1914	fq_nmod_mat_t FLINTN;
1915	convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1916	// rank
1917	long rk= fq_nmod_mat_rref (FLINTN,ctx);
1918	#elif defined(HAVE_NTL)
1919	int p= getCharacteristic ();
1920	if (fac_NTL_char != p)
1921	{
1922	fac_NTL_char= p;
1923	zz_p::init (p);
1924	}
1925	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1926	zz_pE::init (NTLMipo);
1927	mat_zz_pE NTLN= convertFacCFMatrix2NTLmat_zz_pE(N);
1928	long rk= gauss (*NTLN);
1929	#else
1930	factoryError("NTL/FLINT missing: solveSystemFq");
1931	#endif
1932
1933	delete N;
1934	if (rk != M.columns())
1935	{
1936	#if defined(HAVE_NTL) && !defined(HAVE_FLINT)
1937	delete NTLN;
1938	#endif
1939	return CFArray();
1940	}
1941	#ifdef HAVE_FLINT
1942	// convert and clean up
1943	N=convertFq_nmod_mat_t2FacCFMatrix(FLINTN,ctx,alpha);
1944	fq_nmod_mat_clear (FLINTN,ctx);
1945	fq_nmod_ctx_clear(ctx);
1946	#elif defined(HAVE_NTL)
1947	N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1948	delete NTLN;
1949	#endif
1950
1951	CFArray A= readOffSolution (*N, rk);
1952
1953	delete N;
1954	return A;
1955	}
1956	#endif
1957
1958	CFArray
1959	getMonoms (const CanonicalForm& F)
1960	{
1961	if (F.inCoeffDomain())
1962	{
1963	CFArray result= CFArray (1);
1964	result [0]= 1;
1965	return result;
1966	}
1967	if (F.isUnivariate())
1968	{
1969	CFArray result= CFArray (size(F));
1970	int j= 0;
1971	for (CFIterator i= F; i.hasTerms(); i++, j++)
1972	result[j]= power (F.mvar(), i.exp());
1973	return result;
1974	}
1975	int numMon= size (F);
1976	CFArray result= CFArray (numMon);
1977	int j= 0;
1978	CFArray recResult;
1979	Variable x= F.mvar();
1980	CanonicalForm powX;
1981	for (CFIterator i= F; i.hasTerms(); i++)
1982	{
1983	powX= power (x, i.exp());
1984	recResult= getMonoms (i.coeff());
1985	for (int k= 0; k < recResult.size(); k++)
1986	result[j+k]= powX*recResult[k];
1987	j += recResult.size();
1988	}
1989	return result;
1990	}
1991
1992	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
1993	CFArray
1994	evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
1995	{
1996	if (F.inCoeffDomain())
1997	{
1998	CFArray result= CFArray (1);
1999	result [0]= F;
2000	return result;
2001	}
2002	if (F.isUnivariate())
2003	{
2004	ASSERT (evalPoints.length() == 1,
2005	"expected an eval point with only one component");
2006	CFArray result= CFArray (size(F));
2007	int j= 0;
2008	CanonicalForm evalPoint= evalPoints.getLast();
2009	for (CFIterator i= F; i.hasTerms(); i++, j++)
2010	result[j]= power (evalPoint, i.exp());
2011	return result;
2012	}
2013	int numMon= size (F);
2014	CFArray result= CFArray (numMon);
2015	int j= 0;
2016	CanonicalForm evalPoint= evalPoints.getLast();
2017	CFList buf= evalPoints;
2018	buf.removeLast();
2019	CFArray recResult;
2020	CanonicalForm powEvalPoint;
2021	for (CFIterator i= F; i.hasTerms(); i++)
2022	{
2023	powEvalPoint= power (evalPoint, i.exp());
2024	recResult= evaluateMonom (i.coeff(), buf);
2025	for (int k= 0; k < recResult.size(); k++)
2026	result[j+k]= powEvalPoint*recResult[k];
2027	j += recResult.size();
2028	}
2029	return result;
2030	}
2031
2032	CFArray
2033	evaluate (const CFArray& A, const CFList& evalPoints)
2034	{
2035	CFArray result= A.size();
2036	CanonicalForm tmp;
2037	int k;
2038	for (int i= 0; i < A.size(); i++)
2039	{
2040	tmp= A[i];
2041	k= 1;
2042	for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2043	tmp= tmp (j.getItem(), k);
2044	result[i]= tmp;
2045	}
2046	return result;
2047	}
2048
2049	CFList
2050	evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
2051	CanonicalForm& Feval, CanonicalForm& Geval,
2052	const CanonicalForm& LCF, const bool& GF,
2053	const Variable& alpha, bool& fail, CFList& list
2054	)
2055	{
2056	int k= tmax (F.level(), G.level()) - 1;
2057	Variable x= Variable (1);
2058	CFList result;
2059	FFRandom genFF;
2060	GFRandom genGF;
2061	int p= getCharacteristic ();
2062	double bound;
2063	if (alpha != Variable (1))
2064	{
2065	bound= pow ((double) p, (double) degree (getMipo(alpha)));
2066	bound= pow (bound, (double) k);
2067	}
2068	else if (GF)
2069	{
2070	bound= pow ((double) p, (double) getGFDegree());
2071	bound= pow ((double) bound, (double) k);
2072	}
2073	else
2074	bound= pow ((double) p, (double) k);
2075
2076	CanonicalForm random;
2077	int j;
2078	bool zeroOneOccured= false;
2079	bool allEqual= false;
2080	CanonicalForm buf;
2081	do
2082	{
2083	random= 0;
2084	// possible overflow if list.length() does not fit into a int
2085	if (list.length() >= bound)
2086	{
2087	fail= true;
2088	break;
2089	}
2090	for (int i= 0; i < k; i++)
2091	{
2092	if (GF)
2093	{
2094	result.append (genGF.generate());
2095	random += result.getLast()*power (x, i);
2096	}
2097	else if (alpha.level() != 1)
2098	{
2099	AlgExtRandomF genAlgExt (alpha);
2100	result.append (genAlgExt.generate());
2101	random += result.getLast()*power (x, i);
2102	}
2103	else
2104	{
2105	result.append (genFF.generate());
2106	random += result.getLast()*power (x, i);
2107	}
2108	if (result.getLast().isOne() \|\| result.getLast().isZero())
2109	zeroOneOccured= true;
2110	}
2111	if (find (list, random))
2112	{
2113	zeroOneOccured= false;
2114	allEqual= false;
2115	result= CFList();
2116	continue;
2117	}
2118	if (zeroOneOccured)
2119	{
2120	list.append (random);
2121	zeroOneOccured= false;
2122	allEqual= false;
2123	result= CFList();
2124	continue;
2125	}
2126	// no zero at this point
2127	if (k > 1)
2128	{
2129	allEqual= true;
2130	CFIterator iter= random;
2131	buf= iter.coeff();
2132	iter++;
2133	for (; iter.hasTerms(); iter++)
2134	if (buf != iter.coeff())
2135	allEqual= false;
2136	}
2137	if (allEqual)
2138	{
2139	list.append (random);
2140	allEqual= false;
2141	zeroOneOccured= false;
2142	result= CFList();
2143	continue;
2144	}
2145
2146	Feval= F;
2147	Geval= G;
2148	CanonicalForm LCeval= LCF;
2149	j= 1;
2150	for (CFListIterator i= result; i.hasItem(); i++, j++)
2151	{
2152	Feval= Feval (i.getItem(), j);
2153	Geval= Geval (i.getItem(), j);
2154	LCeval= LCeval (i.getItem(), j);
2155	}
2156
2157	if (LCeval.isZero())
2158	{
2159	if (!find (list, random))
2160	list.append (random);
2161	zeroOneOccured= false;
2162	allEqual= false;
2163	result= CFList();
2164	continue;
2165	}
2166
2167	if (list.length() >= bound)
2168	{
2169	fail= true;
2170	break;
2171	}
2172	} while (find (list, random));
2173
2174	return result;
2175	}
2176
2177	/// multiply two lists componentwise
2178	void mult (CFList& L1, const CFList& L2)
2179	{
2180	ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2181
2182	CFListIterator j= L2;
2183	for (CFListIterator i= L1; i.hasItem(); i++, j++)
2184	i.getItem() *= j.getItem();
2185	}
2186
2187	void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2188	CanonicalForm& Beval, const CFList& L)
2189	{
2190	Aeval= A;
2191	Beval= B;
2192	int j= 1;
2193	for (CFListIterator i= L; i.hasItem(); i++, j++)
2194	{
2195	Aeval= Aeval (i.getItem(), j);
2196	Beval= Beval (i.getItem(), j);
2197	}
2198	}
2199
2200	CanonicalForm
2201	monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2202	const CanonicalForm& skeleton, const Variable& alpha,
2203	bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2204	)
2205	{
2206	CanonicalForm A= F;
2207	CanonicalForm B= G;
2208	if (F.isZero() && degree(G) > 0) return G/Lc(G);
2209	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2210	else if (F.isZero() && G.isZero()) return F.genOne();
2211	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
2212	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2213	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2214	if (F == G) return F/Lc(F);
2215
2216	ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2217	ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2218
2219	CFMap M,N;
2220	int best_level= myCompress (A, B, M, N, false);
2221
2222	if (best_level == 0)
2223	return B.genOne();
2224
2225	A= M(A);
2226	B= M(B);
2227
2228	Variable x= Variable (1);
2229
2230	//univariate case
2231	if (A.isUnivariate() && B.isUnivariate())
2232	return N (gcd (A, B));
2233
2234	CanonicalForm skel= M(skeleton);
2235	CanonicalForm cA, cB; // content of A and B
2236	CanonicalForm ppA, ppB; // primitive part of A and B
2237	CanonicalForm gcdcAcB;
2238	cA = uni_content (A);
2239	cB = uni_content (B);
2240	gcdcAcB= gcd (cA, cB);
2241	ppA= A/cA;
2242	ppB= B/cB;
2243
2244	CanonicalForm lcA, lcB; // leading coefficients of A and B
2245	CanonicalForm gcdlcAlcB;
2246	lcA= uni_lcoeff (ppA);
2247	lcB= uni_lcoeff (ppB);
2248
2249	if (fdivides (lcA, lcB))
2250	{
2251	if (fdivides (A, B))
2252	return F/Lc(F);
2253	}
2254	if (fdivides (lcB, lcA))
2255	{
2256	if (fdivides (B, A))
2257	return G/Lc(G);
2258	}
2259
2260	gcdlcAlcB= gcd (lcA, lcB);
2261	int skelSize= size (skel, skel.mvar());
2262
2263	int j= 0;
2264	int biggestSize= 0;
2265
2266	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2267	biggestSize= tmax (biggestSize, size (i.coeff()));
2268
2269	CanonicalForm g, Aeval, Beval;
2270
2271	CFList evalPoints;
2272	bool evalFail= false;
2273	CFList list;
2274	bool GF= false;
2275	CanonicalForm LCA= LC (A);
2276	CanonicalForm tmp;
2277	CFArray gcds= CFArray (biggestSize);
2278	CFList * pEvalPoints= new CFList [biggestSize];
2279	Variable V_buf= alpha, V_buf4= alpha;
2280	CFList source, dest;
2281	CanonicalForm prim_elem, im_prim_elem;
2282	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2283	for (int i= 0; i < biggestSize; i++)
2284	{
2285	if (i == 0)
2286	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2287	list);
2288	else
2289	{
2290	mult (evalPoints, pEvalPoints [0]);
2291	eval (A, B, Aeval, Beval, evalPoints);
2292	}
2293
2294	if (evalFail)
2295	{
2296	if (V_buf.level() != 1)
2297	{
2298	do
2299	{
2300	Variable V_buf3= V_buf;
2301	V_buf= chooseExtension (V_buf);
2302	source= CFList();
2303	dest= CFList();
2304
2305	bool prim_fail= false;
2306	Variable V_buf2;
2307	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2308	if (V_buf4 == alpha && alpha.level() != 1)
2309	prim_elem_alpha= prim_elem;
2310
2311	ASSERT (!prim_fail, "failure in integer factorizer");
2312	if (prim_fail)
2313	; //ERROR
2314	else
2315	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2316
2317	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2318	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2319
2320	if (V_buf4 == alpha && alpha.level() != 1)
2321	im_prim_elem_alpha= im_prim_elem;
2322	else if (alpha.level() != 1)
2323	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2324	prim_elem, im_prim_elem, source, dest);
2325
2326	for (CFListIterator j= list; j.hasItem(); j++)
2327	j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2328	im_prim_elem, source, dest);
2329	for (int k= 0; k < i; k++)
2330	{
2331	for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2332	j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2333	im_prim_elem, source, dest);
2334	gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2335	source, dest);
2336	}
2337
2338	if (alpha.level() != 1)
2339	{
2340	A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2341	B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2342	}
2343	V_buf4= V_buf;
2344	evalFail= false;
2345	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2346	evalFail, list);
2347	} while (evalFail);
2348	}
2349	else
2350	{
2351	CanonicalForm mipo;
2352	int deg= 2;
2353	bool initialized= false;
2354	do
2355	{
2356	mipo= randomIrredpoly (deg, x);
2357	if (initialized)
2358	setMipo (V_buf, mipo);
2359	else
2360	V_buf= rootOf (mipo);
2361	evalFail= false;
2362	initialized= true;
2363	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2364	evalFail, list);
2365	deg++;
2366	} while (evalFail);
2367	V_buf4= V_buf;
2368	}
2369	}
2370
2371	g= gcd (Aeval, Beval);
2372	g /= Lc (g);
2373
2374	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2375	{
2376	delete[] pEvalPoints;
2377	fail= true;
2378	if (alpha.level() != 1 && V_buf != alpha)
2379	prune1 (alpha);
2380	return 0;
2381	}
2382	CFIterator l= skel;
2383	for (CFIterator k= g; k.hasTerms(); k++, l++)
2384	{
2385	if (k.exp() != l.exp())
2386	{
2387	delete[] pEvalPoints;
2388	fail= true;
2389	if (alpha.level() != 1 && V_buf != alpha)
2390	prune1 (alpha);
2391	return 0;
2392	}
2393	}
2394	pEvalPoints[i]= evalPoints;
2395	gcds[i]= g;
2396
2397	tmp= 0;
2398	int j= 0;
2399	for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2400	tmp += k.getItem()*power (x, j);
2401	list.append (tmp);
2402	}
2403
2404	if (Monoms.size() == 0)
2405	Monoms= getMonoms (skel);
2406
2407	coeffMonoms= new CFArray [skelSize];
2408	j= 0;
2409	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2410	coeffMonoms[j]= getMonoms (i.coeff());
2411
2412	CFArray* pL= new CFArray [skelSize];
2413	CFArray* pM= new CFArray [skelSize];
2414	for (int i= 0; i < biggestSize; i++)
2415	{
2416	CFIterator l= gcds [i];
2417	evalPoints= pEvalPoints [i];
2418	for (int k= 0; k < skelSize; k++, l++)
2419	{
2420	if (i == 0)
2421	pL[k]= CFArray (biggestSize);
2422	pL[k] [i]= l.coeff();
2423
2424	if (i == 0)
2425	pM[k]= evaluate (coeffMonoms [k], evalPoints);
2426	}
2427	}
2428
2429	CFArray solution;
2430	CanonicalForm result= 0;
2431	int ind= 0;
2432	CFArray bufArray;
2433	CFMatrix Mat;
2434	for (int k= 0; k < skelSize; k++)
2435	{
2436	if (biggestSize != coeffMonoms[k].size())
2437	{
2438	bufArray= CFArray (coeffMonoms[k].size());
2439	for (int i= 0; i < coeffMonoms[k].size(); i++)
2440	bufArray [i]= pL[k] [i];
2441	solution= solveGeneralVandermonde (pM [k], bufArray);
2442	}
2443	else
2444	solution= solveGeneralVandermonde (pM [k], pL[k]);
2445
2446	if (solution.size() == 0)
2447	{
2448	delete[] pEvalPoints;
2449	delete[] pM;
2450	delete[] pL;
2451	delete[] coeffMonoms;
2452	fail= true;
2453	if (alpha.level() != 1 && V_buf != alpha)
2454	prune1 (alpha);
2455	return 0;
2456	}
2457	for (int l= 0; l < solution.size(); l++)
2458	result += solution[l]*Monoms [ind + l];
2459	ind += solution.size();
2460	}
2461
2462	delete[] pEvalPoints;
2463	delete[] pM;
2464	delete[] pL;
2465	delete[] coeffMonoms;
2466
2467	if (alpha.level() != 1 && V_buf != alpha)
2468	{
2469	CFList u, v;
2470	result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2471	prune1 (alpha);
2472	}
2473
2474	result= N(result);
2475	if (fdivides (result, F) && fdivides (result, G))
2476	return result;
2477	else
2478	{
2479	fail= true;
2480	return 0;
2481	}
2482	}
2483
2484	CanonicalForm
2485	nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2486	const CanonicalForm& skeleton, const Variable& alpha,
2487	bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2488	)
2489	{
2490	CanonicalForm A= F;
2491	CanonicalForm B= G;
2492	if (F.isZero() && degree(G) > 0) return G/Lc(G);
2493	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2494	else if (F.isZero() && G.isZero()) return F.genOne();
2495	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
2496	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2497	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2498	if (F == G) return F/Lc(F);
2499
2500	ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2501	ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2502
2503	CFMap M,N;
2504	int best_level= myCompress (A, B, M, N, false);
2505
2506	if (best_level == 0)
2507	return B.genOne();
2508
2509	A= M(A);
2510	B= M(B);
2511
2512	Variable x= Variable (1);
2513
2514	//univariate case
2515	if (A.isUnivariate() && B.isUnivariate())
2516	return N (gcd (A, B));
2517
2518	CanonicalForm skel= M(skeleton);
2519
2520	CanonicalForm cA, cB; // content of A and B
2521	CanonicalForm ppA, ppB; // primitive part of A and B
2522	CanonicalForm gcdcAcB;
2523	cA = uni_content (A);
2524	cB = uni_content (B);
2525	gcdcAcB= gcd (cA, cB);
2526	ppA= A/cA;
2527	ppB= B/cB;
2528
2529	CanonicalForm lcA, lcB; // leading coefficients of A and B
2530	CanonicalForm gcdlcAlcB;
2531	lcA= uni_lcoeff (ppA);
2532	lcB= uni_lcoeff (ppB);
2533
2534	if (fdivides (lcA, lcB))
2535	{
2536	if (fdivides (A, B))
2537	return F/Lc(F);
2538	}
2539	if (fdivides (lcB, lcA))
2540	{
2541	if (fdivides (B, A))
2542	return G/Lc(G);
2543	}
2544
2545	gcdlcAlcB= gcd (lcA, lcB);
2546	int skelSize= size (skel, skel.mvar());
2547
2548	int j= 0;
2549	int biggestSize= 0;
2550	int bufSize;
2551	int numberUni= 0;
2552	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2553	{
2554	bufSize= size (i.coeff());
2555	biggestSize= tmax (biggestSize, bufSize);
2556	numberUni += bufSize;
2557	}
2558	numberUni--;
2559	numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2560	biggestSize= tmax (biggestSize , numberUni);
2561
2562	numberUni= biggestSize + size (LC(skel)) - 1;
2563	int biggestSize2= tmax (numberUni, biggestSize);
2564
2565	CanonicalForm g, Aeval, Beval;
2566
2567	CFList evalPoints;
2568	CFArray coeffEval;
2569	bool evalFail= false;
2570	CFList list;
2571	bool GF= false;
2572	CanonicalForm LCA= LC (A);
2573	CanonicalForm tmp;
2574	CFArray gcds= CFArray (biggestSize);
2575	CFList * pEvalPoints= new CFList [biggestSize];
2576	Variable V_buf= alpha, V_buf4= alpha;
2577	CFList source, dest;
2578	CanonicalForm prim_elem, im_prim_elem;
2579	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2580	for (int i= 0; i < biggestSize; i++)
2581	{
2582	if (i == 0)
2583	{
2584	if (getCharacteristic() > 3)
2585	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2586	evalFail, list);
2587	else
2588	evalFail= true;
2589
2590	if (evalFail)
2591	{
2592	if (V_buf.level() != 1)
2593	{
2594	do
2595	{
2596	Variable V_buf3= V_buf;
2597	V_buf= chooseExtension (V_buf);
2598	source= CFList();
2599	dest= CFList();
2600
2601	bool prim_fail= false;
2602	Variable V_buf2;
2603	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2604	if (V_buf4 == alpha && alpha.level() != 1)
2605	prim_elem_alpha= prim_elem;
2606
2607	ASSERT (!prim_fail, "failure in integer factorizer");
2608	if (prim_fail)
2609	; //ERROR
2610	else
2611	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2612
2613	DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2614	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2615
2616	if (V_buf4 == alpha && alpha.level() != 1)
2617	im_prim_elem_alpha= im_prim_elem;
2618	else if (alpha.level() != 1)
2619	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2620	prim_elem, im_prim_elem, source, dest);
2621
2622	for (CFListIterator i= list; i.hasItem(); i++)
2623	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2624	im_prim_elem, source, dest);
2625	if (alpha.level() != 1)
2626	{
2627	A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2628	B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2629	}
2630	evalFail= false;
2631	V_buf4= V_buf;
2632	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2633	evalFail, list);
2634	} while (evalFail);
2635	}
2636	else
2637	{
2638	CanonicalForm mipo;
2639	int deg= 2;
2640	bool initialized= false;
2641	do
2642	{
2643	mipo= randomIrredpoly (deg, x);
2644	if (initialized)
2645	setMipo (V_buf, mipo);
2646	else
2647	V_buf= rootOf (mipo);
2648	evalFail= false;
2649	initialized= true;
2650	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2651	evalFail, list);
2652	deg++;
2653	} while (evalFail);
2654	V_buf4= V_buf;
2655	}
2656	}
2657	}
2658	else
2659	{
2660	mult (evalPoints, pEvalPoints[0]);
2661	eval (A, B, Aeval, Beval, evalPoints);
2662	}
2663
2664	g= gcd (Aeval, Beval);
2665	g /= Lc (g);
2666
2667	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2668	{
2669	delete[] pEvalPoints;
2670	fail= true;
2671	if (alpha.level() != 1 && V_buf != alpha)
2672	prune1 (alpha);
2673	return 0;
2674	}
2675	CFIterator l= skel;
2676	for (CFIterator k= g; k.hasTerms(); k++, l++)
2677	{
2678	if (k.exp() != l.exp())
2679	{
2680	delete[] pEvalPoints;
2681	fail= true;
2682	if (alpha.level() != 1 && V_buf != alpha)
2683	prune1 (alpha);
2684	return 0;
2685	}
2686	}
2687	pEvalPoints[i]= evalPoints;
2688	gcds[i]= g;
2689
2690	tmp= 0;
2691	int j= 0;
2692	for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2693	tmp += k.getItem()*power (x, j);
2694	list.append (tmp);
2695	}
2696
2697	if (Monoms.size() == 0)
2698	Monoms= getMonoms (skel);
2699
2700	coeffMonoms= new CFArray [skelSize];
2701
2702	j= 0;
2703	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2704	coeffMonoms[j]= getMonoms (i.coeff());
2705
2706	int minimalColumnsIndex;
2707	if (skelSize > 1)
2708	minimalColumnsIndex= 1;
2709	else
2710	minimalColumnsIndex= 0;
2711	int minimalColumns=-1;
2712
2713	CFArray* pM= new CFArray [skelSize];
2714	CFMatrix Mat;
2715	// find the Matrix with minimal number of columns
2716	for (int i= 0; i < skelSize; i++)
2717	{
2718	pM[i]= CFArray (coeffMonoms[i].size());
2719	if (i == 1)
2720	minimalColumns= coeffMonoms[i].size();
2721	if (i > 1)
2722	{
2723	if (minimalColumns > coeffMonoms[i].size())
2724	{
2725	minimalColumns= coeffMonoms[i].size();
2726	minimalColumnsIndex= i;
2727	}
2728	}
2729	}
2730	CFMatrix* pMat= new CFMatrix [2];
2731	pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2732	pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2733	CFArray* pL= new CFArray [skelSize];
2734	for (int i= 0; i < biggestSize; i++)
2735	{
2736	CFIterator l= gcds [i];
2737	evalPoints= pEvalPoints [i];
2738	for (int k= 0; k < skelSize; k++, l++)
2739	{
2740	if (i == 0)
2741	pL[k]= CFArray (biggestSize);
2742	pL[k] [i]= l.coeff();
2743
2744	if (i == 0 && k != 0 && k != minimalColumnsIndex)
2745	pM[k]= evaluate (coeffMonoms[k], evalPoints);
2746	else if (i == 0 && (k == 0 \|\| k == minimalColumnsIndex))
2747	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2748	if (i > 0 && (k == 0 \|\| k == minimalColumnsIndex))
2749	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2750
2751	if (k == 0)
2752	{
2753	if (pMat[k].rows() >= i + 1)
2754	{
2755	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2756	pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2757	}
2758	}
2759	if (k == minimalColumnsIndex)
2760	{
2761	if (pMat[1].rows() >= i + 1)
2762	{
2763	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2764	pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2765	}
2766	}
2767	}
2768	}
2769
2770	CFArray solution;
2771	CanonicalForm result= 0;
2772	int ind= 1;
2773	int matRows, matColumns;
2774	matRows= pMat[1].rows();
2775	matColumns= pMat[0].columns() - 1;
2776	matColumns += pMat[1].columns();
2777
2778	Mat= CFMatrix (matRows, matColumns);
2779	for (int i= 1; i <= matRows; i++)
2780	for (int j= 1; j <= pMat[1].columns(); j++)
2781	Mat (i, j)= pMat[1] (i, j);
2782
2783	for (int j= pMat[1].columns() + 1; j <= matColumns;
2784	j++, ind++)
2785	{
2786	for (int i= 1; i <= matRows; i++)
2787	Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2788	}
2789
2790	CFArray firstColumn= CFArray (pMat[0].rows());
2791	for (int i= 0; i < pMat[0].rows(); i++)
2792	firstColumn [i]= pMat[0] (i + 1, 1);
2793
2794	CFArray bufArray= pL[minimalColumnsIndex];
2795
2796	for (int i= 0; i < biggestSize; i++)
2797	pL[minimalColumnsIndex] [i] *= firstColumn[i];
2798
2799	CFMatrix bufMat= pMat[1];
2800	pMat[1]= Mat;
2801
2802	if (V_buf.level() != 1)
2803	solution= solveSystemFq (pMat[1],
2804	pL[minimalColumnsIndex], V_buf);
2805	else
2806	solution= solveSystemFp (pMat[1],
2807	pL[minimalColumnsIndex]);
2808
2809	if (solution.size() == 0)
2810	{ //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2811	CFMatrix bufMat0= pMat[0];
2812	delete [] pMat;
2813	pMat= new CFMatrix [skelSize];
2814	pL[minimalColumnsIndex]= bufArray;
2815	CFList* bufpEvalPoints= NULL;
2816	CFArray bufGcds;
2817	if (biggestSize != biggestSize2)
2818	{
2819	bufpEvalPoints= pEvalPoints;
2820	pEvalPoints= new CFList [biggestSize2];
2821	bufGcds= gcds;
2822	gcds= CFArray (biggestSize2);
2823	for (int i= 0; i < biggestSize; i++)
2824	{
2825	pEvalPoints[i]= bufpEvalPoints [i];
2826	gcds[i]= bufGcds[i];
2827	}
2828	for (int i= 0; i < biggestSize2 - biggestSize; i++)
2829	{
2830	mult (evalPoints, pEvalPoints[0]);
2831	eval (A, B, Aeval, Beval, evalPoints);
2832	g= gcd (Aeval, Beval);
2833	g /= Lc (g);
2834	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2835	{
2836	delete[] pEvalPoints;
2837	delete[] pMat;
2838	delete[] pL;
2839	delete[] coeffMonoms;
2840	delete[] pM;
2841	if (bufpEvalPoints != NULL)
2842	delete [] bufpEvalPoints;
2843	fail= true;
2844	if (alpha.level() != 1 && V_buf != alpha)
2845	prune1 (alpha);
2846	return 0;
2847	}
2848	CFIterator l= skel;
2849	for (CFIterator k= g; k.hasTerms(); k++, l++)
2850	{
2851	if (k.exp() != l.exp())
2852	{
2853	delete[] pEvalPoints;
2854	delete[] pMat;
2855	delete[] pL;
2856	delete[] coeffMonoms;
2857	delete[] pM;
2858	if (bufpEvalPoints != NULL)
2859	delete [] bufpEvalPoints;
2860	fail= true;
2861	if (alpha.level() != 1 && V_buf != alpha)
2862	prune1 (alpha);
2863	return 0;
2864	}
2865	}
2866	pEvalPoints[i + biggestSize]= evalPoints;
2867	gcds[i + biggestSize]= g;
2868	}
2869	}
2870	for (int i= 0; i < biggestSize; i++)
2871	{
2872	CFIterator l= gcds [i];
2873	evalPoints= pEvalPoints [i];
2874	for (int k= 1; k < skelSize; k++, l++)
2875	{
2876	if (i == 0)
2877	pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2878	if (k == minimalColumnsIndex)
2879	continue;
2880	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2881	if (pMat[k].rows() >= i + 1)
2882	{
2883	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2884	pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2885	}
2886	}
2887	}
2888	Mat= bufMat0;
2889	pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2890	for (int j= 1; j <= Mat.rows(); j++)
2891	for (int k= 1; k <= Mat.columns(); k++)
2892	pMat [0] (j,k)= Mat (j,k);
2893	Mat= bufMat;
2894	for (int j= 1; j <= Mat.rows(); j++)
2895	for (int k= 1; k <= Mat.columns(); k++)
2896	pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2897	// write old matrix entries into new matrices
2898	for (int i= 0; i < skelSize; i++)
2899	{
2900	bufArray= pL[i];
2901	pL[i]= CFArray (biggestSize2);
2902	for (int j= 0; j < bufArray.size(); j++)
2903	pL[i] [j]= bufArray [j];
2904	}
2905	//write old vector entries into new and add new entries to old matrices
2906	for (int i= 0; i < biggestSize2 - biggestSize; i++)
2907	{
2908	CFIterator l= gcds [i + biggestSize];
2909	evalPoints= pEvalPoints [i + biggestSize];
2910	for (int k= 0; k < skelSize; k++, l++)
2911	{
2912	pL[k] [i + biggestSize]= l.coeff();
2913	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2914	if (pMat[k].rows() >= i + biggestSize + 1)
2915	{
2916	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2917	pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2918	}
2919	}
2920	}
2921	// begin new
2922	for (int i= 0; i < skelSize; i++)
2923	{
2924	if (pL[i].size() > 1)
2925	{
2926	for (int j= 2; j <= pMat[i].rows(); j++)
2927	pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2928	-pL[i] [j - 1];
2929	}
2930	}
2931
2932	matColumns= biggestSize2 - 1;
2933	matRows= 0;
2934	for (int i= 0; i < skelSize; i++)
2935	{
2936	if (V_buf.level() == 1)
2937	(void) gaussianElimFp (pMat[i], pL[i]);
2938	else
2939	(void) gaussianElimFq (pMat[i], pL[i], V_buf);
2940
2941	if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2942	{
2943	delete[] pEvalPoints;
2944	delete[] pMat;
2945	delete[] pL;
2946	delete[] coeffMonoms;
2947	delete[] pM;
2948	if (bufpEvalPoints != NULL)
2949	delete [] bufpEvalPoints;
2950	fail= true;
2951	if (alpha.level() != 1 && V_buf != alpha)
2952	prune1 (alpha);
2953	return 0;
2954	}
2955	matRows += pMat[i].rows() - coeffMonoms[i].size();
2956	}
2957
2958	CFMatrix bufMat;
2959	Mat= CFMatrix (matRows, matColumns);
2960	ind= 0;
2961	bufArray= CFArray (matRows);
2962	CFArray bufArray2;
2963	for (int i= 0; i < skelSize; i++)
2964	{
2965	if (coeffMonoms[i].size() + 1 >= pMat[i].rows() \|\| coeffMonoms[i].size() + 1 >= pMat[i].columns())
2966	{
2967	delete[] pEvalPoints;
2968	delete[] pMat;
2969	delete[] pL;
2970	delete[] coeffMonoms;
2971	delete[] pM;
2972	if (bufpEvalPoints != NULL)
2973	delete [] bufpEvalPoints;
2974	fail= true;
2975	if (alpha.level() != 1 && V_buf != alpha)
2976	prune1 (alpha);
2977	return 0;
2978	}
2979	bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2980	coeffMonoms[i].size() + 1, pMat[i].columns());
2981
2982	for (int j= 1; j <= bufMat.rows(); j++)
2983	for (int k= 1; k <= bufMat.columns(); k++)
2984	Mat (j + ind, k)= bufMat(j, k);
2985	bufArray2= coeffMonoms[i].size();
2986	for (int j= 1; j <= pMat[i].rows(); j++)
2987	{
2988	if (j > coeffMonoms[i].size())
2989	bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2990	else
2991	bufArray2 [j - 1]= pL[i] [j - 1];
2992	}
2993	pL[i]= bufArray2;
2994	ind += bufMat.rows();
2995	pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2996	}
2997
2998	if (V_buf.level() != 1)
2999	solution= solveSystemFq (Mat, bufArray, V_buf);
3000	else
3001	solution= solveSystemFp (Mat, bufArray);
3002
3003	if (solution.size() == 0)
3004	{
3005	delete[] pEvalPoints;
3006	delete[] pMat;
3007	delete[] pL;
3008	delete[] coeffMonoms;
3009	delete[] pM;
3010	if (bufpEvalPoints != NULL)
3011	delete [] bufpEvalPoints;
3012	fail= true;
3013	if (alpha.level() != 1 && V_buf != alpha)
3014	prune1 (alpha);
3015	return 0;
3016	}
3017
3018	ind= 0;
3019	result= 0;
3020	CFArray bufSolution;
3021	for (int i= 0; i < skelSize; i++)
3022	{
3023	bufSolution= readOffSolution (pMat[i], pL[i], solution);
3024	for (int i= 0; i < bufSolution.size(); i++, ind++)
3025	result += Monoms [ind]*bufSolution[i];
3026	}
3027	if (alpha.level() != 1 && V_buf != alpha)
3028	{
3029	CFList u, v;
3030	result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3031	prune1 (alpha);
3032	}
3033	result= N(result);
3034	delete[] pEvalPoints;
3035	delete[] pMat;
3036	delete[] pL;
3037	delete[] coeffMonoms;
3038	delete[] pM;
3039
3040	if (bufpEvalPoints != NULL)
3041	delete [] bufpEvalPoints;
3042	if (fdivides (result, F) && fdivides (result, G))
3043	return result;
3044	else
3045	{
3046	fail= true;
3047	return 0;
3048	}
3049	} // end of deKleine, Monagan & Wittkopf
3050
3051	result += Monoms[0];
3052	int ind2= 0, ind3= 2;
3053	ind= 0;
3054	for (int l= 0; l < minimalColumnsIndex; l++)
3055	ind += coeffMonoms[l].size();
3056	for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3057	l++, ind2++, ind3++)
3058	{
3059	result += solution[l]*Monoms [1 + ind2];
3060	for (int i= 0; i < pMat[0].rows(); i++)
3061	firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3062	}
3063	for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3064	result += solution[l]*Monoms [ind + l];
3065	ind= coeffMonoms[0].size();
3066	for (int k= 1; k < skelSize; k++)
3067	{
3068	if (k == minimalColumnsIndex)
3069	{
3070	ind += coeffMonoms[k].size();
3071	continue;
3072	}
3073	if (k != minimalColumnsIndex)
3074	{
3075	for (int i= 0; i < biggestSize; i++)
3076	pL[k] [i] *= firstColumn [i];
3077	}
3078
3079	if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3080	{
3081	bufArray= CFArray (coeffMonoms[k].size());
3082	for (int i= 0; i < bufArray.size(); i++)
3083	bufArray [i]= pL[k] [i];
3084	solution= solveGeneralVandermonde (pM[k], bufArray);
3085	}
3086	else
3087	solution= solveGeneralVandermonde (pM[k], pL[k]);
3088
3089	if (solution.size() == 0)
3090	{
3091	delete[] pEvalPoints;
3092	delete[] pMat;
3093	delete[] pL;
3094	delete[] coeffMonoms;
3095	delete[] pM;
3096	fail= true;
3097	if (alpha.level() != 1 && V_buf != alpha)
3098	prune1 (alpha);
3099	return 0;
3100	}
3101	if (k != minimalColumnsIndex)
3102	{
3103	for (int l= 0; l < solution.size(); l++)
3104	result += solution[l]*Monoms [ind + l];
3105	ind += solution.size();
3106	}
3107	}
3108
3109	delete[] pEvalPoints;
3110	delete[] pMat;
3111	delete[] pL;
3112	delete[] pM;
3113	delete[] coeffMonoms;
3114
3115	if (alpha.level() != 1 && V_buf != alpha)
3116	{
3117	CFList u, v;
3118	result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3119	prune1 (alpha);
3120	}
3121	result= N(result);
3122
3123	if (fdivides (result, F) && fdivides (result, G))
3124	return result;
3125	else
3126	{
3127	fail= true;
3128	return 0;
3129	}
3130	}
3131
3132	CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
3133	const Variable & alpha, CFList& l, bool& topLevel)
3134	{
3135	CanonicalForm A= F;
3136	CanonicalForm B= G;
3137	if (F.isZero() && degree(G) > 0) return G/Lc(G);
3138	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3139	else if (F.isZero() && G.isZero()) return F.genOne();
3140	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
3141	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3142	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3143	if (F == G) return F/Lc(F);
3144
3145	CFMap M,N;
3146	int best_level= myCompress (A, B, M, N, topLevel);
3147
3148	if (best_level == 0) return B.genOne();
3149
3150	A= M(A);
3151	B= M(B);
3152
3153	Variable x= Variable (1);
3154
3155	//univariate case
3156	if (A.isUnivariate() && B.isUnivariate())
3157	return N (gcd (A, B));
3158
3159	CanonicalForm cA, cB; // content of A and B
3160	CanonicalForm ppA, ppB; // primitive part of A and B
3161	CanonicalForm gcdcAcB;
3162
3163	cA = uni_content (A);
3164	cB = uni_content (B);
3165	gcdcAcB= gcd (cA, cB);
3166	ppA= A/cA;
3167	ppB= B/cB;
3168
3169	CanonicalForm lcA, lcB; // leading coefficients of A and B
3170	CanonicalForm gcdlcAlcB;
3171	lcA= uni_lcoeff (ppA);
3172	lcB= uni_lcoeff (ppB);
3173
3174	if (fdivides (lcA, lcB))
3175	{
3176	if (fdivides (A, B))
3177	return F/Lc(F);
3178	}
3179	if (fdivides (lcB, lcA))
3180	{
3181	if (fdivides (B, A))
3182	return G/Lc(G);
3183	}
3184
3185	gcdlcAlcB= gcd (lcA, lcB);
3186
3187	int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3188	int d0;
3189
3190	if (d == 0)
3191	return N(gcdcAcB);
3192	d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3193
3194	if (d0 < d)
3195	d= d0;
3196
3197	if (d == 0)
3198	return N(gcdcAcB);
3199
3200	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3201	CanonicalForm newtonPoly= 1;
3202	m= gcdlcAlcB;
3203	G_m= 0;
3204	H= 0;
3205	bool fail= false;
3206	topLevel= false;
3207	bool inextension= false;
3208	Variable V_buf= alpha, V_buf4= alpha;
3209	CanonicalForm prim_elem, im_prim_elem;
3210	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3211	CFList source, dest;
3212	do // first do
3213	{
3214	random_element= randomElement (m, V_buf, l, fail);
3215	if (random_element == 0 && !fail)
3216	{
3217	if (!find (l, random_element))
3218	l.append (random_element);
3219	continue;
3220	}
3221	if (fail)
3222	{
3223	source= CFList();
3224	dest= CFList();
3225
3226	Variable V_buf3= V_buf;
3227	V_buf= chooseExtension (V_buf);
3228	bool prim_fail= false;
3229	Variable V_buf2;
3230	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3231	if (V_buf4 == alpha)
3232	prim_elem_alpha= prim_elem;
3233
3234	if (V_buf3 != V_buf4)
3235	{
3236	m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3237	G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3238	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3239	source, dest);
3240	ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3241	ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3242	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3243	dest);
3244	for (CFListIterator i= l; i.hasItem(); i++)
3245	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3246	source, dest);
3247	}
3248
3249	ASSERT (!prim_fail, "failure in integer factorizer");
3250	if (prim_fail)
3251	; //ERROR
3252	else
3253	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3254
3255	if (V_buf4 == alpha)
3256	im_prim_elem_alpha= im_prim_elem;
3257	else
3258	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3259	im_prim_elem, source, dest);
3260
3261	DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3262	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3263	inextension= true;
3264	for (CFListIterator i= l; i.hasItem(); i++)
3265	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3266	im_prim_elem, source, dest);
3267	m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3268	G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3269	newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3270	source, dest);
3271	ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3272	ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3273	gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3274	source, dest);
3275
3276	fail= false;
3277	random_element= randomElement (m, V_buf, l, fail );
3278	DEBOUTLN (cerr, "fail= " << fail);
3279	CFList list;
3280	TIMING_START (gcd_recursion);
3281	G_random_element=
3282	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3283	list, topLevel);
3284	TIMING_END_AND_PRINT (gcd_recursion,
3285	"time for recursive call: ");
3286	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3287	V_buf4= V_buf;
3288	}
3289	else
3290	{
3291	CFList list;
3292	TIMING_START (gcd_recursion);
3293	G_random_element=
3294	sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3295	list, topLevel);
3296	TIMING_END_AND_PRINT (gcd_recursion,
3297	"time for recursive call: ");
3298	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3299	}
3300
3301	if (!G_random_element.inCoeffDomain())
3302	d0= totaldegree (G_random_element, Variable(2),
3303	Variable (G_random_element.level()));
3304	else
3305	d0= 0;
3306
3307	if (d0 == 0)
3308	{
3309	if (inextension)
3310	prune1 (alpha);
3311	return N(gcdcAcB);
3312	}
3313	if (d0 > d)
3314	{
3315	if (!find (l, random_element))
3316	l.append (random_element);
3317	continue;
3318	}
3319
3320	G_random_element=
3321	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3322	* G_random_element;
3323
3324	skeleton= G_random_element;
3325	if (!G_random_element.inCoeffDomain())
3326	d0= totaldegree (G_random_element, Variable(2),
3327	Variable (G_random_element.level()));
3328	else
3329	d0= 0;
3330
3331	if (d0 < d)
3332	{
3333	m= gcdlcAlcB;
3334	newtonPoly= 1;
3335	G_m= 0;
3336	d= d0;
3337	}
3338
3339	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3340	if (uni_lcoeff (H) == gcdlcAlcB)
3341	{
3342	cH= uni_content (H);
3343	ppH= H/cH;
3344	if (inextension)
3345	{
3346	CFList u, v;
3347	//maybe it's better to test if ppH is an element of F(\alpha) before
3348	//mapping down
3349	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3350	{
3351	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3352	ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3353	ppH /= Lc(ppH);
3354	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3355	prune1 (alpha);
3356	return N(gcdcAcB*ppH);
3357	}
3358	}
3359	else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3360	return N(gcdcAcB*ppH);
3361	}
3362	G_m= H;
3363	newtonPoly= newtonPoly*(x - random_element);
3364	m= m*(x - random_element);
3365	if (!find (l, random_element))
3366	l.append (random_element);
3367
3368	if (getCharacteristic () > 3 && size (skeleton) < 100)
3369	{
3370	CFArray Monoms;
3371	CFArray *coeffMonoms;
3372	do //second do
3373	{
3374	random_element= randomElement (m, V_buf, l, fail);
3375	if (random_element == 0 && !fail)
3376	{
3377	if (!find (l, random_element))
3378	l.append (random_element);
3379	continue;
3380	}
3381	if (fail)
3382	{
3383	source= CFList();
3384	dest= CFList();
3385
3386	Variable V_buf3= V_buf;
3387	V_buf= chooseExtension (V_buf);
3388	bool prim_fail= false;
3389	Variable V_buf2;
3390	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3391	if (V_buf4 == alpha)
3392	prim_elem_alpha= prim_elem;
3393
3394	if (V_buf3 != V_buf4)
3395	{
3396	m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3397	G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3398	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3399	source, dest);
3400	ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3401	ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3402	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3403	source, dest);
3404	for (CFListIterator i= l; i.hasItem(); i++)
3405	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3406	source, dest);
3407	}
3408
3409	ASSERT (!prim_fail, "failure in integer factorizer");
3410	if (prim_fail)
3411	; //ERROR
3412	else
3413	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3414
3415	if (V_buf4 == alpha)
3416	im_prim_elem_alpha= im_prim_elem;
3417	else
3418	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3419	prim_elem, im_prim_elem, source, dest);
3420
3421	DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3422	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3423	inextension= true;
3424	for (CFListIterator i= l; i.hasItem(); i++)
3425	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3426	im_prim_elem, source, dest);
3427	m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3428	G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3429	newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3430	source, dest);
3431	ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3432	ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3433
3434	gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3435	source, dest);
3436
3437	fail= false;
3438	random_element= randomElement (m, V_buf, l, fail);
3439	DEBOUTLN (cerr, "fail= " << fail);
3440	CFList list;
3441	TIMING_START (gcd_recursion);
3442
3443	V_buf4= V_buf;
3444
3445	//sparseInterpolation
3446	bool sparseFail= false;
3447	if (LC (skeleton).inCoeffDomain())
3448	G_random_element=
3449	monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3450	skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3451	else
3452	G_random_element=
3453	nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3454	skeleton, V_buf, sparseFail, coeffMonoms,
3455	Monoms);
3456	if (sparseFail)
3457	break;
3458
3459	TIMING_END_AND_PRINT (gcd_recursion,
3460	"time for recursive call: ");
3461	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3462	}
3463	else
3464	{
3465	CFList list;
3466	TIMING_START (gcd_recursion);
3467	bool sparseFail= false;
3468	if (LC (skeleton).inCoeffDomain())
3469	G_random_element=
3470	monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3471	skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3472	else
3473	G_random_element=
3474	nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3475	skeleton, V_buf, sparseFail, coeffMonoms,
3476	Monoms);
3477	if (sparseFail)
3478	break;
3479
3480	TIMING_END_AND_PRINT (gcd_recursion,
3481	"time for recursive call: ");
3482	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3483	}
3484
3485	if (!G_random_element.inCoeffDomain())
3486	d0= totaldegree (G_random_element, Variable(2),
3487	Variable (G_random_element.level()));
3488	else
3489	d0= 0;
3490
3491	if (d0 == 0)
3492	{
3493	if (inextension)
3494	prune1 (alpha);
3495	return N(gcdcAcB);
3496	}
3497	if (d0 > d)
3498	{
3499	if (!find (l, random_element))
3500	l.append (random_element);
3501	continue;
3502	}
3503
3504	G_random_element=
3505	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3506	* G_random_element;
3507
3508	if (!G_random_element.inCoeffDomain())
3509	d0= totaldegree (G_random_element, Variable(2),
3510	Variable (G_random_element.level()));
3511	else
3512	d0= 0;
3513
3514	if (d0 < d)
3515	{
3516	m= gcdlcAlcB;
3517	newtonPoly= 1;
3518	G_m= 0;
3519	d= d0;
3520	}
3521
3522	TIMING_START (newton_interpolation);
3523	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3524	TIMING_END_AND_PRINT (newton_interpolation,
3525	"time for newton interpolation: ");
3526
3527	//termination test
3528	if (uni_lcoeff (H) == gcdlcAlcB)
3529	{
3530	cH= uni_content (H);
3531	ppH= H/cH;
3532	if (inextension)
3533	{
3534	CFList u, v;
3535	//maybe it's better to test if ppH is an element of F(\alpha) before
3536	//mapping down
3537	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3538	{
3539	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3540	ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3541	ppH /= Lc(ppH);
3542	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3543	prune1 (alpha);
3544	return N(gcdcAcB*ppH);
3545	}
3546	}
3547	else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3548	{
3549	return N(gcdcAcB*ppH);
3550	}
3551	}
3552
3553	G_m= H;
3554	newtonPoly= newtonPoly*(x - random_element);
3555	m= m*(x - random_element);
3556	if (!find (l, random_element))
3557	l.append (random_element);
3558
3559	} while (1);
3560	}
3561	} while (1);
3562	}
3563
3564	CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3565	bool& topLevel, CFList& l)
3566	{
3567	CanonicalForm A= F;
3568	CanonicalForm B= G;
3569	if (F.isZero() && degree(G) > 0) return G/Lc(G);
3570	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3571	else if (F.isZero() && G.isZero()) return F.genOne();
3572	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
3573	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3574	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3575	if (F == G) return F/Lc(F);
3576
3577	CFMap M,N;
3578	int best_level= myCompress (A, B, M, N, topLevel);
3579
3580	if (best_level == 0) return B.genOne();
3581
3582	A= M(A);
3583	B= M(B);
3584
3585	Variable x= Variable (1);
3586
3587	//univariate case
3588	if (A.isUnivariate() && B.isUnivariate())
3589	return N (gcd (A, B));
3590
3591	CanonicalForm cA, cB; // content of A and B
3592	CanonicalForm ppA, ppB; // primitive part of A and B
3593	CanonicalForm gcdcAcB;
3594
3595	cA = uni_content (A);
3596	cB = uni_content (B);
3597	gcdcAcB= gcd (cA, cB);
3598	ppA= A/cA;
3599	ppB= B/cB;
3600
3601	CanonicalForm lcA, lcB; // leading coefficients of A and B
3602	CanonicalForm gcdlcAlcB;
3603	lcA= uni_lcoeff (ppA);
3604	lcB= uni_lcoeff (ppB);
3605
3606	if (fdivides (lcA, lcB))
3607	{
3608	if (fdivides (A, B))
3609	return F/Lc(F);
3610	}
3611	if (fdivides (lcB, lcA))
3612	{
3613	if (fdivides (B, A))
3614	return G/Lc(G);
3615	}
3616
3617	gcdlcAlcB= gcd (lcA, lcB);
3618
3619	int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3620	int d0;
3621
3622	if (d == 0)
3623	return N(gcdcAcB);
3624
3625	d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3626
3627	if (d0 < d)
3628	d= d0;
3629
3630	if (d == 0)
3631	return N(gcdcAcB);
3632
3633	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3634	CanonicalForm newtonPoly= 1;
3635	m= gcdlcAlcB;
3636	G_m= 0;
3637	H= 0;
3638	bool fail= false;
3639	topLevel= false;
3640	bool inextension= false;
3641	Variable V_buf, alpha, cleanUp;
3642	CanonicalForm prim_elem, im_prim_elem;
3643	CFList source, dest;
3644	do //first do
3645	{
3646	if (inextension)
3647	random_element= randomElement (m, V_buf, l, fail);
3648	else
3649	random_element= FpRandomElement (m, l, fail);
3650	if (random_element == 0 && !fail)
3651	{
3652	if (!find (l, random_element))
3653	l.append (random_element);
3654	continue;
3655	}
3656
3657	if (!fail && !inextension)
3658	{
3659	CFList list;
3660	TIMING_START (gcd_recursion);
3661	G_random_element=
3662	sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3663	list);
3664	TIMING_END_AND_PRINT (gcd_recursion,
3665	"time for recursive call: ");
3666	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3667	}
3668	else if (!fail && inextension)
3669	{
3670	CFList list;
3671	TIMING_START (gcd_recursion);
3672	G_random_element=
3673	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3674	list, topLevel);
3675	TIMING_END_AND_PRINT (gcd_recursion,
3676	"time for recursive call: ");
3677	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3678	}
3679	else if (fail && !inextension)
3680	{
3681	source= CFList();
3682	dest= CFList();
3683	CFList list;
3684	CanonicalForm mipo;
3685	int deg= 2;
3686	bool initialized= false;
3687	do
3688	{
3689	mipo= randomIrredpoly (deg, x);
3690	if (initialized)
3691	setMipo (alpha, mipo);
3692	else
3693	alpha= rootOf (mipo);
3694	inextension= true;
3695	fail= false;
3696	initialized= true;
3697	random_element= randomElement (m, alpha, l, fail);
3698	deg++;
3699	} while (fail);
3700	cleanUp= alpha;
3701	V_buf= alpha;
3702	list= CFList();
3703	TIMING_START (gcd_recursion);
3704	G_random_element=
3705	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3706	list, topLevel);
3707	TIMING_END_AND_PRINT (gcd_recursion,
3708	"time for recursive call: ");
3709	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3710	}
3711	else if (fail && inextension)
3712	{
3713	source= CFList();
3714	dest= CFList();
3715
3716	Variable V_buf3= V_buf;
3717	V_buf= chooseExtension (V_buf);
3718	bool prim_fail= false;
3719	Variable V_buf2;
3720	CanonicalForm prim_elem, im_prim_elem;
3721	prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3722
3723	if (V_buf3 != alpha)
3724	{
3725	m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3726	G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3727	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3728	dest);
3729	ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3730	ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3731	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3732	dest);
3733	for (CFListIterator i= l; i.hasItem(); i++)
3734	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3735	source, dest);
3736	}
3737
3738	ASSERT (!prim_fail, "failure in integer factorizer");
3739	if (prim_fail)
3740	; //ERROR
3741	else
3742	im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3743
3744	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3745	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3746
3747	for (CFListIterator i= l; i.hasItem(); i++)
3748	i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3749	im_prim_elem, source, dest);
3750	m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3751	G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3752	newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3753	source, dest);
3754	ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3755	ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3756	gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3757	source, dest);
3758	fail= false;
3759	random_element= randomElement (m, V_buf, l, fail );
3760	DEBOUTLN (cerr, "fail= " << fail);
3761	CFList list;
3762	TIMING_START (gcd_recursion);
3763	G_random_element=
3764	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3765	list, topLevel);
3766	TIMING_END_AND_PRINT (gcd_recursion,
3767	"time for recursive call: ");
3768	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3769	alpha= V_buf;
3770	}
3771
3772	if (!G_random_element.inCoeffDomain())
3773	d0= totaldegree (G_random_element, Variable(2),
3774	Variable (G_random_element.level()));
3775	else
3776	d0= 0;
3777
3778	if (d0 == 0)
3779	{
3780	if (inextension)
3781	prune (cleanUp);
3782	return N(gcdcAcB);
3783	}
3784	if (d0 > d)
3785	{
3786	if (!find (l, random_element))
3787	l.append (random_element);
3788	continue;
3789	}
3790
3791	G_random_element=
3792	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3793	* G_random_element;
3794
3795	skeleton= G_random_element;
3796
3797	if (!G_random_element.inCoeffDomain())
3798	d0= totaldegree (G_random_element, Variable(2),
3799	Variable (G_random_element.level()));
3800	else
3801	d0= 0;
3802
3803	if (d0 < d)
3804	{
3805	m= gcdlcAlcB;
3806	newtonPoly= 1;
3807	G_m= 0;
3808	d= d0;
3809	}
3810
3811	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3812
3813	if (uni_lcoeff (H) == gcdlcAlcB)
3814	{
3815	cH= uni_content (H);
3816	ppH= H/cH;
3817	ppH /= Lc (ppH);
3818	DEBOUTLN (cerr, "ppH= " << ppH);
3819
3820	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3821	{
3822	if (inextension)
3823	prune (cleanUp);
3824	return N(gcdcAcB*ppH);
3825	}
3826	}
3827	G_m= H;
3828	newtonPoly= newtonPoly*(x - random_element);
3829	m= m*(x - random_element);
3830	if (!find (l, random_element))
3831	l.append (random_element);
3832
3833	if ((getCharacteristic() > 3 && size (skeleton) < 200))
3834	{
3835	CFArray Monoms;
3836	CFArray* coeffMonoms;
3837
3838	do //second do
3839	{
3840	if (inextension)
3841	random_element= randomElement (m, alpha, l, fail);
3842	else
3843	random_element= FpRandomElement (m, l, fail);
3844	if (random_element == 0 && !fail)
3845	{
3846	if (!find (l, random_element))
3847	l.append (random_element);
3848	continue;
3849	}
3850
3851	bool sparseFail= false;
3852	if (!fail && !inextension)
3853	{
3854	CFList list;
3855	TIMING_START (gcd_recursion);
3856	if (LC (skeleton).inCoeffDomain())
3857	G_random_element=
3858	monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3859	skeleton, x, sparseFail, coeffMonoms,
3860	Monoms);
3861	else
3862	G_random_element=
3863	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3864	skeleton, x, sparseFail,
3865	coeffMonoms, Monoms);
3866	TIMING_END_AND_PRINT (gcd_recursion,
3867	"time for recursive call: ");
3868	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3869	}
3870	else if (!fail && inextension)
3871	{
3872	CFList list;
3873	TIMING_START (gcd_recursion);
3874	if (LC (skeleton).inCoeffDomain())
3875	G_random_element=
3876	monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3877	skeleton, alpha, sparseFail, coeffMonoms,
3878	Monoms);
3879	else
3880	G_random_element=
3881	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3882	skeleton, alpha, sparseFail, coeffMonoms,
3883	Monoms);
3884	TIMING_END_AND_PRINT (gcd_recursion,
3885	"time for recursive call: ");
3886	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3887	}
3888	else if (fail && !inextension)
3889	{
3890	source= CFList();
3891	dest= CFList();
3892	CFList list;
3893	CanonicalForm mipo;
3894	int deg= 2;
3895	bool initialized= false;
3896	do
3897	{
3898	mipo= randomIrredpoly (deg, x);
3899	if (initialized)
3900	setMipo (alpha, mipo);
3901	else
3902	alpha= rootOf (mipo);
3903	inextension= true;
3904	fail= false;
3905	initialized= true;
3906	random_element= randomElement (m, alpha, l, fail);
3907	deg++;
3908	} while (fail);
3909	cleanUp= alpha;
3910	V_buf= alpha;
3911	list= CFList();
3912	TIMING_START (gcd_recursion);
3913	if (LC (skeleton).inCoeffDomain())
3914	G_random_element=
3915	monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3916	skeleton, alpha, sparseFail, coeffMonoms,
3917	Monoms);
3918	else
3919	G_random_element=
3920	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3921	skeleton, alpha, sparseFail, coeffMonoms,
3922	Monoms);
3923	TIMING_END_AND_PRINT (gcd_recursion,
3924	"time for recursive call: ");
3925	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3926	}
3927	else if (fail && inextension)
3928	{
3929	source= CFList();
3930	dest= CFList();
3931
3932	Variable V_buf3= V_buf;
3933	V_buf= chooseExtension (V_buf);
3934	bool prim_fail= false;
3935	Variable V_buf2;
3936	CanonicalForm prim_elem, im_prim_elem;
3937	prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3938
3939	if (V_buf3 != alpha)
3940	{
3941	m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3942	G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3943	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3944	source, dest);
3945	ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3946	ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3947	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3948	source, dest);
3949	for (CFListIterator i= l; i.hasItem(); i++)
3950	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3951	source, dest);
3952	}
3953
3954	ASSERT (!prim_fail, "failure in integer factorizer");
3955	if (prim_fail)
3956	; //ERROR
3957	else
3958	im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3959
3960	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3961	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3962
3963	for (CFListIterator i= l; i.hasItem(); i++)
3964	i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3965	im_prim_elem, source, dest);
3966	m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3967	G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3968	newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3969	source, dest);
3970	ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3971	ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3972	gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3973	source, dest);
3974	fail= false;
3975	random_element= randomElement (m, V_buf, l, fail );
3976	DEBOUTLN (cerr, "fail= " << fail);
3977	CFList list;
3978	TIMING_START (gcd_recursion);
3979	if (LC (skeleton).inCoeffDomain())
3980	G_random_element=
3981	monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3982	skeleton, V_buf, sparseFail, coeffMonoms,
3983	Monoms);
3984	else
3985	G_random_element=
3986	nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3987	skeleton, V_buf, sparseFail,
3988	coeffMonoms, Monoms);
3989	TIMING_END_AND_PRINT (gcd_recursion,
3990	"time for recursive call: ");
3991	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3992	alpha= V_buf;
3993	}
3994
3995	if (sparseFail)
3996	break;
3997
3998	if (!G_random_element.inCoeffDomain())
3999	d0= totaldegree (G_random_element, Variable(2),
4000	Variable (G_random_element.level()));
4001	else
4002	d0= 0;
4003
4004	if (d0 == 0)
4005	{
4006	if (inextension)
4007	prune (cleanUp);
4008	return N(gcdcAcB);
4009	}
4010	if (d0 > d)
4011	{
4012	if (!find (l, random_element))
4013	l.append (random_element);
4014	continue;
4015	}
4016
4017	G_random_element=
4018	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
4019	* G_random_element;
4020
4021	if (!G_random_element.inCoeffDomain())
4022	d0= totaldegree (G_random_element, Variable(2),
4023	Variable (G_random_element.level()));
4024	else
4025	d0= 0;
4026
4027	if (d0 < d)
4028	{
4029	m= gcdlcAlcB;
4030	newtonPoly= 1;
4031	G_m= 0;
4032	d= d0;
4033	}
4034
4035	TIMING_START (newton_interpolation);
4036	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4037	TIMING_END_AND_PRINT (newton_interpolation,
4038	"time for newton interpolation: ");
4039
4040	//termination test
4041	if (uni_lcoeff (H) == gcdlcAlcB)
4042	{
4043	cH= uni_content (H);
4044	ppH= H/cH;
4045	ppH /= Lc (ppH);
4046	DEBOUTLN (cerr, "ppH= " << ppH);
4047	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4048	{
4049	if (inextension)
4050	prune (cleanUp);
4051	return N(gcdcAcB*ppH);
4052	}
4053	}
4054
4055	G_m= H;
4056	newtonPoly= newtonPoly*(x - random_element);
4057	m= m*(x - random_element);
4058	if (!find (l, random_element))
4059	l.append (random_element);
4060
4061	} while (1); //end of second do
4062	}
4063	else
4064	{
4065	if (inextension)
4066	prune (cleanUp);
4067	return N(gcdcAcB*modGCDFp (ppA, ppB));
4068	}
4069	} while (1); //end of first do
4070	}
4071
4072	TIMING_DEFINE_PRINT(modZ_termination)
4073	TIMING_DEFINE_PRINT(modZ_recursion)
4074
4075	/// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
4076	/// Algebra", Algorithm 7.1
4077	CanonicalForm modGCDZ ( const CanonicalForm & FF, const CanonicalForm & GG )
4078	{
4079	CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
4080	int p, i, dp_deg, d_deg=-1;
4081
4082	CanonicalForm cd ( bCommonDen( FF ));
4083	f=cd*FF;
4084	Variable x= Variable (1);
4085	CanonicalForm cf, cg;
4086	cf= icontent (f);
4087	f /= cf;
4088	//cd = bCommonDen( f ); f *=cd;
4089	//f /=vcontent(f,Variable(1));
4090
4091	cd = bCommonDen( GG );
4092	g=cd*GG;
4093	cg= icontent (g);
4094	g /= cg;
4095	//cd = bCommonDen( g ); g *=cd;
4096	//g /=vcontent(g,Variable(1));
4097
4098	CanonicalForm Dn, test= 0;
4099	CanonicalForm lcf, lcg;
4100	lcf= f.lc();
4101	lcg= g.lc();
4102	cl = gcd (f.lc(),g.lc());
4103	CanonicalForm gcdcfcg= gcd (cf, cg);
4104	CanonicalForm fp, gp;
4105	CanonicalForm b= 1;
4106	int minCommonDeg= 0;
4107	for (i= tmax (f.level(), g.level()); i > 0; i--)
4108	{
4109	if (degree (f, i) <= 0 \|\| degree (g, i) <= 0)
4110	continue;
4111	else
4112	{
4113	minCommonDeg= tmin (degree (g, i), degree (f, i));
4114	break;
4115	}
4116	}
4117	if (i == 0)
4118	return gcdcfcg;
4119	for (; i > 0; i--)
4120	{
4121	if (degree (f, i) <= 0 \|\| degree (g, i) <= 0)
4122	continue;
4123	else
4124	minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
4125	}
4126	b= 2tmin (maxNorm (f), maxNorm (g))abs (cl)*
4127	power (CanonicalForm (2), minCommonDeg);
4128	bool equal= false;
4129	i = cf_getNumBigPrimes() - 1;
4130
4131	CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn, cDn;
4132	int maxNumVars= tmax (getNumVars (f), getNumVars (g));
4133	//Off (SW_RATIONAL);
4134	while ( true )
4135	{
4136	p = cf_getBigPrime( i );
4137	i--;
4138	while ( i >= 0 && mod( cl(lc(f)/cl)(lc(g)/cl), p ) == 0 )
4139	{
4140	p = cf_getBigPrime( i );
4141	i--;
4142	}
4143	//printf("try p=%d\n",p);
4144	setCharacteristic( p );
4145	fp= mapinto (f);
4146	gp= mapinto (g);
4147	TIMING_START (modZ_recursion)
4148	#if defined(HAVE_NTL) \|\| defined(HAVE_FLINT)
4149	if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4150	Dp = modGCDFp (fp, gp, cofp, cogp);
4151	else
4152	{
4153	Dp= gcd_poly (fp, gp);
4154	cofp= fp/Dp;
4155	cogp= gp/Dp;
4156	}
4157	#else
4158	Dp= gcd_poly (fp, gp);
4159	cofp= fp/Dp;
4160	cogp= gp/Dp;
4161	#endif
4162	TIMING_END_AND_PRINT (modZ_recursion,
4163	"time for gcd mod p in modular gcd: ");
4164	Dp /=Dp.lc();
4165	Dp *= mapinto (cl);
4166	cofp /= lc (cofp);
4167	cofp *= mapinto (lcf);
4168	cogp /= lc (cogp);
4169	cogp *= mapinto (lcg);
4170	setCharacteristic( 0 );
4171	dp_deg=totaldegree(Dp);
4172	if ( dp_deg == 0 )
4173	{
4174	//printf(" -> 1\n");
4175	return CanonicalForm(gcdcfcg);
4176	}
4177	if ( q.isZero() )
4178	{
4179	D = mapinto( Dp );
4180	cof= mapinto (cofp);
4181	cog= mapinto (cogp);
4182	d_deg=dp_deg;
4183	q = p;
4184	Dn= balance_p (D, p);
4185	cofn= balance_p (cof, p);
4186	cogn= balance_p (cog, p);
4187	}
4188	else
4189	{
4190	if ( dp_deg == d_deg )
4191	{
4192	chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4193	chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4194	chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4195	cof= newCof;
4196	cog= newCog;
4197	newqh= newq/2;
4198	Dn= balance_p (newD, newq, newqh);
4199	cofn= balance_p (newCof, newq, newqh);
4200	cogn= balance_p (newCog, newq, newqh);
4201	if (test != Dn) //balance_p (newD, newq))
4202	test= balance_p (newD, newq);
4203	else
4204	equal= true;
4205	q = newq;
4206	D = newD;
4207	}
4208	else if ( dp_deg < d_deg )
4209	{
4210	//printf(" were all bad, try more\n");
4211	// all previous p's are bad primes
4212	q = p;
4213	D = mapinto( Dp );
4214	cof= mapinto (cof);
4215	cog= mapinto (cog);
4216	d_deg=dp_deg;
4217	test= 0;
4218	equal= false;
4219	Dn= balance_p (D, p);
4220	cofn= balance_p (cof, p);
4221	cogn= balance_p (cog, p);
4222	}
4223	else
4224	{
4225	//printf(" was bad, try more\n");
4226	}
4227	//else dp_deg > d_deg: bad prime
4228	}
4229	if ( i >= 0 )
4230	{
4231	cDn= icontent (Dn);
4232	Dn /= cDn;
4233	cofn /= cl/cDn;
4234	//cofn /= icontent (cofn);
4235	cogn /= cl/cDn;
4236	//cogn /= icontent (cogn);
4237	TIMING_START (modZ_termination);
4238	if ((terminationTest (f,g, cofn, cogn, Dn)) \|\|
4239	((equal \|\| q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4240	{
4241	TIMING_END_AND_PRINT (modZ_termination,
4242	"time for successful termination in modular gcd: ");
4243	//printf(" -> success\n");
4244	return Dn*gcdcfcg;
4245	}
4246	TIMING_END_AND_PRINT (modZ_termination,
4247	"time for unsuccessful termination in modular gcd: ");
4248	equal= false;
4249	//else: try more primes
4250	}
4251	else
4252	{ // try other method
4253	//printf("try other gcd\n");
4254	Off(SW_USE_CHINREM_GCD);
4255	D=gcd_poly( f, g );
4256	On(SW_USE_CHINREM_GCD);
4257	return D*gcdcfcg;
4258	}
4259	}
4260	}
4261	#endif

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

Download in other formats: