Context Navigation

source: git/factory/cfModGcd.cc @ 4e2cc1

spielwiese

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

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

Download in other formats: