Context Navigation

source: git/factory/cf_gcd_smallp.cc @ ed66770

spielwiese

Last change on this file since ed66770 was ea095d, checked in by Martin Lee <martinlee84@…>, 12 years ago
fix: switched off compression
Property mode set to `100644`
File size: 131.0 KB

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

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

Download in other formats: