Context Navigation

source: git/factory/cfModGcd.cc @ b7cc2b

spielwiese

Last change on this file since b7cc2b was b7cc2b, checked in by Hans Schoenemann <hannes@…>, 3 years ago
factory: FindRoot via FLINT
Property mode set to `100644`
File size: 115.5 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	#ifdef HAVE_NTL
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	#ifdef HAVE_NTL // modGCDFq
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	#ifdef HAVE_NTL // mapPrimElem
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	#ifdef HAVE_NTL
1206	CanonicalForm
1207	modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1208	CanonicalForm& coF, CanonicalForm& coG,
1209	bool& topLevel, CFList& l);
1210	#endif
1211
1212	#ifdef HAVE_NTL
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	#ifdef HAVE_NTL // mapPrimElem
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	#ifdef HAVE_NTL //gauss, zz_pE
1786	long
1787	gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
1788	{
1789	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1790	CFMatrix *N;
1791	N= new CFMatrix (M.rows(), M.columns() + 1);
1792
1793	for (int i= 1; i <= M.rows(); i++)
1794	for (int j= 1; j <= M.columns(); j++)
1795	(*N) (i, j)= M (i, j);
1796
1797	int j= 1;
1798	for (int i= 0; i < L.size(); i++, j++)
1799	(*N) (j, M.columns() + 1)= L[i];
1800	int p= getCharacteristic ();
1801	if (fac_NTL_char != p)
1802	{
1803	fac_NTL_char= p;
1804	zz_p::init (p);
1805	}
1806	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1807	zz_pE::init (NTLMipo);
1808	mat_zz_pE NTLN= convertFacCFMatrix2NTLmat_zz_pE(N);
1809	long rk= gauss (*NTLN);
1810
1811	delete N;
1812	N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1813
1814	delete NTLN;
1815
1816	M= (*N) (1, M.rows(), 1, M.columns());
1817	L= CFArray (M.rows());
1818	for (int i= 0; i < M.rows(); i++)
1819	L[i]= (*N) (i + 1, M.columns() + 1);
1820
1821	delete N;
1822	return rk;
1823	}
1824	#endif
1825
1826	CFArray
1827	solveSystemFp (const CFMatrix& M, const CFArray& L)
1828	{
1829	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1830	CFMatrix *N;
1831	N= new CFMatrix (M.rows(), M.columns() + 1);
1832
1833	for (int i= 1; i <= M.rows(); i++)
1834	for (int j= 1; j <= M.columns(); j++)
1835	(*N) (i, j)= M (i, j);
1836
1837	int j= 1;
1838	for (int i= 0; i < L.size(); i++, j++)
1839	(*N) (j, M.columns() + 1)= L[i];
1840
1841	#ifdef HAVE_FLINT
1842	nmod_mat_t FLINTN;
1843	convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1844	long rk= nmod_mat_rref (FLINTN);
1845	#else
1846	int p= getCharacteristic ();
1847	if (fac_NTL_char != p)
1848	{
1849	fac_NTL_char= p;
1850	zz_p::init (p);
1851	}
1852	mat_zz_p NTLN= convertFacCFMatrix2NTLmat_zz_p(N);
1853	long rk= gauss (*NTLN);
1854	#endif
1855	delete N;
1856	if (rk != M.columns())
1857	{
1858	#ifdef HAVE_FLINT
1859	nmod_mat_clear (FLINTN);
1860	#else
1861	delete NTLN;
1862	#endif
1863	return CFArray();
1864	}
1865	#ifdef HAVE_FLINT
1866	N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1867	nmod_mat_clear (FLINTN);
1868	#else
1869	N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1870	delete NTLN;
1871	#endif
1872	CFArray A= readOffSolution (*N, rk);
1873
1874	delete N;
1875	return A;
1876	}
1877
1878	#ifdef HAVE_NTL //gauss, zz_pE
1879	CFArray
1880	solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1881	{
1882	ASSERT (L.size() <= M.rows(), "dimension exceeded");
1883	CFMatrix *N;
1884	N= new CFMatrix (M.rows(), M.columns() + 1);
1885
1886	for (int i= 1; i <= M.rows(); i++)
1887	for (int j= 1; j <= M.columns(); j++)
1888	(*N) (i, j)= M (i, j);
1889	int j= 1;
1890	for (int i= 0; i < L.size(); i++, j++)
1891	(*N) (j, M.columns() + 1)= L[i];
1892	int p= getCharacteristic ();
1893	if (fac_NTL_char != p)
1894	{
1895	fac_NTL_char= p;
1896	zz_p::init (p);
1897	}
1898	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1899	zz_pE::init (NTLMipo);
1900	mat_zz_pE NTLN= convertFacCFMatrix2NTLmat_zz_pE(N);
1901	long rk= gauss (*NTLN);
1902
1903	delete N;
1904	if (rk != M.columns())
1905	{
1906	delete NTLN;
1907	return CFArray();
1908	}
1909	N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1910
1911	delete NTLN;
1912
1913	CFArray A= readOffSolution (*N, rk);
1914
1915	delete N;
1916	return A;
1917	}
1918	#endif
1919	#endif
1920
1921	CFArray
1922	getMonoms (const CanonicalForm& F)
1923	{
1924	if (F.inCoeffDomain())
1925	{
1926	CFArray result= CFArray (1);
1927	result [0]= 1;
1928	return result;
1929	}
1930	if (F.isUnivariate())
1931	{
1932	CFArray result= CFArray (size(F));
1933	int j= 0;
1934	for (CFIterator i= F; i.hasTerms(); i++, j++)
1935	result[j]= power (F.mvar(), i.exp());
1936	return result;
1937	}
1938	int numMon= size (F);
1939	CFArray result= CFArray (numMon);
1940	int j= 0;
1941	CFArray recResult;
1942	Variable x= F.mvar();
1943	CanonicalForm powX;
1944	for (CFIterator i= F; i.hasTerms(); i++)
1945	{
1946	powX= power (x, i.exp());
1947	recResult= getMonoms (i.coeff());
1948	for (int k= 0; k < recResult.size(); k++)
1949	result[j+k]= powX*recResult[k];
1950	j += recResult.size();
1951	}
1952	return result;
1953	}
1954
1955	#ifdef HAVE_NTL
1956	CFArray
1957	evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
1958	{
1959	if (F.inCoeffDomain())
1960	{
1961	CFArray result= CFArray (1);
1962	result [0]= F;
1963	return result;
1964	}
1965	if (F.isUnivariate())
1966	{
1967	ASSERT (evalPoints.length() == 1,
1968	"expected an eval point with only one component");
1969	CFArray result= CFArray (size(F));
1970	int j= 0;
1971	CanonicalForm evalPoint= evalPoints.getLast();
1972	for (CFIterator i= F; i.hasTerms(); i++, j++)
1973	result[j]= power (evalPoint, i.exp());
1974	return result;
1975	}
1976	int numMon= size (F);
1977	CFArray result= CFArray (numMon);
1978	int j= 0;
1979	CanonicalForm evalPoint= evalPoints.getLast();
1980	CFList buf= evalPoints;
1981	buf.removeLast();
1982	CFArray recResult;
1983	CanonicalForm powEvalPoint;
1984	for (CFIterator i= F; i.hasTerms(); i++)
1985	{
1986	powEvalPoint= power (evalPoint, i.exp());
1987	recResult= evaluateMonom (i.coeff(), buf);
1988	for (int k= 0; k < recResult.size(); k++)
1989	result[j+k]= powEvalPoint*recResult[k];
1990	j += recResult.size();
1991	}
1992	return result;
1993	}
1994
1995	CFArray
1996	evaluate (const CFArray& A, const CFList& evalPoints)
1997	{
1998	CFArray result= A.size();
1999	CanonicalForm tmp;
2000	int k;
2001	for (int i= 0; i < A.size(); i++)
2002	{
2003	tmp= A[i];
2004	k= 1;
2005	for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2006	tmp= tmp (j.getItem(), k);
2007	result[i]= tmp;
2008	}
2009	return result;
2010	}
2011
2012	CFList
2013	evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
2014	CanonicalForm& Feval, CanonicalForm& Geval,
2015	const CanonicalForm& LCF, const bool& GF,
2016	const Variable& alpha, bool& fail, CFList& list
2017	)
2018	{
2019	int k= tmax (F.level(), G.level()) - 1;
2020	Variable x= Variable (1);
2021	CFList result;
2022	FFRandom genFF;
2023	GFRandom genGF;
2024	int p= getCharacteristic ();
2025	double bound;
2026	if (alpha != Variable (1))
2027	{
2028	bound= pow ((double) p, (double) degree (getMipo(alpha)));
2029	bound= pow (bound, (double) k);
2030	}
2031	else if (GF)
2032	{
2033	bound= pow ((double) p, (double) getGFDegree());
2034	bound= pow ((double) bound, (double) k);
2035	}
2036	else
2037	bound= pow ((double) p, (double) k);
2038
2039	CanonicalForm random;
2040	int j;
2041	bool zeroOneOccured= false;
2042	bool allEqual= false;
2043	CanonicalForm buf;
2044	do
2045	{
2046	random= 0;
2047	// possible overflow if list.length() does not fit into a int
2048	if (list.length() >= bound)
2049	{
2050	fail= true;
2051	break;
2052	}
2053	for (int i= 0; i < k; i++)
2054	{
2055	if (GF)
2056	{
2057	result.append (genGF.generate());
2058	random += result.getLast()*power (x, i);
2059	}
2060	else if (alpha.level() != 1)
2061	{
2062	AlgExtRandomF genAlgExt (alpha);
2063	result.append (genAlgExt.generate());
2064	random += result.getLast()*power (x, i);
2065	}
2066	else
2067	{
2068	result.append (genFF.generate());
2069	random += result.getLast()*power (x, i);
2070	}
2071	if (result.getLast().isOne() \|\| result.getLast().isZero())
2072	zeroOneOccured= true;
2073	}
2074	if (find (list, random))
2075	{
2076	zeroOneOccured= false;
2077	allEqual= false;
2078	result= CFList();
2079	continue;
2080	}
2081	if (zeroOneOccured)
2082	{
2083	list.append (random);
2084	zeroOneOccured= false;
2085	allEqual= false;
2086	result= CFList();
2087	continue;
2088	}
2089	// no zero at this point
2090	if (k > 1)
2091	{
2092	allEqual= true;
2093	CFIterator iter= random;
2094	buf= iter.coeff();
2095	iter++;
2096	for (; iter.hasTerms(); iter++)
2097	if (buf != iter.coeff())
2098	allEqual= false;
2099	}
2100	if (allEqual)
2101	{
2102	list.append (random);
2103	allEqual= false;
2104	zeroOneOccured= false;
2105	result= CFList();
2106	continue;
2107	}
2108
2109	Feval= F;
2110	Geval= G;
2111	CanonicalForm LCeval= LCF;
2112	j= 1;
2113	for (CFListIterator i= result; i.hasItem(); i++, j++)
2114	{
2115	Feval= Feval (i.getItem(), j);
2116	Geval= Geval (i.getItem(), j);
2117	LCeval= LCeval (i.getItem(), j);
2118	}
2119
2120	if (LCeval.isZero())
2121	{
2122	if (!find (list, random))
2123	list.append (random);
2124	zeroOneOccured= false;
2125	allEqual= false;
2126	result= CFList();
2127	continue;
2128	}
2129
2130	if (list.length() >= bound)
2131	{
2132	fail= true;
2133	break;
2134	}
2135	} while (find (list, random));
2136
2137	return result;
2138	}
2139
2140	/// multiply two lists componentwise
2141	void mult (CFList& L1, const CFList& L2)
2142	{
2143	ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2144
2145	CFListIterator j= L2;
2146	for (CFListIterator i= L1; i.hasItem(); i++, j++)
2147	i.getItem() *= j.getItem();
2148	}
2149
2150	void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2151	CanonicalForm& Beval, const CFList& L)
2152	{
2153	Aeval= A;
2154	Beval= B;
2155	int j= 1;
2156	for (CFListIterator i= L; i.hasItem(); i++, j++)
2157	{
2158	Aeval= Aeval (i.getItem(), j);
2159	Beval= Beval (i.getItem(), j);
2160	}
2161	}
2162
2163	CanonicalForm
2164	monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2165	const CanonicalForm& skeleton, const Variable& alpha,
2166	bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2167	)
2168	{
2169	CanonicalForm A= F;
2170	CanonicalForm B= G;
2171	if (F.isZero() && degree(G) > 0) return G/Lc(G);
2172	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2173	else if (F.isZero() && G.isZero()) return F.genOne();
2174	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
2175	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2176	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2177	if (F == G) return F/Lc(F);
2178
2179	ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2180	ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2181
2182	CFMap M,N;
2183	int best_level= myCompress (A, B, M, N, false);
2184
2185	if (best_level == 0)
2186	return B.genOne();
2187
2188	A= M(A);
2189	B= M(B);
2190
2191	Variable x= Variable (1);
2192
2193	//univariate case
2194	if (A.isUnivariate() && B.isUnivariate())
2195	return N (gcd (A, B));
2196
2197	CanonicalForm skel= M(skeleton);
2198	CanonicalForm cA, cB; // content of A and B
2199	CanonicalForm ppA, ppB; // primitive part of A and B
2200	CanonicalForm gcdcAcB;
2201	cA = uni_content (A);
2202	cB = uni_content (B);
2203	gcdcAcB= gcd (cA, cB);
2204	ppA= A/cA;
2205	ppB= B/cB;
2206
2207	CanonicalForm lcA, lcB; // leading coefficients of A and B
2208	CanonicalForm gcdlcAlcB;
2209	lcA= uni_lcoeff (ppA);
2210	lcB= uni_lcoeff (ppB);
2211
2212	if (fdivides (lcA, lcB))
2213	{
2214	if (fdivides (A, B))
2215	return F/Lc(F);
2216	}
2217	if (fdivides (lcB, lcA))
2218	{
2219	if (fdivides (B, A))
2220	return G/Lc(G);
2221	}
2222
2223	gcdlcAlcB= gcd (lcA, lcB);
2224	int skelSize= size (skel, skel.mvar());
2225
2226	int j= 0;
2227	int biggestSize= 0;
2228
2229	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2230	biggestSize= tmax (biggestSize, size (i.coeff()));
2231
2232	CanonicalForm g, Aeval, Beval;
2233
2234	CFList evalPoints;
2235	bool evalFail= false;
2236	CFList list;
2237	bool GF= false;
2238	CanonicalForm LCA= LC (A);
2239	CanonicalForm tmp;
2240	CFArray gcds= CFArray (biggestSize);
2241	CFList * pEvalPoints= new CFList [biggestSize];
2242	Variable V_buf= alpha, V_buf4= alpha;
2243	CFList source, dest;
2244	CanonicalForm prim_elem, im_prim_elem;
2245	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2246	for (int i= 0; i < biggestSize; i++)
2247	{
2248	if (i == 0)
2249	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2250	list);
2251	else
2252	{
2253	mult (evalPoints, pEvalPoints [0]);
2254	eval (A, B, Aeval, Beval, evalPoints);
2255	}
2256
2257	if (evalFail)
2258	{
2259	if (V_buf.level() != 1)
2260	{
2261	do
2262	{
2263	Variable V_buf3= V_buf;
2264	V_buf= chooseExtension (V_buf);
2265	source= CFList();
2266	dest= CFList();
2267
2268	bool prim_fail= false;
2269	Variable V_buf2;
2270	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2271	if (V_buf4 == alpha && alpha.level() != 1)
2272	prim_elem_alpha= prim_elem;
2273
2274	ASSERT (!prim_fail, "failure in integer factorizer");
2275	if (prim_fail)
2276	; //ERROR
2277	else
2278	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2279
2280	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2281	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2282
2283	if (V_buf4 == alpha && alpha.level() != 1)
2284	im_prim_elem_alpha= im_prim_elem;
2285	else if (alpha.level() != 1)
2286	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2287	prim_elem, im_prim_elem, source, dest);
2288
2289	for (CFListIterator j= list; j.hasItem(); j++)
2290	j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2291	im_prim_elem, source, dest);
2292	for (int k= 0; k < i; k++)
2293	{
2294	for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2295	j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2296	im_prim_elem, source, dest);
2297	gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2298	source, dest);
2299	}
2300
2301	if (alpha.level() != 1)
2302	{
2303	A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2304	B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2305	}
2306	V_buf4= V_buf;
2307	evalFail= false;
2308	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2309	evalFail, list);
2310	} while (evalFail);
2311	}
2312	else
2313	{
2314	CanonicalForm mipo;
2315	int deg= 2;
2316	bool initialized= false;
2317	do
2318	{
2319	mipo= randomIrredpoly (deg, x);
2320	if (initialized)
2321	setMipo (V_buf, mipo);
2322	else
2323	V_buf= rootOf (mipo);
2324	evalFail= false;
2325	initialized= true;
2326	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2327	evalFail, list);
2328	deg++;
2329	} while (evalFail);
2330	V_buf4= V_buf;
2331	}
2332	}
2333
2334	g= gcd (Aeval, Beval);
2335	g /= Lc (g);
2336
2337	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2338	{
2339	delete[] pEvalPoints;
2340	fail= true;
2341	if (alpha.level() != 1 && V_buf != alpha)
2342	prune1 (alpha);
2343	return 0;
2344	}
2345	CFIterator l= skel;
2346	for (CFIterator k= g; k.hasTerms(); k++, l++)
2347	{
2348	if (k.exp() != l.exp())
2349	{
2350	delete[] pEvalPoints;
2351	fail= true;
2352	if (alpha.level() != 1 && V_buf != alpha)
2353	prune1 (alpha);
2354	return 0;
2355	}
2356	}
2357	pEvalPoints[i]= evalPoints;
2358	gcds[i]= g;
2359
2360	tmp= 0;
2361	int j= 0;
2362	for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2363	tmp += k.getItem()*power (x, j);
2364	list.append (tmp);
2365	}
2366
2367	if (Monoms.size() == 0)
2368	Monoms= getMonoms (skel);
2369
2370	coeffMonoms= new CFArray [skelSize];
2371	j= 0;
2372	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2373	coeffMonoms[j]= getMonoms (i.coeff());
2374
2375	CFArray* pL= new CFArray [skelSize];
2376	CFArray* pM= new CFArray [skelSize];
2377	for (int i= 0; i < biggestSize; i++)
2378	{
2379	CFIterator l= gcds [i];
2380	evalPoints= pEvalPoints [i];
2381	for (int k= 0; k < skelSize; k++, l++)
2382	{
2383	if (i == 0)
2384	pL[k]= CFArray (biggestSize);
2385	pL[k] [i]= l.coeff();
2386
2387	if (i == 0)
2388	pM[k]= evaluate (coeffMonoms [k], evalPoints);
2389	}
2390	}
2391
2392	CFArray solution;
2393	CanonicalForm result= 0;
2394	int ind= 0;
2395	CFArray bufArray;
2396	CFMatrix Mat;
2397	for (int k= 0; k < skelSize; k++)
2398	{
2399	if (biggestSize != coeffMonoms[k].size())
2400	{
2401	bufArray= CFArray (coeffMonoms[k].size());
2402	for (int i= 0; i < coeffMonoms[k].size(); i++)
2403	bufArray [i]= pL[k] [i];
2404	solution= solveGeneralVandermonde (pM [k], bufArray);
2405	}
2406	else
2407	solution= solveGeneralVandermonde (pM [k], pL[k]);
2408
2409	if (solution.size() == 0)
2410	{
2411	delete[] pEvalPoints;
2412	delete[] pM;
2413	delete[] pL;
2414	delete[] coeffMonoms;
2415	fail= true;
2416	if (alpha.level() != 1 && V_buf != alpha)
2417	prune1 (alpha);
2418	return 0;
2419	}
2420	for (int l= 0; l < solution.size(); l++)
2421	result += solution[l]*Monoms [ind + l];
2422	ind += solution.size();
2423	}
2424
2425	delete[] pEvalPoints;
2426	delete[] pM;
2427	delete[] pL;
2428	delete[] coeffMonoms;
2429
2430	if (alpha.level() != 1 && V_buf != alpha)
2431	{
2432	CFList u, v;
2433	result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2434	prune1 (alpha);
2435	}
2436
2437	result= N(result);
2438	if (fdivides (result, F) && fdivides (result, G))
2439	return result;
2440	else
2441	{
2442	fail= true;
2443	return 0;
2444	}
2445	}
2446
2447	CanonicalForm
2448	nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2449	const CanonicalForm& skeleton, const Variable& alpha,
2450	bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2451	)
2452	{
2453	CanonicalForm A= F;
2454	CanonicalForm B= G;
2455	if (F.isZero() && degree(G) > 0) return G/Lc(G);
2456	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2457	else if (F.isZero() && G.isZero()) return F.genOne();
2458	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
2459	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2460	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2461	if (F == G) return F/Lc(F);
2462
2463	ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2464	ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2465
2466	CFMap M,N;
2467	int best_level= myCompress (A, B, M, N, false);
2468
2469	if (best_level == 0)
2470	return B.genOne();
2471
2472	A= M(A);
2473	B= M(B);
2474
2475	Variable x= Variable (1);
2476
2477	//univariate case
2478	if (A.isUnivariate() && B.isUnivariate())
2479	return N (gcd (A, B));
2480
2481	CanonicalForm skel= M(skeleton);
2482
2483	CanonicalForm cA, cB; // content of A and B
2484	CanonicalForm ppA, ppB; // primitive part of A and B
2485	CanonicalForm gcdcAcB;
2486	cA = uni_content (A);
2487	cB = uni_content (B);
2488	gcdcAcB= gcd (cA, cB);
2489	ppA= A/cA;
2490	ppB= B/cB;
2491
2492	CanonicalForm lcA, lcB; // leading coefficients of A and B
2493	CanonicalForm gcdlcAlcB;
2494	lcA= uni_lcoeff (ppA);
2495	lcB= uni_lcoeff (ppB);
2496
2497	if (fdivides (lcA, lcB))
2498	{
2499	if (fdivides (A, B))
2500	return F/Lc(F);
2501	}
2502	if (fdivides (lcB, lcA))
2503	{
2504	if (fdivides (B, A))
2505	return G/Lc(G);
2506	}
2507
2508	gcdlcAlcB= gcd (lcA, lcB);
2509	int skelSize= size (skel, skel.mvar());
2510
2511	int j= 0;
2512	int biggestSize= 0;
2513	int bufSize;
2514	int numberUni= 0;
2515	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2516	{
2517	bufSize= size (i.coeff());
2518	biggestSize= tmax (biggestSize, bufSize);
2519	numberUni += bufSize;
2520	}
2521	numberUni--;
2522	numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2523	biggestSize= tmax (biggestSize , numberUni);
2524
2525	numberUni= biggestSize + size (LC(skel)) - 1;
2526	int biggestSize2= tmax (numberUni, biggestSize);
2527
2528	CanonicalForm g, Aeval, Beval;
2529
2530	CFList evalPoints;
2531	CFArray coeffEval;
2532	bool evalFail= false;
2533	CFList list;
2534	bool GF= false;
2535	CanonicalForm LCA= LC (A);
2536	CanonicalForm tmp;
2537	CFArray gcds= CFArray (biggestSize);
2538	CFList * pEvalPoints= new CFList [biggestSize];
2539	Variable V_buf= alpha, V_buf4= alpha;
2540	CFList source, dest;
2541	CanonicalForm prim_elem, im_prim_elem;
2542	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2543	for (int i= 0; i < biggestSize; i++)
2544	{
2545	if (i == 0)
2546	{
2547	if (getCharacteristic() > 3)
2548	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2549	evalFail, list);
2550	else
2551	evalFail= true;
2552
2553	if (evalFail)
2554	{
2555	if (V_buf.level() != 1)
2556	{
2557	do
2558	{
2559	Variable V_buf3= V_buf;
2560	V_buf= chooseExtension (V_buf);
2561	source= CFList();
2562	dest= CFList();
2563
2564	bool prim_fail= false;
2565	Variable V_buf2;
2566	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2567	if (V_buf4 == alpha && alpha.level() != 1)
2568	prim_elem_alpha= prim_elem;
2569
2570	ASSERT (!prim_fail, "failure in integer factorizer");
2571	if (prim_fail)
2572	; //ERROR
2573	else
2574	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2575
2576	DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2577	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2578
2579	if (V_buf4 == alpha && alpha.level() != 1)
2580	im_prim_elem_alpha= im_prim_elem;
2581	else if (alpha.level() != 1)
2582	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2583	prim_elem, im_prim_elem, source, dest);
2584
2585	for (CFListIterator i= list; i.hasItem(); i++)
2586	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2587	im_prim_elem, source, dest);
2588	if (alpha.level() != 1)
2589	{
2590	A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2591	B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2592	}
2593	evalFail= false;
2594	V_buf4= V_buf;
2595	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2596	evalFail, list);
2597	} while (evalFail);
2598	}
2599	else
2600	{
2601	CanonicalForm mipo;
2602	int deg= 2;
2603	bool initialized= false;
2604	do
2605	{
2606	mipo= randomIrredpoly (deg, x);
2607	if (initialized)
2608	setMipo (V_buf, mipo);
2609	else
2610	V_buf= rootOf (mipo);
2611	evalFail= false;
2612	initialized= true;
2613	evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2614	evalFail, list);
2615	deg++;
2616	} while (evalFail);
2617	V_buf4= V_buf;
2618	}
2619	}
2620	}
2621	else
2622	{
2623	mult (evalPoints, pEvalPoints[0]);
2624	eval (A, B, Aeval, Beval, evalPoints);
2625	}
2626
2627	g= gcd (Aeval, Beval);
2628	g /= Lc (g);
2629
2630	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2631	{
2632	delete[] pEvalPoints;
2633	fail= true;
2634	if (alpha.level() != 1 && V_buf != alpha)
2635	prune1 (alpha);
2636	return 0;
2637	}
2638	CFIterator l= skel;
2639	for (CFIterator k= g; k.hasTerms(); k++, l++)
2640	{
2641	if (k.exp() != l.exp())
2642	{
2643	delete[] pEvalPoints;
2644	fail= true;
2645	if (alpha.level() != 1 && V_buf != alpha)
2646	prune1 (alpha);
2647	return 0;
2648	}
2649	}
2650	pEvalPoints[i]= evalPoints;
2651	gcds[i]= g;
2652
2653	tmp= 0;
2654	int j= 0;
2655	for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2656	tmp += k.getItem()*power (x, j);
2657	list.append (tmp);
2658	}
2659
2660	if (Monoms.size() == 0)
2661	Monoms= getMonoms (skel);
2662
2663	coeffMonoms= new CFArray [skelSize];
2664
2665	j= 0;
2666	for (CFIterator i= skel; i.hasTerms(); i++, j++)
2667	coeffMonoms[j]= getMonoms (i.coeff());
2668
2669	int minimalColumnsIndex;
2670	if (skelSize > 1)
2671	minimalColumnsIndex= 1;
2672	else
2673	minimalColumnsIndex= 0;
2674	int minimalColumns=-1;
2675
2676	CFArray* pM= new CFArray [skelSize];
2677	CFMatrix Mat;
2678	// find the Matrix with minimal number of columns
2679	for (int i= 0; i < skelSize; i++)
2680	{
2681	pM[i]= CFArray (coeffMonoms[i].size());
2682	if (i == 1)
2683	minimalColumns= coeffMonoms[i].size();
2684	if (i > 1)
2685	{
2686	if (minimalColumns > coeffMonoms[i].size())
2687	{
2688	minimalColumns= coeffMonoms[i].size();
2689	minimalColumnsIndex= i;
2690	}
2691	}
2692	}
2693	CFMatrix* pMat= new CFMatrix [2];
2694	pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2695	pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2696	CFArray* pL= new CFArray [skelSize];
2697	for (int i= 0; i < biggestSize; i++)
2698	{
2699	CFIterator l= gcds [i];
2700	evalPoints= pEvalPoints [i];
2701	for (int k= 0; k < skelSize; k++, l++)
2702	{
2703	if (i == 0)
2704	pL[k]= CFArray (biggestSize);
2705	pL[k] [i]= l.coeff();
2706
2707	if (i == 0 && k != 0 && k != minimalColumnsIndex)
2708	pM[k]= evaluate (coeffMonoms[k], evalPoints);
2709	else if (i == 0 && (k == 0 \|\| k == minimalColumnsIndex))
2710	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2711	if (i > 0 && (k == 0 \|\| k == minimalColumnsIndex))
2712	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2713
2714	if (k == 0)
2715	{
2716	if (pMat[k].rows() >= i + 1)
2717	{
2718	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2719	pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2720	}
2721	}
2722	if (k == minimalColumnsIndex)
2723	{
2724	if (pMat[1].rows() >= i + 1)
2725	{
2726	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2727	pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2728	}
2729	}
2730	}
2731	}
2732
2733	CFArray solution;
2734	CanonicalForm result= 0;
2735	int ind= 1;
2736	int matRows, matColumns;
2737	matRows= pMat[1].rows();
2738	matColumns= pMat[0].columns() - 1;
2739	matColumns += pMat[1].columns();
2740
2741	Mat= CFMatrix (matRows, matColumns);
2742	for (int i= 1; i <= matRows; i++)
2743	for (int j= 1; j <= pMat[1].columns(); j++)
2744	Mat (i, j)= pMat[1] (i, j);
2745
2746	for (int j= pMat[1].columns() + 1; j <= matColumns;
2747	j++, ind++)
2748	{
2749	for (int i= 1; i <= matRows; i++)
2750	Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2751	}
2752
2753	CFArray firstColumn= CFArray (pMat[0].rows());
2754	for (int i= 0; i < pMat[0].rows(); i++)
2755	firstColumn [i]= pMat[0] (i + 1, 1);
2756
2757	CFArray bufArray= pL[minimalColumnsIndex];
2758
2759	for (int i= 0; i < biggestSize; i++)
2760	pL[minimalColumnsIndex] [i] *= firstColumn[i];
2761
2762	CFMatrix bufMat= pMat[1];
2763	pMat[1]= Mat;
2764
2765	if (V_buf.level() != 1)
2766	solution= solveSystemFq (pMat[1],
2767	pL[minimalColumnsIndex], V_buf);
2768	else
2769	solution= solveSystemFp (pMat[1],
2770	pL[minimalColumnsIndex]);
2771
2772	if (solution.size() == 0)
2773	{ //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2774	CFMatrix bufMat0= pMat[0];
2775	delete [] pMat;
2776	pMat= new CFMatrix [skelSize];
2777	pL[minimalColumnsIndex]= bufArray;
2778	CFList* bufpEvalPoints= NULL;
2779	CFArray bufGcds;
2780	if (biggestSize != biggestSize2)
2781	{
2782	bufpEvalPoints= pEvalPoints;
2783	pEvalPoints= new CFList [biggestSize2];
2784	bufGcds= gcds;
2785	gcds= CFArray (biggestSize2);
2786	for (int i= 0; i < biggestSize; i++)
2787	{
2788	pEvalPoints[i]= bufpEvalPoints [i];
2789	gcds[i]= bufGcds[i];
2790	}
2791	for (int i= 0; i < biggestSize2 - biggestSize; i++)
2792	{
2793	mult (evalPoints, pEvalPoints[0]);
2794	eval (A, B, Aeval, Beval, evalPoints);
2795	g= gcd (Aeval, Beval);
2796	g /= Lc (g);
2797	if (degree (g) != degree (skel) \|\| (skelSize != size (g)))
2798	{
2799	delete[] pEvalPoints;
2800	delete[] pMat;
2801	delete[] pL;
2802	delete[] coeffMonoms;
2803	delete[] pM;
2804	if (bufpEvalPoints != NULL)
2805	delete [] bufpEvalPoints;
2806	fail= true;
2807	if (alpha.level() != 1 && V_buf != alpha)
2808	prune1 (alpha);
2809	return 0;
2810	}
2811	CFIterator l= skel;
2812	for (CFIterator k= g; k.hasTerms(); k++, l++)
2813	{
2814	if (k.exp() != l.exp())
2815	{
2816	delete[] pEvalPoints;
2817	delete[] pMat;
2818	delete[] pL;
2819	delete[] coeffMonoms;
2820	delete[] pM;
2821	if (bufpEvalPoints != NULL)
2822	delete [] bufpEvalPoints;
2823	fail= true;
2824	if (alpha.level() != 1 && V_buf != alpha)
2825	prune1 (alpha);
2826	return 0;
2827	}
2828	}
2829	pEvalPoints[i + biggestSize]= evalPoints;
2830	gcds[i + biggestSize]= g;
2831	}
2832	}
2833	for (int i= 0; i < biggestSize; i++)
2834	{
2835	CFIterator l= gcds [i];
2836	evalPoints= pEvalPoints [i];
2837	for (int k= 1; k < skelSize; k++, l++)
2838	{
2839	if (i == 0)
2840	pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2841	if (k == minimalColumnsIndex)
2842	continue;
2843	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2844	if (pMat[k].rows() >= i + 1)
2845	{
2846	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2847	pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2848	}
2849	}
2850	}
2851	Mat= bufMat0;
2852	pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2853	for (int j= 1; j <= Mat.rows(); j++)
2854	for (int k= 1; k <= Mat.columns(); k++)
2855	pMat [0] (j,k)= Mat (j,k);
2856	Mat= bufMat;
2857	for (int j= 1; j <= Mat.rows(); j++)
2858	for (int k= 1; k <= Mat.columns(); k++)
2859	pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2860	// write old matrix entries into new matrices
2861	for (int i= 0; i < skelSize; i++)
2862	{
2863	bufArray= pL[i];
2864	pL[i]= CFArray (biggestSize2);
2865	for (int j= 0; j < bufArray.size(); j++)
2866	pL[i] [j]= bufArray [j];
2867	}
2868	//write old vector entries into new and add new entries to old matrices
2869	for (int i= 0; i < biggestSize2 - biggestSize; i++)
2870	{
2871	CFIterator l= gcds [i + biggestSize];
2872	evalPoints= pEvalPoints [i + biggestSize];
2873	for (int k= 0; k < skelSize; k++, l++)
2874	{
2875	pL[k] [i + biggestSize]= l.coeff();
2876	coeffEval= evaluate (coeffMonoms[k], evalPoints);
2877	if (pMat[k].rows() >= i + biggestSize + 1)
2878	{
2879	for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2880	pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2881	}
2882	}
2883	}
2884	// begin new
2885	for (int i= 0; i < skelSize; i++)
2886	{
2887	if (pL[i].size() > 1)
2888	{
2889	for (int j= 2; j <= pMat[i].rows(); j++)
2890	pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2891	-pL[i] [j - 1];
2892	}
2893	}
2894
2895	matColumns= biggestSize2 - 1;
2896	matRows= 0;
2897	for (int i= 0; i < skelSize; i++)
2898	{
2899	if (V_buf.level() == 1)
2900	(void) gaussianElimFp (pMat[i], pL[i]);
2901	else
2902	(void) gaussianElimFq (pMat[i], pL[i], V_buf);
2903
2904	if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2905	{
2906	delete[] pEvalPoints;
2907	delete[] pMat;
2908	delete[] pL;
2909	delete[] coeffMonoms;
2910	delete[] pM;
2911	if (bufpEvalPoints != NULL)
2912	delete [] bufpEvalPoints;
2913	fail= true;
2914	if (alpha.level() != 1 && V_buf != alpha)
2915	prune1 (alpha);
2916	return 0;
2917	}
2918	matRows += pMat[i].rows() - coeffMonoms[i].size();
2919	}
2920
2921	CFMatrix bufMat;
2922	Mat= CFMatrix (matRows, matColumns);
2923	ind= 0;
2924	bufArray= CFArray (matRows);
2925	CFArray bufArray2;
2926	for (int i= 0; i < skelSize; i++)
2927	{
2928	if (coeffMonoms[i].size() + 1 >= pMat[i].rows() \|\| coeffMonoms[i].size() + 1 >= pMat[i].columns())
2929	{
2930	delete[] pEvalPoints;
2931	delete[] pMat;
2932	delete[] pL;
2933	delete[] coeffMonoms;
2934	delete[] pM;
2935	if (bufpEvalPoints != NULL)
2936	delete [] bufpEvalPoints;
2937	fail= true;
2938	if (alpha.level() != 1 && V_buf != alpha)
2939	prune1 (alpha);
2940	return 0;
2941	}
2942	bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2943	coeffMonoms[i].size() + 1, pMat[i].columns());
2944
2945	for (int j= 1; j <= bufMat.rows(); j++)
2946	for (int k= 1; k <= bufMat.columns(); k++)
2947	Mat (j + ind, k)= bufMat(j, k);
2948	bufArray2= coeffMonoms[i].size();
2949	for (int j= 1; j <= pMat[i].rows(); j++)
2950	{
2951	if (j > coeffMonoms[i].size())
2952	bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2953	else
2954	bufArray2 [j - 1]= pL[i] [j - 1];
2955	}
2956	pL[i]= bufArray2;
2957	ind += bufMat.rows();
2958	pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2959	}
2960
2961	if (V_buf.level() != 1)
2962	solution= solveSystemFq (Mat, bufArray, V_buf);
2963	else
2964	solution= solveSystemFp (Mat, bufArray);
2965
2966	if (solution.size() == 0)
2967	{
2968	delete[] pEvalPoints;
2969	delete[] pMat;
2970	delete[] pL;
2971	delete[] coeffMonoms;
2972	delete[] pM;
2973	if (bufpEvalPoints != NULL)
2974	delete [] bufpEvalPoints;
2975	fail= true;
2976	if (alpha.level() != 1 && V_buf != alpha)
2977	prune1 (alpha);
2978	return 0;
2979	}
2980
2981	ind= 0;
2982	result= 0;
2983	CFArray bufSolution;
2984	for (int i= 0; i < skelSize; i++)
2985	{
2986	bufSolution= readOffSolution (pMat[i], pL[i], solution);
2987	for (int i= 0; i < bufSolution.size(); i++, ind++)
2988	result += Monoms [ind]*bufSolution[i];
2989	}
2990	if (alpha.level() != 1 && V_buf != alpha)
2991	{
2992	CFList u, v;
2993	result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2994	prune1 (alpha);
2995	}
2996	result= N(result);
2997	delete[] pEvalPoints;
2998	delete[] pMat;
2999	delete[] pL;
3000	delete[] coeffMonoms;
3001	delete[] pM;
3002
3003	if (bufpEvalPoints != NULL)
3004	delete [] bufpEvalPoints;
3005	if (fdivides (result, F) && fdivides (result, G))
3006	return result;
3007	else
3008	{
3009	fail= true;
3010	return 0;
3011	}
3012	} // end of deKleine, Monagan & Wittkopf
3013
3014	result += Monoms[0];
3015	int ind2= 0, ind3= 2;
3016	ind= 0;
3017	for (int l= 0; l < minimalColumnsIndex; l++)
3018	ind += coeffMonoms[l].size();
3019	for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3020	l++, ind2++, ind3++)
3021	{
3022	result += solution[l]*Monoms [1 + ind2];
3023	for (int i= 0; i < pMat[0].rows(); i++)
3024	firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3025	}
3026	for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3027	result += solution[l]*Monoms [ind + l];
3028	ind= coeffMonoms[0].size();
3029	for (int k= 1; k < skelSize; k++)
3030	{
3031	if (k == minimalColumnsIndex)
3032	{
3033	ind += coeffMonoms[k].size();
3034	continue;
3035	}
3036	if (k != minimalColumnsIndex)
3037	{
3038	for (int i= 0; i < biggestSize; i++)
3039	pL[k] [i] *= firstColumn [i];
3040	}
3041
3042	if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3043	{
3044	bufArray= CFArray (coeffMonoms[k].size());
3045	for (int i= 0; i < bufArray.size(); i++)
3046	bufArray [i]= pL[k] [i];
3047	solution= solveGeneralVandermonde (pM[k], bufArray);
3048	}
3049	else
3050	solution= solveGeneralVandermonde (pM[k], pL[k]);
3051
3052	if (solution.size() == 0)
3053	{
3054	delete[] pEvalPoints;
3055	delete[] pMat;
3056	delete[] pL;
3057	delete[] coeffMonoms;
3058	delete[] pM;
3059	fail= true;
3060	if (alpha.level() != 1 && V_buf != alpha)
3061	prune1 (alpha);
3062	return 0;
3063	}
3064	if (k != minimalColumnsIndex)
3065	{
3066	for (int l= 0; l < solution.size(); l++)
3067	result += solution[l]*Monoms [ind + l];
3068	ind += solution.size();
3069	}
3070	}
3071
3072	delete[] pEvalPoints;
3073	delete[] pMat;
3074	delete[] pL;
3075	delete[] pM;
3076	delete[] coeffMonoms;
3077
3078	if (alpha.level() != 1 && V_buf != alpha)
3079	{
3080	CFList u, v;
3081	result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3082	prune1 (alpha);
3083	}
3084	result= N(result);
3085
3086	if (fdivides (result, F) && fdivides (result, G))
3087	return result;
3088	else
3089	{
3090	fail= true;
3091	return 0;
3092	}
3093	}
3094
3095	CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
3096	const Variable & alpha, CFList& l, bool& topLevel)
3097	{
3098	CanonicalForm A= F;
3099	CanonicalForm B= G;
3100	if (F.isZero() && degree(G) > 0) return G/Lc(G);
3101	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3102	else if (F.isZero() && G.isZero()) return F.genOne();
3103	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
3104	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3105	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3106	if (F == G) return F/Lc(F);
3107
3108	CFMap M,N;
3109	int best_level= myCompress (A, B, M, N, topLevel);
3110
3111	if (best_level == 0) return B.genOne();
3112
3113	A= M(A);
3114	B= M(B);
3115
3116	Variable x= Variable (1);
3117
3118	//univariate case
3119	if (A.isUnivariate() && B.isUnivariate())
3120	return N (gcd (A, B));
3121
3122	CanonicalForm cA, cB; // content of A and B
3123	CanonicalForm ppA, ppB; // primitive part of A and B
3124	CanonicalForm gcdcAcB;
3125
3126	cA = uni_content (A);
3127	cB = uni_content (B);
3128	gcdcAcB= gcd (cA, cB);
3129	ppA= A/cA;
3130	ppB= B/cB;
3131
3132	CanonicalForm lcA, lcB; // leading coefficients of A and B
3133	CanonicalForm gcdlcAlcB;
3134	lcA= uni_lcoeff (ppA);
3135	lcB= uni_lcoeff (ppB);
3136
3137	if (fdivides (lcA, lcB))
3138	{
3139	if (fdivides (A, B))
3140	return F/Lc(F);
3141	}
3142	if (fdivides (lcB, lcA))
3143	{
3144	if (fdivides (B, A))
3145	return G/Lc(G);
3146	}
3147
3148	gcdlcAlcB= gcd (lcA, lcB);
3149
3150	int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3151	int d0;
3152
3153	if (d == 0)
3154	return N(gcdcAcB);
3155	d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3156
3157	if (d0 < d)
3158	d= d0;
3159
3160	if (d == 0)
3161	return N(gcdcAcB);
3162
3163	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3164	CanonicalForm newtonPoly= 1;
3165	m= gcdlcAlcB;
3166	G_m= 0;
3167	H= 0;
3168	bool fail= false;
3169	topLevel= false;
3170	bool inextension= false;
3171	Variable V_buf= alpha, V_buf4= alpha;
3172	CanonicalForm prim_elem, im_prim_elem;
3173	CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3174	CFList source, dest;
3175	do // first do
3176	{
3177	random_element= randomElement (m, V_buf, l, fail);
3178	if (random_element == 0 && !fail)
3179	{
3180	if (!find (l, random_element))
3181	l.append (random_element);
3182	continue;
3183	}
3184	if (fail)
3185	{
3186	source= CFList();
3187	dest= CFList();
3188
3189	Variable V_buf3= V_buf;
3190	V_buf= chooseExtension (V_buf);
3191	bool prim_fail= false;
3192	Variable V_buf2;
3193	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3194	if (V_buf4 == alpha)
3195	prim_elem_alpha= prim_elem;
3196
3197	if (V_buf3 != V_buf4)
3198	{
3199	m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3200	G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3201	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3202	source, dest);
3203	ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3204	ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3205	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3206	dest);
3207	for (CFListIterator i= l; i.hasItem(); i++)
3208	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3209	source, dest);
3210	}
3211
3212	ASSERT (!prim_fail, "failure in integer factorizer");
3213	if (prim_fail)
3214	; //ERROR
3215	else
3216	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3217
3218	if (V_buf4 == alpha)
3219	im_prim_elem_alpha= im_prim_elem;
3220	else
3221	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3222	im_prim_elem, source, dest);
3223
3224	DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3225	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3226	inextension= true;
3227	for (CFListIterator i= l; i.hasItem(); i++)
3228	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3229	im_prim_elem, source, dest);
3230	m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3231	G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3232	newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3233	source, dest);
3234	ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3235	ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3236	gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3237	source, dest);
3238
3239	fail= false;
3240	random_element= randomElement (m, V_buf, l, fail );
3241	DEBOUTLN (cerr, "fail= " << fail);
3242	CFList list;
3243	TIMING_START (gcd_recursion);
3244	G_random_element=
3245	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3246	list, topLevel);
3247	TIMING_END_AND_PRINT (gcd_recursion,
3248	"time for recursive call: ");
3249	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3250	V_buf4= V_buf;
3251	}
3252	else
3253	{
3254	CFList list;
3255	TIMING_START (gcd_recursion);
3256	G_random_element=
3257	sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3258	list, topLevel);
3259	TIMING_END_AND_PRINT (gcd_recursion,
3260	"time for recursive call: ");
3261	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3262	}
3263
3264	if (!G_random_element.inCoeffDomain())
3265	d0= totaldegree (G_random_element, Variable(2),
3266	Variable (G_random_element.level()));
3267	else
3268	d0= 0;
3269
3270	if (d0 == 0)
3271	{
3272	if (inextension)
3273	prune1 (alpha);
3274	return N(gcdcAcB);
3275	}
3276	if (d0 > d)
3277	{
3278	if (!find (l, random_element))
3279	l.append (random_element);
3280	continue;
3281	}
3282
3283	G_random_element=
3284	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3285	* G_random_element;
3286
3287	skeleton= G_random_element;
3288	if (!G_random_element.inCoeffDomain())
3289	d0= totaldegree (G_random_element, Variable(2),
3290	Variable (G_random_element.level()));
3291	else
3292	d0= 0;
3293
3294	if (d0 < d)
3295	{
3296	m= gcdlcAlcB;
3297	newtonPoly= 1;
3298	G_m= 0;
3299	d= d0;
3300	}
3301
3302	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3303	if (uni_lcoeff (H) == gcdlcAlcB)
3304	{
3305	cH= uni_content (H);
3306	ppH= H/cH;
3307	if (inextension)
3308	{
3309	CFList u, v;
3310	//maybe it's better to test if ppH is an element of F(\alpha) before
3311	//mapping down
3312	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3313	{
3314	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3315	ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3316	ppH /= Lc(ppH);
3317	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3318	prune1 (alpha);
3319	return N(gcdcAcB*ppH);
3320	}
3321	}
3322	else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3323	return N(gcdcAcB*ppH);
3324	}
3325	G_m= H;
3326	newtonPoly= newtonPoly*(x - random_element);
3327	m= m*(x - random_element);
3328	if (!find (l, random_element))
3329	l.append (random_element);
3330
3331	if (getCharacteristic () > 3 && size (skeleton) < 100)
3332	{
3333	CFArray Monoms;
3334	CFArray *coeffMonoms;
3335	do //second do
3336	{
3337	random_element= randomElement (m, V_buf, l, fail);
3338	if (random_element == 0 && !fail)
3339	{
3340	if (!find (l, random_element))
3341	l.append (random_element);
3342	continue;
3343	}
3344	if (fail)
3345	{
3346	source= CFList();
3347	dest= CFList();
3348
3349	Variable V_buf3= V_buf;
3350	V_buf= chooseExtension (V_buf);
3351	bool prim_fail= false;
3352	Variable V_buf2;
3353	prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3354	if (V_buf4 == alpha)
3355	prim_elem_alpha= prim_elem;
3356
3357	if (V_buf3 != V_buf4)
3358	{
3359	m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3360	G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3361	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3362	source, dest);
3363	ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3364	ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3365	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3366	source, dest);
3367	for (CFListIterator i= l; i.hasItem(); i++)
3368	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3369	source, dest);
3370	}
3371
3372	ASSERT (!prim_fail, "failure in integer factorizer");
3373	if (prim_fail)
3374	; //ERROR
3375	else
3376	im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3377
3378	if (V_buf4 == alpha)
3379	im_prim_elem_alpha= im_prim_elem;
3380	else
3381	im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3382	prim_elem, im_prim_elem, source, dest);
3383
3384	DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3385	DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3386	inextension= true;
3387	for (CFListIterator i= l; i.hasItem(); i++)
3388	i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3389	im_prim_elem, source, dest);
3390	m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3391	G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3392	newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3393	source, dest);
3394	ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3395	ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3396
3397	gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3398	source, dest);
3399
3400	fail= false;
3401	random_element= randomElement (m, V_buf, l, fail);
3402	DEBOUTLN (cerr, "fail= " << fail);
3403	CFList list;
3404	TIMING_START (gcd_recursion);
3405
3406	V_buf4= V_buf;
3407
3408	//sparseInterpolation
3409	bool sparseFail= false;
3410	if (LC (skeleton).inCoeffDomain())
3411	G_random_element=
3412	monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3413	skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3414	else
3415	G_random_element=
3416	nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3417	skeleton, V_buf, sparseFail, coeffMonoms,
3418	Monoms);
3419	if (sparseFail)
3420	break;
3421
3422	TIMING_END_AND_PRINT (gcd_recursion,
3423	"time for recursive call: ");
3424	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3425	}
3426	else
3427	{
3428	CFList list;
3429	TIMING_START (gcd_recursion);
3430	bool sparseFail= false;
3431	if (LC (skeleton).inCoeffDomain())
3432	G_random_element=
3433	monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3434	skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3435	else
3436	G_random_element=
3437	nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3438	skeleton, V_buf, sparseFail, coeffMonoms,
3439	Monoms);
3440	if (sparseFail)
3441	break;
3442
3443	TIMING_END_AND_PRINT (gcd_recursion,
3444	"time for recursive call: ");
3445	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3446	}
3447
3448	if (!G_random_element.inCoeffDomain())
3449	d0= totaldegree (G_random_element, Variable(2),
3450	Variable (G_random_element.level()));
3451	else
3452	d0= 0;
3453
3454	if (d0 == 0)
3455	{
3456	if (inextension)
3457	prune1 (alpha);
3458	return N(gcdcAcB);
3459	}
3460	if (d0 > d)
3461	{
3462	if (!find (l, random_element))
3463	l.append (random_element);
3464	continue;
3465	}
3466
3467	G_random_element=
3468	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3469	* G_random_element;
3470
3471	if (!G_random_element.inCoeffDomain())
3472	d0= totaldegree (G_random_element, Variable(2),
3473	Variable (G_random_element.level()));
3474	else
3475	d0= 0;
3476
3477	if (d0 < d)
3478	{
3479	m= gcdlcAlcB;
3480	newtonPoly= 1;
3481	G_m= 0;
3482	d= d0;
3483	}
3484
3485	TIMING_START (newton_interpolation);
3486	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3487	TIMING_END_AND_PRINT (newton_interpolation,
3488	"time for newton interpolation: ");
3489
3490	//termination test
3491	if (uni_lcoeff (H) == gcdlcAlcB)
3492	{
3493	cH= uni_content (H);
3494	ppH= H/cH;
3495	if (inextension)
3496	{
3497	CFList u, v;
3498	//maybe it's better to test if ppH is an element of F(\alpha) before
3499	//mapping down
3500	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3501	{
3502	DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3503	ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3504	ppH /= Lc(ppH);
3505	DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3506	prune1 (alpha);
3507	return N(gcdcAcB*ppH);
3508	}
3509	}
3510	else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3511	{
3512	return N(gcdcAcB*ppH);
3513	}
3514	}
3515
3516	G_m= H;
3517	newtonPoly= newtonPoly*(x - random_element);
3518	m= m*(x - random_element);
3519	if (!find (l, random_element))
3520	l.append (random_element);
3521
3522	} while (1);
3523	}
3524	} while (1);
3525	}
3526
3527	CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3528	bool& topLevel, CFList& l)
3529	{
3530	CanonicalForm A= F;
3531	CanonicalForm B= G;
3532	if (F.isZero() && degree(G) > 0) return G/Lc(G);
3533	else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3534	else if (F.isZero() && G.isZero()) return F.genOne();
3535	if (F.inBaseDomain() \|\| G.inBaseDomain()) return F.genOne();
3536	if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3537	if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3538	if (F == G) return F/Lc(F);
3539
3540	CFMap M,N;
3541	int best_level= myCompress (A, B, M, N, topLevel);
3542
3543	if (best_level == 0) return B.genOne();
3544
3545	A= M(A);
3546	B= M(B);
3547
3548	Variable x= Variable (1);
3549
3550	//univariate case
3551	if (A.isUnivariate() && B.isUnivariate())
3552	return N (gcd (A, B));
3553
3554	CanonicalForm cA, cB; // content of A and B
3555	CanonicalForm ppA, ppB; // primitive part of A and B
3556	CanonicalForm gcdcAcB;
3557
3558	cA = uni_content (A);
3559	cB = uni_content (B);
3560	gcdcAcB= gcd (cA, cB);
3561	ppA= A/cA;
3562	ppB= B/cB;
3563
3564	CanonicalForm lcA, lcB; // leading coefficients of A and B
3565	CanonicalForm gcdlcAlcB;
3566	lcA= uni_lcoeff (ppA);
3567	lcB= uni_lcoeff (ppB);
3568
3569	if (fdivides (lcA, lcB))
3570	{
3571	if (fdivides (A, B))
3572	return F/Lc(F);
3573	}
3574	if (fdivides (lcB, lcA))
3575	{
3576	if (fdivides (B, A))
3577	return G/Lc(G);
3578	}
3579
3580	gcdlcAlcB= gcd (lcA, lcB);
3581
3582	int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3583	int d0;
3584
3585	if (d == 0)
3586	return N(gcdcAcB);
3587
3588	d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3589
3590	if (d0 < d)
3591	d= d0;
3592
3593	if (d == 0)
3594	return N(gcdcAcB);
3595
3596	CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3597	CanonicalForm newtonPoly= 1;
3598	m= gcdlcAlcB;
3599	G_m= 0;
3600	H= 0;
3601	bool fail= false;
3602	topLevel= false;
3603	bool inextension= false;
3604	Variable V_buf, alpha, cleanUp;
3605	CanonicalForm prim_elem, im_prim_elem;
3606	CFList source, dest;
3607	do //first do
3608	{
3609	if (inextension)
3610	random_element= randomElement (m, V_buf, l, fail);
3611	else
3612	random_element= FpRandomElement (m, l, fail);
3613	if (random_element == 0 && !fail)
3614	{
3615	if (!find (l, random_element))
3616	l.append (random_element);
3617	continue;
3618	}
3619
3620	if (!fail && !inextension)
3621	{
3622	CFList list;
3623	TIMING_START (gcd_recursion);
3624	G_random_element=
3625	sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3626	list);
3627	TIMING_END_AND_PRINT (gcd_recursion,
3628	"time for recursive call: ");
3629	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3630	}
3631	else if (!fail && inextension)
3632	{
3633	CFList list;
3634	TIMING_START (gcd_recursion);
3635	G_random_element=
3636	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3637	list, topLevel);
3638	TIMING_END_AND_PRINT (gcd_recursion,
3639	"time for recursive call: ");
3640	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3641	}
3642	else if (fail && !inextension)
3643	{
3644	source= CFList();
3645	dest= CFList();
3646	CFList list;
3647	CanonicalForm mipo;
3648	int deg= 2;
3649	bool initialized= false;
3650	do
3651	{
3652	mipo= randomIrredpoly (deg, x);
3653	if (initialized)
3654	setMipo (alpha, mipo);
3655	else
3656	alpha= rootOf (mipo);
3657	inextension= true;
3658	fail= false;
3659	initialized= true;
3660	random_element= randomElement (m, alpha, l, fail);
3661	deg++;
3662	} while (fail);
3663	cleanUp= alpha;
3664	V_buf= alpha;
3665	list= CFList();
3666	TIMING_START (gcd_recursion);
3667	G_random_element=
3668	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3669	list, topLevel);
3670	TIMING_END_AND_PRINT (gcd_recursion,
3671	"time for recursive call: ");
3672	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3673	}
3674	else if (fail && inextension)
3675	{
3676	source= CFList();
3677	dest= CFList();
3678
3679	Variable V_buf3= V_buf;
3680	V_buf= chooseExtension (V_buf);
3681	bool prim_fail= false;
3682	Variable V_buf2;
3683	CanonicalForm prim_elem, im_prim_elem;
3684	prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3685
3686	if (V_buf3 != alpha)
3687	{
3688	m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3689	G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3690	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3691	dest);
3692	ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3693	ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3694	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3695	dest);
3696	for (CFListIterator i= l; i.hasItem(); i++)
3697	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3698	source, dest);
3699	}
3700
3701	ASSERT (!prim_fail, "failure in integer factorizer");
3702	if (prim_fail)
3703	; //ERROR
3704	else
3705	im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3706
3707	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3708	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3709
3710	for (CFListIterator i= l; i.hasItem(); i++)
3711	i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3712	im_prim_elem, source, dest);
3713	m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3714	G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3715	newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3716	source, dest);
3717	ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3718	ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3719	gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3720	source, dest);
3721	fail= false;
3722	random_element= randomElement (m, V_buf, l, fail );
3723	DEBOUTLN (cerr, "fail= " << fail);
3724	CFList list;
3725	TIMING_START (gcd_recursion);
3726	G_random_element=
3727	sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3728	list, topLevel);
3729	TIMING_END_AND_PRINT (gcd_recursion,
3730	"time for recursive call: ");
3731	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3732	alpha= V_buf;
3733	}
3734
3735	if (!G_random_element.inCoeffDomain())
3736	d0= totaldegree (G_random_element, Variable(2),
3737	Variable (G_random_element.level()));
3738	else
3739	d0= 0;
3740
3741	if (d0 == 0)
3742	{
3743	if (inextension)
3744	prune (cleanUp);
3745	return N(gcdcAcB);
3746	}
3747	if (d0 > d)
3748	{
3749	if (!find (l, random_element))
3750	l.append (random_element);
3751	continue;
3752	}
3753
3754	G_random_element=
3755	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3756	* G_random_element;
3757
3758	skeleton= G_random_element;
3759
3760	if (!G_random_element.inCoeffDomain())
3761	d0= totaldegree (G_random_element, Variable(2),
3762	Variable (G_random_element.level()));
3763	else
3764	d0= 0;
3765
3766	if (d0 < d)
3767	{
3768	m= gcdlcAlcB;
3769	newtonPoly= 1;
3770	G_m= 0;
3771	d= d0;
3772	}
3773
3774	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3775
3776	if (uni_lcoeff (H) == gcdlcAlcB)
3777	{
3778	cH= uni_content (H);
3779	ppH= H/cH;
3780	ppH /= Lc (ppH);
3781	DEBOUTLN (cerr, "ppH= " << ppH);
3782
3783	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3784	{
3785	if (inextension)
3786	prune (cleanUp);
3787	return N(gcdcAcB*ppH);
3788	}
3789	}
3790	G_m= H;
3791	newtonPoly= newtonPoly*(x - random_element);
3792	m= m*(x - random_element);
3793	if (!find (l, random_element))
3794	l.append (random_element);
3795
3796	if ((getCharacteristic() > 3 && size (skeleton) < 200))
3797	{
3798	CFArray Monoms;
3799	CFArray* coeffMonoms;
3800
3801	do //second do
3802	{
3803	if (inextension)
3804	random_element= randomElement (m, alpha, l, fail);
3805	else
3806	random_element= FpRandomElement (m, l, fail);
3807	if (random_element == 0 && !fail)
3808	{
3809	if (!find (l, random_element))
3810	l.append (random_element);
3811	continue;
3812	}
3813
3814	bool sparseFail= false;
3815	if (!fail && !inextension)
3816	{
3817	CFList list;
3818	TIMING_START (gcd_recursion);
3819	if (LC (skeleton).inCoeffDomain())
3820	G_random_element=
3821	monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3822	skeleton, x, sparseFail, coeffMonoms,
3823	Monoms);
3824	else
3825	G_random_element=
3826	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3827	skeleton, x, sparseFail,
3828	coeffMonoms, Monoms);
3829	TIMING_END_AND_PRINT (gcd_recursion,
3830	"time for recursive call: ");
3831	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3832	}
3833	else if (!fail && inextension)
3834	{
3835	CFList list;
3836	TIMING_START (gcd_recursion);
3837	if (LC (skeleton).inCoeffDomain())
3838	G_random_element=
3839	monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3840	skeleton, alpha, sparseFail, coeffMonoms,
3841	Monoms);
3842	else
3843	G_random_element=
3844	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3845	skeleton, alpha, sparseFail, coeffMonoms,
3846	Monoms);
3847	TIMING_END_AND_PRINT (gcd_recursion,
3848	"time for recursive call: ");
3849	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3850	}
3851	else if (fail && !inextension)
3852	{
3853	source= CFList();
3854	dest= CFList();
3855	CFList list;
3856	CanonicalForm mipo;
3857	int deg= 2;
3858	bool initialized= false;
3859	do
3860	{
3861	mipo= randomIrredpoly (deg, x);
3862	if (initialized)
3863	setMipo (alpha, mipo);
3864	else
3865	alpha= rootOf (mipo);
3866	inextension= true;
3867	fail= false;
3868	initialized= true;
3869	random_element= randomElement (m, alpha, l, fail);
3870	deg++;
3871	} while (fail);
3872	cleanUp= alpha;
3873	V_buf= alpha;
3874	list= CFList();
3875	TIMING_START (gcd_recursion);
3876	if (LC (skeleton).inCoeffDomain())
3877	G_random_element=
3878	monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3879	skeleton, alpha, sparseFail, coeffMonoms,
3880	Monoms);
3881	else
3882	G_random_element=
3883	nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3884	skeleton, alpha, sparseFail, coeffMonoms,
3885	Monoms);
3886	TIMING_END_AND_PRINT (gcd_recursion,
3887	"time for recursive call: ");
3888	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3889	}
3890	else if (fail && inextension)
3891	{
3892	source= CFList();
3893	dest= CFList();
3894
3895	Variable V_buf3= V_buf;
3896	V_buf= chooseExtension (V_buf);
3897	bool prim_fail= false;
3898	Variable V_buf2;
3899	CanonicalForm prim_elem, im_prim_elem;
3900	prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3901
3902	if (V_buf3 != alpha)
3903	{
3904	m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3905	G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3906	newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3907	source, dest);
3908	ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3909	ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3910	gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3911	source, dest);
3912	for (CFListIterator i= l; i.hasItem(); i++)
3913	i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3914	source, dest);
3915	}
3916
3917	ASSERT (!prim_fail, "failure in integer factorizer");
3918	if (prim_fail)
3919	; //ERROR
3920	else
3921	im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3922
3923	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3924	DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3925
3926	for (CFListIterator i= l; i.hasItem(); i++)
3927	i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3928	im_prim_elem, source, dest);
3929	m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3930	G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3931	newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3932	source, dest);
3933	ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3934	ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3935	gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3936	source, dest);
3937	fail= false;
3938	random_element= randomElement (m, V_buf, l, fail );
3939	DEBOUTLN (cerr, "fail= " << fail);
3940	CFList list;
3941	TIMING_START (gcd_recursion);
3942	if (LC (skeleton).inCoeffDomain())
3943	G_random_element=
3944	monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3945	skeleton, V_buf, sparseFail, coeffMonoms,
3946	Monoms);
3947	else
3948	G_random_element=
3949	nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3950	skeleton, V_buf, sparseFail,
3951	coeffMonoms, Monoms);
3952	TIMING_END_AND_PRINT (gcd_recursion,
3953	"time for recursive call: ");
3954	DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3955	alpha= V_buf;
3956	}
3957
3958	if (sparseFail)
3959	break;
3960
3961	if (!G_random_element.inCoeffDomain())
3962	d0= totaldegree (G_random_element, Variable(2),
3963	Variable (G_random_element.level()));
3964	else
3965	d0= 0;
3966
3967	if (d0 == 0)
3968	{
3969	if (inextension)
3970	prune (cleanUp);
3971	return N(gcdcAcB);
3972	}
3973	if (d0 > d)
3974	{
3975	if (!find (l, random_element))
3976	l.append (random_element);
3977	continue;
3978	}
3979
3980	G_random_element=
3981	(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3982	* G_random_element;
3983
3984	if (!G_random_element.inCoeffDomain())
3985	d0= totaldegree (G_random_element, Variable(2),
3986	Variable (G_random_element.level()));
3987	else
3988	d0= 0;
3989
3990	if (d0 < d)
3991	{
3992	m= gcdlcAlcB;
3993	newtonPoly= 1;
3994	G_m= 0;
3995	d= d0;
3996	}
3997
3998	TIMING_START (newton_interpolation);
3999	H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4000	TIMING_END_AND_PRINT (newton_interpolation,
4001	"time for newton interpolation: ");
4002
4003	//termination test
4004	if (uni_lcoeff (H) == gcdlcAlcB)
4005	{
4006	cH= uni_content (H);
4007	ppH= H/cH;
4008	ppH /= Lc (ppH);
4009	DEBOUTLN (cerr, "ppH= " << ppH);
4010	if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4011	{
4012	if (inextension)
4013	prune (cleanUp);
4014	return N(gcdcAcB*ppH);
4015	}
4016	}
4017
4018	G_m= H;
4019	newtonPoly= newtonPoly*(x - random_element);
4020	m= m*(x - random_element);
4021	if (!find (l, random_element))
4022	l.append (random_element);
4023
4024	} while (1); //end of second do
4025	}
4026	else
4027	{
4028	if (inextension)
4029	prune (cleanUp);
4030	return N(gcdcAcB*modGCDFp (ppA, ppB));
4031	}
4032	} while (1); //end of first do
4033	}
4034
4035	TIMING_DEFINE_PRINT(modZ_termination)
4036	TIMING_DEFINE_PRINT(modZ_recursion)
4037
4038	/// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
4039	/// Algebra", Algorithm 7.1
4040	CanonicalForm modGCDZ ( const CanonicalForm & FF, const CanonicalForm & GG )
4041	{
4042	CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
4043	int p, i, dp_deg, d_deg=-1;
4044
4045	CanonicalForm cd ( bCommonDen( FF ));
4046	f=cd*FF;
4047	Variable x= Variable (1);
4048	CanonicalForm cf, cg;
4049	cf= icontent (f);
4050	f /= cf;
4051	//cd = bCommonDen( f ); f *=cd;
4052	//f /=vcontent(f,Variable(1));
4053
4054	cd = bCommonDen( GG );
4055	g=cd*GG;
4056	cg= icontent (g);
4057	g /= cg;
4058	//cd = bCommonDen( g ); g *=cd;
4059	//g /=vcontent(g,Variable(1));
4060
4061	CanonicalForm Dn, test= 0;
4062	CanonicalForm lcf, lcg;
4063	lcf= f.lc();
4064	lcg= g.lc();
4065	cl = gcd (f.lc(),g.lc());
4066	CanonicalForm gcdcfcg= gcd (cf, cg);
4067	CanonicalForm fp, gp;
4068	CanonicalForm b= 1;
4069	int minCommonDeg= 0;
4070	for (i= tmax (f.level(), g.level()); i > 0; i--)
4071	{
4072	if (degree (f, i) <= 0 \|\| degree (g, i) <= 0)
4073	continue;
4074	else
4075	{
4076	minCommonDeg= tmin (degree (g, i), degree (f, i));
4077	break;
4078	}
4079	}
4080	if (i == 0)
4081	return gcdcfcg;
4082	for (; i > 0; i--)
4083	{
4084	if (degree (f, i) <= 0 \|\| degree (g, i) <= 0)
4085	continue;
4086	else
4087	minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
4088	}
4089	b= 2tmin (maxNorm (f), maxNorm (g))abs (cl)*
4090	power (CanonicalForm (2), minCommonDeg);
4091	bool equal= false;
4092	i = cf_getNumBigPrimes() - 1;
4093
4094	CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn, cDn;
4095	int maxNumVars= tmax (getNumVars (f), getNumVars (g));
4096	//Off (SW_RATIONAL);
4097	while ( true )
4098	{
4099	p = cf_getBigPrime( i );
4100	i--;
4101	while ( i >= 0 && mod( cl(lc(f)/cl)(lc(g)/cl), p ) == 0 )
4102	{
4103	p = cf_getBigPrime( i );
4104	i--;
4105	}
4106	//printf("try p=%d\n",p);
4107	setCharacteristic( p );
4108	fp= mapinto (f);
4109	gp= mapinto (g);
4110	TIMING_START (modZ_recursion)
4111	#ifdef HAVE_NTL
4112	if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4113	Dp = modGCDFp (fp, gp, cofp, cogp);
4114	else
4115	{
4116	Dp= gcd_poly (fp, gp);
4117	cofp= fp/Dp;
4118	cogp= gp/Dp;
4119	}
4120	#else
4121	Dp= gcd_poly (fp, gp);
4122	cofp= fp/Dp;
4123	cogp= gp/Dp;
4124	#endif
4125	TIMING_END_AND_PRINT (modZ_recursion,
4126	"time for gcd mod p in modular gcd: ");
4127	Dp /=Dp.lc();
4128	Dp *= mapinto (cl);
4129	cofp /= lc (cofp);
4130	cofp *= mapinto (lcf);
4131	cogp /= lc (cogp);
4132	cogp *= mapinto (lcg);
4133	setCharacteristic( 0 );
4134	dp_deg=totaldegree(Dp);
4135	if ( dp_deg == 0 )
4136	{
4137	//printf(" -> 1\n");
4138	return CanonicalForm(gcdcfcg);
4139	}
4140	if ( q.isZero() )
4141	{
4142	D = mapinto( Dp );
4143	cof= mapinto (cofp);
4144	cog= mapinto (cogp);
4145	d_deg=dp_deg;
4146	q = p;
4147	Dn= balance_p (D, p);
4148	cofn= balance_p (cof, p);
4149	cogn= balance_p (cog, p);
4150	}
4151	else
4152	{
4153	if ( dp_deg == d_deg )
4154	{
4155	chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4156	chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4157	chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4158	cof= newCof;
4159	cog= newCog;
4160	newqh= newq/2;
4161	Dn= balance_p (newD, newq, newqh);
4162	cofn= balance_p (newCof, newq, newqh);
4163	cogn= balance_p (newCog, newq, newqh);
4164	if (test != Dn) //balance_p (newD, newq))
4165	test= balance_p (newD, newq);
4166	else
4167	equal= true;
4168	q = newq;
4169	D = newD;
4170	}
4171	else if ( dp_deg < d_deg )
4172	{
4173	//printf(" were all bad, try more\n");
4174	// all previous p's are bad primes
4175	q = p;
4176	D = mapinto( Dp );
4177	cof= mapinto (cof);
4178	cog= mapinto (cog);
4179	d_deg=dp_deg;
4180	test= 0;
4181	equal= false;
4182	Dn= balance_p (D, p);
4183	cofn= balance_p (cof, p);
4184	cogn= balance_p (cog, p);
4185	}
4186	else
4187	{
4188	//printf(" was bad, try more\n");
4189	}
4190	//else dp_deg > d_deg: bad prime
4191	}
4192	if ( i >= 0 )
4193	{
4194	cDn= icontent (Dn);
4195	Dn /= cDn;
4196	cofn /= cl/cDn;
4197	//cofn /= icontent (cofn);
4198	cogn /= cl/cDn;
4199	//cogn /= icontent (cogn);
4200	TIMING_START (modZ_termination);
4201	if ((terminationTest (f,g, cofn, cogn, Dn)) \|\|
4202	((equal \|\| q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4203	{
4204	TIMING_END_AND_PRINT (modZ_termination,
4205	"time for successful termination in modular gcd: ");
4206	//printf(" -> success\n");
4207	return Dn*gcdcfcg;
4208	}
4209	TIMING_END_AND_PRINT (modZ_termination,
4210	"time for unsuccessful termination in modular gcd: ");
4211	equal= false;
4212	//else: try more primes
4213	}
4214	else
4215	{ // try other method
4216	//printf("try other gcd\n");
4217	Off(SW_USE_CHINREM_GCD);
4218	D=gcd_poly( f, g );
4219	On(SW_USE_CHINREM_GCD);
4220	return D*gcdcfcg;
4221	}
4222	}
4223	}
4224
4225
4226	#endif

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

Download in other formats: