Context Navigation

source: git/factory/facMul.cc @ f63d2e8

spielwiese

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

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

Download in other formats: