src/polynomial/elem/cl_UP_number.h

   1 // Univariate Polynomials over some subring of the numbers.
   2
   3 #include "cl_SV_number.h"
   4 #include "cl_number.h"
   5 #include "cl_integer.h"
   6 #include "cl_abort.h"
   7
   8 // Assume a ring is a number ring.
   9   inline cl_heap_number_ring* TheNumberRing (const cl_ring& R)
  10   { return (cl_heap_number_ring*) R.heappointer; }
  11
  12 // Normalize a vector: remove leading zero coefficients.
  13 // The result vector is known to have length len > 0.
  14 static inline void num_normalize (cl_number_ring_ops<cl_number>& ops, cl_SV_number& result, uintL len)
  15 {
  16         if (ops.zerop(result[len-1])) {
  17                 len--;
  18                 while (len > 0) {
  19                         if (!ops.zerop(result[len-1]))
  20                                 break;
  21                         len--;
  22                 }
  23                 var cl_SV_number newresult = cl_SV_number(cl_make_heap_SV_number_uninit(len));
  24                 for (var sintL i = len-1; i >= 0; i--)
  25                         init1(cl_number, newresult[i]) (result[i]);
  26                 result = newresult;
  27         }
  28 }
  29
  30 static void num_fprint (cl_heap_univpoly_ring* UPR, cl_ostream stream, const _cl_UP& x)
  31 {{
  32         DeclarePoly(cl_SV_number,x);
  33         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
  34         var sintL xlen = x.length();
  35         if (xlen == 0)
  36                 fprint(stream, "0");
  37         else {
  38                 var const cl_ring& R = UPR->basering();
  39                 var cl_string varname = get_varname(UPR);
  40                 for (var sintL i = xlen-1; i >= 0; i--)
  41                         if (!ops.zerop(x[i])) {
  42                                 if (i < xlen-1)
  43                                         fprint(stream, " + ");
  44                                 fprint(stream, cl_ring_element(R,x[i]));
  45                                 if (i > 0) {
  46                                         fprint(stream, "*");
  47                                         fprint(stream, varname);
  48                                         if (i != 1) {
  49                                                 fprint(stream, "^");
  50                                                 fprintdecimal(stream, i);
  51                                         }
  52                                 }
  53                         }
  54         }
  55 }}
  56
  57 static cl_boolean num_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  58 {{
  59         DeclarePoly(cl_SV_number,x);
  60         DeclarePoly(cl_SV_number,y);
  61         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
  62         var sintL xlen = x.length();
  63         var sintL ylen = y.length();
  64         if (!(xlen == ylen))
  65                 return cl_false;
  66         for (var sintL i = xlen-1; i >= 0; i--)
  67                 if (!ops.equal(x[i],y[i]))
  68                         return cl_false;
  69         return cl_true;
  70 }}
  71
  72 static const _cl_UP num_zero (cl_heap_univpoly_ring* UPR)
  73 {
  74         return _cl_UP(UPR, cl_null_SV_number);
  75 }
  76
  77 static cl_boolean num_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  78 {
  79         unused UPR;
  80  {      DeclarePoly(cl_SV_number,x);
  81         var sintL xlen = x.length();
  82         if (xlen == 0)
  83                 return cl_true;
  84         else
  85                 return cl_false;
  86 }}
  87
  88 static const _cl_UP num_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  89 {{
  90         DeclarePoly(cl_SV_number,x);
  91         DeclarePoly(cl_SV_number,y);
  92         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
  93         var sintL xlen = x.length();
  94         var sintL ylen = y.length();
  95         if (xlen == 0)
  96                 return _cl_UP(UPR, y);
  97         if (ylen == 0)
  98                 return _cl_UP(UPR, x);
  99         // Now xlen > 0, ylen > 0.
 100         if (xlen > ylen) {
 101                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 102                 var sintL i;
 103                 for (i = xlen-1; i >= ylen; i--)
 104                         init1(cl_number, result[i]) (x[i]);
 105                 for (i = ylen-1; i >= 0; i--)
 106                         init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
 107                 return _cl_UP(UPR, result);
 108         }
 109         if (xlen < ylen) {
 110                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 111                 var sintL i;
 112                 for (i = ylen-1; i >= xlen; i--)
 113                         init1(cl_number, result[i]) (y[i]);
 114                 for (i = xlen-1; i >= 0; i--)
 115                         init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
 116                 return _cl_UP(UPR, result);
 117         }
 118         // Now xlen = ylen > 0. Add and normalize simultaneously.
 119         for (var sintL i = xlen-1; i >= 0; i--) {
 120                 var cl_number hicoeff = ops.plus(x[i],y[i]);
 121                 if (!ops.zerop(hicoeff)) {
 122                         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(i+1));
 123                         init1(cl_number, result[i]) (hicoeff);
 124                         for (i-- ; i >= 0; i--)
 125                                 init1(cl_number, result[i]) (ops.plus(x[i],y[i]));
 126                         return _cl_UP(UPR, result);
 127                 }
 128         }
 129         return _cl_UP(UPR, cl_null_SV_number);
 130 }}
 131
 132 static const _cl_UP num_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 133 {{
 134         DeclarePoly(cl_SV_number,x);
 135         DeclarePoly(cl_SV_number,y);
 136         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 137         var sintL xlen = x.length();
 138         var sintL ylen = y.length();
 139         if (xlen == 0)
 140                 return _cl_UP(UPR, y);
 141         if (ylen == 0)
 142                 return _cl_UP(UPR, x);
 143         // Now xlen > 0, ylen > 0.
 144         if (xlen > ylen) {
 145                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 146                 var sintL i;
 147                 for (i = xlen-1; i >= ylen; i--)
 148                         init1(cl_number, result[i]) (x[i]);
 149                 for (i = ylen-1; i >= 0; i--)
 150                         init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 151                 return _cl_UP(UPR, result);
 152         }
 153         if (xlen < ylen) {
 154                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 155                 var sintL i;
 156                 for (i = ylen-1; i >= xlen; i--)
 157                         init1(cl_number, result[i]) (ops.uminus(y[i]));
 158                 for (i = xlen-1; i >= 0; i--)
 159                         init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 160                 return _cl_UP(UPR, result);
 161         }
 162         // Now xlen = ylen > 0. Add and normalize simultaneously.
 163         for (var sintL i = xlen-1; i >= 0; i--) {
 164                 var cl_number hicoeff = ops.minus(x[i],y[i]);
 165                 if (!ops.zerop(hicoeff)) {
 166                         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(i+1));
 167                         init1(cl_number, result[i]) (hicoeff);
 168                         for (i-- ; i >= 0; i--)
 169                                 init1(cl_number, result[i]) (ops.minus(x[i],y[i]));
 170                         return _cl_UP(UPR, result);
 171                 }
 172         }
 173         return _cl_UP(UPR, cl_null_SV_number);
 174 }}
 175
 176 static const _cl_UP num_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 177 {{
 178         DeclarePoly(cl_SV_number,x);
 179         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 180         var sintL xlen = x.length();
 181         if (xlen == 0)
 182                 return _cl_UP(UPR, x);
 183         // Now xlen > 0.
 184         // Negate. No normalization necessary, since the degree doesn't change.
 185         var sintL i = xlen-1;
 186         var cl_number hicoeff = ops.uminus(x[i]);
 187         if (ops.zerop(hicoeff)) cl_abort();
 188         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 189         init1(cl_number, result[i]) (hicoeff);
 190         for (i-- ; i >= 0; i--)
 191                 init1(cl_number, result[i]) (ops.uminus(x[i]));
 192         return _cl_UP(UPR, result);
 193 }}
 194
 195 static const _cl_UP num_one (cl_heap_univpoly_ring* UPR)
 196 {
 197         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(1));
 198         init1(cl_number, result[0]) (1);
 199         return _cl_UP(UPR, result);
 200 }
 201
 202 static const _cl_UP num_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
 203 {
 204         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(1));
 205         init1(cl_number, result[0]) (x);
 206         return _cl_UP(UPR, result);
 207 }
 208
 209 static const _cl_UP num_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 210 {{
 211         DeclarePoly(cl_SV_number,x);
 212         DeclarePoly(cl_SV_number,y);
 213         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 214         var sintL xlen = x.length();
 215         var sintL ylen = y.length();
 216         if (xlen == 0)
 217                 return _cl_UP(UPR, x);
 218         if (ylen == 0)
 219                 return _cl_UP(UPR, y);
 220         // Multiply.
 221         var sintL len = xlen + ylen - 1;
 222         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(len));
 223         if (xlen < ylen) {
 224                 {
 225                         var sintL i = xlen-1;
 226                         var cl_number xi = x[i];
 227                         for (sintL j = ylen-1; j >= 0; j--)
 228                                 init1(cl_number, result[i+j]) (ops.mul(xi,y[j]));
 229                 }
 230                 for (sintL i = xlen-2; i >= 0; i--) {
 231                         var cl_number xi = x[i];
 232                         for (sintL j = ylen-1; j > 0; j--)
 233                                 result[i+j] = ops.plus(result[i+j],ops.mul(xi,y[j]));
 234                         /* j=0 */ init1(cl_number, result[i]) (ops.mul(xi,y[0]));
 235                 }
 236         } else {
 237                 {
 238                         var sintL j = ylen-1;
 239                         var cl_number yj = y[j];
 240                         for (sintL i = xlen-1; i >= 0; i--)
 241                                 init1(cl_number, result[i+j]) (ops.mul(x[i],yj));
 242                 }
 243                 for (sintL j = ylen-2; j >= 0; j--) {
 244                         var cl_number yj = y[j];
 245                         for (sintL i = xlen-1; i > 0; i--)
 246                                 result[i+j] = ops.plus(result[i+j],ops.mul(x[i],yj));
 247                         /* i=0 */ init1(cl_number, result[j]) (ops.mul(x[0],yj));
 248                 }
 249         }
 250         // Normalize (not necessary in integral domains).
 251         //num_normalize(ops,result,len);
 252         if (ops.zerop(result[len-1])) cl_abort();
 253         return _cl_UP(UPR, result);
 254 }}
 255
 256 static const _cl_UP num_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 257 {{
 258         DeclarePoly(cl_SV_number,x);
 259         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 260         var sintL xlen = x.length();
 261         if (xlen == 0)
 262                 return cl_UP(UPR, x);
 263         var sintL len = 2*xlen-1;
 264         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(len));
 265         if (xlen > 1) {
 266                 // Loop through all 0 <= j < i <= xlen-1.
 267                 {
 268                         var sintL i = xlen-1;
 269                         var cl_number xi = x[i];
 270                         for (sintL j = i-1; j >= 0; j--)
 271                                 init1(cl_number, result[i+j]) (ops.mul(xi,x[j]));
 272                 }
 273                 {for (sintL i = xlen-2; i >= 1; i--) {
 274                         var cl_number xi = x[i];
 275                         for (sintL j = i-1; j >= 1; j--)
 276                                 result[i+j] = ops.plus(result[i+j],ops.mul(xi,x[j]));
 277                         /* j=0 */ init1(cl_number, result[i]) (ops.mul(xi,x[0]));
 278                 }}
 279                 // Double.
 280                 {for (sintL i = len-2; i >= 1; i--)
 281                         result[i] = ops.plus(result[i],result[i]);
 282                 }
 283                 // Add squares.
 284                 init1(cl_number, result[2*(xlen-1)]) (ops.square(x[xlen-1]));
 285                 for (sintL i = xlen-2; i >= 1; i--)
 286                         result[2*i] = ops.plus(result[2*i],ops.square(x[i]));
 287         }
 288         init1(cl_number, result[0]) (ops.square(x[0]));
 289         // Normalize (not necessary in integral domains).
 290         //num_normalize(ops,result,len);
 291         if (ops.zerop(result[len-1])) cl_abort();
 292         return _cl_UP(UPR, result);
 293 }}
 294
 295 static const _cl_UP num_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
 296 {
 297         var _cl_UP a = x;
 298         var cl_I b = y;
 299         while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
 300         var _cl_UP c = a;
 301         until (b == 1)
 302           { b = b >> 1;
 303             a = UPR->_square(a);
 304             if (oddp(b)) { c = UPR->_mul(a,c); }
 305           }
 306         return c;
 307 }
 308
 309 static const _cl_UP num_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
 310 {
 311         if (!(UPR->basering() == x.ring())) cl_abort();
 312  {
 313         DeclarePoly(cl_number,x);
 314         DeclarePoly(cl_SV_number,y);
 315         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 316         var sintL ylen = y.length();
 317         if (ylen == 0)
 318                 return _cl_UP(UPR, y);
 319         if (ops.zerop(x))
 320                 return _cl_UP(UPR, cl_null_SV_number);
 321         // Now ylen > 0.
 322         // No normalization necessary, since the degree doesn't change.
 323         var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 324         for (sintL i = ylen-1; i >= 0; i--)
 325                 init1(cl_number, result[i]) (ops.mul(x,y[i]));
 326         return _cl_UP(UPR, result);
 327 }}
 328
 329 static sintL num_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 330 {
 331         unused UPR;
 332  {      DeclarePoly(cl_SV_number,x);
 333         return (sintL) x.length() - 1;
 334 }}
 335
 336 static const _cl_UP num_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 337 {
 338         if (!(UPR->basering() == x.ring())) cl_abort();
 339  {      DeclarePoly(cl_number,x);
 340         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 341         if (ops.zerop(x))
 342                 return _cl_UP(UPR, cl_null_SV_number);
 343         else {
 344                 var sintL len = e+1;
 345                 var cl_SV_number result = cl_SV_number(len);
 346                 result[e] = x;
 347                 return _cl_UP(UPR, result);
 348         }
 349 }}
 350
 351 static const cl_ring_element num_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 352 {{
 353         DeclarePoly(cl_SV_number,x);
 354         var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
 355         if (index < x.length())
 356                 return cl_ring_element(R, x[index]);
 357         else
 358                 return R->zero();
 359 }}
 360
 361 static const _cl_UP num_create (cl_heap_univpoly_ring* UPR, sintL deg)
 362 {
 363         if (deg < 0)
 364                 return _cl_UP(UPR, cl_null_SV_number);
 365         else {
 366                 var sintL len = deg+1;
 367                 return _cl_UP(UPR, cl_SV_number(len));
 368         }
 369 }
 370
 371 static void num_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 372 {{
 373         DeclareMutablePoly(cl_SV_number,x);
 374         if (!(UPR->basering() == y.ring())) cl_abort();
 375   {     DeclarePoly(cl_number,y);
 376         if (!(index < x.length())) cl_abort();
 377         x[index] = y;
 378 }}}
 379
 380 static void num_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 381 {{
 382         DeclareMutablePoly(cl_SV_number,x); // NB: x is modified by reference!
 383         var cl_number_ring_ops<cl_number>& ops = *TheNumberRing(UPR->basering())->ops;
 384         var uintL len = x.length();
 385         if (len > 0)
 386                 num_normalize(ops,x,len);
 387 }}
 388
 389 static const cl_ring_element num_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 390 {{
 391         // Method:
 392         // If x = 0, return 0.
 393         // If y = 0, return x[0].
 394         // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
 395         DeclarePoly(cl_SV_number,x);
 396         if (!(UPR->basering() == y.ring())) cl_abort();
 397   {     DeclarePoly(cl_number,y);
 398         var cl_heap_number_ring* R = TheNumberRing(UPR->basering());
 399         var cl_number_ring_ops<cl_number>& ops = *R->ops;
 400         var uintL len = x.length();
 401         if (len==0)
 402                 return R->zero();
 403         if (ops.zerop(y))
 404                 return cl_ring_element(R, x[0]);
 405         var sintL i = len-1;
 406         var cl_number z = x[i];
 407         for ( ; --i >= 0; )
 408                 z = ops.plus(ops.mul(z,y),x[i]);
 409         return cl_ring_element(R, z);
 410 }}}
 411
 412 static cl_univpoly_setops num_setops = {
 413         num_fprint,
 414         num_equal
 415 };
 416
 417 static cl_univpoly_addops num_addops = {
 418         num_zero,
 419         num_zerop,
 420         num_plus,
 421         num_minus,
 422         num_uminus
 423 };
 424
 425 static cl_univpoly_mulops num_mulops = {
 426         num_one,
 427         num_canonhom,
 428         num_mul,
 429         num_square,
 430         num_exptpos
 431 };
 432
 433 static cl_univpoly_modulops num_modulops = {
 434         num_scalmul
 435 };
 436
 437 static cl_univpoly_polyops num_polyops = {
 438         num_degree,
 439         num_monomial,
 440         num_coeff,
 441         num_create,
 442         num_set_coeff,
 443         num_finalize,
 444         num_eval
 445 };
 446
 447 class cl_heap_num_univpoly_ring : public cl_heap_univpoly_ring {
 448         SUBCLASS_cl_heap_univpoly_ring()
 449 public:
 450         // Constructor.
 451         cl_heap_num_univpoly_ring (const cl_ring& r)
 452                 : cl_heap_univpoly_ring (r, &num_setops, &num_addops, &num_mulops, &num_modulops, &num_polyops) {}
 453 };