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/abort.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 cl_boolean 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 cl_false;
  68         for (var sintL i = xlen-1; i >= 0; i--)
  69                 if (!R->_equal(x[i],y[i]))
  70                         return cl_false;
  71         return cl_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 cl_boolean 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 cl_true;
  86         else
  87                 return cl_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_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 135 {{
 136         DeclarePoly(cl_SV_ringelt,x);
 137         DeclarePoly(cl_SV_ringelt,y);
 138         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 148                 var sintL i;
 149                 for (i = xlen-1; i >= ylen; i--)
 150                         init1(_cl_ring_element, result[i]) (x[i]);
 151                 for (i = ylen-1; i >= 0; i--)
 152                         init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 153                 return _cl_UP(UPR, result);
 154         }
 155         if (xlen < ylen) {
 156                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 157                 var sintL i;
 158                 for (i = ylen-1; i >= xlen; i--)
 159                         init1(_cl_ring_element, result[i]) (R->_uminus(y[i]));
 160                 for (i = xlen-1; i >= 0; i--)
 161                         init1(_cl_ring_element, result[i]) (R->_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_ring_element hicoeff = R->_minus(x[i],y[i]);
 167                 if (!R->_zerop(hicoeff)) {
 168                         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(i+1));
 169                         init1(_cl_ring_element, result[i]) (hicoeff);
 170                         for (i-- ; i >= 0; i--)
 171                                 init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 172                         return _cl_UP(UPR, result);
 173                 }
 174         }
 175         return _cl_UP(UPR, cl_null_SV_ringelt);
 176 }}
 177
 178 static const _cl_UP gen_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 179 {{
 180         DeclarePoly(cl_SV_ringelt,x);
 181         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ring_element hicoeff = R->_uminus(x[i]);
 189         if (R->_zerop(hicoeff)) cl_abort();
 190         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 191         init1(_cl_ring_element, result[i]) (hicoeff);
 192         for (i-- ; i >= 0; i--)
 193                 init1(_cl_ring_element, result[i]) (R->_uminus(x[i]));
 194         return _cl_UP(UPR, result);
 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])) cl_abort();
 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])) cl_abort();
 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())) cl_abort();
 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])) cl_abort();
 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 const _cl_UP gen_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 341 {
 342         if (!(UPR->basering() == x.ring())) cl_abort();
 343         var cl_heap_ring* R = TheRing(UPR->basering());
 344         if (R->_zerop(x))
 345                 return _cl_UP(UPR, cl_null_SV_ringelt);
 346         else {
 347                 var sintL len = e+1;
 348                 var cl_SV_ringelt result = cl_SV_ringelt(len);
 349                 result[e] = x;
 350                 return _cl_UP(UPR, result);
 351         }
 352 }
 353
 354 static const cl_ring_element gen_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 355 {{
 356         DeclarePoly(cl_SV_ringelt,x);
 357         var cl_heap_ring* R = TheRing(UPR->basering());
 358         if (index < x.length())
 359                 return cl_ring_element(R, x[index]);
 360         else
 361                 return R->zero();
 362 }}
 363
 364 static const _cl_UP gen_create (cl_heap_univpoly_ring* UPR, sintL deg)
 365 {
 366         if (deg < 0)
 367                 return _cl_UP(UPR, cl_null_SV_ringelt);
 368         else {
 369                 var sintL len = deg+1;
 370                 return _cl_UP(UPR, cl_SV_ringelt(len));
 371         }
 372 }
 373
 374 static void gen_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 375 {{
 376         DeclareMutablePoly(cl_SV_ringelt,x);
 377         if (!(UPR->basering() == y.ring())) cl_abort();
 378         if (!(index < x.length())) cl_abort();
 379         x[index] = y;
 380 }}
 381
 382 static void gen_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 383 {{
 384         DeclareMutablePoly(cl_SV_ringelt,x); // NB: x is modified by reference!
 385         var cl_heap_ring* R = TheRing(UPR->basering());
 386         var uintL len = x.length();
 387         if (len > 0)
 388                 gen_normalize(R,x,len);
 389 }}
 390
 391 static const cl_ring_element gen_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_ringelt,x);
 398         var cl_heap_ring* R = TheRing(UPR->basering());
 399         if (!(y.ring() == R)) cl_abort();
 400         var uintL len = x.length();
 401         if (len==0)
 402                 return R->zero();
 403         if (R->_zerop(y))
 404                 return cl_ring_element(R, x[0]);
 405         var sintL i = len-1;
 406         var _cl_ring_element z = x[i];
 407         for ( ; --i >= 0; )
 408                 z = R->_plus(R->_mul(z,y),x[i]);
 409         return cl_ring_element(R, z);
 410 }}
 411
 412 static cl_univpoly_setops gen_setops = {
 413         gen_fprint,
 414         gen_equal
 415 };
 416
 417 static cl_univpoly_addops gen_addops = {
 418         gen_zero,
 419         gen_zerop,
 420         gen_plus,
 421         gen_minus,
 422         gen_uminus
 423 };
 424
 425 static cl_univpoly_mulops gen_mulops = {
 426         gen_one,
 427         gen_canonhom,
 428         gen_mul,
 429         gen_square,
 430         gen_exptpos
 431 };
 432
 433 static cl_univpoly_modulops gen_modulops = {
 434         gen_scalmul
 435 };
 436
 437 static cl_univpoly_polyops gen_polyops = {
 438         gen_degree,
 439         gen_monomial,
 440         gen_coeff,
 441         gen_create,
 442         gen_set_coeff,
 443         gen_finalize,
 444         gen_eval
 445 };
 446
 447 class cl_heap_gen_univpoly_ring : public cl_heap_univpoly_ring {
 448         SUBCLASS_cl_heap_univpoly_ring()
 449 public:
 450         // Constructor.
 451         cl_heap_gen_univpoly_ring (const cl_ring& r)
 452                 : cl_heap_univpoly_ring (r, &gen_setops, &gen_addops, &gen_mulops, &gen_modulops, &gen_polyops) {}
 453 };
 454
 455 }  // namespace cln