#### /cln-1.3.2/src/polynomial/elem/cl_UP_MI.h

#
C Header | 506 lines | 444 code | 35 blank | 27 comment | 96 complexity | cea3179322a68b56733904fe49200059 MD5 | raw file
``````  1// Univariate Polynomials over a ring of modular integers.
2
3#include "cln/GV_modinteger.h"
4#include "cln/modinteger.h"
5#include "cln/exception.h"
6
7namespace cln {
8
9// Assume a ring is a modint ring.
10  inline cl_heap_modint_ring* TheModintRing (const cl_ring& R)
11  { return (cl_heap_modint_ring*) R.heappointer; }
12
13// Normalize a vector: remove leading zero coefficients.
14// The result vector is known to have length len > 0.
15static inline void modint_normalize (cl_heap_modint_ring* R, cl_GV_MI& result, uintL len)
16{
17	if (R->_zerop(result[len-1])) {
18		len--;
19		while (len > 0) {
20			if (!R->_zerop(result[len-1]))
21				break;
22			len--;
23		}
24		var cl_GV_MI newresult = cl_GV_MI(len,R);
25		#if 0
26		for (var sintL i = len-1; i >= 0; i--)
27			newresult[i] = result[i];
28		#else
29		cl_GV_MI::copy_elements(result,0,newresult,0,len);
30		#endif
31		result = newresult;
32	}
33}
34
35static void modint_fprint (cl_heap_univpoly_ring* UPR, std::ostream& stream, const _cl_UP& x)
36{{
37	DeclarePoly(cl_GV_MI,x);
38	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
39	var sintL xlen = x.size();
40	if (xlen == 0)
41		fprint(stream, "0");
42	else {
43		var cl_string varname = get_varname(UPR);
44		for (var sintL i = xlen-1; i >= 0; i--)
45			if (!R->_zerop(x[i])) {
46				if (i < xlen-1)
47					fprint(stream, " + ");
48				fprint(stream, "(");
49				R->_fprint(stream, x[i]);
50				fprint(stream, ")");
51				if (i > 0) {
52					fprint(stream, "*");
53					fprint(stream, varname);
54					if (i != 1) {
55						fprint(stream, "^");
56						fprintdecimal(stream, i);
57					}
58				}
59			}
60	}
61}}
62
63static bool modint_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
64{{
65	DeclarePoly(cl_GV_MI,x);
66	DeclarePoly(cl_GV_MI,y);
67	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
68	var sintL xlen = x.size();
69	var sintL ylen = y.size();
70	if (!(xlen == ylen))
71		return false;
72	for (var sintL i = xlen-1; i >= 0; i--)
73		if (!R->_equal(x[i],y[i]))
74			return false;
75	return true;
76}}
77
78static const _cl_UP modint_zero (cl_heap_univpoly_ring* UPR)
79{
80	return _cl_UP(UPR, cl_null_GV_I);
81}
82
83static bool modint_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
84{
85	unused UPR;
86 {	DeclarePoly(cl_GV_MI,x);
87	var sintL xlen = x.size();
88	if (xlen == 0)
89		return true;
90	else
91		return false;
92}}
93
94static const _cl_UP modint_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
95{{
96	DeclarePoly(cl_GV_MI,x);
97	DeclarePoly(cl_GV_MI,y);
98	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
99	var sintL xlen = x.size();
100	var sintL ylen = y.size();
101	if (xlen == 0)
102		return _cl_UP(UPR, y);
103	if (ylen == 0)
104		return _cl_UP(UPR, x);
105	// Now xlen > 0, ylen > 0.
106	if (xlen > ylen) {
107		var cl_GV_MI result = cl_GV_MI(xlen,R);
108		var sintL i;
109		#if 0
110		for (i = xlen-1; i >= ylen; i--)
111			result[i] = x[i];
112		#else
113		cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
114		#endif
115		for (i = ylen-1; i >= 0; i--)
116			result[i] = R->_plus(x[i],y[i]);
117		return _cl_UP(UPR, result);
118	}
119	if (xlen < ylen) {
120		var cl_GV_MI result = cl_GV_MI(ylen,R);
121		var sintL i;
122		#if 0
123		for (i = ylen-1; i >= xlen; i--)
124			result[i] = y[i];
125		#else
126		cl_GV_MI::copy_elements(y,xlen,result,xlen,ylen-xlen);
127		#endif
128		for (i = xlen-1; i >= 0; i--)
129			result[i] = R->_plus(x[i],y[i]);
130		return _cl_UP(UPR, result);
131	}
132	// Now xlen = ylen > 0. Add and normalize simultaneously.
133	for (var sintL i = xlen-1; i >= 0; i--) {
134		var _cl_MI hicoeff = R->_plus(x[i],y[i]);
135		if (!R->_zerop(hicoeff)) {
136			var cl_GV_MI result = cl_GV_MI(i+1,R);
137			result[i] = hicoeff;
138			for (i-- ; i >= 0; i--)
139				result[i] = R->_plus(x[i],y[i]);
140			return _cl_UP(UPR, result);
141		}
142	}
143	return _cl_UP(UPR, cl_null_GV_I);
144}}
145
146static const _cl_UP modint_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
147{{
148	DeclarePoly(cl_GV_MI,x);
149	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
150	var sintL xlen = x.size();
151	if (xlen == 0)
152		return _cl_UP(UPR, x);
153	// Now xlen > 0.
154	// Negate. No normalization necessary, since the degree doesn't change.
155	var sintL i = xlen-1;
156	var _cl_MI hicoeff = R->_uminus(x[i]);
157	if (R->_zerop(hicoeff)) throw runtime_exception();
158	var cl_GV_MI result = cl_GV_MI(xlen,R);
159	result[i] = hicoeff;
160	for (i-- ; i >= 0; i--)
161		result[i] = R->_uminus(x[i]);
162	return _cl_UP(UPR, result);
163}}
164
165static const _cl_UP modint_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
166{{
167	DeclarePoly(cl_GV_MI,x);
168	DeclarePoly(cl_GV_MI,y);
169	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
170	var sintL xlen = x.size();
171	var sintL ylen = y.size();
172	if (ylen == 0)
173		return _cl_UP(UPR, x);
174	if (xlen == 0)
175		return modint_uminus(UPR, _cl_UP(UPR, y));
176	// Now xlen > 0, ylen > 0.
177	if (xlen > ylen) {
178		var cl_GV_MI result = cl_GV_MI(xlen,R);
179		var sintL i;
180		#if 0
181		for (i = xlen-1; i >= ylen; i--)
182			result[i] = x[i];
183		#else
184		cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
185		#endif
186		for (i = ylen-1; i >= 0; i--)
187			result[i] = R->_minus(x[i],y[i]);
188		return _cl_UP(UPR, result);
189	}
190	if (xlen < ylen) {
191		var cl_GV_MI result = cl_GV_MI(ylen,R);
192		var sintL i;
193		for (i = ylen-1; i >= xlen; i--)
194			result[i] = R->_uminus(y[i]);
195		for (i = xlen-1; i >= 0; i--)
196			result[i] = R->_minus(x[i],y[i]);
197		return _cl_UP(UPR, result);
198	}
199	// Now xlen = ylen > 0. Add and normalize simultaneously.
200	for (var sintL i = xlen-1; i >= 0; i--) {
201		var _cl_MI hicoeff = R->_minus(x[i],y[i]);
202		if (!R->_zerop(hicoeff)) {
203			var cl_GV_MI result = cl_GV_MI(i+1,R);
204			result[i] = hicoeff;
205			for (i-- ; i >= 0; i--)
206				result[i] = R->_minus(x[i],y[i]);
207			return _cl_UP(UPR, result);
208		}
209	}
210	return _cl_UP(UPR, cl_null_GV_I);
211}}
212
213static const _cl_UP modint_one (cl_heap_univpoly_ring* UPR)
214{
215	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
216	var cl_GV_MI result = cl_GV_MI(1,R);
217	result[0] = R->_one();
218	return _cl_UP(UPR, result);
219}
220
221static const _cl_UP modint_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
222{
223	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
224	var cl_GV_MI result = cl_GV_MI(1,R);
225	result[0] = R->_canonhom(x);
226	return _cl_UP(UPR, result);
227}
228
229static const _cl_UP modint_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
230{{
231	DeclarePoly(cl_GV_MI,x);
232	DeclarePoly(cl_GV_MI,y);
233	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
234	var sintL xlen = x.size();
235	var sintL ylen = y.size();
236	if (xlen == 0)
237		return _cl_UP(UPR, x);
238	if (ylen == 0)
239		return _cl_UP(UPR, y);
240	// Multiply.
241	var sintL len = xlen + ylen - 1;
242	var cl_GV_MI result = cl_GV_MI(len,R);
243	if (xlen < ylen) {
244		{
245			var sintL i = xlen-1;
246			var _cl_MI xi = x[i];
247			for (sintL j = ylen-1; j >= 0; j--)
248				result[i+j] = R->_mul(xi,y[j]);
249		}
250		for (sintL i = xlen-2; i >= 0; i--) {
251			var _cl_MI xi = x[i];
252			for (sintL j = ylen-1; j > 0; j--)
253				result[i+j] = R->_plus(result[i+j],R->_mul(xi,y[j]));
254			/* j=0 */ result[i] = R->_mul(xi,y[0]);
255		}
256	} else {
257		{
258			var sintL j = ylen-1;
259			var _cl_MI yj = y[j];
260			for (sintL i = xlen-1; i >= 0; i--)
261				result[i+j] = R->_mul(x[i],yj);
262		}
263		for (sintL j = ylen-2; j >= 0; j--) {
264			var _cl_MI yj = y[j];
265			for (sintL i = xlen-1; i > 0; i--)
266				result[i+j] = R->_plus(result[i+j],R->_mul(x[i],yj));
267			/* i=0 */ result[j] = R->_mul(x[0],yj);
268		}
269	}
270	// Normalize (not necessary in integral domains).
271	//modint_normalize(R,result,len);
272	if (R->_zerop(result[len-1])) throw runtime_exception();
273	return _cl_UP(UPR, result);
274}}
275
276static const _cl_UP modint_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
277{{
278	DeclarePoly(cl_GV_MI,x);
279	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
280	var sintL xlen = x.size();
281	if (xlen == 0)
282		return cl_UP(UPR, x);
283	var sintL len = 2*xlen-1;
284	var cl_GV_MI result = cl_GV_MI(len,R);
285	if (xlen > 1) {
286		// Loop through all 0 <= j < i <= xlen-1.
287		{
288			var sintL i = xlen-1;
289			var _cl_MI xi = x[i];
290			for (sintL j = i-1; j >= 0; j--)
291				result[i+j] = R->_mul(xi,x[j]);
292		}
293		{for (sintL i = xlen-2; i >= 1; i--) {
294			var _cl_MI xi = x[i];
295			for (sintL j = i-1; j >= 1; j--)
296				result[i+j] = R->_plus(result[i+j],R->_mul(xi,x[j]));
297			/* j=0 */ result[i] = R->_mul(xi,x[0]);
298		}}
299		// Double.
300		{for (sintL i = len-2; i >= 1; i--)
301			result[i] = R->_plus(result[i],result[i]);
302		}
304		result[2*(xlen-1)] = R->_square(x[xlen-1]);
305		for (sintL i = xlen-2; i >= 1; i--)
306			result[2*i] = R->_plus(result[2*i],R->_square(x[i]));
307	}
308	result[0] = R->_square(x[0]);
309	// Normalize (not necessary in integral domains).
310	//modint_normalize(R,result,len);
311	if (R->_zerop(result[len-1])) throw runtime_exception();
312	return _cl_UP(UPR, result);
313}}
314
315static const _cl_UP modint_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
316{
317	var _cl_UP a = x;
318	var cl_I b = y;
319	while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
320	var _cl_UP c = a;
321	until (b == 1)
322	  { b = b >> 1;
323	    a = UPR->_square(a);
324	    if (oddp(b)) { c = UPR->_mul(a,c); }
325	  }
326	return c;
327}
328
329static const _cl_UP modint_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
330{
331	if (!(UPR->basering() == x.ring())) throw runtime_exception();
332 {
333	DeclarePoly(_cl_MI,x);
334	DeclarePoly(cl_GV_MI,y);
335	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
336	var sintL ylen = y.size();
337	if (ylen == 0)
338		return _cl_UP(UPR, y);
339	if (R->_zerop(x))
340		return _cl_UP(UPR, cl_null_GV_I);
341	// Now ylen > 0.
342	// No normalization necessary, since the degree doesn't change.
343	var cl_GV_MI result = cl_GV_MI(ylen,R);
344	for (sintL i = ylen-1; i >= 0; i--)
345		result[i] = R->_mul(x,y[i]);
346	return _cl_UP(UPR, result);
347}}
348
349static sintL modint_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
350{
351	unused UPR;
352 {	DeclarePoly(cl_GV_MI,x);
353	return (sintL) x.size() - 1;
354}}
355
356static sintL modint_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
357{{
358	DeclarePoly(cl_GV_MI,x);
359	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
360	var sintL xlen = x.size();
361	for (sintL i = 0; i < xlen; i++) {
362		if (!R->_zerop(x[i]))
363			return i;
364	}
365	return -1;
366}}
367
368static const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
369{
370	if (!(UPR->basering() == x.ring())) throw runtime_exception();
371 {	DeclarePoly(_cl_MI,x);
372	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
373	if (R->_zerop(x))
374		return _cl_UP(UPR, cl_null_GV_I);
375	else {
376		var sintL len = e+1;
377		var cl_GV_MI result = cl_GV_MI(len,R);
378		result[e] = x;
379		return _cl_UP(UPR, result);
380	}
381}}
382
383static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
384{{
385	DeclarePoly(cl_GV_MI,x);
386	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
387	if (index < x.size())
388		return cl_MI(R, x[index]);
389	else
390		return R->zero();
391}}
392
393static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
394{
395	if (deg < 0)
396		return _cl_UP(UPR, cl_null_GV_I);
397	else {
398		var sintL len = deg+1;
399		var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
400		return _cl_UP(UPR, cl_GV_MI(len,R));
401	}
402}
403
404static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
405{{
406	DeclareMutablePoly(cl_GV_MI,x);
407	if (!(UPR->basering() == y.ring())) throw runtime_exception();
408  {	DeclarePoly(_cl_MI,y);
409	if (!(index < x.size())) throw runtime_exception();
410	x[index] = y;
411}}}
412
413static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
414{{
415	DeclareMutablePoly(cl_GV_MI,x); // NB: x is modified by reference!
416	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
417	var uintL len = x.size();
418	if (len > 0)
419		modint_normalize(R,x,len);
420}}
421
422static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
423{{
424	// Method:
425	// If x = 0, return 0.
426	// If y = 0, return x[0].
427	// Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
428	DeclarePoly(cl_GV_MI,x);
429	if (!(UPR->basering() == y.ring())) throw runtime_exception();
430  {	DeclarePoly(_cl_MI,y);
431	var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
432	var uintL len = x.size();
433	if (len==0)
434		return R->zero();
435	if (R->_zerop(y))
436		return cl_MI(R, x[0]);
437	var sintL i = len-1;
438	var _cl_MI z = x[i];
439	for ( ; --i >= 0; )
440		z = R->_plus(R->_mul(z,y),x[i]);
441	return cl_MI(R, z);
442}}}
443
444static cl_univpoly_setops modint_setops = {
445	modint_fprint,
446	modint_equal
447};
448
450	modint_zero,
451	modint_zerop,
452	modint_plus,
453	modint_minus,
454	modint_uminus
455};
456
457static cl_univpoly_mulops modint_mulops = {
458	modint_one,
459	modint_canonhom,
460	modint_mul,
461	modint_square,
462	modint_exptpos
463};
464
465static cl_univpoly_modulops modint_modulops = {
466	modint_scalmul
467};
468
469static cl_univpoly_polyops modint_polyops = {
470	modint_degree,
471	modint_ldegree,
472	modint_monomial,
473	modint_coeff,
474	modint_create,
475	modint_set_coeff,
476	modint_finalize,
477	modint_eval
478};
479
480class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
481	SUBCLASS_cl_heap_univpoly_ring()
482public:
483	// Constructor.
484	cl_heap_modint_univpoly_ring (const cl_ring& r);
485	// Destructor.
486	~cl_heap_modint_univpoly_ring () {}
487};
488
489static void cl_heap_modint_univpoly_ring_destructor (cl_heap* pointer)
490{
491	(*(cl_heap_modint_univpoly_ring*)pointer).~cl_heap_modint_univpoly_ring();
492}
493
494cl_class cl_class_modint_univpoly_ring = {
495	cl_heap_modint_univpoly_ring_destructor,
496	cl_class_flags_univpoly_ring
497};
498
499// Constructor.
500inline cl_heap_modint_univpoly_ring::cl_heap_modint_univpoly_ring (const cl_ring& r)
501	: cl_heap_univpoly_ring (r, &modint_setops, &modint_addops, &modint_mulops, &modint_modulops, &modint_polyops)
502{
503	type = &cl_class_modint_univpoly_ring;
504}
505
506}  // namespace cln
``````