Context Navigation

source: git/factory/facMul.cc @ c22249

spielwiese

Last change on this file since c22249 was c22249, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: use NTL/FLINT in *=
Property mode set to `100644`
File size: 91.3 KB

Line
1	/*****************************************************************************\
2	* Computer Algebra System SINGULAR
3	\*****************************************************************************/
4	/** @file facMul.cc
5	*
6	* This file implements functions for fast multiplication and division with
7	* remainder
8	*
9	* @author Martin Lee
10	*
11	**/
12	/*****************************************************************************/
13
14	#include "debug.h"
15	#ifdef HAVE_CONFIG_H
16	#include "config.h"
17	#endif /* HAVE_CONFIG_H */
18
19	#include "canonicalform.h"
20	#include "facMul.h"
21	#include "algext.h"
22	#include "cf_util.h"
23	#include "templates/ftmpl_functions.h"
24
25	#ifdef HAVE_NTL
26	#include <NTL/lzz_pEX.h>
27	#include "NTLconvert.h"
28
29	#ifdef HAVE_FLINT
30	#include "FLINTconvert.h"
31	#endif
32
33	// univariate polys
34
35	#ifdef HAVE_FLINT
36	void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
37	{
38	int degAy= degree (A);
39	fmpz_poly_init2 (result, d*(degAy + 1));
40	_fmpz_poly_set_length (result, d*(degAy + 1));
41	CFIterator j;
42	for (CFIterator i= A; i.hasTerms(); i++)
43	{
44	if (i.coeff().inBaseDomain())
45	convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
46	else
47	for (j= i.coeff(); j.hasTerms(); j++)
48	convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
49	j.coeff());
50	}
51	_fmpz_poly_normalise(result);
52	}
53
54
55	CanonicalForm
56	reverseSubstQa (const fmpz_poly_t F, int d, const Variable& x,
57	const Variable& alpha, const CanonicalForm& den)
58	{
59	CanonicalForm result= 0;
60	int i= 0;
61	int degf= fmpz_poly_degree (F);
62	int k= 0;
63	int degfSubK;
64	int repLength;
65	fmpq_poly_t buf;
66	fmpq_poly_t mipo;
67	convertFacCF2Fmpq_poly_t (mipo, getMipo(alpha));
68	while (degf >= k)
69	{
70	degfSubK= degf - k;
71	if (degfSubK >= d)
72	repLength= d;
73	else
74	repLength= degfSubK + 1;
75
76	fmpq_poly_init2 (buf, repLength);
77	_fmpq_poly_set_length (buf, repLength);
78	_fmpz_vec_set (buf->coeffs, F->coeffs + k, repLength);
79	_fmpq_poly_normalise (buf);
80	fmpq_poly_rem (buf, buf, mipo);
81
82	result += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, i);
83	fmpq_poly_clear (buf);
84	i++;
85	k= d*i;
86	}
87	fmpq_poly_clear (mipo);
88	result /= den;
89	return result;
90	}
91
92	CanonicalForm
93	mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
94	const Variable& alpha)
95	{
96	CanonicalForm A= F;
97	CanonicalForm B= G;
98
99	CanonicalForm denA= bCommonDen (A);
100	CanonicalForm denB= bCommonDen (B);
101
102	A *= denA;
103	B *= denB;
104	int degAa= degree (A, alpha);
105	int degBa= degree (B, alpha);
106	int d= degAa + 1 + degBa;
107
108	fmpz_poly_t FLINTA,FLINTB;
109	kronSub (FLINTA, A, d);
110	kronSub (FLINTB, B, d);
111
112	fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
113
114	denA *= denB;
115	A= reverseSubstQa (FLINTA, d, F.mvar(), alpha, denA);
116
117	fmpz_poly_clear (FLINTA);
118	fmpz_poly_clear (FLINTB);
119	return A;
120	}
121
122	CanonicalForm
123	mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
124	{
125	CanonicalForm A= F;
126	CanonicalForm B= G;
127
128	CanonicalForm denA= bCommonDen (A);
129	CanonicalForm denB= bCommonDen (B);
130
131	A *= denA;
132	B *= denB;
133	fmpz_poly_t FLINTA,FLINTB;
134	convertFacCF2Fmpz_poly_t (FLINTA, A);
135	convertFacCF2Fmpz_poly_t (FLINTB, B);
136	fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
137	denA *= denB;
138	A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
139	A /= denA;
140	fmpz_poly_clear (FLINTA);
141	fmpz_poly_clear (FLINTB);
142
143	return A;
144	}
145
146	/*CanonicalForm
147	mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
148	{
149	CanonicalForm A= F;
150	CanonicalForm B= G;
151
152	fmpq_poly_t FLINTA,FLINTB;
153	convertFacCF2Fmpq_poly_t (FLINTA, A);
154	convertFacCF2Fmpq_poly_t (FLINTB, B);
155
156	fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
157	A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
158	fmpq_poly_clear (FLINTA);
159	fmpq_poly_clear (FLINTB);
160	return A;
161	}*/
162
163	CanonicalForm
164	divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
165	{
166	CanonicalForm A= F;
167	CanonicalForm B= G;
168
169	fmpq_poly_t FLINTA,FLINTB;
170	convertFacCF2Fmpq_poly_t (FLINTA, A);
171	convertFacCF2Fmpq_poly_t (FLINTB, B);
172
173	fmpq_poly_div (FLINTA, FLINTA, FLINTB);
174	A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
175
176	fmpq_poly_clear (FLINTA);
177	fmpq_poly_clear (FLINTB);
178	return A;
179	}
180
181	CanonicalForm
182	modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
183	{
184	CanonicalForm A= F;
185	CanonicalForm B= G;
186
187	fmpq_poly_t FLINTA,FLINTB;
188	convertFacCF2Fmpq_poly_t (FLINTA, A);
189	convertFacCF2Fmpq_poly_t (FLINTB, B);
190
191	fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
192	A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
193
194	fmpq_poly_clear (FLINTA);
195	fmpq_poly_clear (FLINTB);
196	return A;
197	}
198
199	CanonicalForm
200	mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
201	const Variable& alpha, int m)
202	{
203	CanonicalForm A= F;
204	CanonicalForm B= G;
205
206	CanonicalForm denA= bCommonDen (A);
207	CanonicalForm denB= bCommonDen (B);
208
209	A *= denA;
210	B *= denB;
211
212	int degAa= degree (A, alpha);
213	int degBa= degree (B, alpha);
214	int d= degAa + 1 + degBa;
215
216	fmpz_poly_t FLINTA,FLINTB;
217	kronSub (FLINTA, A, d);
218	kronSub (FLINTB, B, d);
219
220	int k= d*m;
221	fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
222
223	denA *= denB;
224	A= reverseSubstQa (FLINTA, d, F.mvar(), alpha, denA);
225	fmpz_poly_clear (FLINTA);
226	fmpz_poly_clear (FLINTB);
227	return A;
228	}
229
230	CanonicalForm
231	mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
232	{
233	if (F.inCoeffDomain() \|\| G.inCoeffDomain())
234	return mod (F*G, power (Variable (1), m));
235	Variable alpha;
236	if (hasFirstAlgVar (F, alpha) \|\| hasFirstAlgVar (G, alpha))
237	return mulFLINTQaTrunc (F, G, alpha, m);
238
239	CanonicalForm A= F;
240	CanonicalForm B= G;
241
242	CanonicalForm denA= bCommonDen (A);
243	CanonicalForm denB= bCommonDen (B);
244
245	A *= denA;
246	B *= denB;
247	fmpz_poly_t FLINTA,FLINTB;
248	convertFacCF2Fmpz_poly_t (FLINTA, A);
249	convertFacCF2Fmpz_poly_t (FLINTB, B);
250	fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
251	denA *= denB;
252	A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
253	A /= denA;
254	fmpz_poly_clear (FLINTA);
255	fmpz_poly_clear (FLINTB);
256
257	return A;
258	}
259
260	CanonicalForm uniReverse (const CanonicalForm& F, int d)
261	{
262	if (d == 0)
263	return F;
264	if (F.inCoeffDomain())
265	return F*power (Variable (1),d);
266	Variable x= Variable (1);
267	CanonicalForm result= 0;
268	CFIterator i= F;
269	while (d - i.exp() < 0)
270	i++;
271
272	for (; i.hasTerms() && (d - i.exp() >= 0); i++)
273	result += i.coeff()*power (x, d - i.exp());
274	return result;
275	}
276
277	CanonicalForm
278	newtonInverse (const CanonicalForm& F, const int n)
279	{
280	int l= ilog2(n);
281
282	CanonicalForm g= F [0];
283
284	ASSERT (!g.isZero(), "expected a unit");
285
286	if (!g.isOne())
287	g = 1/g;
288	Variable x= Variable (1);
289	CanonicalForm result;
290	int exp= 0;
291	if (n & 1)
292	{
293	result= g;
294	exp= 1;
295	}
296	CanonicalForm h;
297
298	for (int i= 1; i <= l; i++)
299	{
300	h= mulNTL (g, mod (F, power (x, (1 << i))));
301	h= mod (h, power (x, (1 << i)) - 1);
302	h= div (h, power (x, (1 << (i - 1))));
303	g -= power (x, (1 << (i - 1)))*
304	mulFLINTQTrunc (g, h, 1 << (i-1));
305
306	if (n & (1 << i))
307	{
308	if (exp)
309	{
310	h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
311	h= mod (h, power (x, exp + (1 << i)) - 1);
312	h= div (h, power (x, exp));
313	result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
314	exp += (1 << i);
315	}
316	else
317	{
318	exp= (1 << i);
319	result= g;
320	}
321	}
322	}
323
324	return result;
325	}
326
327	void
328	newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
329	CanonicalForm& R)
330	{
331	CanonicalForm A= F;
332	CanonicalForm B= G;
333	Variable x= Variable (1);
334	int degA= degree (A, x);
335	int degB= degree (B, x);
336	int m= degA - degB;
337
338	if (m < 0)
339	{
340	R= A;
341	Q= 0;
342	return;
343	}
344
345	if (degB <= 1)
346	divrem (A, B, Q, R);
347	else
348	{
349	R= uniReverse (A, degA);
350
351	CanonicalForm revB= uniReverse (B, degB);
352	CanonicalForm buf= revB;
353	revB= newtonInverse (revB, m + 1);
354	Q= mulFLINTQTrunc (R, revB, m + 1);
355	Q= uniReverse (Q, m);
356
357	R= A - mulNTL (Q, B);
358	}
359	}
360
361	void
362	newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
363	{
364	CanonicalForm A= F;
365	CanonicalForm B= G;
366	Variable x= Variable (1);
367	int degA= degree (A, x);
368	int degB= degree (B, x);
369	int m= degA - degB;
370
371	if (m < 0)
372	{
373	Q= 0;
374	return;
375	}
376
377	if (degB <= 1)
378	Q= div (A, B);
379	else
380	{
381	CanonicalForm R= uniReverse (A, degA);
382
383	CanonicalForm revB= uniReverse (B, degB);
384	revB= newtonInverse (revB, m + 1);
385	Q= mulFLINTQTrunc (R, revB, m + 1);
386	Q= uniReverse (Q, m);
387	}
388	}
389
390	#endif
391
392	CanonicalForm
393	mulNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
394	{
395	if (CFFactory::gettype() == GaloisFieldDomain)
396	return F*G;
397	if (getCharacteristic() == 0)
398	{
399	Variable alpha;
400	if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
401	(hasFirstAlgVar (F, alpha) \|\| hasFirstAlgVar (G, alpha)))
402	{
403	if (b.getp() != 0)
404	{
405	CanonicalForm mipo= getMipo (alpha);
406	bool is_rat= isOn (SW_RATIONAL);
407	if (!is_rat)
408	On (SW_RATIONAL);
409	mipo *=bCommonDen (mipo);
410	if (!is_rat)
411	Off (SW_RATIONAL);
412	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
413	fmpz_t FLINTp;
414	fmpz_mod_poly_t FLINTmipo;
415	fq_ctx_t fq_con;
416	fq_poly_t FLINTF, FLINTG;
417
418	fmpz_init (FLINTp);
419
420	convertCF2Fmpz (FLINTp, b.getpk());
421
422	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, mipo, FLINTp);
423
424	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
425
426	convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
427	convertFacCF2Fq_poly_t (FLINTG, G, fq_con);
428
429	fq_poly_mul (FLINTF, FLINTF, FLINTG, fq_con);
430
431	CanonicalForm result= convertFq_poly_t2FacCF (FLINTF, F.mvar(),
432	alpha, fq_con);
433
434	fmpz_clear (FLINTp);
435	fmpz_mod_poly_clear (FLINTmipo);
436	fq_poly_clear (FLINTF, fq_con);
437	fq_poly_clear (FLINTG, fq_con);
438	fq_ctx_clear (fq_con);
439	return b (result);
440	#else
441	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
442	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (mipo));
443	ZZ_pE::init (NTLmipo);
444	ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
445	ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
446	mul (NTLf, NTLf, NTLg);
447
448	return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
449	#endif
450	}
451	#ifdef HAVE_FLINT
452	CanonicalForm result= mulFLINTQa (F, G, alpha);
453	return result;
454	#else
455	return F*G;
456	#endif
457	}
458	else if (!F.inCoeffDomain() && !G.inCoeffDomain())
459	{
460	#ifdef HAVE_FLINT
461	if (b.getp() != 0)
462	{
463	fmpz_t FLINTpk;
464	fmpz_init (FLINTpk);
465	convertCF2Fmpz (FLINTpk, b.getpk());
466	fmpz_mod_poly_t FLINTF, FLINTG;
467	convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
468	convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
469	fmpz_mod_poly_mul (FLINTF, FLINTF, FLINTG);
470	CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF, F.mvar(),b);
471	fmpz_mod_poly_clear (FLINTG);
472	fmpz_mod_poly_clear (FLINTF);
473	fmpz_clear (FLINTpk);
474	return result;
475	}
476	return mulFLINTQ (F, G);
477	#else
478	if (b.getp() != 0)
479	{
480	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
481	ZZX ZZf= convertFacCF2NTLZZX (F);
482	ZZX ZZg= convertFacCF2NTLZZX (G);
483	ZZ_pX NTLf= to_ZZ_pX (ZZf);
484	ZZ_pX NTLg= to_ZZ_pX (ZZg);
485	mul (NTLf, NTLf, NTLg);
486	return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
487	}
488	return F*G;
489	#endif
490	}
491	if (b.getp() != 0)
492	{
493	if (!F.inBaseDomain() && !G.inBaseDomain())
494	{
495	if (hasFirstAlgVar (G, alpha) \|\| hasFirstAlgVar (F, alpha))
496	{
497	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
498	fmpz_t FLINTp;
499	fmpz_mod_poly_t FLINTmipo;
500	fq_ctx_t fq_con;
501
502	fmpz_init (FLINTp);
503	convertCF2Fmpz (FLINTp, b.getpk());
504
505	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
506
507	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
508
509	CanonicalForm result;
510
511	if (F.inCoeffDomain() && !G.inCoeffDomain())
512	{
513	fq_poly_t FLINTG;
514	fmpz_poly_t FLINTF;
515	convertFacCF2Fmpz_poly_t (FLINTF, F);
516	convertFacCF2Fq_poly_t (FLINTG, G, fq_con);
517
518	fq_poly_scalar_mul_fq (FLINTG, FLINTG, FLINTF, fq_con);
519
520	result= convertFq_poly_t2FacCF (FLINTG, G.mvar(), alpha, fq_con);
521	fmpz_poly_clear (FLINTF);
522	fq_poly_clear (FLINTG, fq_con);
523	}
524	else if (!F.inCoeffDomain() && G.inCoeffDomain())
525	{
526	fq_poly_t FLINTF;
527	fmpz_poly_t FLINTG;
528
529	convertFacCF2Fmpz_poly_t (FLINTG, G);
530	convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
531
532	fq_poly_scalar_mul_fq (FLINTF, FLINTF, FLINTG, fq_con);
533
534	result= convertFq_poly_t2FacCF (FLINTF, F.mvar(), alpha, fq_con);
535	fmpz_poly_clear (FLINTG);
536	fq_poly_clear (FLINTF, fq_con);
537	}
538	else
539	{
540	fq_t FLINTF, FLINTG;
541
542	convertFacCF2Fq_t (FLINTF, F, fq_con);
543	convertFacCF2Fq_t (FLINTG, G, fq_con);
544
545	fq_mul (FLINTF, FLINTF, FLINTG, fq_con);
546
547	result= convertFq_t2FacCF (FLINTF, alpha);
548	fq_clear (FLINTF, fq_con);
549	fq_clear (FLINTG, fq_con);
550	}
551
552	fmpz_clear (FLINTp);
553	fmpz_mod_poly_clear (FLINTmipo);
554	fq_ctx_clear (fq_con);
555
556	return b (result);
557	#else
558	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
559	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
560	ZZ_pE::init (NTLmipo);
561
562	if (F.inCoeffDomain() && !G.inCoeffDomain())
563	{
564	ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
565	ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
566	mul (NTLg, to_ZZ_pE (NTLf), NTLg);
567	return b (convertNTLZZ_pEX2CF (NTLg, G.mvar(), alpha));
568	}
569	else if (!F.inCoeffDomain() && G.inCoeffDomain())
570	{
571	ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
572	ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
573	mul (NTLf, NTLf, to_ZZ_pE (NTLg));
574	return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
575	}
576	else
577	{
578	ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
579	ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
580	ZZ_pE result;
581	mul (result, to_ZZ_pE (NTLg), to_ZZ_pE (NTLf));
582	return b (convertNTLZZpX2CF (rep (result), alpha));
583	}
584	#endif
585	}
586	}
587	return b (F*G);
588	}
589	return F*G;
590	}
591	else if (F.inCoeffDomain() \|\| G.inCoeffDomain())
592	return F*G;
593	ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
594	ASSERT (F.level() == G.level(), "expected polys of same level");
595	if (fac_NTL_char != getCharacteristic())
596	{
597	fac_NTL_char= getCharacteristic();
598	zz_p::init (getCharacteristic());
599	}
600	Variable alpha;
601	CanonicalForm result;
602	if (hasFirstAlgVar (F, alpha) \|\| hasFirstAlgVar (G, alpha))
603	{
604	if (!getReduce (alpha))
605	{
606	result= 0;
607	for (CFIterator i= F; i.hasTerms(); i++)
608	result += i.coeff()Gpower (F.mvar(),i.exp());
609	return result;
610	}
611	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
612	nmod_poly_t FLINTmipo;
613	fq_nmod_ctx_t fq_con;
614
615	nmod_poly_init (FLINTmipo, getCharacteristic());
616	convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
617
618	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
619
620	fq_nmod_poly_t FLINTF, FLINTG;
621	convertFacCF2Fq_nmod_poly_t (FLINTF, F, fq_con);
622	convertFacCF2Fq_nmod_poly_t (FLINTG, G, fq_con);
623
624	fq_nmod_poly_mul (FLINTF, FLINTF, FLINTG, fq_con);
625
626	result= convertFq_nmod_poly_t2FacCF (FLINTF, F.mvar(), alpha, fq_con);
627
628	fq_nmod_poly_clear (FLINTF, fq_con);
629	fq_nmod_poly_clear (FLINTG, fq_con);
630	nmod_poly_clear (FLINTmipo);
631	fq_nmod_ctx_clear (fq_con);
632	#else
633	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
634	zz_pE::init (NTLMipo);
635	zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
636	zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
637	mul (NTLF, NTLF, NTLG);
638	result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
639	#endif
640	}
641	else
642	{
643	#ifdef HAVE_FLINT
644	nmod_poly_t FLINTF, FLINTG;
645	convertFacCF2nmod_poly_t (FLINTF, F);
646	convertFacCF2nmod_poly_t (FLINTG, G);
647	nmod_poly_mul (FLINTF, FLINTF, FLINTG);
648	result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
649	nmod_poly_clear (FLINTF);
650	nmod_poly_clear (FLINTG);
651	#else
652	zz_pX NTLF= convertFacCF2NTLzzpX (F);
653	zz_pX NTLG= convertFacCF2NTLzzpX (G);
654	mul (NTLF, NTLF, NTLG);
655	result= convertNTLzzpX2CF(NTLF, F.mvar());
656	#endif
657	}
658	return result;
659	}
660
661	CanonicalForm
662	modNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
663	{
664	if (CFFactory::gettype() == GaloisFieldDomain)
665	return mod (F, G);
666	if (F.inCoeffDomain() && G.isUnivariate() && !G.inCoeffDomain())
667	{
668	if (b.getp() != 0)
669	return b(F);
670	return F;
671	}
672	else if (F.inCoeffDomain() && G.inCoeffDomain())
673	{
674	if (b.getp() != 0)
675	return b(F%G);
676	return mod (F, G);
677	}
678	else if (F.isUnivariate() && G.inCoeffDomain())
679	{
680	if (b.getp() != 0)
681	return b(F%G);
682	return mod (F,G);
683	}
684
685	if (getCharacteristic() == 0)
686	{
687	Variable alpha;
688	if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
689	{
690	#ifdef HAVE_FLINT
691	if (b.getp() != 0)
692	{
693	fmpz_t FLINTpk;
694	fmpz_init (FLINTpk);
695	convertCF2Fmpz (FLINTpk, b.getpk());
696	fmpz_mod_poly_t FLINTF, FLINTG;
697	convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
698	convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
699	fmpz_mod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
700	CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF,F.mvar(),b);
701	fmpz_mod_poly_clear (FLINTG);
702	fmpz_mod_poly_clear (FLINTF);
703	fmpz_clear (FLINTpk);
704	return result;
705	}
706	return modFLINTQ (F, G);
707	#else
708	if (b.getp() != 0)
709	{
710	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
711	ZZX ZZf= convertFacCF2NTLZZX (F);
712	ZZX ZZg= convertFacCF2NTLZZX (G);
713	ZZ_pX NTLf= to_ZZ_pX (ZZf);
714	ZZ_pX NTLg= to_ZZ_pX (ZZg);
715	rem (NTLf, NTLf, NTLg);
716	return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
717	}
718	return mod (F, G);
719	#endif
720	}
721	else
722	{
723	if (b.getp() != 0)
724	{
725	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
726	fmpz_t FLINTp;
727	fmpz_mod_poly_t FLINTmipo;
728	fq_ctx_t fq_con;
729	fq_poly_t FLINTF, FLINTG;
730
731	fmpz_init (FLINTp);
732
733	convertCF2Fmpz (FLINTp, b.getpk());
734
735	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
736
737	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
738
739	convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
740	convertFacCF2Fq_poly_t (FLINTG, G, fq_con);
741
742	fq_poly_rem (FLINTF, FLINTF, FLINTG, fq_con);
743
744	CanonicalForm result= convertFq_poly_t2FacCF (FLINTF, F.mvar(),
745	alpha, fq_con);
746
747	fmpz_clear (FLINTp);
748	fmpz_mod_poly_clear (FLINTmipo);
749	fq_poly_clear (FLINTF, fq_con);
750	fq_poly_clear (FLINTG, fq_con);
751	fq_ctx_clear (fq_con);
752
753	return b(result);
754	#else
755	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
756	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
757	ZZ_pE::init (NTLmipo);
758	ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
759	ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
760	rem (NTLf, NTLf, NTLg);
761	return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
762	#endif
763	}
764	#ifdef HAVE_FLINT
765	CanonicalForm Q, R;
766	newtonDivrem (F, G, Q, R);
767	return R;
768	#else
769	return mod (F,G);
770	#endif
771	}
772	}
773
774	ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
775	ASSERT (F.level() == G.level(), "expected polys of same level");
776	if (fac_NTL_char != getCharacteristic())
777	{
778	fac_NTL_char= getCharacteristic();
779	zz_p::init (getCharacteristic());
780	}
781	Variable alpha;
782	CanonicalForm result;
783	if (hasFirstAlgVar (F, alpha) \|\| hasFirstAlgVar (G, alpha))
784	{
785	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
786	nmod_poly_t FLINTmipo;
787	fq_nmod_ctx_t fq_con;
788
789	nmod_poly_init (FLINTmipo, getCharacteristic());
790	convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
791
792	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
793
794	fq_nmod_poly_t FLINTF, FLINTG;
795	convertFacCF2Fq_nmod_poly_t (FLINTF, F, fq_con);
796	convertFacCF2Fq_nmod_poly_t (FLINTG, G, fq_con);
797
798	fq_nmod_poly_rem (FLINTF, FLINTF, FLINTG, fq_con);
799
800	result= convertFq_nmod_poly_t2FacCF (FLINTF, F.mvar(), alpha, fq_con);
801
802	fq_nmod_poly_clear (FLINTF, fq_con);
803	fq_nmod_poly_clear (FLINTG, fq_con);
804	nmod_poly_clear (FLINTmipo);
805	fq_nmod_ctx_clear (fq_con);
806	#else
807	zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
808	zz_pE::init (NTLMipo);
809	zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
810	zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
811	rem (NTLF, NTLF, NTLG);
812	result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
813	#endif
814	}
815	else
816	{
817	#ifdef HAVE_FLINT
818	nmod_poly_t FLINTF, FLINTG;
819	convertFacCF2nmod_poly_t (FLINTF, F);
820	convertFacCF2nmod_poly_t (FLINTG, G);
821	nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
822	result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
823	nmod_poly_clear (FLINTF);
824	nmod_poly_clear (FLINTG);
825	#else
826	zz_pX NTLF= convertFacCF2NTLzzpX (F);
827	zz_pX NTLG= convertFacCF2NTLzzpX (G);
828	rem (NTLF, NTLF, NTLG);
829	result= convertNTLzzpX2CF(NTLF, F.mvar());
830	#endif
831	}
832	return result;
833	}
834
835	CanonicalForm
836	divNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
837	{
838	if (CFFactory::gettype() == GaloisFieldDomain)
839	return div (F, G);
840	if (F.inCoeffDomain() && G.isUnivariate() && !G.inCoeffDomain())
841	{
842	return 0;
843	}
844	else if (F.inCoeffDomain() && G.inCoeffDomain())
845	{
846	if (b.getp() != 0)
847	{
848	if (!F.inBaseDomain() \|\| !G.inBaseDomain())
849	{
850	Variable alpha;
851	hasFirstAlgVar (F, alpha);
852	hasFirstAlgVar (G, alpha);
853	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
854	fmpz_t FLINTp;
855	fmpz_mod_poly_t FLINTmipo;
856	fq_ctx_t fq_con;
857	fq_t FLINTF, FLINTG;
858
859	fmpz_init (FLINTp);
860	convertCF2Fmpz (FLINTp, b.getpk());
861
862	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
863
864	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
865
866	convertFacCF2Fq_t (FLINTF, F, fq_con);
867	convertFacCF2Fq_t (FLINTG, G, fq_con);
868
869	fq_inv (FLINTG, FLINTG, fq_con);
870	fq_mul (FLINTF, FLINTF, FLINTG, fq_con);
871
872	CanonicalForm result= convertFq_t2FacCF (FLINTF, alpha);
873
874	fmpz_clear (FLINTp);
875	fmpz_mod_poly_clear (FLINTmipo);
876	fq_clear (FLINTF, fq_con);
877	fq_clear (FLINTG, fq_con);
878	fq_ctx_clear (fq_con);
879	return b (result);
880	#else
881	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
882	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
883	ZZ_pE::init (NTLmipo);
884	ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
885	ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
886	ZZ_pE result;
887	div (result, to_ZZ_pE (NTLf), to_ZZ_pE (NTLg));
888	return b (convertNTLZZpX2CF (rep (result), alpha));
889	#endif
890	}
891	return b(div (F,G));
892	}
893	return div (F, G);
894	}
895	else if (F.isUnivariate() && G.inCoeffDomain())
896	{
897	if (b.getp() != 0)
898	{
899	if (!G.inBaseDomain())
900	{
901	Variable alpha;
902	hasFirstAlgVar (G, alpha);
903	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
904	fmpz_t FLINTp;
905	fmpz_mod_poly_t FLINTmipo;
906	fq_ctx_t fq_con;
907	fq_poly_t FLINTF;
908	fq_t FLINTG;
909
910	fmpz_init (FLINTp);
911	convertCF2Fmpz (FLINTp, b.getpk());
912
913	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
914
915	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
916
917	convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
918	convertFacCF2Fq_t (FLINTG, G, fq_con);
919
920	fq_inv (FLINTG, FLINTG, fq_con);
921	fq_poly_scalar_mul_fq (FLINTF, FLINTF, FLINTG, fq_con);
922
923	CanonicalForm result= convertFq_poly_t2FacCF (FLINTF, F.mvar(),
924	alpha, fq_con);
925
926	fmpz_clear (FLINTp);
927	fmpz_mod_poly_clear (FLINTmipo);
928	fq_poly_clear (FLINTF, fq_con);
929	fq_clear (FLINTG, fq_con);
930	fq_ctx_clear (fq_con);
931	return b (result);
932	#else
933	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
934	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
935	ZZ_pE::init (NTLmipo);
936	ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
937	ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
938	div (NTLf, NTLf, to_ZZ_pE (NTLg));
939	return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
940	#endif
941	}
942	return b(div (F,G));
943	}
944	return div (F, G);
945	}
946
947	if (getCharacteristic() == 0)
948	{
949
950	Variable alpha;
951	if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
952	{
953	#ifdef HAVE_FLINT
954	if (b.getp() != 0)
955	{
956	fmpz_t FLINTpk;
957	fmpz_init (FLINTpk);
958	convertCF2Fmpz (FLINTpk, b.getpk());
959	fmpz_mod_poly_t FLINTF, FLINTG;
960	convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
961	convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
962	fmpz_mod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG);
963	CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF,F.mvar(),b);
964	fmpz_mod_poly_clear (FLINTG);
965	fmpz_mod_poly_clear (FLINTF);
966	fmpz_clear (FLINTpk);
967	return result;
968	}
969	return divFLINTQ (F,G);
970	#else
971	if (b.getp() != 0)
972	{
973	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
974	ZZX ZZf= convertFacCF2NTLZZX (F);
975	ZZX ZZg= convertFacCF2NTLZZX (G);
976	ZZ_pX NTLf= to_ZZ_pX (ZZf);
977	ZZ_pX NTLg= to_ZZ_pX (ZZg);
978	div (NTLf, NTLf, NTLg);
979	return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
980	}
981	return div (F, G);
982	#endif
983	}
984	else
985	{
986	if (b.getp() != 0)
987	{
988	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
989	fmpz_t FLINTp;
990	fmpz_mod_poly_t FLINTmipo;
991	fq_ctx_t fq_con;
992	fq_poly_t FLINTF, FLINTG;
993
994	fmpz_init (FLINTp);
995	convertCF2Fmpz (FLINTp, b.getpk());
996
997	convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
998
999	fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
1000
1001	convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
1002	convertFacCF2Fq_poly_t (FLINTG, G, fq_con);
1003
1004	fq_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG, fq_con);
1005
1006	CanonicalForm result= convertFq_poly_t2FacCF (FLINTF, F.mvar(),
1007	alpha, fq_con);
1008
1009	fmpz_clear (FLINTp);
1010	fmpz_mod_poly_clear (FLINTmipo);
1011	fq_ctx_clear (fq_con);
1012	fq_poly_clear (FLINTF, fq_con);
1013	fq_poly_clear (FLINTG, fq_con);
1014	return b (result);
1015	#else
1016	ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
1017	ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
1018	ZZ_pE::init (NTLmipo);
1019	ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
1020	ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
1021	div (NTLf, NTLf, NTLg);
1022	return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
1023	#endif
1024	}
1025	#ifdef HAVE_FLINT
1026	CanonicalForm Q;
1027	newtonDiv (F, G, Q);
1028	return Q;
1029	#else
1030	return div (F,G);
1031	#endif
1032	}
1033	}
1034
1035	ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
1036	ASSERT (F.level() == G.level(), "expected polys of same level");
1037	if (fac_NTL_char != getCharacteristic())
1038	{
1039	fac_NTL_char= getCharacteristic();
1040	zz_p::init (getCharacteristic());
1041	}
1042	Variable alpha;
1043	CanonicalForm result;
1044	if (hasFirstAlgVar (F, alpha) \|\| hasFirstAlgVar (G, alpha))
1045	{
1046	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
1047	nmod_poly_t FLINTmipo;
1048	fq_nmod_ctx_t fq_con;
1049
1050	nmod_poly_init (FLINTmipo, getCharacteristic());
1051	convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
1052
1053	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
1054
1055	fq_nmod_poly_t FLINTF, FLINTG;
1056	convertFacCF2Fq_nmod_poly_t (FLINTF, F, fq_con);
1057	convertFacCF2Fq_nmod_poly_t (FLINTG, G, fq_con);
1058
1059	fq_nmod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG, fq_con);
1060
1061	result= convertFq_nmod_poly_t2FacCF (FLINTF, F.mvar(), alpha, fq_con);
1062
1063	fq_nmod_poly_clear (FLINTF, fq_con);
1064	fq_nmod_poly_clear (FLINTG, fq_con);
1065	nmod_poly_clear (FLINTmipo);
1066	fq_nmod_ctx_clear (fq_con);
1067	#else
1068	zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
1069	zz_pE::init (NTLMipo);
1070	zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
1071	zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
1072	div (NTLF, NTLF, NTLG);
1073	result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
1074	#endif
1075	}
1076	else
1077	{
1078	#ifdef HAVE_FLINT
1079	nmod_poly_t FLINTF, FLINTG;
1080	convertFacCF2nmod_poly_t (FLINTF, F);
1081	convertFacCF2nmod_poly_t (FLINTG, G);
1082	nmod_poly_div (FLINTF, FLINTF, FLINTG);
1083	result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
1084	nmod_poly_clear (FLINTF);
1085	nmod_poly_clear (FLINTG);
1086	#else
1087	zz_pX NTLF= convertFacCF2NTLzzpX (F);
1088	zz_pX NTLG= convertFacCF2NTLzzpX (G);
1089	div (NTLF, NTLF, NTLG);
1090	result= convertNTLzzpX2CF(NTLF, F.mvar());
1091	#endif
1092	}
1093	return result;
1094	}
1095
1096	// end univariate polys
1097	//*************************
1098	// bivariate polys
1099
1100	#ifdef HAVE_FLINT
1101	void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
1102	{
1103	int degAy= degree (A);
1104	nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
1105	result->length= d*(degAy + 1);
1106	flint_mpn_zero (result->coeffs, d*(degAy+1));
1107
1108	nmod_poly_t buf;
1109
1110	int k;
1111	for (CFIterator i= A; i.hasTerms(); i++)
1112	{
1113	convertFacCF2nmod_poly_t (buf, i.coeff());
1114	k= i.exp()*d;
1115	flint_mpn_copyi (result->coeffs+k, buf->coeffs, nmod_poly_length(buf));
1116
1117	nmod_poly_clear (buf);
1118	}
1119	_nmod_poly_normalise (result);
1120	}
1121
1122	#if (__FLINT_VERSION_MINOR >= 4)
1123	void
1124	kronSubFq (fq_nmod_poly_t result, const CanonicalForm& A, int d,
1125	const fq_nmod_ctx_t fq_con)
1126	{
1127	int degAy= degree (A);
1128	fq_nmod_poly_init2 (result, d*(degAy + 1), fq_con);
1129	_fq_nmod_poly_set_length (result, d*(degAy + 1), fq_con);
1130	_fq_nmod_vec_zero (result->coeffs, d*(degAy + 1), fq_con);
1131
1132	fq_nmod_poly_t buf1;
1133
1134	nmod_poly_t buf2;
1135
1136	int k;
1137
1138	for (CFIterator i= A; i.hasTerms(); i++)
1139	{
1140	if (i.coeff().inCoeffDomain())
1141	{
1142	convertFacCF2nmod_poly_t (buf2, i.coeff());
1143	fq_nmod_poly_init2 (buf1, 1, fq_con);
1144	fq_nmod_poly_set_coeff (buf1, 0, buf2, fq_con);
1145	nmod_poly_clear (buf2);
1146	}
1147	else
1148	convertFacCF2Fq_nmod_poly_t (buf1, i.coeff(), fq_con);
1149
1150	k= i.exp()*d;
1151	_fq_nmod_vec_set (result->coeffs+k, buf1->coeffs,
1152	fq_nmod_poly_length (buf1, fq_con), fq_con);
1153
1154	fq_nmod_poly_clear (buf1, fq_con);
1155	}
1156
1157	_fq_nmod_poly_normalise (result, fq_con);
1158	}
1159	#endif
1160
1161	/*void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
1162	{
1163	int degAy= degree (A);
1164	fmpq_poly_init2 (result, d1*(degAy + 1));
1165
1166	fmpq_poly_t buf;
1167	fmpq_t coeff;
1168
1169	int k, l, bufRepLength;
1170	CFIterator j;
1171	for (CFIterator i= A; i.hasTerms(); i++)
1172	{
1173	if (i.coeff().inCoeffDomain())
1174	{
1175	k= d1*i.exp();
1176	convertFacCF2Fmpq_poly_t (buf, i.coeff());
1177	bufRepLength= (int) fmpq_poly_length(buf);
1178	for (l= 0; l < bufRepLength; l++)
1179	{
1180	fmpq_poly_get_coeff_fmpq (coeff, buf, l);
1181	fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
1182	}
1183	fmpq_poly_clear (buf);
1184	}
1185	else
1186	{
1187	for (j= i.coeff(); j.hasTerms(); j++)
1188	{
1189	k= d1*i.exp();
1190	k += d2*j.exp();
1191	convertFacCF2Fmpq_poly_t (buf, j.coeff());
1192	bufRepLength= (int) fmpq_poly_length(buf);
1193	for (l= 0; l < bufRepLength; l++)
1194	{
1195	fmpq_poly_get_coeff_fmpq (coeff, buf, l);
1196	fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
1197	}
1198	fmpq_poly_clear (buf);
1199	}
1200	}
1201	}
1202	fmpq_clear (coeff);
1203	_fmpq_poly_normalise (result);
1204	}*/
1205
1206	void kronSubQa (fmpz_poly_t result, const CanonicalForm& A, int d1, int d2)
1207	{
1208	int degAy= degree (A);
1209	fmpz_poly_init2 (result, d1*(degAy + 1));
1210	_fmpz_poly_set_length (result, d1*(degAy + 1));
1211
1212	fmpz_poly_t buf;
1213
1214	int k;
1215	CFIterator j;
1216	for (CFIterator i= A; i.hasTerms(); i++)
1217	{
1218	if (i.coeff().inCoeffDomain())
1219	{
1220	k= d1*i.exp();
1221	convertFacCF2Fmpz_poly_t (buf, i.coeff());
1222	_fmpz_vec_set (result->coeffs + k, buf->coeffs, buf->length);
1223	fmpz_poly_clear (buf);
1224	}
1225	else
1226	{
1227	for (j= i.coeff(); j.hasTerms(); j++)
1228	{
1229	k= d1*i.exp();
1230	k += d2*j.exp();
1231	convertFacCF2Fmpz_poly_t (buf, j.coeff());
1232	_fmpz_vec_set (result->coeffs + k, buf->coeffs, buf->length);
1233	fmpz_poly_clear (buf);
1234	}
1235	}
1236	}
1237	_fmpz_poly_normalise (result);
1238	}
1239
1240	void
1241	kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
1242	int d)
1243	{
1244	int degAy= degree (A);
1245	mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1246	nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
1247	nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
1248
1249	nmod_poly_t buf;
1250
1251	int k, kk, j, bufRepLength;
1252	for (CFIterator i= A; i.hasTerms(); i++)
1253	{
1254	convertFacCF2nmod_poly_t (buf, i.coeff());
1255
1256	k= i.exp()*d;
1257	kk= (degAy - i.exp())*d;
1258	bufRepLength= (int) nmod_poly_length (buf);
1259	for (j= 0; j < bufRepLength; j++)
1260	{
1261	nmod_poly_set_coeff_ui (subA1, j + k,
1262	n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
1263	nmod_poly_get_coeff_ui (buf, j),
1264	getCharacteristic()
1265	)
1266	);
1267	nmod_poly_set_coeff_ui (subA2, j + kk,
1268	n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
1269	nmod_poly_get_coeff_ui (buf, j),
1270	getCharacteristic()
1271	)
1272	);
1273	}
1274	nmod_poly_clear (buf);
1275	}
1276	_nmod_poly_normalise (subA1);
1277	_nmod_poly_normalise (subA2);
1278	}
1279
1280	#if (__FLINT_VERSION_MINOR >= 4)
1281	void
1282	kronSubReciproFq (fq_nmod_poly_t subA1, fq_nmod_poly_t subA2,
1283	const CanonicalForm& A, int d, const fq_nmod_ctx_t fq_con)
1284	{
1285	int degAy= degree (A);
1286	fq_nmod_poly_init2 (subA1, d*(degAy + 2), fq_con);
1287	fq_nmod_poly_init2 (subA2, d*(degAy + 2), fq_con);
1288
1289	_fq_nmod_poly_set_length (subA1, d*(degAy + 2), fq_con);
1290	_fq_nmod_vec_zero (subA1->coeffs, d*(degAy + 2), fq_con);
1291
1292	_fq_nmod_poly_set_length (subA2, d*(degAy + 2), fq_con);
1293	_fq_nmod_vec_zero (subA2->coeffs, d*(degAy + 2), fq_con);
1294
1295	fq_nmod_poly_t buf1;
1296
1297	nmod_poly_t buf2;
1298
1299	int k, kk;
1300	for (CFIterator i= A; i.hasTerms(); i++)
1301	{
1302	if (i.coeff().inCoeffDomain())
1303	{
1304	convertFacCF2nmod_poly_t (buf2, i.coeff());
1305	fq_nmod_poly_init2 (buf1, 1, fq_con);
1306	fq_nmod_poly_set_coeff (buf1, 0, buf2, fq_con);
1307	nmod_poly_clear (buf2);
1308	}
1309	else
1310	convertFacCF2Fq_nmod_poly_t (buf1, i.coeff(), fq_con);
1311
1312	k= i.exp()*d;
1313	kk= (degAy - i.exp())*d;
1314	_fq_nmod_vec_add (subA1->coeffs+k, subA1->coeffs+k, buf1->coeffs,
1315	fq_nmod_poly_length(buf1, fq_con), fq_con);
1316	_fq_nmod_vec_add (subA2->coeffs+kk, subA2->coeffs+kk, buf1->coeffs,
1317	fq_nmod_poly_length(buf1, fq_con), fq_con);
1318
1319	fq_nmod_poly_clear (buf1, fq_con);
1320	}
1321	_fq_nmod_poly_normalise (subA1, fq_con);
1322	_fq_nmod_poly_normalise (subA2, fq_con);
1323	}
1324	#endif
1325
1326	void
1327	kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
1328	int d)
1329	{
1330	int degAy= degree (A);
1331	fmpz_poly_init2 (subA1, d*(degAy + 2));
1332	fmpz_poly_init2 (subA2, d*(degAy + 2));
1333
1334	fmpz_poly_t buf;
1335
1336	int k, kk;
1337	for (CFIterator i= A; i.hasTerms(); i++)
1338	{
1339	convertFacCF2Fmpz_poly_t (buf, i.coeff());
1340
1341	k= i.exp()*d;
1342	kk= (degAy - i.exp())*d;
1343	_fmpz_vec_add (subA1->coeffs+k, subA1->coeffs + k, buf->coeffs, buf->length);
1344	_fmpz_vec_add (subA2->coeffs+kk, subA2->coeffs + kk, buf->coeffs, buf->length);
1345	fmpz_poly_clear (buf);
1346	}
1347
1348	_fmpz_poly_normalise (subA1);
1349	_fmpz_poly_normalise (subA2);
1350	}
1351
1352	CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
1353	{
1354	Variable y= Variable (2);
1355	Variable x= Variable (1);
1356
1357	fmpz_poly_t buf;
1358	CanonicalForm result= 0;
1359	int i= 0;
1360	int degf= fmpz_poly_degree(F);
1361	int k= 0;
1362	int degfSubK, repLength;
1363	while (degf >= k)
1364	{
1365	degfSubK= degf - k;
1366	if (degfSubK >= d)
1367	repLength= d;
1368	else
1369	repLength= degfSubK + 1;
1370
1371	fmpz_poly_init2 (buf, repLength);
1372	_fmpz_poly_set_length (buf, repLength);
1373	_fmpz_vec_set (buf->coeffs, F->coeffs+k, repLength);
1374	_fmpz_poly_normalise (buf);
1375
1376	result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
1377	i++;
1378	k= d*i;
1379	fmpz_poly_clear (buf);
1380	}
1381
1382	return result;
1383	}
1384
1385	/*CanonicalForm
1386	reverseSubstQa (const fmpq_poly_t F, int d1, int d2, const Variable& alpha,
1387	const fmpq_poly_t mipo)
1388	{
1389	Variable y= Variable (2);
1390	Variable x= Variable (1);
1391
1392	fmpq_poly_t f;
1393	fmpq_poly_init (f);
1394	fmpq_poly_set (f, F);
1395
1396	fmpq_poly_t buf;
1397	CanonicalForm result= 0, result2;
1398	int i= 0;
1399	int degf= fmpq_poly_degree(f);
1400	int k= 0;
1401	int degfSubK;
1402	int repLength;
1403	fmpq_t coeff;
1404	while (degf >= k)
1405	{
1406	degfSubK= degf - k;
1407	if (degfSubK >= d1)
1408	repLength= d1;
1409	else
1410	repLength= degfSubK + 1;
1411
1412	fmpq_init (coeff);
1413	int j= 0;
1414	int l;
1415	result2= 0;
1416	while (j*d2 < repLength)
1417	{
1418	fmpq_poly_init2 (buf, d2);
1419	for (l= 0; l < d2; l++)
1420	{
1421	fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
1422	fmpq_poly_set_coeff_fmpq (buf, l, coeff);
1423	}
1424	_fmpq_poly_normalise (buf);
1425	fmpq_poly_rem (buf, buf, mipo);
1426	result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1427	j++;
1428	fmpq_poly_clear (buf);
1429	}
1430	if (repLength - jd2 != 0 && jd2 - repLength < d2)
1431	{
1432	j--;
1433	repLength -= j*d2;
1434	fmpq_poly_init2 (buf, repLength);
1435	j++;
1436	for (l= 0; l < repLength; l++)
1437	{
1438	fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
1439	fmpq_poly_set_coeff_fmpq (buf, l, coeff);
1440	}
1441	_fmpq_poly_normalise (buf);
1442	fmpq_poly_rem (buf, buf, mipo);
1443	result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1444	fmpq_poly_clear (buf);
1445	}
1446	fmpq_clear (coeff);
1447
1448	result += result2*power (y, i);
1449	i++;
1450	k= d1*i;
1451	}
1452
1453	fmpq_poly_clear (f);
1454	return result;
1455	}*/
1456
1457	CanonicalForm
1458	reverseSubstQa (const fmpz_poly_t F, int d1, int d2, const Variable& alpha,
1459	const fmpq_poly_t mipo)
1460	{
1461	Variable y= Variable (2);
1462	Variable x= Variable (1);
1463
1464	fmpq_poly_t buf;
1465	CanonicalForm result= 0, result2;
1466	int i= 0;
1467	int degf= fmpz_poly_degree(F);
1468	int k= 0;
1469	int degfSubK;
1470	int repLength;
1471	while (degf >= k)
1472	{
1473	degfSubK= degf - k;
1474	if (degfSubK >= d1)
1475	repLength= d1;
1476	else
1477	repLength= degfSubK + 1;
1478
1479	int j= 0;
1480	result2= 0;
1481	while (j*d2 < repLength)
1482	{
1483	fmpq_poly_init2 (buf, d2);
1484	_fmpq_poly_set_length (buf, d2);
1485	_fmpz_vec_set (buf->coeffs, F->coeffs + k + j*d2, d2);
1486	_fmpq_poly_normalise (buf);
1487	fmpq_poly_rem (buf, buf, mipo);
1488	result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1489	j++;
1490	fmpq_poly_clear (buf);
1491	}
1492	if (repLength - jd2 != 0 && jd2 - repLength < d2)
1493	{
1494	j--;
1495	repLength -= j*d2;
1496	fmpq_poly_init2 (buf, repLength);
1497	_fmpq_poly_set_length (buf, repLength);
1498	j++;
1499	_fmpz_vec_set (buf->coeffs, F->coeffs + k + j*d2, repLength);
1500	_fmpq_poly_normalise (buf);
1501	fmpq_poly_rem (buf, buf, mipo);
1502	result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1503	fmpq_poly_clear (buf);
1504	}
1505
1506	result += result2*power (y, i);
1507	i++;
1508	k= d1*i;
1509	}
1510
1511	return result;
1512	}
1513
1514	CanonicalForm
1515	reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
1516	{
1517	Variable y= Variable (2);
1518	Variable x= Variable (1);
1519
1520	nmod_poly_t f, g;
1521	mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1522	nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1523	nmod_poly_init_preinv (g, getCharacteristic(), ninv);
1524	nmod_poly_set (f, F);
1525	nmod_poly_set (g, G);
1526	int degf= nmod_poly_degree(f);
1527	int degg= nmod_poly_degree(g);
1528
1529
1530	nmod_poly_t buf1,buf2, buf3;
1531
1532	if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
1533	nmod_poly_fit_length (f,(long)d*(k+1));
1534
1535	CanonicalForm result= 0;
1536	int i= 0;
1537	int lf= 0;
1538	int lg= d*k;
1539	int degfSubLf= degf;
1540	int deggSubLg= degg-lg;
1541	int repLengthBuf2, repLengthBuf1, ind, tmp;
1542	while (degf >= lf \|\| lg >= 0)
1543	{
1544	if (degfSubLf >= d)
1545	repLengthBuf1= d;
1546	else if (degfSubLf < 0)
1547	repLengthBuf1= 0;
1548	else
1549	repLengthBuf1= degfSubLf + 1;
1550	nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
1551
1552	for (ind= 0; ind < repLengthBuf1; ind++)
1553	nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
1554	_nmod_poly_normalise (buf1);
1555
1556	repLengthBuf1= nmod_poly_length (buf1);
1557
1558	if (deggSubLg >= d - 1)
1559	repLengthBuf2= d - 1;
1560	else if (deggSubLg < 0)
1561	repLengthBuf2= 0;
1562	else
1563	repLengthBuf2= deggSubLg + 1;
1564
1565	nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
1566	for (ind= 0; ind < repLengthBuf2; ind++)
1567	nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
1568
1569	_nmod_poly_normalise (buf2);
1570	repLengthBuf2= nmod_poly_length (buf2);
1571
1572	nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
1573	for (ind= 0; ind < repLengthBuf1; ind++)
1574	nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
1575	for (ind= repLengthBuf1; ind < d; ind++)
1576	nmod_poly_set_coeff_ui (buf3, ind, 0);
1577	for (ind= 0; ind < repLengthBuf2; ind++)
1578	nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
1579	_nmod_poly_normalise (buf3);
1580
1581	result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
1582	i++;
1583
1584
1585	lf= i*d;
1586	degfSubLf= degf - lf;
1587
1588	lg= d*(k-i);
1589	deggSubLg= degg - lg;
1590
1591	if (lg >= 0 && deggSubLg > 0)
1592	{
1593	if (repLengthBuf2 > degfSubLf + 1)
1594	degfSubLf= repLengthBuf2 - 1;
1595	tmp= tmin (repLengthBuf1, deggSubLg + 1);
1596	for (ind= 0; ind < tmp; ind++)
1597	nmod_poly_set_coeff_ui (g, ind + lg,
1598	n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
1599	nmod_poly_get_coeff_ui (buf1, ind),
1600	getCharacteristic()
1601	)
1602	);
1603	}
1604	if (lg < 0)
1605	{
1606	nmod_poly_clear (buf1);
1607	nmod_poly_clear (buf2);
1608	nmod_poly_clear (buf3);
1609	break;
1610	}
1611	if (degfSubLf >= 0)
1612	{
1613	for (ind= 0; ind < repLengthBuf2; ind++)
1614	nmod_poly_set_coeff_ui (f, ind + lf,
1615	n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
1616	nmod_poly_get_coeff_ui (buf2, ind),
1617	getCharacteristic()
1618	)
1619	);
1620	}
1621	nmod_poly_clear (buf1);
1622	nmod_poly_clear (buf2);
1623	nmod_poly_clear (buf3);
1624	}
1625
1626	nmod_poly_clear (f);
1627	nmod_poly_clear (g);
1628
1629	return result;
1630	}
1631
1632	#if (__FLINT_VERSION_MINOR >= 4)
1633	CanonicalForm
1634	reverseSubstReciproFq (const fq_nmod_poly_t F, const fq_nmod_poly_t G, int d,
1635	int k, const Variable& alpha, const fq_nmod_ctx_t fq_con)
1636	{
1637	Variable y= Variable (2);
1638	Variable x= Variable (1);
1639
1640	fq_nmod_poly_t f, g;
1641	int degf= fq_nmod_poly_degree(F, fq_con);
1642	int degg= fq_nmod_poly_degree(G, fq_con);
1643
1644	fq_nmod_poly_t buf1,buf2, buf3;
1645
1646	fq_nmod_poly_init (f, fq_con);
1647	fq_nmod_poly_init (g, fq_con);
1648	fq_nmod_poly_set (f, F, fq_con);
1649	fq_nmod_poly_set (g, G, fq_con);
1650	if (fq_nmod_poly_length (f, fq_con) < (long) d*(k + 1)) //zero padding
1651	fq_nmod_poly_fit_length (f, (long) d*(k + 1), fq_con);
1652
1653	CanonicalForm result= 0;
1654	int i= 0;
1655	int lf= 0;
1656	int lg= d*k;
1657	int degfSubLf= degf;
1658	int deggSubLg= degg-lg;
1659	int repLengthBuf2, repLengthBuf1, tmp;
1660	while (degf >= lf \|\| lg >= 0)
1661	{
1662	if (degfSubLf >= d)
1663	repLengthBuf1= d;
1664	else if (degfSubLf < 0)
1665	repLengthBuf1= 0;
1666	else
1667	repLengthBuf1= degfSubLf + 1;
1668	fq_nmod_poly_init2 (buf1, repLengthBuf1, fq_con);
1669	_fq_nmod_poly_set_length (buf1, repLengthBuf1, fq_con);
1670
1671	_fq_nmod_vec_set (buf1->coeffs, f->coeffs + lf, repLengthBuf1, fq_con);
1672	_fq_nmod_poly_normalise (buf1, fq_con);
1673
1674	repLengthBuf1= fq_nmod_poly_length (buf1, fq_con);
1675
1676	if (deggSubLg >= d - 1)
1677	repLengthBuf2= d - 1;
1678	else if (deggSubLg < 0)
1679	repLengthBuf2= 0;
1680	else
1681	repLengthBuf2= deggSubLg + 1;
1682
1683	fq_nmod_poly_init2 (buf2, repLengthBuf2, fq_con);
1684	_fq_nmod_poly_set_length (buf2, repLengthBuf2, fq_con);
1685	_fq_nmod_vec_set (buf2->coeffs, g->coeffs + lg, repLengthBuf2, fq_con);
1686
1687	_fq_nmod_poly_normalise (buf2, fq_con);
1688	repLengthBuf2= fq_nmod_poly_length (buf2, fq_con);
1689
1690	fq_nmod_poly_init2 (buf3, repLengthBuf2 + d, fq_con);
1691	_fq_nmod_poly_set_length (buf3, repLengthBuf2 + d, fq_con);
1692	_fq_nmod_vec_set (buf3->coeffs, buf1->coeffs, repLengthBuf1, fq_con);
1693	_fq_nmod_vec_set (buf3->coeffs + d, buf2->coeffs, repLengthBuf2, fq_con);
1694
1695	_fq_nmod_poly_normalise (buf3, fq_con);
1696
1697	result += convertFq_nmod_poly_t2FacCF (buf3, x, alpha, fq_con)*power (y, i);
1698	i++;
1699
1700
1701	lf= i*d;
1702	degfSubLf= degf - lf;
1703
1704	lg= d*(k - i);
1705	deggSubLg= degg - lg;
1706
1707	if (lg >= 0 && deggSubLg > 0)
1708	{
1709	if (repLengthBuf2 > degfSubLf + 1)
1710	degfSubLf= repLengthBuf2 - 1;
1711	tmp= tmin (repLengthBuf1, deggSubLg + 1);
1712	_fq_nmod_vec_sub (g->coeffs + lg, g->coeffs + lg, buf1-> coeffs,
1713	tmp, fq_con);
1714	}
1715	if (lg < 0)
1716	{
1717	fq_nmod_poly_clear (buf1, fq_con);
1718	fq_nmod_poly_clear (buf2, fq_con);
1719	fq_nmod_poly_clear (buf3, fq_con);
1720	break;
1721	}
1722	if (degfSubLf >= 0)
1723	_fq_nmod_vec_sub (f->coeffs + lf, f->coeffs + lf, buf2->coeffs,
1724	repLengthBuf2, fq_con);
1725	fq_nmod_poly_clear (buf1, fq_con);
1726	fq_nmod_poly_clear (buf2, fq_con);
1727	fq_nmod_poly_clear (buf3, fq_con);
1728	}
1729
1730	fq_nmod_poly_clear (f, fq_con);
1731	fq_nmod_poly_clear (g, fq_con);
1732
1733	return result;
1734	}
1735	#endif
1736
1737	CanonicalForm
1738	reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
1739	{
1740	Variable y= Variable (2);
1741	Variable x= Variable (1);
1742
1743	fmpz_poly_t f, g;
1744	fmpz_poly_init (f);
1745	fmpz_poly_init (g);
1746	fmpz_poly_set (f, F);
1747	fmpz_poly_set (g, G);
1748	int degf= fmpz_poly_degree(f);
1749	int degg= fmpz_poly_degree(g);
1750
1751	fmpz_poly_t buf1,buf2, buf3;
1752
1753	if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
1754	fmpz_poly_fit_length (f,(long)d*(k+1));
1755
1756	CanonicalForm result= 0;
1757	int i= 0;
1758	int lf= 0;
1759	int lg= d*k;
1760	int degfSubLf= degf;
1761	int deggSubLg= degg-lg;
1762	int repLengthBuf2, repLengthBuf1, ind, tmp;
1763	fmpz_t tmp1, tmp2;
1764	while (degf >= lf \|\| lg >= 0)
1765	{
1766	if (degfSubLf >= d)
1767	repLengthBuf1= d;
1768	else if (degfSubLf < 0)
1769	repLengthBuf1= 0;
1770	else
1771	repLengthBuf1= degfSubLf + 1;
1772
1773	fmpz_poly_init2 (buf1, repLengthBuf1);
1774
1775	for (ind= 0; ind < repLengthBuf1; ind++)
1776	{
1777	fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1778	fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
1779	}
1780	_fmpz_poly_normalise (buf1);
1781
1782	repLengthBuf1= fmpz_poly_length (buf1);
1783
1784	if (deggSubLg >= d - 1)
1785	repLengthBuf2= d - 1;
1786	else if (deggSubLg < 0)
1787	repLengthBuf2= 0;
1788	else
1789	repLengthBuf2= deggSubLg + 1;
1790
1791	fmpz_poly_init2 (buf2, repLengthBuf2);
1792
1793	for (ind= 0; ind < repLengthBuf2; ind++)
1794	{
1795	fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1796	fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
1797	}
1798
1799	_fmpz_poly_normalise (buf2);
1800	repLengthBuf2= fmpz_poly_length (buf2);
1801
1802	fmpz_poly_init2 (buf3, repLengthBuf2 + d);
1803	for (ind= 0; ind < repLengthBuf1; ind++)
1804	{
1805	fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);
1806	fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
1807	}
1808	for (ind= repLengthBuf1; ind < d; ind++)
1809	fmpz_poly_set_coeff_ui (buf3, ind, 0);
1810	for (ind= 0; ind < repLengthBuf2; ind++)
1811	{
1812	fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
1813	fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
1814	}
1815	_fmpz_poly_normalise (buf3);
1816
1817	result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
1818	i++;
1819
1820
1821	lf= i*d;
1822	degfSubLf= degf - lf;
1823
1824	lg= d*(k-i);
1825	deggSubLg= degg - lg;
1826
1827	if (lg >= 0 && deggSubLg > 0)
1828	{
1829	if (repLengthBuf2 > degfSubLf + 1)
1830	degfSubLf= repLengthBuf2 - 1;
1831	tmp= tmin (repLengthBuf1, deggSubLg + 1);
1832	for (ind= 0; ind < tmp; ind++)
1833	{
1834	fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1835	fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
1836	fmpz_sub (tmp1, tmp1, tmp2);
1837	fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
1838	}
1839	}
1840	if (lg < 0)
1841	{
1842	fmpz_poly_clear (buf1);
1843	fmpz_poly_clear (buf2);
1844	fmpz_poly_clear (buf3);
1845	break;
1846	}
1847	if (degfSubLf >= 0)
1848	{
1849	for (ind= 0; ind < repLengthBuf2; ind++)
1850	{
1851	fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1852	fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
1853	fmpz_sub (tmp1, tmp1, tmp2);
1854	fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
1855	}
1856	}
1857	fmpz_poly_clear (buf1);
1858	fmpz_poly_clear (buf2);
1859	fmpz_poly_clear (buf3);
1860	}
1861
1862	fmpz_poly_clear (f);
1863	fmpz_poly_clear (g);
1864	fmpz_clear (tmp1);
1865	fmpz_clear (tmp2);
1866
1867	return result;
1868	}
1869
1870	#if (__FLINT_VERSION_MINOR >= 4)
1871	CanonicalForm
1872	reverseSubstFq (const fq_nmod_poly_t F, int d, const Variable& alpha,
1873	const fq_nmod_ctx_t fq_con)
1874	{
1875	Variable y= Variable (2);
1876	Variable x= Variable (1);
1877
1878	fq_nmod_poly_t buf;
1879	CanonicalForm result= 0;
1880	int i= 0;
1881	int degf= fq_nmod_poly_degree(F, fq_con);
1882	int k= 0;
1883	int degfSubK, repLength;
1884	while (degf >= k)
1885	{
1886	degfSubK= degf - k;
1887	if (degfSubK >= d)
1888	repLength= d;
1889	else
1890	repLength= degfSubK + 1;
1891
1892	fq_nmod_poly_init2 (buf, repLength, fq_con);
1893	_fq_nmod_poly_set_length (buf, repLength, fq_con);
1894	_fq_nmod_vec_set (buf->coeffs, F->coeffs+k, repLength, fq_con);
1895	_fq_nmod_poly_normalise (buf, fq_con);
1896
1897	result += convertFq_nmod_poly_t2FacCF (buf, x, alpha, fq_con)*power (y, i);
1898	i++;
1899	k= d*i;
1900	fq_nmod_poly_clear (buf, fq_con);
1901	}
1902
1903	return result;
1904	}
1905	#endif
1906
1907	CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
1908	{
1909	Variable y= Variable (2);
1910	Variable x= Variable (1);
1911
1912	mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1913
1914	nmod_poly_t buf;
1915	CanonicalForm result= 0;
1916	int i= 0;
1917	int degf= nmod_poly_degree(F);
1918	int k= 0;
1919	int degfSubK, repLength, j;
1920	while (degf >= k)
1921	{
1922	degfSubK= degf - k;
1923	if (degfSubK >= d)
1924	repLength= d;
1925	else
1926	repLength= degfSubK + 1;
1927
1928	nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
1929	for (j= 0; j < repLength; j++)
1930	nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (F, j + k));
1931	_nmod_poly_normalise (buf);
1932
1933	result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
1934	i++;
1935	k= d*i;
1936	nmod_poly_clear (buf);
1937	}
1938
1939	return result;
1940	}
1941
1942	CanonicalForm
1943	mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1944	CanonicalForm& M)
1945	{
1946	int d1= degree (F, 1) + degree (G, 1) + 1;
1947	d1 /= 2;
1948	d1 += 1;
1949
1950	nmod_poly_t F1, F2;
1951	kronSubReciproFp (F1, F2, F, d1);
1952
1953	nmod_poly_t G1, G2;
1954	kronSubReciproFp (G1, G2, G, d1);
1955
1956	int k= d1*degree (M);
1957	nmod_poly_mullow (F1, F1, G1, (long) k);
1958
1959	int degtailF= degree (tailcoeff (F), 1);;
1960	int degtailG= degree (tailcoeff (G), 1);
1961	int taildegF= taildegree (F);
1962	int taildegG= taildegree (G);
1963
1964	int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
1965	+ d1*(2+taildegF + taildegG);
1966	nmod_poly_mulhigh (F2, F2, G2, b);
1967	nmod_poly_shift_right (F2, F2, b);
1968	int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
1969
1970
1971	CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
1972
1973	nmod_poly_clear (F1);
1974	nmod_poly_clear (F2);
1975	nmod_poly_clear (G1);
1976	nmod_poly_clear (G2);
1977	return result;
1978	}
1979
1980	CanonicalForm
1981	mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1982	CanonicalForm& M)
1983	{
1984	CanonicalForm A= F;
1985	CanonicalForm B= G;
1986
1987	int degAx= degree (A, 1);
1988	int degAy= degree (A, 2);
1989	int degBx= degree (B, 1);
1990	int degBy= degree (B, 2);
1991	int d1= degAx + 1 + degBx;
1992	int d2= tmax (degAy, degBy);
1993
1994	if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1995	return mulMod2FLINTFpReci (A, B, M);
1996
1997	nmod_poly_t FLINTA, FLINTB;
1998	kronSubFp (FLINTA, A, d1);
1999	kronSubFp (FLINTB, B, d1);
2000
2001	int k= d1*degree (M);
2002	nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2003
2004	A= reverseSubstFp (FLINTA, d1);
2005
2006	nmod_poly_clear (FLINTA);
2007	nmod_poly_clear (FLINTB);
2008	return A;
2009	}
2010
2011	#if (__FLINT_VERSION_MINOR >= 4)
2012	CanonicalForm
2013	mulMod2FLINTFqReci (const CanonicalForm& F, const CanonicalForm& G, const
2014	CanonicalForm& M, const Variable& alpha,
2015	const fq_nmod_ctx_t fq_con)
2016	{
2017	int d1= degree (F, 1) + degree (G, 1) + 1;
2018	d1 /= 2;
2019	d1 += 1;
2020
2021	fq_nmod_poly_t F1, F2;
2022	kronSubReciproFq (F1, F2, F, d1, fq_con);
2023
2024	fq_nmod_poly_t G1, G2;
2025	kronSubReciproFq (G1, G2, G, d1, fq_con);
2026
2027	int k= d1*degree (M);
2028	fq_nmod_poly_mullow (F1, F1, G1, (long) k, fq_con);
2029
2030	int degtailF= degree (tailcoeff (F), 1);
2031	int degtailG= degree (tailcoeff (G), 1);
2032	int taildegF= taildegree (F);
2033	int taildegG= taildegree (G);
2034
2035	int b= k + degtailF + degtailG - d1*(2+taildegF + taildegG);
2036
2037	fq_nmod_poly_reverse (F2, F2, fq_nmod_poly_length (F2, fq_con), fq_con);
2038	fq_nmod_poly_reverse (G2, G2, fq_nmod_poly_length (G2, fq_con), fq_con);
2039	fq_nmod_poly_mullow (F2, F2, G2, b+1, fq_con);
2040	fq_nmod_poly_reverse (F2, F2, b+1, fq_con);
2041
2042	int d2= tmax (fq_nmod_poly_degree (F2, fq_con)/d1,
2043	fq_nmod_poly_degree (F1, fq_con)/d1);
2044
2045	CanonicalForm result= reverseSubstReciproFq (F1, F2, d1, d2, alpha, fq_con);
2046
2047	fq_nmod_poly_clear (F1, fq_con);
2048	fq_nmod_poly_clear (F2, fq_con);
2049	fq_nmod_poly_clear (G1, fq_con);
2050	fq_nmod_poly_clear (G2, fq_con);
2051	return result;
2052	}
2053
2054	CanonicalForm
2055	mulMod2FLINTFq (const CanonicalForm& F, const CanonicalForm& G, const
2056	CanonicalForm& M, const Variable& alpha,
2057	const fq_nmod_ctx_t fq_con)
2058	{
2059	CanonicalForm A= F;
2060	CanonicalForm B= G;
2061
2062	int degAx= degree (A, 1);
2063	int degAy= degree (A, 2);
2064	int degBx= degree (B, 1);
2065	int degBy= degree (B, 2);
2066	int d1= degAx + 1 + degBx;
2067	int d2= tmax (degAy, degBy);
2068
2069	if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2070	return mulMod2FLINTFqReci (A, B, M, alpha, fq_con);
2071
2072	fq_nmod_poly_t FLINTA, FLINTB;
2073	kronSubFq (FLINTA, A, d1, fq_con);
2074	kronSubFq (FLINTB, B, d1, fq_con);
2075
2076	int k= d1*degree (M);
2077	fq_nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k, fq_con);
2078
2079	A= reverseSubstFq (FLINTA, d1, alpha, fq_con);
2080
2081	fq_nmod_poly_clear (FLINTA, fq_con);
2082	fq_nmod_poly_clear (FLINTB, fq_con);
2083	return A;
2084	}
2085	#endif
2086
2087	CanonicalForm
2088	mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
2089	CanonicalForm& M)
2090	{
2091	int d1= degree (F, 1) + degree (G, 1) + 1;
2092	d1 /= 2;
2093	d1 += 1;
2094
2095	fmpz_poly_t F1, F2;
2096	kronSubReciproQ (F1, F2, F, d1);
2097
2098	fmpz_poly_t G1, G2;
2099	kronSubReciproQ (G1, G2, G, d1);
2100
2101	int k= d1*degree (M);
2102	fmpz_poly_mullow (F1, F1, G1, (long) k);
2103
2104	int degtailF= degree (tailcoeff (F), 1);;
2105	int degtailG= degree (tailcoeff (G), 1);
2106	int taildegF= taildegree (F);
2107	int taildegG= taildegree (G);
2108
2109	int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
2110	+ d1*(2+taildegF + taildegG);
2111	fmpz_poly_mulhigh_n (F2, F2, G2, b);
2112	fmpz_poly_shift_right (F2, F2, b);
2113	int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
2114
2115	CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
2116
2117	fmpz_poly_clear (F1);
2118	fmpz_poly_clear (F2);
2119	fmpz_poly_clear (G1);
2120	fmpz_poly_clear (G2);
2121	return result;
2122	}
2123
2124	CanonicalForm
2125	mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
2126	CanonicalForm& M)
2127	{
2128	CanonicalForm A= F;
2129	CanonicalForm B= G;
2130
2131	int degAx= degree (A, 1);
2132	int degBx= degree (B, 1);
2133	int d1= degAx + 1 + degBx;
2134
2135	CanonicalForm f= bCommonDen (F);
2136	CanonicalForm g= bCommonDen (G);
2137	A *= f;
2138	B *= g;
2139
2140	fmpz_poly_t FLINTA, FLINTB;
2141	kronSub (FLINTA, A, d1);
2142	kronSub (FLINTB, B, d1);
2143	int k= d1*degree (M);
2144
2145	fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2146	A= reverseSubstQ (FLINTA, d1);
2147	fmpz_poly_clear (FLINTA);
2148	fmpz_poly_clear (FLINTB);
2149	return A/(f*g);
2150	}
2151
2152	/*CanonicalForm
2153	mulMod2FLINTQa (const CanonicalForm& F, const CanonicalForm& G,
2154	const CanonicalForm& M)
2155	{
2156	Variable a;
2157	if (!hasFirstAlgVar (F,a) && !hasFirstAlgVar (G, a))
2158	return mulMod2FLINTQ (F, G, M);
2159	CanonicalForm A= F;
2160
2161	int degFx= degree (F, 1);
2162	int degFa= degree (F, a);
2163	int degGx= degree (G, 1);
2164	int degGa= degree (G, a);
2165
2166	int d2= degFa+degGa+1;
2167	int d1= degFx + 1 + degGx;
2168	d1 *= d2;
2169
2170	fmpq_poly_t FLINTF, FLINTG;
2171	kronSubQa (FLINTF, F, d1, d2);
2172	kronSubQa (FLINTG, G, d1, d2);
2173
2174	fmpq_poly_mullow (FLINTF, FLINTF, FLINTG, d1*degree (M));
2175
2176	fmpq_poly_t mipo;
2177	convertFacCF2Fmpq_poly_t (mipo, getMipo (a));
2178	CanonicalForm result= reverseSubstQa (FLINTF, d1, d2, a, mipo);
2179	fmpq_poly_clear (FLINTF);
2180	fmpq_poly_clear (FLINTG);
2181	return result;
2182	}*/
2183
2184	CanonicalForm
2185	mulMod2FLINTQa (const CanonicalForm& F, const CanonicalForm& G,
2186	const CanonicalForm& M)
2187	{
2188	Variable a;
2189	if (!hasFirstAlgVar (F,a) && !hasFirstAlgVar (G, a))
2190	return mulMod2FLINTQ (F, G, M);
2191	CanonicalForm A= F, B= G;
2192
2193	int degFx= degree (F, 1);
2194	int degFa= degree (F, a);
2195	int degGx= degree (G, 1);
2196	int degGa= degree (G, a);
2197
2198	int d2= degFa+degGa+1;
2199	int d1= degFx + 1 + degGx;
2200	d1 *= d2;
2201
2202	CanonicalForm f= bCommonDen (F);
2203	CanonicalForm g= bCommonDen (G);
2204	A *= f;
2205	B *= g;
2206
2207	fmpz_poly_t FLINTF, FLINTG;
2208	kronSubQa (FLINTF, A, d1, d2);
2209	kronSubQa (FLINTG, B, d1, d2);
2210
2211	fmpz_poly_mullow (FLINTF, FLINTF, FLINTG, d1*degree (M));
2212
2213	fmpq_poly_t mipo;
2214	convertFacCF2Fmpq_poly_t (mipo, getMipo (a));
2215	A= reverseSubstQa (FLINTF, d1, d2, a, mipo);
2216	fmpz_poly_clear (FLINTF);
2217	fmpz_poly_clear (FLINTG);
2218	return A/(f*g);
2219	}
2220
2221	#endif
2222
2223	#ifndef HAVE_FLINT
2224	zz_pX kronSubFp (const CanonicalForm& A, int d)
2225	{
2226	int degAy= degree (A);
2227	zz_pX result;
2228	result.rep.SetLength (d*(degAy + 1));
2229
2230	zz_p *resultp;
2231	resultp= result.rep.elts();
2232	zz_pX buf;
2233	zz_p *bufp;
2234	int j, k, bufRepLength;
2235
2236	for (CFIterator i= A; i.hasTerms(); i++)
2237	{
2238	if (i.coeff().inCoeffDomain())
2239	buf= convertFacCF2NTLzzpX (i.coeff());
2240	else
2241	buf= convertFacCF2NTLzzpX (i.coeff());
2242
2243	k= i.exp()*d;
2244	bufp= buf.rep.elts();
2245	bufRepLength= (int) buf.rep.length();
2246	for (j= 0; j < bufRepLength; j++)
2247	resultp [j + k]= bufp [j];
2248	}
2249	result.normalize();
2250
2251	return result;
2252	}
2253	#endif
2254
2255	#if (!(HAVE_FLINT && __FLINT_VERSION_MINOR >= 4))
2256	zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
2257	{
2258	int degAy= degree (A);
2259	zz_pEX result;
2260	result.rep.SetLength (d*(degAy + 1));
2261
2262	Variable v;
2263	zz_pE *resultp;
2264	resultp= result.rep.elts();
2265	zz_pEX buf1;
2266	zz_pE *buf1p;
2267	zz_pX buf2;
2268	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2269	int j, k, buf1RepLength;
2270
2271	for (CFIterator i= A; i.hasTerms(); i++)
2272	{
2273	if (i.coeff().inCoeffDomain())
2274	{
2275	buf2= convertFacCF2NTLzzpX (i.coeff());
2276	buf1= to_zz_pEX (to_zz_pE (buf2));
2277	}
2278	else
2279	buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
2280
2281	k= i.exp()*d;
2282	buf1p= buf1.rep.elts();
2283	buf1RepLength= (int) buf1.rep.length();
2284	for (j= 0; j < buf1RepLength; j++)
2285	resultp [j + k]= buf1p [j];
2286	}
2287	result.normalize();
2288
2289	return result;
2290	}
2291
2292	void
2293	kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
2294	const Variable& alpha)
2295	{
2296	int degAy= degree (A);
2297	subA1.rep.SetLength ((long) d*(degAy + 2));
2298	subA2.rep.SetLength ((long) d*(degAy + 2));
2299
2300	Variable v;
2301	zz_pE *subA1p;
2302	zz_pE *subA2p;
2303	subA1p= subA1.rep.elts();
2304	subA2p= subA2.rep.elts();
2305	zz_pEX buf;
2306	zz_pE *bufp;
2307	zz_pX buf2;
2308	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2309	int j, k, kk, bufRepLength;
2310
2311	for (CFIterator i= A; i.hasTerms(); i++)
2312	{
2313	if (i.coeff().inCoeffDomain())
2314	{
2315	buf2= convertFacCF2NTLzzpX (i.coeff());
2316	buf= to_zz_pEX (to_zz_pE (buf2));
2317	}
2318	else
2319	buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
2320
2321	k= i.exp()*d;
2322	kk= (degAy - i.exp())*d;
2323	bufp= buf.rep.elts();
2324	bufRepLength= (int) buf.rep.length();
2325	for (j= 0; j < bufRepLength; j++)
2326	{
2327	subA1p [j + k] += bufp [j];
2328	subA2p [j + kk] += bufp [j];
2329	}
2330	}
2331	subA1.normalize();
2332	subA2.normalize();
2333	}
2334	#endif
2335
2336	#ifndef HAVE_FLINT
2337	void
2338	kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
2339	{
2340	int degAy= degree (A);
2341	subA1.rep.SetLength ((long) d*(degAy + 2));
2342	subA2.rep.SetLength ((long) d*(degAy + 2));
2343
2344	zz_p *subA1p;
2345	zz_p *subA2p;
2346	subA1p= subA1.rep.elts();
2347	subA2p= subA2.rep.elts();
2348	zz_pX buf;
2349	zz_p *bufp;
2350	int j, k, kk, bufRepLength;
2351
2352	for (CFIterator i= A; i.hasTerms(); i++)
2353	{
2354	buf= convertFacCF2NTLzzpX (i.coeff());
2355
2356	k= i.exp()*d;
2357	kk= (degAy - i.exp())*d;
2358	bufp= buf.rep.elts();
2359	bufRepLength= (int) buf.rep.length();
2360	for (j= 0; j < bufRepLength; j++)
2361	{
2362	subA1p [j + k] += bufp [j];
2363	subA2p [j + kk] += bufp [j];
2364	}
2365	}
2366	subA1.normalize();
2367	subA2.normalize();
2368	}
2369	#endif
2370
2371	#if (!(HAVE_FLINT && __FLINT_VERSION_MINOR >= 4))
2372	CanonicalForm
2373	reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
2374	const Variable& alpha)
2375	{
2376	Variable y= Variable (2);
2377	Variable x= Variable (1);
2378
2379	zz_pEX f= F;
2380	zz_pEX g= G;
2381	int degf= deg(f);
2382	int degg= deg(g);
2383
2384	zz_pEX buf1;
2385	zz_pEX buf2;
2386	zz_pEX buf3;
2387
2388	zz_pE *buf1p;
2389	zz_pE *buf2p;
2390	zz_pE *buf3p;
2391	if (f.rep.length() < (long) d*(k+1)) //zero padding
2392	f.rep.SetLength ((long)d*(k+1));
2393
2394	zz_pE *gp= g.rep.elts();
2395	zz_pE *fp= f.rep.elts();
2396	CanonicalForm result= 0;
2397	int i= 0;
2398	int lf= 0;
2399	int lg= d*k;
2400	int degfSubLf= degf;
2401	int deggSubLg= degg-lg;
2402	int repLengthBuf2, repLengthBuf1, ind, tmp;
2403	zz_pE zzpEZero= zz_pE();
2404
2405	while (degf >= lf \|\| lg >= 0)
2406	{
2407	if (degfSubLf >= d)
2408	repLengthBuf1= d;
2409	else if (degfSubLf < 0)
2410	repLengthBuf1= 0;
2411	else
2412	repLengthBuf1= degfSubLf + 1;
2413	buf1.rep.SetLength((long) repLengthBuf1);
2414
2415	buf1p= buf1.rep.elts();
2416	for (ind= 0; ind < repLengthBuf1; ind++)
2417	buf1p [ind]= fp [ind + lf];
2418	buf1.normalize();
2419
2420	repLengthBuf1= buf1.rep.length();
2421
2422	if (deggSubLg >= d - 1)
2423	repLengthBuf2= d - 1;
2424	else if (deggSubLg < 0)
2425	repLengthBuf2= 0;
2426	else
2427	repLengthBuf2= deggSubLg + 1;
2428
2429	buf2.rep.SetLength ((long) repLengthBuf2);
2430	buf2p= buf2.rep.elts();
2431	for (ind= 0; ind < repLengthBuf2; ind++)
2432	buf2p [ind]= gp [ind + lg];
2433	buf2.normalize();
2434
2435	repLengthBuf2= buf2.rep.length();
2436
2437	buf3.rep.SetLength((long) repLengthBuf2 + d);
2438	buf3p= buf3.rep.elts();
2439	buf2p= buf2.rep.elts();
2440	buf1p= buf1.rep.elts();
2441	for (ind= 0; ind < repLengthBuf1; ind++)
2442	buf3p [ind]= buf1p [ind];
2443	for (ind= repLengthBuf1; ind < d; ind++)
2444	buf3p [ind]= zzpEZero;
2445	for (ind= 0; ind < repLengthBuf2; ind++)
2446	buf3p [ind + d]= buf2p [ind];
2447	buf3.normalize();
2448
2449	result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
2450	i++;
2451
2452
2453	lf= i*d;
2454	degfSubLf= degf - lf;
2455
2456	lg= d*(k-i);
2457	deggSubLg= degg - lg;
2458
2459	buf1p= buf1.rep.elts();
2460
2461	if (lg >= 0 && deggSubLg > 0)
2462	{
2463	if (repLengthBuf2 > degfSubLf + 1)
2464	degfSubLf= repLengthBuf2 - 1;
2465	tmp= tmin (repLengthBuf1, deggSubLg + 1);
2466	for (ind= 0; ind < tmp; ind++)
2467	gp [ind + lg] -= buf1p [ind];
2468	}
2469
2470	if (lg < 0)
2471	break;
2472
2473	buf2p= buf2.rep.elts();
2474	if (degfSubLf >= 0)
2475	{
2476	for (ind= 0; ind < repLengthBuf2; ind++)
2477	fp [ind + lf] -= buf2p [ind];
2478	}
2479	}
2480
2481	return result;
2482	}
2483	#endif
2484
2485	#ifndef HAVE_FLINT
2486	CanonicalForm
2487	reverseSubstReciproFp (const zz_pX& F, const zz_pX& G, int d, int k)
2488	{
2489	Variable y= Variable (2);
2490	Variable x= Variable (1);
2491
2492	zz_pX f= F;
2493	zz_pX g= G;
2494	int degf= deg(f);
2495	int degg= deg(g);
2496
2497	zz_pX buf1;
2498	zz_pX buf2;
2499	zz_pX buf3;
2500
2501	zz_p *buf1p;
2502	zz_p *buf2p;
2503	zz_p *buf3p;
2504
2505	if (f.rep.length() < (long) d*(k+1)) //zero padding
2506	f.rep.SetLength ((long)d*(k+1));
2507
2508	zz_p *gp= g.rep.elts();
2509	zz_p *fp= f.rep.elts();
2510	CanonicalForm result= 0;
2511	int i= 0;
2512	int lf= 0;
2513	int lg= d*k;
2514	int degfSubLf= degf;
2515	int deggSubLg= degg-lg;
2516	int repLengthBuf2, repLengthBuf1, ind, tmp;
2517	zz_p zzpZero= zz_p();
2518	while (degf >= lf \|\| lg >= 0)
2519	{
2520	if (degfSubLf >= d)
2521	repLengthBuf1= d;
2522	else if (degfSubLf < 0)
2523	repLengthBuf1= 0;
2524	else
2525	repLengthBuf1= degfSubLf + 1;
2526	buf1.rep.SetLength((long) repLengthBuf1);
2527
2528	buf1p= buf1.rep.elts();
2529	for (ind= 0; ind < repLengthBuf1; ind++)
2530	buf1p [ind]= fp [ind + lf];
2531	buf1.normalize();
2532
2533	repLengthBuf1= buf1.rep.length();
2534
2535	if (deggSubLg >= d - 1)
2536	repLengthBuf2= d - 1;
2537	else if (deggSubLg < 0)
2538	repLengthBuf2= 0;
2539	else
2540	repLengthBuf2= deggSubLg + 1;
2541
2542	buf2.rep.SetLength ((long) repLengthBuf2);
2543	buf2p= buf2.rep.elts();
2544	for (ind= 0; ind < repLengthBuf2; ind++)
2545	buf2p [ind]= gp [ind + lg];
2546
2547	buf2.normalize();
2548
2549	repLengthBuf2= buf2.rep.length();
2550
2551
2552	buf3.rep.SetLength((long) repLengthBuf2 + d);
2553	buf3p= buf3.rep.elts();
2554	buf2p= buf2.rep.elts();
2555	buf1p= buf1.rep.elts();
2556	for (ind= 0; ind < repLengthBuf1; ind++)
2557	buf3p [ind]= buf1p [ind];
2558	for (ind= repLengthBuf1; ind < d; ind++)
2559	buf3p [ind]= zzpZero;
2560	for (ind= 0; ind < repLengthBuf2; ind++)
2561	buf3p [ind + d]= buf2p [ind];
2562	buf3.normalize();
2563
2564	result += convertNTLzzpX2CF (buf3, x)*power (y, i);
2565	i++;
2566
2567
2568	lf= i*d;
2569	degfSubLf= degf - lf;
2570
2571	lg= d*(k-i);
2572	deggSubLg= degg - lg;
2573
2574	buf1p= buf1.rep.elts();
2575
2576	if (lg >= 0 && deggSubLg > 0)
2577	{
2578	if (repLengthBuf2 > degfSubLf + 1)
2579	degfSubLf= repLengthBuf2 - 1;
2580	tmp= tmin (repLengthBuf1, deggSubLg + 1);
2581	for (ind= 0; ind < tmp; ind++)
2582	gp [ind + lg] -= buf1p [ind];
2583	}
2584	if (lg < 0)
2585	break;
2586
2587	buf2p= buf2.rep.elts();
2588	if (degfSubLf >= 0)
2589	{
2590	for (ind= 0; ind < repLengthBuf2; ind++)
2591	fp [ind + lf] -= buf2p [ind];
2592	}
2593	}
2594
2595	return result;
2596	}
2597	#endif
2598
2599	#if (!(HAVE_FLINT && __FLINT_VERSION_MINOR >= 4))
2600	CanonicalForm reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
2601	{
2602	Variable y= Variable (2);
2603	Variable x= Variable (1);
2604
2605	zz_pEX f= F;
2606	zz_pE *fp= f.rep.elts();
2607
2608	zz_pEX buf;
2609	zz_pE *bufp;
2610	CanonicalForm result= 0;
2611	int i= 0;
2612	int degf= deg(f);
2613	int k= 0;
2614	int degfSubK, repLength, j;
2615	while (degf >= k)
2616	{
2617	degfSubK= degf - k;
2618	if (degfSubK >= d)
2619	repLength= d;
2620	else
2621	repLength= degfSubK + 1;
2622
2623	buf.rep.SetLength ((long) repLength);
2624	bufp= buf.rep.elts();
2625	for (j= 0; j < repLength; j++)
2626	bufp [j]= fp [j + k];
2627	buf.normalize();
2628
2629	result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
2630	i++;
2631	k= d*i;
2632	}
2633
2634	return result;
2635	}
2636	#endif
2637
2638	#ifndef HAVE_FLINT
2639	CanonicalForm reverseSubstFp (const zz_pX& F, int d)
2640	{
2641	Variable y= Variable (2);
2642	Variable x= Variable (1);
2643
2644	zz_pX f= F;
2645	zz_p *fp= f.rep.elts();
2646
2647	zz_pX buf;
2648	zz_p *bufp;
2649	CanonicalForm result= 0;
2650	int i= 0;
2651	int degf= deg(f);
2652	int k= 0;
2653	int degfSubK, repLength, j;
2654	while (degf >= k)
2655	{
2656	degfSubK= degf - k;
2657	if (degfSubK >= d)
2658	repLength= d;
2659	else
2660	repLength= degfSubK + 1;
2661
2662	buf.rep.SetLength ((long) repLength);
2663	bufp= buf.rep.elts();
2664	for (j= 0; j < repLength; j++)
2665	bufp [j]= fp [j + k];
2666	buf.normalize();
2667
2668	result += convertNTLzzpX2CF (buf, x)*power (y, i);
2669	i++;
2670	k= d*i;
2671	}
2672
2673	return result;
2674	}
2675
2676	// assumes input to be reduced mod M and to be an element of Fp
2677	CanonicalForm
2678	mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
2679	CanonicalForm& M)
2680	{
2681	int d1= degree (F, 1) + degree (G, 1) + 1;
2682	d1 /= 2;
2683	d1 += 1;
2684
2685	zz_pX F1, F2;
2686	kronSubReciproFp (F1, F2, F, d1);
2687	zz_pX G1, G2;
2688	kronSubReciproFp (G1, G2, G, d1);
2689
2690	int k= d1*degree (M);
2691	MulTrunc (F1, F1, G1, (long) k);
2692
2693	int degtailF= degree (tailcoeff (F), 1);
2694	int degtailG= degree (tailcoeff (G), 1);
2695	int taildegF= taildegree (F);
2696	int taildegG= taildegree (G);
2697	int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
2698
2699	reverse (F2, F2);
2700	reverse (G2, G2);
2701	MulTrunc (F2, F2, G2, b + 1);
2702	reverse (F2, F2, b);
2703
2704	int d2= tmax (deg (F2)/d1, deg (F1)/d1);
2705	return reverseSubstReciproFp (F1, F2, d1, d2);
2706	}
2707
2708	//Kronecker substitution
2709	CanonicalForm
2710	mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
2711	CanonicalForm& M)
2712	{
2713	CanonicalForm A= F;
2714	CanonicalForm B= G;
2715
2716	int degAx= degree (A, 1);
2717	int degAy= degree (A, 2);
2718	int degBx= degree (B, 1);
2719	int degBy= degree (B, 2);
2720	int d1= degAx + 1 + degBx;
2721	int d2= tmax (degAy, degBy);
2722
2723	if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2724	return mulMod2NTLFpReci (A, B, M);
2725
2726	zz_pX NTLA= kronSubFp (A, d1);
2727	zz_pX NTLB= kronSubFp (B, d1);
2728
2729	int k= d1*degree (M);
2730	MulTrunc (NTLA, NTLA, NTLB, (long) k);
2731
2732	A= reverseSubstFp (NTLA, d1);
2733
2734	return A;
2735	}
2736	#endif
2737
2738	#if (!(HAVE_FLINT && __FLINT_VERSION_MINOR >= 4))
2739	// assumes input to be reduced mod M and to be an element of Fq not Fp
2740	CanonicalForm
2741	mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
2742	CanonicalForm& M, const Variable& alpha)
2743	{
2744	int d1= degree (F, 1) + degree (G, 1) + 1;
2745	d1 /= 2;
2746	d1 += 1;
2747
2748	zz_pEX F1, F2;
2749	kronSubReciproFq (F1, F2, F, d1, alpha);
2750	zz_pEX G1, G2;
2751	kronSubReciproFq (G1, G2, G, d1, alpha);
2752
2753	int k= d1*degree (M);
2754	MulTrunc (F1, F1, G1, (long) k);
2755
2756	int degtailF= degree (tailcoeff (F), 1);
2757	int degtailG= degree (tailcoeff (G), 1);
2758	int taildegF= taildegree (F);
2759	int taildegG= taildegree (G);
2760	int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
2761
2762	reverse (F2, F2);
2763	reverse (G2, G2);
2764	MulTrunc (F2, F2, G2, b + 1);
2765	reverse (F2, F2, b);
2766
2767	int d2= tmax (deg (F2)/d1, deg (F1)/d1);
2768	return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
2769	}
2770	#endif
2771
2772	#ifdef HAVE_FLINT
2773	CanonicalForm
2774	mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
2775	CanonicalForm& M);
2776	#endif
2777
2778	CanonicalForm
2779	mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
2780	CanonicalForm& M)
2781	{
2782	Variable alpha;
2783	CanonicalForm A= F;
2784	CanonicalForm B= G;
2785
2786	if (hasFirstAlgVar (A, alpha) \|\| hasFirstAlgVar (B, alpha))
2787	{
2788	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
2789	nmod_poly_t FLINTmipo;
2790	convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
2791
2792	fq_nmod_ctx_t fq_con;
2793	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
2794
2795	A= mulMod2FLINTFq (A, B, M, alpha, fq_con);
2796	nmod_poly_clear (FLINTmipo);
2797	fq_nmod_ctx_clear (fq_con);
2798	#else
2799	int degAx= degree (A, 1);
2800	int degAy= degree (A, 2);
2801	int degBx= degree (B, 1);
2802	int degBy= degree (B, 2);
2803	int d1= degAx + degBx + 1;
2804	int d2= tmax (degAy, degBy);
2805	if (fac_NTL_char != getCharacteristic())
2806	{
2807	fac_NTL_char= getCharacteristic();
2808	zz_p::init (getCharacteristic());
2809	}
2810	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2811	zz_pE::init (NTLMipo);
2812
2813	int degMipo= degree (getMipo (alpha));
2814	if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
2815	(2*degAy > degree (M)))
2816	return mulMod2NTLFqReci (A, B, M, alpha);
2817
2818	zz_pEX NTLA= kronSubFq (A, d1, alpha);
2819	zz_pEX NTLB= kronSubFq (B, d1, alpha);
2820
2821	int k= d1*degree (M);
2822
2823	MulTrunc (NTLA, NTLA, NTLB, (long) k);
2824
2825	A= reverseSubstFq (NTLA, d1, alpha);
2826	#endif
2827	}
2828	else
2829	{
2830	#ifdef HAVE_FLINT
2831	A= mulMod2FLINTFp (A, B, M);
2832	#else
2833	A= mulMod2NTLFp (A, B, M);
2834	#endif
2835	}
2836	return A;
2837	}
2838
2839	CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
2840	const CanonicalForm& M)
2841	{
2842	if (A.isZero() \|\| B.isZero())
2843	return 0;
2844
2845	ASSERT (M.isUnivariate(), "M must be univariate");
2846
2847	CanonicalForm F= mod (A, M);
2848	CanonicalForm G= mod (B, M);
2849	if (F.inCoeffDomain())
2850	return G*F;
2851	if (G.inCoeffDomain())
2852	return F*G;
2853
2854	Variable y= M.mvar();
2855	int degF= degree (F, y);
2856	int degG= degree (G, y);
2857
2858	if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
2859	(F.level() == G.level()))
2860	{
2861	CanonicalForm result= mulNTL (F, G);
2862	return mod (result, M);
2863	}
2864	else if (degF <= 1 && degG <= 1)
2865	{
2866	CanonicalForm result= F*G;
2867	return mod (result, M);
2868	}
2869
2870	int sizeF= size (F);
2871	int sizeG= size (G);
2872
2873	int fallBackToNaive= 50;
2874	if (sizeF < fallBackToNaive \|\| sizeG < fallBackToNaive)
2875	{
2876	if (sizeF < sizeG)
2877	return mod (G*F, M);
2878	else
2879	return mod (F*G, M);
2880	}
2881
2882	#ifdef HAVE_FLINT
2883	if (getCharacteristic() == 0)
2884	return mulMod2FLINTQa (F, G, M);
2885	#endif
2886
2887	if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
2888	(((degF-degG) < 50 && degF > degG) \|\| ((degG-degF) < 50 && degF <= degG)))
2889	return mulMod2NTLFq (F, G, M);
2890
2891	int m= (int) ceil (degree (M)/2.0);
2892	if (degF >= m \|\| degG >= m)
2893	{
2894	CanonicalForm MLo= power (y, m);
2895	CanonicalForm MHi= power (y, degree (M) - m);
2896	CanonicalForm F0= mod (F, MLo);
2897	CanonicalForm F1= div (F, MLo);
2898	CanonicalForm G0= mod (G, MLo);
2899	CanonicalForm G1= div (G, MLo);
2900	CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
2901	CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
2902	CanonicalForm F0G0= mulMod2 (F0, G0, M);
2903	return F0G0 + MLo*(F0G1 + F1G0);
2904	}
2905	else
2906	{
2907	m= (int) ceil (tmax (degF, degG)/2.0);
2908	CanonicalForm yToM= power (y, m);
2909	CanonicalForm F0= mod (F, yToM);
2910	CanonicalForm F1= div (F, yToM);
2911	CanonicalForm G0= mod (G, yToM);
2912	CanonicalForm G1= div (G, yToM);
2913	CanonicalForm H00= mulMod2 (F0, G0, M);
2914	CanonicalForm H11= mulMod2 (F1, G1, M);
2915	CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
2916	return H11yToMyToM + (H01 - H11 - H00)*yToM + H00;
2917	}
2918	DEBOUTLN (cerr, "fatal end in mulMod2");
2919	}
2920
2921	// end bivariate polys
2922	//**********************
2923	// multivariate polys
2924
2925	CanonicalForm mod (const CanonicalForm& F, const CFList& M)
2926	{
2927	CanonicalForm A= F;
2928	for (CFListIterator i= M; i.hasItem(); i++)
2929	A= mod (A, i.getItem());
2930	return A;
2931	}
2932
2933	CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
2934	const CFList& MOD)
2935	{
2936	if (A.isZero() \|\| B.isZero())
2937	return 0;
2938
2939	if (MOD.length() == 1)
2940	return mulMod2 (A, B, MOD.getLast());
2941
2942	CanonicalForm M= MOD.getLast();
2943	CanonicalForm F= mod (A, M);
2944	CanonicalForm G= mod (B, M);
2945	if (F.inCoeffDomain())
2946	return G*F;
2947	if (G.inCoeffDomain())
2948	return F*G;
2949
2950	int sizeF= size (F);
2951	int sizeG= size (G);
2952
2953	if (sizeF / MOD.length() < 100 \|\| sizeG / MOD.length() < 100)
2954	{
2955	if (sizeF < sizeG)
2956	return mod (G*F, MOD);
2957	else
2958	return mod (F*G, MOD);
2959	}
2960
2961	Variable y= M.mvar();
2962	int degF= degree (F, y);
2963	int degG= degree (G, y);
2964
2965	if ((degF <= 1 && F.level() <= M.level()) &&
2966	(degG <= 1 && G.level() <= M.level()))
2967	{
2968	CFList buf= MOD;
2969	buf.removeLast();
2970	if (degF == 1 && degG == 1)
2971	{
2972	CanonicalForm F0= mod (F, y);
2973	CanonicalForm F1= div (F, y);
2974	CanonicalForm G0= mod (G, y);
2975	CanonicalForm G1= div (G, y);
2976	if (degree (M) > 2)
2977	{
2978	CanonicalForm H00= mulMod (F0, G0, buf);
2979	CanonicalForm H11= mulMod (F1, G1, buf);
2980	CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
2981	return H11yy + (H01 - H00 - H11)*y + H00;
2982	}
2983	else //here degree (M) == 2
2984	{
2985	buf.append (y);
2986	CanonicalForm F0G1= mulMod (F0, G1, buf);
2987	CanonicalForm F1G0= mulMod (F1, G0, buf);
2988	CanonicalForm F0G0= mulMod (F0, G0, MOD);
2989	CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
2990	return result;
2991	}
2992	}
2993	else if (degF == 1 && degG == 0)
2994	return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
2995	else if (degF == 0 && degG == 1)
2996	return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
2997	else
2998	return mulMod (F, G, buf);
2999	}
3000	int m= (int) ceil (degree (M)/2.0);
3001	if (degF >= m \|\| degG >= m)
3002	{
3003	CanonicalForm MLo= power (y, m);
3004	CanonicalForm MHi= power (y, degree (M) - m);
3005	CanonicalForm F0= mod (F, MLo);
3006	CanonicalForm F1= div (F, MLo);
3007	CanonicalForm G0= mod (G, MLo);
3008	CanonicalForm G1= div (G, MLo);
3009	CFList buf= MOD;
3010	buf.removeLast();
3011	buf.append (MHi);
3012	CanonicalForm F0G1= mulMod (F0, G1, buf);
3013	CanonicalForm F1G0= mulMod (F1, G0, buf);
3014	CanonicalForm F0G0= mulMod (F0, G0, MOD);
3015	return F0G0 + MLo*(F0G1 + F1G0);
3016	}
3017	else
3018	{
3019	m= (int) ceil (tmin (degF, degG)/2.0);
3020	CanonicalForm yToM= power (y, m);
3021	CanonicalForm F0= mod (F, yToM);
3022	CanonicalForm F1= div (F, yToM);
3023	CanonicalForm G0= mod (G, yToM);
3024	CanonicalForm G1= div (G, yToM);
3025	CanonicalForm H00= mulMod (F0, G0, MOD);
3026	CanonicalForm H11= mulMod (F1, G1, MOD);
3027	CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
3028	return H11yToMyToM + (H01 - H11 - H00)*yToM + H00;
3029	}
3030	DEBOUTLN (cerr, "fatal end in mulMod");
3031	}
3032
3033	CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
3034	{
3035	if (L.isEmpty())
3036	return 1;
3037	int l= L.length();
3038	if (l == 1)
3039	return mod (L.getFirst(), M);
3040	else if (l == 2) {
3041	CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
3042	return result;
3043	}
3044	else
3045	{
3046	l /= 2;
3047	CFList tmp1, tmp2;
3048	CFListIterator i= L;
3049	CanonicalForm buf1, buf2;
3050	for (int j= 1; j <= l; j++, i++)
3051	tmp1.append (i.getItem());
3052	tmp2= Difference (L, tmp1);
3053	buf1= prodMod (tmp1, M);
3054	buf2= prodMod (tmp2, M);
3055	CanonicalForm result= mulMod2 (buf1, buf2, M);
3056	return result;
3057	}
3058	}
3059
3060	CanonicalForm prodMod (const CFList& L, const CFList& M)
3061	{
3062	if (L.isEmpty())
3063	return 1;
3064	else if (L.length() == 1)
3065	return L.getFirst();
3066	else if (L.length() == 2)
3067	return mulMod (L.getFirst(), L.getLast(), M);
3068	else
3069	{
3070	int l= L.length()/2;
3071	CFListIterator i= L;
3072	CFList tmp1, tmp2;
3073	CanonicalForm buf1, buf2;
3074	for (int j= 1; j <= l; j++, i++)
3075	tmp1.append (i.getItem());
3076	tmp2= Difference (L, tmp1);
3077	buf1= prodMod (tmp1, M);
3078	buf2= prodMod (tmp2, M);
3079	return mulMod (buf1, buf2, M);
3080	}
3081	}
3082
3083	// end multivariate polys
3084	//***************************
3085	// division
3086
3087	CanonicalForm reverse (const CanonicalForm& F, int d)
3088	{
3089	if (d == 0)
3090	return F;
3091	CanonicalForm A= F;
3092	Variable y= Variable (2);
3093	Variable x= Variable (1);
3094	if (degree (A, x) > 0)
3095	{
3096	A= swapvar (A, x, y);
3097	CanonicalForm result= 0;
3098	CFIterator i= A;
3099	while (d - i.exp() < 0)
3100	i++;
3101
3102	for (; i.hasTerms() && (d - i.exp() >= 0); i++)
3103	result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
3104	return result;
3105	}
3106	else
3107	return A*power (x, d);
3108	}
3109
3110	CanonicalForm
3111	newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
3112	{
3113	int l= ilog2(n);
3114
3115	CanonicalForm g= mod (F, M)[0] [0];
3116
3117	ASSERT (!g.isZero(), "expected a unit");
3118
3119	Variable alpha;
3120
3121	if (!g.isOne())
3122	g = 1/g;
3123	Variable x= Variable (1);
3124	CanonicalForm result;
3125	int exp= 0;
3126	if (n & 1)
3127	{
3128	result= g;
3129	exp= 1;
3130	}
3131	CanonicalForm h;
3132
3133	for (int i= 1; i <= l; i++)
3134	{
3135	h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
3136	h= mod (h, power (x, (1 << i)) - 1);
3137	h= div (h, power (x, (1 << (i - 1))));
3138	h= mod (h, M);
3139	g -= power (x, (1 << (i - 1)))*
3140	mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
3141
3142	if (n & (1 << i))
3143	{
3144	if (exp)
3145	{
3146	h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
3147	h= mod (h, power (x, exp + (1 << i)) - 1);
3148	h= div (h, power (x, exp));
3149	h= mod (h, M);
3150	result -= power(x, exp)*mod (mulMod2 (g, h, M),
3151	power (x, (1 << i)));
3152	exp += (1 << i);
3153	}
3154	else
3155	{
3156	exp= (1 << i);
3157	result= g;
3158	}
3159	}
3160	}
3161
3162	return result;
3163	}
3164
3165	CanonicalForm
3166	newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
3167	M)
3168	{
3169	ASSERT (getCharacteristic() > 0, "positive characteristic expected");
3170
3171	CanonicalForm A= mod (F, M);
3172	CanonicalForm B= mod (G, M);
3173
3174	Variable x= Variable (1);
3175	int degA= degree (A, x);
3176	int degB= degree (B, x);
3177	int m= degA - degB;
3178	if (m < 0)
3179	return 0;
3180
3181	Variable v;
3182	CanonicalForm Q;
3183	if (degB < 1 \|\| CFFactory::gettype() == GaloisFieldDomain)
3184	{
3185	CanonicalForm R;
3186	divrem2 (A, B, Q, R, M);
3187	}
3188	else
3189	{
3190	if (hasFirstAlgVar (A, v) \|\| hasFirstAlgVar (B, v))
3191	{
3192	CanonicalForm R= reverse (A, degA);
3193	CanonicalForm revB= reverse (B, degB);
3194	revB= newtonInverse (revB, m + 1, M);
3195	Q= mulMod2 (R, revB, M);
3196	Q= mod (Q, power (x, m + 1));
3197	Q= reverse (Q, m);
3198	}
3199	else
3200	{
3201	Variable y= Variable (2);
3202	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
3203	nmod_poly_t FLINTmipo;
3204	fq_nmod_ctx_t fq_con;
3205
3206	nmod_poly_init (FLINTmipo, getCharacteristic());
3207	convertFacCF2nmod_poly_t (FLINTmipo, M);
3208
3209	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3210
3211
3212	fq_nmod_poly_t FLINTA, FLINTB;
3213	convertFacCF2Fq_nmod_poly_t (FLINTA, swapvar (A, x, y), fq_con);
3214	convertFacCF2Fq_nmod_poly_t (FLINTB, swapvar (B, x, y), fq_con);
3215
3216	fq_nmod_poly_divrem (FLINTA, FLINTB, FLINTA, FLINTB, fq_con);
3217
3218	Q= convertFq_nmod_poly_t2FacCF (FLINTA, x, y, fq_con);
3219
3220	fq_nmod_poly_clear (FLINTA, fq_con);
3221	fq_nmod_poly_clear (FLINTB, fq_con);
3222	nmod_poly_clear (FLINTmipo);
3223	fq_nmod_ctx_clear (fq_con);
3224	#else
3225	bool zz_pEbak= zz_pE::initialized();
3226	zz_pEBak bak;
3227	if (zz_pEbak)
3228	bak.save();
3229	zz_pX mipo= convertFacCF2NTLzzpX (M);
3230	zz_pEX NTLA, NTLB;
3231	NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
3232	NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
3233	div (NTLA, NTLA, NTLB);
3234	Q= convertNTLzz_pEX2CF (NTLA, x, y);
3235	if (zz_pEbak)
3236	bak.restore();
3237	#endif
3238	}
3239	}
3240
3241	return Q;
3242	}
3243
3244	void
3245	newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3246	CanonicalForm& R, const CanonicalForm& M)
3247	{
3248	CanonicalForm A= mod (F, M);
3249	CanonicalForm B= mod (G, M);
3250	Variable x= Variable (1);
3251	int degA= degree (A, x);
3252	int degB= degree (B, x);
3253	int m= degA - degB;
3254
3255	if (m < 0)
3256	{
3257	R= A;
3258	Q= 0;
3259	return;
3260	}
3261
3262	Variable v;
3263	if (degB <= 1 \|\| CFFactory::gettype() == GaloisFieldDomain)
3264	{
3265	divrem2 (A, B, Q, R, M);
3266	}
3267	else
3268	{
3269	if (hasFirstAlgVar (A, v) \|\| hasFirstAlgVar (B, v))
3270	{
3271	R= reverse (A, degA);
3272
3273	CanonicalForm revB= reverse (B, degB);
3274	revB= newtonInverse (revB, m + 1, M);
3275	Q= mulMod2 (R, revB, M);
3276
3277	Q= mod (Q, power (x, m + 1));
3278	Q= reverse (Q, m);
3279
3280	R= A - mulMod2 (Q, B, M);
3281	}
3282	else
3283	{
3284	Variable y= Variable (2);
3285	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
3286	nmod_poly_t FLINTmipo;
3287	fq_nmod_ctx_t fq_con;
3288
3289	nmod_poly_init (FLINTmipo, getCharacteristic());
3290	convertFacCF2nmod_poly_t (FLINTmipo, M);
3291
3292	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3293
3294	fq_nmod_poly_t FLINTA, FLINTB;
3295	convertFacCF2Fq_nmod_poly_t (FLINTA, swapvar (A, x, y), fq_con);
3296	convertFacCF2Fq_nmod_poly_t (FLINTB, swapvar (B, x, y), fq_con);
3297
3298	fq_nmod_poly_divrem (FLINTA, FLINTB, FLINTA, FLINTB, fq_con);
3299
3300	Q= convertFq_nmod_poly_t2FacCF (FLINTA, x, y, fq_con);
3301	R= convertFq_nmod_poly_t2FacCF (FLINTB, x, y, fq_con);
3302
3303	fq_nmod_poly_clear (FLINTA, fq_con);
3304	fq_nmod_poly_clear (FLINTB, fq_con);
3305	nmod_poly_clear (FLINTmipo);
3306	fq_nmod_ctx_clear (fq_con);
3307	#else
3308	zz_pX mipo= convertFacCF2NTLzzpX (M);
3309	zz_pEX NTLA, NTLB;
3310	NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
3311	NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
3312	zz_pEX NTLQ, NTLR;
3313	DivRem (NTLQ, NTLR, NTLA, NTLB);
3314	Q= convertNTLzz_pEX2CF (NTLQ, x, y);
3315	R= convertNTLzz_pEX2CF (NTLR, x, y);
3316	#endif
3317	}
3318	}
3319	}
3320
3321	static inline
3322	CFList split (const CanonicalForm& F, const int m, const Variable& x)
3323	{
3324	CanonicalForm A= F;
3325	CanonicalForm buf= 0;
3326	bool swap= false;
3327	if (degree (A, x) <= 0)
3328	return CFList(A);
3329	else if (x.level() != A.level())
3330	{
3331	swap= true;
3332	A= swapvar (A, x, A.mvar());
3333	}
3334
3335	int j= (int) floor ((double) degree (A)/ m);
3336	CFList result;
3337	CFIterator i= A;
3338	for (; j >= 0; j--)
3339	{
3340	while (i.hasTerms() && i.exp() - j*m >= 0)
3341	{
3342	if (swap)
3343	buf += i.coeff()power (A.mvar(), i.exp() - jm);
3344	else
3345	buf += i.coeff()power (x, i.exp() - jm);
3346	i++;
3347	}
3348	if (swap)
3349	result.append (swapvar (buf, x, F.mvar()));
3350	else
3351	result.append (buf);
3352	buf= 0;
3353	}
3354	return result;
3355	}
3356
3357	static inline
3358	void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3359	CanonicalForm& R, const CFList& M);
3360
3361	static inline
3362	void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3363	CanonicalForm& R, const CFList& M)
3364	{
3365	CanonicalForm A= mod (F, M);
3366	CanonicalForm B= mod (G, M);
3367	Variable x= Variable (1);
3368	int degB= degree (B, x);
3369	int degA= degree (A, x);
3370	if (degA < degB)
3371	{
3372	Q= 0;
3373	R= A;
3374	return;
3375	}
3376	if (degB < 1)
3377	{
3378	divrem (A, B, Q, R);
3379	Q= mod (Q, M);
3380	R= mod (R, M);
3381	return;
3382	}
3383	int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
3384	ASSERT (4m >= degA, "expected degree (F, 1) < 2degree (G, 1)");
3385	CFList splitA= split (A, m, x);
3386	if (splitA.length() == 3)
3387	splitA.insert (0);
3388	if (splitA.length() == 2)
3389	{
3390	splitA.insert (0);
3391	splitA.insert (0);
3392	}
3393	if (splitA.length() == 1)
3394	{
3395	splitA.insert (0);
3396	splitA.insert (0);
3397	splitA.insert (0);
3398	}
3399
3400	CanonicalForm xToM= power (x, m);
3401
3402	CFListIterator i= splitA;
3403	CanonicalForm H= i.getItem();
3404	i++;
3405	H *= xToM;
3406	H += i.getItem();
3407	i++;
3408	H *= xToM;
3409	H += i.getItem();
3410	i++;
3411
3412	divrem32 (H, B, Q, R, M);
3413
3414	CFList splitR= split (R, m, x);
3415	if (splitR.length() == 1)
3416	splitR.insert (0);
3417
3418	H= splitR.getFirst();
3419	H *= xToM;
3420	H += splitR.getLast();
3421	H *= xToM;
3422	H += i.getItem();
3423
3424	CanonicalForm bufQ;
3425	divrem32 (H, B, bufQ, R, M);
3426
3427	Q *= xToM;
3428	Q += bufQ;
3429	return;
3430	}
3431
3432	static inline
3433	void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3434	CanonicalForm& R, const CFList& M)
3435	{
3436	CanonicalForm A= mod (F, M);
3437	CanonicalForm B= mod (G, M);
3438	Variable x= Variable (1);
3439	int degB= degree (B, x);
3440	int degA= degree (A, x);
3441	if (degA < degB)
3442	{
3443	Q= 0;
3444	R= A;
3445	return;
3446	}
3447	if (degB < 1)
3448	{
3449	divrem (A, B, Q, R);
3450	Q= mod (Q, M);
3451	R= mod (R, M);
3452	return;
3453	}
3454	int m= (int) ceil ((double) (degB + 1)/ 2.0);
3455	ASSERT (3m > degA, "expected degree (F, 1) < 3degree (G, 1)");
3456	CFList splitA= split (A, m, x);
3457	CFList splitB= split (B, m, x);
3458
3459	if (splitA.length() == 2)
3460	{
3461	splitA.insert (0);
3462	}
3463	if (splitA.length() == 1)
3464	{
3465	splitA.insert (0);
3466	splitA.insert (0);
3467	}
3468	CanonicalForm xToM= power (x, m);
3469
3470	CanonicalForm H;
3471	CFListIterator i= splitA;
3472	i++;
3473
3474	if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
3475	{
3476	H= splitA.getFirst()*xToM + i.getItem();
3477	divrem21 (H, splitB.getFirst(), Q, R, M);
3478	}
3479	else
3480	{
3481	R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
3482	splitB.getFirst()*xToM;
3483	Q= xToM - 1;
3484	}
3485
3486	H= mulMod (Q, splitB.getLast(), M);
3487
3488	R= R*xToM + splitA.getLast() - H;
3489
3490	while (degree (R, x) >= degB)
3491	{
3492	xToM= power (x, degree (R, x) - degB);
3493	Q += LC (R, x)*xToM;
3494	R -= mulMod (LC (R, x), B, M)*xToM;
3495	Q= mod (Q, M);
3496	R= mod (R, M);
3497	}
3498
3499	return;
3500	}
3501
3502	void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3503	CanonicalForm& R, const CanonicalForm& M)
3504	{
3505	CanonicalForm A= mod (F, M);
3506	CanonicalForm B= mod (G, M);
3507
3508	if (B.inCoeffDomain())
3509	{
3510	divrem (A, B, Q, R);
3511	return;
3512	}
3513	if (A.inCoeffDomain() && !B.inCoeffDomain())
3514	{
3515	Q= 0;
3516	R= A;
3517	return;
3518	}
3519
3520	if (B.level() < A.level())
3521	{
3522	divrem (A, B, Q, R);
3523	return;
3524	}
3525	if (A.level() > B.level())
3526	{
3527	R= A;
3528	Q= 0;
3529	return;
3530	}
3531	if (B.level() == 1 && B.isUnivariate())
3532	{
3533	divrem (A, B, Q, R);
3534	return;
3535	}
3536
3537	Variable x= Variable (1);
3538	int degB= degree (B, x);
3539	if (degB > degree (A, x))
3540	{
3541	Q= 0;
3542	R= A;
3543	return;
3544	}
3545
3546	CFList splitA= split (A, degB, x);
3547
3548	CanonicalForm xToDegB= power (x, degB);
3549	CanonicalForm H, bufQ;
3550	Q= 0;
3551	CFListIterator i= splitA;
3552	H= i.getItem()*xToDegB;
3553	i++;
3554	H += i.getItem();
3555	CFList buf;
3556	while (i.hasItem())
3557	{
3558	buf= CFList (M);
3559	divrem21 (H, B, bufQ, R, buf);
3560	i++;
3561	if (i.hasItem())
3562	H= R*xToDegB + i.getItem();
3563	Q *= xToDegB;
3564	Q += bufQ;
3565	}
3566	return;
3567	}
3568
3569	void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3570	CanonicalForm& R, const CFList& MOD)
3571	{
3572	CanonicalForm A= mod (F, MOD);
3573	CanonicalForm B= mod (G, MOD);
3574	Variable x= Variable (1);
3575	int degB= degree (B, x);
3576	if (degB > degree (A, x))
3577	{
3578	Q= 0;
3579	R= A;
3580	return;
3581	}
3582
3583	if (degB <= 0)
3584	{
3585	divrem (A, B, Q, R);
3586	Q= mod (Q, MOD);
3587	R= mod (R, MOD);
3588	return;
3589	}
3590	CFList splitA= split (A, degB, x);
3591
3592	CanonicalForm xToDegB= power (x, degB);
3593	CanonicalForm H, bufQ;
3594	Q= 0;
3595	CFListIterator i= splitA;
3596	H= i.getItem()*xToDegB;
3597	i++;
3598	H += i.getItem();
3599	while (i.hasItem())
3600	{
3601	divrem21 (H, B, bufQ, R, MOD);
3602	i++;
3603	if (i.hasItem())
3604	H= R*xToDegB + i.getItem();
3605	Q *= xToDegB;
3606	Q += bufQ;
3607	}
3608	return;
3609	}
3610
3611	bool
3612	uniFdivides (const CanonicalForm& A, const CanonicalForm& B)
3613	{
3614	if (B.isZero())
3615	return true;
3616	if (A.isZero())
3617	return false;
3618	if (CFFactory::gettype() == GaloisFieldDomain)
3619	return fdivides (A, B);
3620	int p= getCharacteristic();
3621	if (A.inCoeffDomain() \|\| B.inCoeffDomain())
3622	{
3623	if (A.inCoeffDomain())
3624	return true;
3625	else
3626	return false;
3627	}
3628	if (p > 0)
3629	{
3630	if (fac_NTL_char != p)
3631	{
3632	fac_NTL_char= p;
3633	zz_p::init (p);
3634	}
3635	Variable alpha;
3636	if (hasFirstAlgVar (A, alpha) \|\| hasFirstAlgVar (B, alpha))
3637	{
3638	#if (HAVE_FLINT && __FLINT_VERSION_MINOR >= 4)
3639	nmod_poly_t FLINTmipo;
3640	fq_nmod_ctx_t fq_con;
3641
3642	nmod_poly_init (FLINTmipo, getCharacteristic());
3643	convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
3644
3645	fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3646
3647	fq_nmod_poly_t FLINTA, FLINTB;
3648	convertFacCF2Fq_nmod_poly_t (FLINTA, A, fq_con);
3649	convertFacCF2Fq_nmod_poly_t (FLINTB, B, fq_con);
3650	int result= fq_nmod_poly_divides (FLINTA, FLINTB, FLINTA, fq_con);
3651	fq_nmod_poly_clear (FLINTA, fq_con);
3652	fq_nmod_poly_clear (FLINTB, fq_con);
3653	nmod_poly_clear (FLINTmipo);
3654	fq_nmod_ctx_clear (fq_con);
3655	return result;
3656	#else
3657	zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
3658	zz_pE::init (NTLMipo);
3659	zz_pEX NTLA= convertFacCF2NTLzz_pEX (A, NTLMipo);
3660	zz_pEX NTLB= convertFacCF2NTLzz_pEX (B, NTLMipo);
3661	return divide (NTLB, NTLA);
3662	#endif
3663	}
3664	#ifdef HAVE_FLINT
3665	nmod_poly_t FLINTA, FLINTB;
3666	convertFacCF2nmod_poly_t (FLINTA, A);
3667	convertFacCF2nmod_poly_t (FLINTB, B);
3668	nmod_poly_divrem (FLINTB, FLINTA, FLINTB, FLINTA);
3669	bool result= nmod_poly_is_zero (FLINTA);
3670	nmod_poly_clear (FLINTA);
3671	nmod_poly_clear (FLINTB);
3672	return result;
3673	#else
3674	zz_pX NTLA= convertFacCF2NTLzzpX (A);
3675	zz_pX NTLB= convertFacCF2NTLzzpX (B);
3676	return divide (NTLB, NTLA);
3677	#endif
3678	}
3679	#ifdef HAVE_FLINT
3680	Variable alpha;
3681	bool isRat= isOn (SW_RATIONAL);
3682	if (!isRat)
3683	On (SW_RATIONAL);
3684	if (!hasFirstAlgVar (A, alpha) && !hasFirstAlgVar (B, alpha))
3685	{
3686	fmpq_poly_t FLINTA,FLINTB;
3687	convertFacCF2Fmpq_poly_t (FLINTA, A);
3688	convertFacCF2Fmpq_poly_t (FLINTB, B);
3689	fmpq_poly_rem (FLINTA, FLINTB, FLINTA);
3690	bool result= fmpq_poly_is_zero (FLINTA);
3691	fmpq_poly_clear (FLINTA);
3692	fmpq_poly_clear (FLINTB);
3693	if (!isRat)
3694	Off (SW_RATIONAL);
3695	return result;
3696	}
3697	CanonicalForm Q, R;
3698	Variable x= Variable (1);
3699	Variable y= Variable (2);
3700	newtonDivrem (swapvar (B, y, x), swapvar (A, y, x), Q, R);
3701	if (!isRat)
3702	Off (SW_RATIONAL);
3703	return R.isZero();
3704	#else
3705	bool isRat= isOn (SW_RATIONAL);
3706	if (!isRat)
3707	On (SW_RATIONAL);
3708	bool result= fdivides (A, B);
3709	if (!isRat)
3710	Off (SW_RATIONAL);
3711	return result; //maybe NTL?
3712	#endif
3713	}
3714
3715	// end division
3716
3717	#endif

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

Download in other formats: