1 // Univariate Polynomials over some subring of the numbers.
3 #include "cln/SV_number.h"
4 #include "cln/number.h"
5 #include "cln/integer.h"
10 // Assume a ring is a number ring.
11 inline cl_heap_number_ring* TheNumberRing (const cl_ring& R)
12 { return (cl_heap_number_ring*) R.heappointer; }
14 // Normalize a vector: remove leading zero coefficients.
15 // The result vector is known to have length len > 0.
16 static inline void num_normalize (cl_number_ring_ops<cl_number>& ops, cl_SV_number& result, uintL len)
18 if (ops.zerop(result[len-1])) {
21 if (!ops.zerop(result[len-1]))
25 var cl_SV_number newresult = cl_SV_number(cl_make_heap_SV_number_uninit(len));
26 for (var sintL i = len-1; i >= 0; i--)
27 init1(cl_number, newresult[i]) (result[i]);
32 static void num_fprint (cl_heap_univpoly_ring* UPR, std::ostream& stream, const _cl_UP& x)
34 DeclarePoly(cl_SV_number,x);
35 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
36 var sintL xlen = x.length();
40 var const cl_ring& R = UPR->basering();
41 var cl_string varname = get_varname(UPR);
42 for (var sintL i = xlen-1; i >= 0; i--)
43 if (!ops.zerop(x[i])) {
45 fprint(stream, " + ");
46 fprint(stream, cl_ring_element(R,x[i]));
49 fprint(stream, varname);
52 fprintdecimal(stream, i);
59 static cl_boolean num_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
61 DeclarePoly(cl_SV_number,x);
62 DeclarePoly(cl_SV_number,y);
63 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
64 var sintL xlen = x.length();
65 var sintL ylen = y.length();
68 for (var sintL i = xlen-1; i >= 0; i--)
69 if (!ops.equal(x[i],y[i]))
74 static const _cl_UP num_zero (cl_heap_univpoly_ring* UPR)
76 return _cl_UP(UPR, cl_null_SV_number);
79 static cl_boolean num_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
82 { DeclarePoly(cl_SV_number,x);
83 var sintL xlen = x.length();
90 static const _cl_UP num_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
92 DeclarePoly(cl_SV_number,x);
93 DeclarePoly(cl_SV_number,y);
94 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
95 var sintL xlen = x.length();
96 var sintL ylen = y.length();
98 return _cl_UP(UPR, y);
100 return _cl_UP(UPR, x);
101 // Now xlen > 0, ylen > 0.
103 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
105 for (i = xlen-1; i >= ylen; i--)
106 init1(cl_number, result[i]) (x[i]);
107 for (i = ylen-1; i >= 0; i--)
108 init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
109 return _cl_UP(UPR, result);
112 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
114 for (i = ylen-1; i >= xlen; i--)
115 init1(cl_number, result[i]) (y[i]);
116 for (i = xlen-1; i >= 0; i--)
117 init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
118 return _cl_UP(UPR, result);
120 // Now xlen = ylen > 0. Add and normalize simultaneously.
121 for (var sintL i = xlen-1; i >= 0; i--) {
122 var cl_number hicoeff = ops.plus(x[i],y[i]);
123 if (!ops.zerop(hicoeff)) {
124 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(i+1));
125 init1(cl_number, result[i]) (hicoeff);
126 for (i-- ; i >= 0; i--)
127 init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
128 return _cl_UP(UPR, result);
131 return _cl_UP(UPR, cl_null_SV_number);
134 static const _cl_UP num_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
136 DeclarePoly(cl_SV_number,x);
137 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
138 var sintL xlen = x.length();
140 return _cl_UP(UPR, x);
142 // Negate. No normalization necessary, since the degree doesn't change.
143 var sintL i = xlen-1;
144 var cl_number hicoeff = ops.uminus(x[i]);
145 if (ops.zerop(hicoeff)) cl_abort();
146 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
147 init1(cl_number, result[i]) (hicoeff);
148 for (i-- ; i >= 0; i--)
149 init1(cl_number, result[i]) (ops.uminus(x[i]));
150 return _cl_UP(UPR, result);
153 static const _cl_UP num_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
155 DeclarePoly(cl_SV_number,x);
156 DeclarePoly(cl_SV_number,y);
157 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
158 var sintL xlen = x.length();
159 var sintL ylen = y.length();
161 return _cl_UP(UPR, x);
163 return num_uminus(UPR, _cl_UP(UPR, y));
164 // Now xlen > 0, ylen > 0.
166 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
168 for (i = xlen-1; i >= ylen; i--)
169 init1(cl_number, result[i]) (x[i]);
170 for (i = ylen-1; i >= 0; i--)
171 init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
172 return _cl_UP(UPR, result);
175 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
177 for (i = ylen-1; i >= xlen; i--)
178 init1(cl_number, result[i]) (ops.uminus(y[i]));
179 for (i = xlen-1; i >= 0; i--)
180 init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
181 return _cl_UP(UPR, result);
183 // Now xlen = ylen > 0. Add and normalize simultaneously.
184 for (var sintL i = xlen-1; i >= 0; i--) {
185 var cl_number hicoeff = ops.minus(x[i],y[i]);
186 if (!ops.zerop(hicoeff)) {
187 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(i+1));
188 init1(cl_number, result[i]) (hicoeff);
189 for (i-- ; i >= 0; i--)
190 init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
191 return _cl_UP(UPR, result);
194 return _cl_UP(UPR, cl_null_SV_number);
197 static const _cl_UP num_one (cl_heap_univpoly_ring* UPR)
199 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(1));
200 init1(cl_number, result[0]) (1);
201 return _cl_UP(UPR, result);
204 static const _cl_UP num_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
206 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(1));
207 init1(cl_number, result[0]) (x);
208 return _cl_UP(UPR, result);
211 static const _cl_UP num_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
213 DeclarePoly(cl_SV_number,x);
214 DeclarePoly(cl_SV_number,y);
215 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
216 var sintL xlen = x.length();
217 var sintL ylen = y.length();
219 return _cl_UP(UPR, x);
221 return _cl_UP(UPR, y);
223 var sintL len = xlen + ylen - 1;
224 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(len));
227 var sintL i = xlen-1;
228 var cl_number xi = x[i];
229 for (sintL j = ylen-1; j >= 0; j--)
230 init1(cl_number, result[i+j]) (ops.mul(xi,y[j]));
232 for (sintL i = xlen-2; i >= 0; i--) {
233 var cl_number xi = x[i];
234 for (sintL j = ylen-1; j > 0; j--)
235 result[i+j] = ops.plus(result[i+j],ops.mul(xi,y[j]));
236 /* j=0 */ init1(cl_number, result[i]) (ops.mul(xi,y[0]));
240 var sintL j = ylen-1;
241 var cl_number yj = y[j];
242 for (sintL i = xlen-1; i >= 0; i--)
243 init1(cl_number, result[i+j]) (ops.mul(x[i],yj));
245 for (sintL j = ylen-2; j >= 0; j--) {
246 var cl_number yj = y[j];
247 for (sintL i = xlen-1; i > 0; i--)
248 result[i+j] = ops.plus(result[i+j],ops.mul(x[i],yj));
249 /* i=0 */ init1(cl_number, result[j]) (ops.mul(x[0],yj));
252 // Normalize (not necessary in integral domains).
253 //num_normalize(ops,result,len);
254 if (ops.zerop(result[len-1])) cl_abort();
255 return _cl_UP(UPR, result);
258 static const _cl_UP num_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
260 DeclarePoly(cl_SV_number,x);
261 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
262 var sintL xlen = x.length();
264 return cl_UP(UPR, x);
265 var sintL len = 2*xlen-1;
266 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(len));
268 // Loop through all 0 <= j < i <= xlen-1.
270 var sintL i = xlen-1;
271 var cl_number xi = x[i];
272 for (sintL j = i-1; j >= 0; j--)
273 init1(cl_number, result[i+j]) (ops.mul(xi,x[j]));
275 {for (sintL i = xlen-2; i >= 1; i--) {
276 var cl_number xi = x[i];
277 for (sintL j = i-1; j >= 1; j--)
278 result[i+j] = ops.plus(result[i+j],ops.mul(xi,x[j]));
279 /* j=0 */ init1(cl_number, result[i]) (ops.mul(xi,x[0]));
282 {for (sintL i = len-2; i >= 1; i--)
283 result[i] = ops.plus(result[i],result[i]);
286 init1(cl_number, result[2*(xlen-1)]) (ops.square(x[xlen-1]));
287 for (sintL i = xlen-2; i >= 1; i--)
288 result[2*i] = ops.plus(result[2*i],ops.square(x[i]));
290 init1(cl_number, result[0]) (ops.square(x[0]));
291 // Normalize (not necessary in integral domains).
292 //num_normalize(ops,result,len);
293 if (ops.zerop(result[len-1])) cl_abort();
294 return _cl_UP(UPR, result);
297 static const _cl_UP num_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
301 while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
306 if (oddp(b)) { c = UPR->_mul(a,c); }
311 static const _cl_UP num_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
313 if (!(UPR->basering() == x.ring())) cl_abort();
315 DeclarePoly(cl_number,x);
316 DeclarePoly(cl_SV_number,y);
317 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
318 var sintL ylen = y.length();
320 return _cl_UP(UPR, y);
322 return _cl_UP(UPR, cl_null_SV_number);
324 // No normalization necessary, since the degree doesn't change.
325 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
326 for (sintL i = ylen-1; i >= 0; i--)
327 init1(cl_number, result[i]) (ops.mul(x,y[i]));
328 return _cl_UP(UPR, result);
331 static sintL num_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
334 { DeclarePoly(cl_SV_number,x);
335 return (sintL) x.length() - 1;
338 static sintL num_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
340 DeclarePoly(cl_SV_number,x);
341 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
342 var sintL xlen = x.length();
343 for (sintL i = 0; i < xlen; i++) {
344 if (!ops.zerop(x[i]))
350 static const _cl_UP num_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
352 if (!(UPR->basering() == x.ring())) cl_abort();
353 { DeclarePoly(cl_number,x);
354 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
356 return _cl_UP(UPR, cl_null_SV_number);
359 var cl_SV_number result = cl_SV_number(len);
361 return _cl_UP(UPR, result);
365 static const cl_ring_element num_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
367 DeclarePoly(cl_SV_number,x);
368 var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
369 if (index < x.length())
370 return cl_ring_element(R, x[index]);
375 static const _cl_UP num_create (cl_heap_univpoly_ring* UPR, sintL deg)
378 return _cl_UP(UPR, cl_null_SV_number);
380 var sintL len = deg+1;
381 return _cl_UP(UPR, cl_SV_number(len));
385 static void num_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
387 DeclareMutablePoly(cl_SV_number,x);
388 if (!(UPR->basering() == y.ring())) cl_abort();
389 { DeclarePoly(cl_number,y);
390 if (!(index < x.length())) cl_abort();
394 static void num_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
396 DeclareMutablePoly(cl_SV_number,x); // NB: x is modified by reference!
397 var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
398 var uintL len = x.length();
400 num_normalize(ops,x,len);
403 static const cl_ring_element num_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
406 // If x = 0, return 0.
407 // If y = 0, return x[0].
408 // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
409 DeclarePoly(cl_SV_number,x);
410 if (!(UPR->basering() == y.ring())) cl_abort();
411 { DeclarePoly(cl_number,y);
412 var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
413 var cl_number_ring_ops<cl_number>& ops = *R->ops;
414 var uintL len = x.length();
418 return cl_ring_element(R, x[0]);
420 var cl_number z = x[i];
422 z = ops.plus(ops.mul(z,y),x[i]);
423 return cl_ring_element(R, z);
426 static cl_univpoly_setops num_setops = {
431 static cl_univpoly_addops num_addops = {
439 static cl_univpoly_mulops num_mulops = {
447 static cl_univpoly_modulops num_modulops = {
451 static cl_univpoly_polyops num_polyops = {
462 class cl_heap_num_univpoly_ring : public cl_heap_univpoly_ring {
463 SUBCLASS_cl_heap_univpoly_ring()
466 cl_heap_num_univpoly_ring (const cl_ring& r);
468 ~cl_heap_num_univpoly_ring () {}
471 static void cl_heap_num_univpoly_ring_destructor (cl_heap* pointer)
473 (*(cl_heap_num_univpoly_ring*)pointer).~cl_heap_num_univpoly_ring();
476 cl_class cl_class_num_univpoly_ring = {
477 cl_heap_num_univpoly_ring_destructor,
478 cl_class_flags_univpoly_ring
482 inline cl_heap_num_univpoly_ring::cl_heap_num_univpoly_ring (const cl_ring& r)
483 : cl_heap_univpoly_ring (r, &num_setops, &num_addops, &num_mulops, &num_modulops, &num_polyops)
485 type = &cl_class_num_univpoly_ring;