src/polynomial/elem/cl_UP_gen.h

   1 // Univariate Polynomials over a general ring.
   2
   3 #include "cln/SV_ringelt.h"
   4 #include "cln/integer.h"
   5 #include "cln/exception.h"
   6
   7 namespace cln {
   8
   9 // Assume a ring is a ring.
  10   inline cl_heap_ring* TheRing (const cl_ring& R)
  11   { return (cl_heap_ring*) R.heappointer; }
  12
  13 // Normalize a vector: remove leading zero coefficients.
  14 // The result vector is known to have length len > 0.
  15 static inline void gen_normalize (cl_heap_ring* R, cl_SV_ringelt& 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_SV_ringelt newresult = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(len));
  25                 for (var sintL i = len-1; i >= 0; i--)
  26                         init1(_cl_ring_element, newresult[i]) (result[i]);
  27                 result = newresult;
  28         }
  29 }
  30
  31 static void gen_fprint (cl_heap_univpoly_ring* UPR, std::ostream& stream, const _cl_UP& x)
  32 {{
  33         DeclarePoly(cl_SV_ringelt,x);
  34         var cl_heap_ring* R = TheRing(UPR->basering());
  35         var sintL xlen = x.length();
  36         if (xlen == 0)
  37                 fprint(stream, "0");
  38         else {
  39                 var cl_string varname = get_varname(UPR);
  40                 for (var sintL i = xlen-1; i >= 0; i--)
  41                         if (!R->_zerop(x[i])) {
  42                                 if (i < xlen-1)
  43                                         fprint(stream, " + ");
  44                                 fprint(stream, "(");
  45                                 R->_fprint(stream, x[i]);
  46                                 fprint(stream, ")");
  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 bool gen_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  60 {{
  61         DeclarePoly(cl_SV_ringelt,x);
  62         DeclarePoly(cl_SV_ringelt,y);
  63         var cl_heap_ring* R = TheRing(UPR->basering());
  64         var sintL xlen = x.length();
  65         var sintL ylen = y.length();
  66         if (!(xlen == ylen))
  67                 return false;
  68         for (var sintL i = xlen-1; i >= 0; i--)
  69                 if (!R->_equal(x[i],y[i]))
  70                         return false;
  71         return true;
  72 }}
  73
  74 static const _cl_UP gen_zero (cl_heap_univpoly_ring* UPR)
  75 {
  76         return _cl_UP(UPR, cl_null_SV_ringelt);
  77 }
  78
  79 static bool gen_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  80 {
  81         unused UPR;
  82  {      DeclarePoly(cl_SV_ringelt,x);
  83         var sintL xlen = x.length();
  84         if (xlen == 0)
  85                 return true;
  86         else
  87                 return false;
  88 }}
  89
  90 static const _cl_UP gen_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  91 {{
  92         DeclarePoly(cl_SV_ringelt,x);
  93         DeclarePoly(cl_SV_ringelt,y);
  94         var cl_heap_ring* R = TheRing(UPR->basering());
  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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 104                 var sintL i;
 105                 for (i = xlen-1; i >= ylen; i--)
 106                         init1(_cl_ring_element, result[i]) (x[i]);
 107                 for (i = ylen-1; i >= 0; i--)
 108                         init1(_cl_ring_element, result[i]) (R->_plus(x[i],y[i]));
 109                 return _cl_UP(UPR, result);
 110         }
 111         if (xlen < ylen) {
 112                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 113                 var sintL i;
 114                 for (i = ylen-1; i >= xlen; i--)
 115                         init1(_cl_ring_element, result[i]) (y[i]);
 116                 for (i = xlen-1; i >= 0; i--)
 117                         init1(_cl_ring_element, result[i]) (R->_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_ring_element hicoeff = R->_plus(x[i],y[i]);
 123                 if (!R->_zerop(hicoeff)) {
 124                         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(i+1));
 125                         init1(_cl_ring_element, result[i]) (hicoeff);
 126                         for (i-- ; i >= 0; i--)
 127                                 init1(_cl_ring_element, result[i]) (R->_plus(x[i],y[i]));
 128                         return _cl_UP(UPR, result);
 129                 }
 130         }
 131         return _cl_UP(UPR, cl_null_SV_ringelt);
 132 }}
 133
 134 static const _cl_UP gen_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 135 {{
 136         DeclarePoly(cl_SV_ringelt,x);
 137         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ring_element hicoeff = R->_uminus(x[i]);
 145         if (R->_zerop(hicoeff)) throw runtime_exception();
 146         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 147         init1(_cl_ring_element, result[i]) (hicoeff);
 148         for (i-- ; i >= 0; i--)
 149                 init1(_cl_ring_element, result[i]) (R->_uminus(x[i]));
 150         return _cl_UP(UPR, result);
 151 }}
 152
 153 static const _cl_UP gen_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 154 {{
 155         DeclarePoly(cl_SV_ringelt,x);
 156         DeclarePoly(cl_SV_ringelt,y);
 157         var cl_heap_ring* R = TheRing(UPR->basering());
 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 gen_uminus(UPR,_cl_UP(UPR, y));
 164         // Now xlen > 0, ylen > 0.
 165         if (xlen > ylen) {
 166                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 167                 var sintL i;
 168                 for (i = xlen-1; i >= ylen; i--)
 169                         init1(_cl_ring_element, result[i]) (x[i]);
 170                 for (i = ylen-1; i >= 0; i--)
 171                         init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 172                 return _cl_UP(UPR, result);
 173         }
 174         if (xlen < ylen) {
 175                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 176                 var sintL i;
 177                 for (i = ylen-1; i >= xlen; i--)
 178                         init1(_cl_ring_element, result[i]) (R->_uminus(y[i]));
 179                 for (i = xlen-1; i >= 0; i--)
 180                         init1(_cl_ring_element, result[i]) (R->_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_ring_element hicoeff = R->_minus(x[i],y[i]);
 186                 if (!R->_zerop(hicoeff)) {
 187                         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(i+1));
 188                         init1(_cl_ring_element, result[i]) (hicoeff);
 189                         for (i-- ; i >= 0; i--)
 190                                 init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 191                         return _cl_UP(UPR, result);
 192                 }
 193         }
 194         return _cl_UP(UPR, cl_null_SV_ringelt);
 195 }}
 196
 197 static const _cl_UP gen_one (cl_heap_univpoly_ring* UPR)
 198 {
 199         var cl_heap_ring* R = TheRing(UPR->basering());
 200         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(1));
 201         init1(_cl_ring_element, result[0]) (R->_one());
 202         return _cl_UP(UPR, result);
 203 }
 204
 205 static const _cl_UP gen_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
 206 {
 207         var cl_heap_ring* R = TheRing(UPR->basering());
 208         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(1));
 209         init1(_cl_ring_element, result[0]) (R->_canonhom(x));
 210         return _cl_UP(UPR, result);
 211 }
 212
 213 static const _cl_UP gen_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 214 {{
 215         DeclarePoly(cl_SV_ringelt,x);
 216         DeclarePoly(cl_SV_ringelt,y);
 217         var cl_heap_ring* R = TheRing(UPR->basering());
 218         var sintL xlen = x.length();
 219         var sintL ylen = y.length();
 220         if (xlen == 0)
 221                 return _cl_UP(UPR, x);
 222         if (ylen == 0)
 223                 return _cl_UP(UPR, y);
 224         // Multiply.
 225         var sintL len = xlen + ylen - 1;
 226         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(len));
 227         if (xlen < ylen) {
 228                 {
 229                         var sintL i = xlen-1;
 230                         var _cl_ring_element xi = x[i];
 231                         for (sintL j = ylen-1; j >= 0; j--)
 232                                 init1(_cl_ring_element, result[i+j]) (R->_mul(xi,y[j]));
 233                 }
 234                 for (sintL i = xlen-2; i >= 0; i--) {
 235                         var _cl_ring_element xi = x[i];
 236                         for (sintL j = ylen-1; j > 0; j--)
 237                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,y[j]));
 238                         /* j=0 */ init1(_cl_ring_element, result[i]) (R->_mul(xi,y[0]));
 239                 }
 240         } else {
 241                 {
 242                         var sintL j = ylen-1;
 243                         var _cl_ring_element yj = y[j];
 244                         for (sintL i = xlen-1; i >= 0; i--)
 245                                 init1(_cl_ring_element, result[i+j]) (R->_mul(x[i],yj));
 246                 }
 247                 for (sintL j = ylen-2; j >= 0; j--) {
 248                         var _cl_ring_element yj = y[j];
 249                         for (sintL i = xlen-1; i > 0; i--)
 250                                 result[i+j] = R->_plus(result[i+j],R->_mul(x[i],yj));
 251                         /* i=0 */ init1(_cl_ring_element, result[j]) (R->_mul(x[0],yj));
 252                 }
 253         }
 254         // Normalize (not necessary in integral domains).
 255         //gen_normalize(R,result,len);
 256         if (R->_zerop(result[len-1])) throw runtime_exception();
 257         return _cl_UP(UPR, result);
 258 }}
 259
 260 static const _cl_UP gen_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 261 {{
 262         DeclarePoly(cl_SV_ringelt,x);
 263         var cl_heap_ring* R = TheRing(UPR->basering());
 264         var sintL xlen = x.length();
 265         if (xlen == 0)
 266                 return cl_UP(UPR, x);
 267         var sintL len = 2*xlen-1;
 268         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(len));
 269         if (xlen > 1) {
 270                 // Loop through all 0 <= j < i <= xlen-1.
 271                 {
 272                         var sintL i = xlen-1;
 273                         var _cl_ring_element xi = x[i];
 274                         for (sintL j = i-1; j >= 0; j--)
 275                                 init1(_cl_ring_element, result[i+j]) (R->_mul(xi,x[j]));
 276                 }
 277                 {for (sintL i = xlen-2; i >= 1; i--) {
 278                         var _cl_ring_element xi = x[i];
 279                         for (sintL j = i-1; j >= 1; j--)
 280                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,x[j]));
 281                         /* j=0 */ init1(_cl_ring_element, result[i]) (R->_mul(xi,x[0]));
 282                 }}
 283                 // Double.
 284                 {for (sintL i = len-2; i >= 1; i--)
 285                         result[i] = R->_plus(result[i],result[i]);
 286                 }
 287                 // Add squares.
 288                 init1(_cl_ring_element, result[2*(xlen-1)]) (R->_square(x[xlen-1]));
 289                 for (sintL i = xlen-2; i >= 1; i--)
 290                         result[2*i] = R->_plus(result[2*i],R->_square(x[i]));
 291         }
 292         init1(_cl_ring_element, result[0]) (R->_square(x[0]));
 293         // Normalize (not necessary in integral domains).
 294         //gen_normalize(R,result,len);
 295         if (R->_zerop(result[len-1])) throw runtime_exception();
 296         return _cl_UP(UPR, result);
 297 }}
 298
 299 static const _cl_UP gen_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
 300 {
 301         var _cl_UP a = x;
 302         var cl_I b = y;
 303         while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
 304         var _cl_UP c = a;
 305         until (b == 1)
 306           { b = b >> 1;
 307             a = UPR->_square(a);
 308             if (oddp(b)) { c = UPR->_mul(a,c); }
 309           }
 310         return c;
 311 }
 312
 313 static const _cl_UP gen_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
 314 {
 315         if (!(UPR->basering() == x.ring())) throw runtime_exception();
 316  {
 317         DeclarePoly(cl_SV_ringelt,y);
 318         var cl_heap_ring* R = TheRing(UPR->basering());
 319         var sintL ylen = y.length();
 320         if (ylen == 0)
 321                 return _cl_UP(UPR, y);
 322         if (R->zerop(x))
 323                 return _cl_UP(UPR, cl_null_SV_ringelt);
 324         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 325         for (sintL i = ylen-1; i >= 0; i--)
 326                 init1(_cl_ring_element, result[i]) (R->_mul(x,y[i]));
 327         // Normalize (not necessary in integral domains).
 328         //gen_normalize(R,result,ylen);
 329         if (R->_zerop(result[ylen-1])) throw runtime_exception();
 330         return _cl_UP(UPR, result);
 331 }}
 332
 333 static sintL gen_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 334 {
 335         unused UPR;
 336  {      DeclarePoly(cl_SV_ringelt,x);
 337         return (sintL) x.length() - 1;
 338 }}
 339
 340 static sintL gen_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 341 {{      DeclarePoly(cl_SV_ringelt,x);
 342         var cl_heap_ring* R = TheRing(UPR->basering());
 343         var sintL xlen = x.length();
 344         for (sintL i = 0; i < xlen; i++) {
 345                 if (!R->_zerop(x[i]))
 346                         return i;
 347         }
 348         return -1;
 349 }}
 350
 351 static const _cl_UP gen_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 352 {
 353         if (!(UPR->basering() == x.ring())) throw runtime_exception();
 354         var cl_heap_ring* R = TheRing(UPR->basering());
 355         if (R->_zerop(x))
 356                 return _cl_UP(UPR, cl_null_SV_ringelt);
 357         else {
 358                 var sintL len = e+1;
 359                 var cl_SV_ringelt result = cl_SV_ringelt(len);
 360                 result[e] = x;
 361                 return _cl_UP(UPR, result);
 362         }
 363 }
 364
 365 static const cl_ring_element gen_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 366 {{
 367         DeclarePoly(cl_SV_ringelt,x);
 368         var cl_heap_ring* R = TheRing(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 gen_create (cl_heap_univpoly_ring* UPR, sintL deg)
 376 {
 377         if (deg < 0)
 378                 return _cl_UP(UPR, cl_null_SV_ringelt);
 379         else {
 380                 var sintL len = deg+1;
 381                 return _cl_UP(UPR, cl_SV_ringelt(len));
 382         }
 383 }
 384
 385 static void gen_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 386 {{
 387         DeclareMutablePoly(cl_SV_ringelt,x);
 388         if (!(UPR->basering() == y.ring())) throw runtime_exception();
 389         if (!(index < x.length())) throw runtime_exception();
 390         x[index] = y;
 391 }}
 392
 393 static void gen_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 394 {{
 395         DeclareMutablePoly(cl_SV_ringelt,x); // NB: x is modified by reference!
 396         var cl_heap_ring* R = TheRing(UPR->basering());
 397         var uintL len = x.length();
 398         if (len > 0)
 399                 gen_normalize(R,x,len);
 400 }}
 401
 402 static const cl_ring_element gen_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 403 {{
 404         // Method:
 405         // If x = 0, return 0.
 406         // If y = 0, return x[0].
 407         // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
 408         DeclarePoly(cl_SV_ringelt,x);
 409         var cl_heap_ring* R = TheRing(UPR->basering());
 410         if (!(y.ring() == R)) throw runtime_exception();
 411         var uintL len = x.length();
 412         if (len==0)
 413                 return R->zero();
 414         if (R->_zerop(y))
 415                 return cl_ring_element(R, x[0]);
 416         var sintL i = len-1;
 417         var _cl_ring_element z = x[i];
 418         for ( ; --i >= 0; )
 419                 z = R->_plus(R->_mul(z,y),x[i]);
 420         return cl_ring_element(R, z);
 421 }}
 422
 423 static cl_univpoly_setops gen_setops = {
 424         gen_fprint,
 425         gen_equal
 426 };
 427
 428 static cl_univpoly_addops gen_addops = {
 429         gen_zero,
 430         gen_zerop,
 431         gen_plus,
 432         gen_minus,
 433         gen_uminus
 434 };
 435
 436 static cl_univpoly_mulops gen_mulops = {
 437         gen_one,
 438         gen_canonhom,
 439         gen_mul,
 440         gen_square,
 441         gen_exptpos
 442 };
 443
 444 static cl_univpoly_modulops gen_modulops = {
 445         gen_scalmul
 446 };
 447
 448 static cl_univpoly_polyops gen_polyops = {
 449         gen_degree,
 450         gen_ldegree,
 451         gen_monomial,
 452         gen_coeff,
 453         gen_create,
 454         gen_set_coeff,
 455         gen_finalize,
 456         gen_eval
 457 };
 458
 459 class cl_heap_gen_univpoly_ring : public cl_heap_univpoly_ring {
 460         SUBCLASS_cl_heap_univpoly_ring()
 461 public:
 462         // Constructor.
 463         cl_heap_gen_univpoly_ring (const cl_ring& r);
 464         // Destructor
 465         ~cl_heap_gen_univpoly_ring () {}
 466 };
 467
 468 static void cl_heap_gen_univpoly_ring_destructor (cl_heap* pointer)
 469 {
 470         (*(cl_heap_gen_univpoly_ring*)pointer).~cl_heap_gen_univpoly_ring();
 471 }
 472
 473 cl_class cl_class_gen_univpoly_ring = {
 474         cl_heap_gen_univpoly_ring_destructor,
 475         cl_class_flags_univpoly_ring
 476 };
 477
 478 // Constructor.
 479 inline cl_heap_gen_univpoly_ring::cl_heap_gen_univpoly_ring (const cl_ring& r)
 480         : cl_heap_univpoly_ring (r, &gen_setops, &gen_addops, &gen_mulops, &gen_modulops, &gen_polyops)
 481 {
 482         type = &cl_class_gen_univpoly_ring;
 483 }
 484
 485 }  // namespace cln