1 // Univariate Polynomials over a ring of modular integers.
3 #include "cln/GV_modinteger.h"
4 #include "cln/modinteger.h"
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; }
13 // Normalize a vector: remove leading zero coefficients.
14 // The result vector is known to have length len > 0.
15 static inline void modint_normalize (cl_heap_modint_ring* R, cl_GV_MI& result, uintL len)
17 if (R->_zerop(result[len-1])) {
20 if (!R->_zerop(result[len-1]))
24 var cl_GV_MI newresult = cl_GV_MI(len,R);
26 for (var sintL i = len-1; i >= 0; i--)
27 newresult[i] = result[i];
29 cl_GV_MI::copy_elements(result,0,newresult,0,len);
35 static void modint_fprint (cl_heap_univpoly_ring* UPR, std::ostream& stream, const _cl_UP& x)
37 DeclarePoly(cl_GV_MI,x);
38 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
39 var sintL xlen = x.length();
43 var cl_string varname = get_varname(UPR);
44 for (var sintL i = xlen-1; i >= 0; i--)
45 if (!R->_zerop(x[i])) {
47 fprint(stream, " + ");
49 R->_fprint(stream, x[i]);
53 fprint(stream, varname);
56 fprintdecimal(stream, i);
63 static cl_boolean modint_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
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.length();
69 var sintL ylen = y.length();
72 for (var sintL i = xlen-1; i >= 0; i--)
73 if (!R->_equal(x[i],y[i]))
78 static const _cl_UP modint_zero (cl_heap_univpoly_ring* UPR)
80 return _cl_UP(UPR, cl_null_GV_I);
83 static cl_boolean modint_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
86 { DeclarePoly(cl_GV_MI,x);
87 var sintL xlen = x.length();
94 static const _cl_UP modint_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
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.length();
100 var sintL ylen = y.length();
102 return _cl_UP(UPR, y);
104 return _cl_UP(UPR, x);
105 // Now xlen > 0, ylen > 0.
107 var cl_GV_MI result = cl_GV_MI(xlen,R);
110 for (i = xlen-1; i >= ylen; i--)
113 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
115 for (i = ylen-1; i >= 0; i--)
116 result[i] = R->_plus(x[i],y[i]);
117 return _cl_UP(UPR, result);
120 var cl_GV_MI result = cl_GV_MI(ylen,R);
123 for (i = ylen-1; i >= xlen; i--)
126 cl_GV_MI::copy_elements(y,xlen,result,xlen,ylen-xlen);
128 for (i = xlen-1; i >= 0; i--)
129 result[i] = R->_plus(x[i],y[i]);
130 return _cl_UP(UPR, result);
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);
138 for (i-- ; i >= 0; i--)
139 result[i] = R->_plus(x[i],y[i]);
140 return _cl_UP(UPR, result);
143 return _cl_UP(UPR, cl_null_GV_I);
146 static const _cl_UP modint_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
148 DeclarePoly(cl_GV_MI,x);
149 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
150 var sintL xlen = x.length();
152 return _cl_UP(UPR, x);
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)) cl_abort();
158 var cl_GV_MI result = cl_GV_MI(xlen,R);
160 for (i-- ; i >= 0; i--)
161 result[i] = R->_uminus(x[i]);
162 return _cl_UP(UPR, result);
165 static const _cl_UP modint_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
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.length();
171 var sintL ylen = y.length();
173 return _cl_UP(UPR, x);
175 return modint_uminus(UPR, _cl_UP(UPR, y));
176 // Now xlen > 0, ylen > 0.
178 var cl_GV_MI result = cl_GV_MI(xlen,R);
181 for (i = xlen-1; i >= ylen; i--)
184 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
186 for (i = ylen-1; i >= 0; i--)
187 result[i] = R->_minus(x[i],y[i]);
188 return _cl_UP(UPR, result);
191 var cl_GV_MI result = cl_GV_MI(ylen,R);
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);
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);
205 for (i-- ; i >= 0; i--)
206 result[i] = R->_minus(x[i],y[i]);
207 return _cl_UP(UPR, result);
210 return _cl_UP(UPR, cl_null_GV_I);
213 static const _cl_UP modint_one (cl_heap_univpoly_ring* UPR)
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);
221 static const _cl_UP modint_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
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);
229 static const _cl_UP modint_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
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.length();
235 var sintL ylen = y.length();
237 return _cl_UP(UPR, x);
239 return _cl_UP(UPR, y);
241 var sintL len = xlen + ylen - 1;
242 var cl_GV_MI result = cl_GV_MI(len,R);
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]);
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]);
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);
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);
270 // Normalize (not necessary in integral domains).
271 //modint_normalize(R,result,len);
272 if (R->_zerop(result[len-1])) cl_abort();
273 return _cl_UP(UPR, result);
276 static const _cl_UP modint_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
278 DeclarePoly(cl_GV_MI,x);
279 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
280 var sintL xlen = x.length();
282 return cl_UP(UPR, x);
283 var sintL len = 2*xlen-1;
284 var cl_GV_MI result = cl_GV_MI(len,R);
286 // Loop through all 0 <= j < i <= xlen-1.
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]);
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]);
300 {for (sintL i = len-2; i >= 1; i--)
301 result[i] = R->_plus(result[i],result[i]);
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]));
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])) cl_abort();
312 return _cl_UP(UPR, result);
315 static const _cl_UP modint_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
319 while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
324 if (oddp(b)) { c = UPR->_mul(a,c); }
329 static const _cl_UP modint_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
331 if (!(UPR->basering() == x.ring())) cl_abort();
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.length();
338 return _cl_UP(UPR, y);
340 return _cl_UP(UPR, cl_null_GV_I);
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);
349 static sintL modint_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
352 { DeclarePoly(cl_GV_MI,x);
353 return (sintL) x.length() - 1;
356 static sintL modint_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
358 DeclarePoly(cl_GV_MI,x);
359 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
360 var sintL xlen = x.length();
361 for (sintL i = 0; i < xlen; i++) {
362 if (!R->_zerop(x[i]))
368 static const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
370 if (!(UPR->basering() == x.ring())) cl_abort();
371 { DeclarePoly(_cl_MI,x);
372 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
374 return _cl_UP(UPR, cl_null_GV_I);
377 var cl_GV_MI result = cl_GV_MI(len,R);
379 return _cl_UP(UPR, result);
383 static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
385 DeclarePoly(cl_GV_MI,x);
386 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
387 if (index < x.length())
388 return cl_MI(R, x[index]);
393 static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
396 return _cl_UP(UPR, cl_null_GV_I);
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));
404 static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
406 DeclareMutablePoly(cl_GV_MI,x);
407 if (!(UPR->basering() == y.ring())) cl_abort();
408 { DeclarePoly(_cl_MI,y);
409 if (!(index < x.length())) cl_abort();
413 static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
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.length();
419 modint_normalize(R,x,len);
422 static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
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())) cl_abort();
430 { DeclarePoly(_cl_MI,y);
431 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
432 var uintL len = x.length();
436 return cl_MI(R, x[0]);
440 z = R->_plus(R->_mul(z,y),x[i]);
444 static cl_univpoly_setops modint_setops = {
449 static cl_univpoly_addops modint_addops = {
457 static cl_univpoly_mulops modint_mulops = {
465 static cl_univpoly_modulops modint_modulops = {
469 static cl_univpoly_polyops modint_polyops = {
480 class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
481 SUBCLASS_cl_heap_univpoly_ring()
484 cl_heap_modint_univpoly_ring (const cl_ring& r);
486 ~cl_heap_modint_univpoly_ring () {}
489 static void cl_heap_modint_univpoly_ring_destructor (cl_heap* pointer)
491 (*(cl_heap_modint_univpoly_ring*)pointer).~cl_heap_modint_univpoly_ring();
494 cl_class cl_class_modint_univpoly_ring = {
495 cl_heap_modint_univpoly_ring_destructor,
496 cl_class_flags_univpoly_ring
500 inline 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)
503 type = &cl_class_modint_univpoly_ring;