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, std::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_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 135 {{
 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();
 139         if (xlen == 0)
 140                 return _cl_UP(UPR, x);
 141         // Now xlen > 0.
 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);
 151 }}
 152
 153 static const _cl_UP num_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 154 {{
 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();
 160         if (ylen == 0)
 161                 return _cl_UP(UPR, x);
 162         if (xlen == 0)
 163                 return num_uminus(UPR, _cl_UP(UPR, y));
 164         // Now xlen > 0, ylen > 0.
 165         if (xlen > ylen) {
 166                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(xlen));
 167                 var sintL i;
 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);
 173         }
 174         if (xlen < ylen) {
 175                 var cl_SV_number result = cl_SV_number(cl_make_heap_SV_number_uninit(ylen));
 176                 var sintL i;
 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);
 182         }
 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);
 192                 }
 193         }
 194         return _cl_UP(UPR, cl_null_SV_number);
 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 sintL num_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 339 {{
 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]))
 345                         return i;
 346         }
 347         return -1;
 348 }}
 349
 350 static const _cl_UP num_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 351 {
 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;
 355         if (ops.zerop(x))
 356                 return _cl_UP(UPR, cl_null_SV_number);
 357         else {
 358                 var sintL len = e+1;
 359                 var cl_SV_number result = cl_SV_number(len);
 360                 result[e] = x;
 361                 return _cl_UP(UPR, result);
 362         }
 363 }}
 364
 365 static const cl_ring_element num_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 366 {{
 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]);
 371         else
 372                 return R->zero();
 373 }}
 374
 375 static const _cl_UP num_create (cl_heap_univpoly_ring* UPR, sintL deg)
 376 {
 377         if (deg < 0)
 378                 return _cl_UP(UPR, cl_null_SV_number);
 379         else {
 380                 var sintL len = deg+1;
 381                 return _cl_UP(UPR, cl_SV_number(len));
 382         }
 383 }
 384
 385 static void num_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 386 {{
 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();
 391         x[index] = y;
 392 }}}
 393
 394 static void num_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 395 {{
 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();
 399         if (len > 0)
 400                 num_normalize(ops,x,len);
 401 }}
 402
 403 static const cl_ring_element num_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 404 {{
 405         // Method:
 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();
 415         if (len==0)
 416                 return R->zero();
 417         if (ops.zerop(y))
 418                 return cl_ring_element(R, x[0]);
 419         var sintL i = len-1;
 420         var cl_number z = x[i];
 421         for ( ; --i >= 0; )
 422                 z = ops.plus(ops.mul(z,y),x[i]);
 423         return cl_ring_element(R, z);
 424 }}}
 425
 426 static cl_univpoly_setops num_setops = {
 427         num_fprint,
 428         num_equal
 429 };
 430
 431 static cl_univpoly_addops num_addops = {
 432         num_zero,
 433         num_zerop,
 434         num_plus,
 435         num_minus,
 436         num_uminus
 437 };
 438
 439 static cl_univpoly_mulops num_mulops = {
 440         num_one,
 441         num_canonhom,
 442         num_mul,
 443         num_square,
 444         num_exptpos
 445 };
 446
 447 static cl_univpoly_modulops num_modulops = {
 448         num_scalmul
 449 };
 450
 451 static cl_univpoly_polyops num_polyops = {
 452         num_degree,
 453         num_ldegree,
 454         num_monomial,
 455         num_coeff,
 456         num_create,
 457         num_set_coeff,
 458         num_finalize,
 459         num_eval
 460 };
 461
 462 class cl_heap_num_univpoly_ring : public cl_heap_univpoly_ring {
 463         SUBCLASS_cl_heap_univpoly_ring()
 464 public:
 465         // Constructor.
 466         cl_heap_num_univpoly_ring (const cl_ring& r);
 467         // Destructor.
 468         ~cl_heap_num_univpoly_ring () {}
 469 };
 470
 471 static void cl_heap_num_univpoly_ring_destructor (cl_heap* pointer)
 472 {
 473         (*(cl_heap_num_univpoly_ring*)pointer).~cl_heap_num_univpoly_ring();
 474 }
 475
 476 cl_class cl_class_num_univpoly_ring = {
 477         cl_heap_num_univpoly_ring_destructor,
 478         0
 479 };
 480
 481 // Constructor.
 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)
 484 {
 485         type = &cl_class_num_univpoly_ring;
 486 }
 487
 488 }  // namespace cln