Context Navigation

source: git/factory/cf_gcd_smallp.cc @ d1dc39

spielwiese

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

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

Download in other formats: