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_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
148 DeclarePoly(cl_GV_MI,x);
149 DeclarePoly(cl_GV_MI,y);
150 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
151 var sintL xlen = x.length();
152 var sintL ylen = y.length();
154 return _cl_UP(UPR, y);
156 return _cl_UP(UPR, x);
157 // Now xlen > 0, ylen > 0.
159 var cl_GV_MI result = cl_GV_MI(xlen,R);
162 for (i = xlen-1; i >= ylen; i--)
165 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
167 for (i = ylen-1; i >= 0; i--)
168 result[i] = R->_minus(x[i],y[i]);
169 return _cl_UP(UPR, result);
172 var cl_GV_MI result = cl_GV_MI(ylen,R);
174 for (i = ylen-1; i >= xlen; i--)
175 result[i] = R->_uminus(y[i]);
176 for (i = xlen-1; i >= 0; i--)
177 result[i] = R->_minus(x[i],y[i]);
178 return _cl_UP(UPR, result);
180 // Now xlen = ylen > 0. Add and normalize simultaneously.
181 for (var sintL i = xlen-1; i >= 0; i--) {
182 var _cl_MI hicoeff = R->_minus(x[i],y[i]);
183 if (!R->_zerop(hicoeff)) {
184 var cl_GV_MI result = cl_GV_MI(i+1,R);
186 for (i-- ; i >= 0; i--)
187 result[i] = R->_minus(x[i],y[i]);
188 return _cl_UP(UPR, result);
191 return _cl_UP(UPR, cl_null_GV_I);
194 static const _cl_UP modint_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
196 DeclarePoly(cl_GV_MI,x);
197 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
198 var sintL xlen = x.length();
200 return _cl_UP(UPR, x);
202 // Negate. No normalization necessary, since the degree doesn't change.
203 var sintL i = xlen-1;
204 var _cl_MI hicoeff = R->_uminus(x[i]);
205 if (R->_zerop(hicoeff)) cl_abort();
206 var cl_GV_MI result = cl_GV_MI(xlen,R);
208 for (i-- ; i >= 0; i--)
209 result[i] = R->_uminus(x[i]);
210 return _cl_UP(UPR, result);
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 const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
358 if (!(UPR->basering() == x.ring())) cl_abort();
359 { DeclarePoly(_cl_MI,x);
360 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
362 return _cl_UP(UPR, cl_null_GV_I);
365 var cl_GV_MI result = cl_GV_MI(len,R);
367 return _cl_UP(UPR, result);
371 static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
373 DeclarePoly(cl_GV_MI,x);
374 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
375 if (index < x.length())
376 return cl_MI(R, x[index]);
381 static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
384 return _cl_UP(UPR, cl_null_GV_I);
386 var sintL len = deg+1;
387 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
388 return _cl_UP(UPR, cl_GV_MI(len,R));
392 static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
394 DeclareMutablePoly(cl_GV_MI,x);
395 if (!(UPR->basering() == y.ring())) cl_abort();
396 { DeclarePoly(_cl_MI,y);
397 if (!(index < x.length())) cl_abort();
401 static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
403 DeclareMutablePoly(cl_GV_MI,x); // NB: x is modified by reference!
404 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
405 var uintL len = x.length();
407 modint_normalize(R,x,len);
410 static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
413 // If x = 0, return 0.
414 // If y = 0, return x[0].
415 // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
416 DeclarePoly(cl_GV_MI,x);
417 if (!(UPR->basering() == y.ring())) cl_abort();
418 { DeclarePoly(_cl_MI,y);
419 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
420 var uintL len = x.length();
424 return cl_MI(R, x[0]);
428 z = R->_plus(R->_mul(z,y),x[i]);
432 static cl_univpoly_setops modint_setops = {
437 static cl_univpoly_addops modint_addops = {
445 static cl_univpoly_mulops modint_mulops = {
453 static cl_univpoly_modulops modint_modulops = {
457 static cl_univpoly_polyops modint_polyops = {
467 class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
468 SUBCLASS_cl_heap_univpoly_ring()
471 cl_heap_modint_univpoly_ring (const cl_ring& r)
472 : cl_heap_univpoly_ring (r, &modint_setops, &modint_addops, &modint_mulops, &modint_modulops, &modint_polyops) {}