Context Navigation

source: git/factory/facMul.cc @ 0ebf2c

spielwiese

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

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

Download in other formats: