src/polynomial/elem/cl_UP_number.h

   1 // Univariate Polynomials over some subring of the numbers.
   2
   3 #include "cln/SV_number.h"
   4 #include "cln/number.h"
   5 #include "cln/integer.h"
   6 #include "cln/abort.h"
   7
   8 namespace cln {
   9
  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; }
  13
  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)
  17 {
  18         if (ops.zerop(result[len-1])) {
  19                 len--;
  20                 while (len > 0) {
  21                         if (!ops.zerop(result[len-1]))
  22                                 break;
  23                         len--;
  24                 }
  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]);
  28                 result = newresult;
  29         }
  30 }
  31
  32 static void num_fprint (cl_heap_univpoly_ring* UPR, cl_ostream stream, const _cl_UP& x)
  33 {{
  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();
  37         if (xlen == 0)
  38                 fprint(stream, "0");
  39         else {
  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])) {
  44                                 if (i < xlen-1)
  45                                         fprint(stream, " + ");
  46                                 fprint(stream, cl_ring_element(R,x[i]));
  47                                 if (i > 0) {
  48                                         fprint(stream, "*");
  49                                         fprint(stream, varname);
  50                                         if (i != 1) {
  51                                                 fprint(stream, "^");
  52                                                 fprintdecimal(stream, i);
  53                                         }
  54                                 }
  55                         }
  56         }
  57 }}
  58
  59 static cl_boolean num_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  60 {{
  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();
  66         if (!(xlen == ylen))
  67                 return cl_false;
  68         for (var sintL i = xlen-1; i >= 0; i--)
  69                 if (!ops.equal(x[i],y[i]))
  70                         return cl_false;
  71         return cl_true;
  72 }}
  73
  74 static const _cl_UP num_zero (cl_heap_univpoly_ring* UPR)
  75 {
  76         return _cl_UP(UPR, cl_null_SV_number);
  77 }
  78
  79 static cl_boolean num_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  80 {
  81         unused UPR;
  82  {      DeclarePoly(cl_SV_number,x);
  83         var sintL xlen = x.length();
  84         if (xlen == 0)
  85                 return cl_true;
  86         else
  87                 return cl_false;
  88 }}
  89
  90 static const _cl_UP num_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  91 {{
  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();
  97         if (xlen == 0)
  98                 return _cl_UP(UPR, y);
  99         if (ylen == 0)
 100                 return _cl_UP(UPR, x);
 101         // Now xlen > 0, ylen > 0.
 102         if (xlen > ylen) {
 103                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 104                 var sintL i;
 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);
 110         }
 111         if (xlen < ylen) {
 112                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 113                 var sintL i;
 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);
 119         }
 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);
 129                 }
 130         }
 131         return _cl_UP(UPR, cl_null_SV_number);
 132 }}
 133
 134 static const _cl_UP num_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 135 {{
 136         DeclarePoly(cl_SV_number,x);
 137         DeclarePoly(cl_SV_number,y);
 138         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 139         var sintL xlen = x.length();
 140         var sintL ylen = y.length();
 141         if (xlen == 0)
 142                 return _cl_UP(UPR, y);
 143         if (ylen == 0)
 144                 return _cl_UP(UPR, x);
 145         // Now xlen > 0, ylen > 0.
 146         if (xlen > ylen) {
 147                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 148                 var sintL i;
 149                 for (i = xlen-1; i >= ylen; i--)
 150                         init1(cl_number, result[i]) (x[i]);
 151                 for (i = ylen-1; i >= 0; i--)
 152                         init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 153                 return _cl_UP(UPR, result);
 154         }
 155         if (xlen < ylen) {
 156                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 157                 var sintL i;
 158                 for (i = ylen-1; i >= xlen; i--)
 159                         init1(cl_number, result[i]) (ops.uminus(y[i]));
 160                 for (i = xlen-1; i >= 0; i--)
 161                         init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 162                 return _cl_UP(UPR, result);
 163         }
 164         // Now xlen = ylen > 0. Add and normalize simultaneously.
 165         for (var sintL i = xlen-1; i >= 0; i--) {
 166                 var cl_number hicoeff = ops.minus(x[i],y[i]);
 167                 if (!ops.zerop(hicoeff)) {
 168                         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(i+1));
 169                         init1(cl_number, result[i]) (hicoeff);
 170                         for (i-- ; i >= 0; i--)
 171                                 init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 172                         return _cl_UP(UPR, result);
 173                 }
 174         }
 175         return _cl_UP(UPR, cl_null_SV_number);
 176 }}
 177
 178 static const _cl_UP num_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 179 {{
 180         DeclarePoly(cl_SV_number,x);
 181         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 182         var sintL xlen = x.length();
 183         if (xlen == 0)
 184                 return _cl_UP(UPR, x);
 185         // Now xlen > 0.
 186         // Negate. No normalization necessary, since the degree doesn't change.
 187         var sintL i = xlen-1;
 188         var cl_number hicoeff = ops.uminus(x[i]);
 189         if (ops.zerop(hicoeff)) cl_abort();
 190         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 191         init1(cl_number, result[i]) (hicoeff);
 192         for (i-- ; i >= 0; i--)
 193                 init1(cl_number, result[i]) (ops.uminus(x[i]));
 194         return _cl_UP(UPR, result);
 195 }}
 196
 197 static const _cl_UP num_one (cl_heap_univpoly_ring* UPR)
 198 {
 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);
 202 }
 203
 204 static const _cl_UP num_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
 205 {
 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);
 209 }
 210
 211 static const _cl_UP num_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 212 {{
 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();
 218         if (xlen == 0)
 219                 return _cl_UP(UPR, x);
 220         if (ylen == 0)
 221                 return _cl_UP(UPR, y);
 222         // Multiply.
 223         var sintL len = xlen + ylen - 1;
 224         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(len));
 225         if (xlen < ylen) {
 226                 {
 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]));
 231                 }
 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]));
 237                 }
 238         } else {
 239                 {
 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));
 244                 }
 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));
 250                 }
 251         }
 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);
 256 }}
 257
 258 static const _cl_UP num_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 259 {{
 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();
 263         if (xlen == 0)
 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));
 267         if (xlen > 1) {
 268                 // Loop through all 0 <= j < i <= xlen-1.
 269                 {
 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]));
 274                 }
 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]));
 280                 }}
 281                 // Double.
 282                 {for (sintL i = len-2; i >= 1; i--)
 283                         result[i] = ops.plus(result[i],result[i]);
 284                 }
 285                 // Add squares.
 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]));
 289         }
 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);
 295 }}
 296
 297 static const _cl_UP num_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
 298 {
 299         var _cl_UP a = x;
 300         var cl_I b = y;
 301         while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
 302         var _cl_UP c = a;
 303         until (b == 1)
 304           { b = b >> 1;
 305             a = UPR->_square(a);
 306             if (oddp(b)) { c = UPR->_mul(a,c); }
 307           }
 308         return c;
 309 }
 310
 311 static const _cl_UP num_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
 312 {
 313         if (!(UPR->basering() == x.ring())) cl_abort();
 314  {
 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();
 319         if (ylen == 0)
 320                 return _cl_UP(UPR, y);
 321         if (ops.zerop(x))
 322                 return _cl_UP(UPR, cl_null_SV_number);
 323         // Now ylen > 0.
 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);
 329 }}
 330
 331 static sintL num_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 332 {
 333         unused UPR;
 334  {      DeclarePoly(cl_SV_number,x);
 335         return (sintL) x.length() - 1;
 336 }}
 337
 338 static const _cl_UP num_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 339 {
 340         if (!(UPR->basering() == x.ring())) cl_abort();
 341  {      DeclarePoly(cl_number,x);
 342         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 343         if (ops.zerop(x))
 344                 return _cl_UP(UPR, cl_null_SV_number);
 345         else {
 346                 var sintL len = e+1;
 347                 var cl_SV_number result = cl_SV_number(len);
 348                 result[e] = x;
 349                 return _cl_UP(UPR, result);
 350         }
 351 }}
 352
 353 static const cl_ring_element num_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 354 {{
 355         DeclarePoly(cl_SV_number,x);
 356         var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
 357         if (index < x.length())
 358                 return cl_ring_element(R, x[index]);
 359         else
 360                 return R->zero();
 361 }}
 362
 363 static const _cl_UP num_create (cl_heap_univpoly_ring* UPR, sintL deg)
 364 {
 365         if (deg < 0)
 366                 return _cl_UP(UPR, cl_null_SV_number);
 367         else {
 368                 var sintL len = deg+1;
 369                 return _cl_UP(UPR, cl_SV_number(len));
 370         }
 371 }
 372
 373 static void num_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 374 {{
 375         DeclareMutablePoly(cl_SV_number,x);
 376         if (!(UPR->basering() == y.ring())) cl_abort();
 377   {     DeclarePoly(cl_number,y);
 378         if (!(index < x.length())) cl_abort();
 379         x[index] = y;
 380 }}}
 381
 382 static void num_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 383 {{
 384         DeclareMutablePoly(cl_SV_number,x); // NB: x is modified by reference!
 385         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 386         var uintL len = x.length();
 387         if (len > 0)
 388                 num_normalize(ops,x,len);
 389 }}
 390
 391 static const cl_ring_element num_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 392 {{
 393         // Method:
 394         // If x = 0, return 0.
 395         // If y = 0, return x[0].
 396         // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
 397         DeclarePoly(cl_SV_number,x);
 398         if (!(UPR->basering() == y.ring())) cl_abort();
 399   {     DeclarePoly(cl_number,y);
 400         var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
 401         var cl_number_ring_ops<cl_number>& ops = *R->ops;
 402         var uintL len = x.length();
 403         if (len==0)
 404                 return R->zero();
 405         if (ops.zerop(y))
 406                 return cl_ring_element(R, x[0]);
 407         var sintL i = len-1;
 408         var cl_number z = x[i];
 409         for ( ; --i >= 0; )
 410                 z = ops.plus(ops.mul(z,y),x[i]);
 411         return cl_ring_element(R, z);
 412 }}}
 413
 414 static cl_univpoly_setops num_setops = {
 415         num_fprint,
 416         num_equal
 417 };
 418
 419 static cl_univpoly_addops num_addops = {
 420         num_zero,
 421         num_zerop,
 422         num_plus,
 423         num_minus,
 424         num_uminus
 425 };
 426
 427 static cl_univpoly_mulops num_mulops = {
 428         num_one,
 429         num_canonhom,
 430         num_mul,
 431         num_square,
 432         num_exptpos
 433 };
 434
 435 static cl_univpoly_modulops num_modulops = {
 436         num_scalmul
 437 };
 438
 439 static cl_univpoly_polyops num_polyops = {
 440         num_degree,
 441         num_monomial,
 442         num_coeff,
 443         num_create,
 444         num_set_coeff,
 445         num_finalize,
 446         num_eval
 447 };
 448
 449 class cl_heap_num_univpoly_ring : public cl_heap_univpoly_ring {
 450         SUBCLASS_cl_heap_univpoly_ring()
 451 public:
 452         // Constructor.
 453         cl_heap_num_univpoly_ring (const cl_ring& r)
 454                 : cl_heap_univpoly_ring (r, &num_setops, &num_addops, &num_mulops, &num_modulops, &num_polyops) {}
 455 };
 456
 457 }  // namespace cln