/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
Possible License(s): GPL-2.0
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 }
303 // Add squares.
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
449static cl_univpoly_addops modint_addops = {
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