Context Navigation

source: git/factory/cf_gcd_smallp.cc @ 597783

spielwiese

Last change on this file since 597783 was 597783, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: compute cofactors of gcd and use them to test divisibility chg: do not use bounds in gcd over Fq
Property mode set to `100644`
File size: 130.5 KB

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

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

Download in other formats: