Context Navigation

source: git/factory/cf_gcd_smallp.cc @ 1e4b53

spielwiese

Last change on this file since 1e4b53 was 1e4b53, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: first test LC's before doing a full termination test in GCD
Property mode set to `100644`
File size: 131.2 KB

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

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

Download in other formats: