src/polynomial/elem/cl_UP_gen.h

   1 // Univariate Polynomials over a general ring.
   2
   3 #include "cl_SV_ringelt.h"
   4 #include "cl_integer.h"
   5 #include "cl_abort.h"
   6
   7 // Assume a ring is a ring.
   8   inline cl_heap_ring* TheRing (const cl_ring& R)
   9   { return (cl_heap_ring*) R.heappointer; }
  10
  11 // Normalize a vector: remove leading zero coefficients.
  12 // The result vector is known to have length len > 0.
  13 static inline void gen_normalize (cl_heap_ring* R, cl_SV_ringelt& result, uintL len)
  14 {
  15         if (R->_zerop(result[len-1])) {
  16                 len--;
  17                 while (len > 0) {
  18                         if (!R->_zerop(result[len-1]))
  19                                 break;
  20                         len--;
  21                 }
  22                 var cl_SV_ringelt newresult = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(len));
  23                 for (var sintL i = len-1; i >= 0; i--)
  24                         init1(_cl_ring_element, newresult[i]) (result[i]);
  25                 result = newresult;
  26         }
  27 }
  28
  29 static void gen_fprint (cl_heap_univpoly_ring* UPR, cl_ostream stream, const _cl_UP& x)
  30 {{
  31         DeclarePoly(cl_SV_ringelt,x);
  32         var cl_heap_ring* R = TheRing(UPR->basering());
  33         var sintL xlen = x.length();
  34         if (xlen == 0)
  35                 fprint(stream, "0");
  36         else {
  37                 var cl_string varname = get_varname(UPR);
  38                 for (var sintL i = xlen-1; i >= 0; i--)
  39                         if (!R->_zerop(x[i])) {
  40                                 if (i < xlen-1)
  41                                         fprint(stream, " + ");
  42                                 fprint(stream, "(");
  43                                 R->_fprint(stream, x[i]);
  44                                 fprint(stream, ")");
  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 gen_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  58 {{
  59         DeclarePoly(cl_SV_ringelt,x);
  60         DeclarePoly(cl_SV_ringelt,y);
  61         var cl_heap_ring* R = TheRing(UPR->basering());
  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 (!R->_equal(x[i],y[i]))
  68                         return cl_false;
  69         return cl_true;
  70 }}
  71
  72 static const _cl_UP gen_zero (cl_heap_univpoly_ring* UPR)
  73 {
  74         return _cl_UP(UPR, cl_null_SV_ringelt);
  75 }
  76
  77 static cl_boolean gen_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  78 {
  79         unused UPR;
  80  {      DeclarePoly(cl_SV_ringelt,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 gen_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  89 {{
  90         DeclarePoly(cl_SV_ringelt,x);
  91         DeclarePoly(cl_SV_ringelt,y);
  92         var cl_heap_ring* R = TheRing(UPR->basering());
  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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 102                 var sintL i;
 103                 for (i = xlen-1; i >= ylen; i--)
 104                         init1(_cl_ring_element, result[i]) (x[i]);
 105                 for (i = ylen-1; i >= 0; i--)
 106                         init1(_cl_ring_element, result[i]) (R->_plus(x[i],y[i]));
 107                 return _cl_UP(UPR, result);
 108         }
 109         if (xlen < ylen) {
 110                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 111                 var sintL i;
 112                 for (i = ylen-1; i >= xlen; i--)
 113                         init1(_cl_ring_element, result[i]) (y[i]);
 114                 for (i = xlen-1; i >= 0; i--)
 115                         init1(_cl_ring_element, result[i]) (R->_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_ring_element hicoeff = R->_plus(x[i],y[i]);
 121                 if (!R->_zerop(hicoeff)) {
 122                         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(i+1));
 123                         init1(_cl_ring_element, result[i]) (hicoeff);
 124                         for (i-- ; i >= 0; i--)
 125                                 init1(_cl_ring_element, result[i]) (R->_plus(x[i],y[i]));
 126                         return _cl_UP(UPR, result);
 127                 }
 128         }
 129         return _cl_UP(UPR, cl_null_SV_ringelt);
 130 }}
 131
 132 static const _cl_UP gen_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 133 {{
 134         DeclarePoly(cl_SV_ringelt,x);
 135         DeclarePoly(cl_SV_ringelt,y);
 136         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 146                 var sintL i;
 147                 for (i = xlen-1; i >= ylen; i--)
 148                         init1(_cl_ring_element, result[i]) (x[i]);
 149                 for (i = ylen-1; i >= 0; i--)
 150                         init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 151                 return _cl_UP(UPR, result);
 152         }
 153         if (xlen < ylen) {
 154                 var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 155                 var sintL i;
 156                 for (i = ylen-1; i >= xlen; i--)
 157                         init1(_cl_ring_element, result[i]) (R->_uminus(y[i]));
 158                 for (i = xlen-1; i >= 0; i--)
 159                         init1(_cl_ring_element, result[i]) (R->_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_ring_element hicoeff = R->_minus(x[i],y[i]);
 165                 if (!R->_zerop(hicoeff)) {
 166                         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(i+1));
 167                         init1(_cl_ring_element, result[i]) (hicoeff);
 168                         for (i-- ; i >= 0; i--)
 169                                 init1(_cl_ring_element, result[i]) (R->_minus(x[i],y[i]));
 170                         return _cl_UP(UPR, result);
 171                 }
 172         }
 173         return _cl_UP(UPR, cl_null_SV_ringelt);
 174 }}
 175
 176 static const _cl_UP gen_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 177 {{
 178         DeclarePoly(cl_SV_ringelt,x);
 179         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ring_element hicoeff = R->_uminus(x[i]);
 187         if (R->_zerop(hicoeff)) cl_abort();
 188         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(xlen));
 189         init1(_cl_ring_element, result[i]) (hicoeff);
 190         for (i-- ; i >= 0; i--)
 191                 init1(_cl_ring_element, result[i]) (R->_uminus(x[i]));
 192         return _cl_UP(UPR, result);
 193 }}
 194
 195 static const _cl_UP gen_one (cl_heap_univpoly_ring* UPR)
 196 {
 197         var cl_heap_ring* R = TheRing(UPR->basering());
 198         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(1));
 199         init1(_cl_ring_element, result[0]) (R->_one());
 200         return _cl_UP(UPR, result);
 201 }
 202
 203 static const _cl_UP gen_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
 204 {
 205         var cl_heap_ring* R = TheRing(UPR->basering());
 206         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(1));
 207         init1(_cl_ring_element, result[0]) (R->_canonhom(x));
 208         return _cl_UP(UPR, result);
 209 }
 210
 211 static const _cl_UP gen_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 212 {{
 213         DeclarePoly(cl_SV_ringelt,x);
 214         DeclarePoly(cl_SV_ringelt,y);
 215         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(len));
 225         if (xlen < ylen) {
 226                 {
 227                         var sintL i = xlen-1;
 228                         var _cl_ring_element xi = x[i];
 229                         for (sintL j = ylen-1; j >= 0; j--)
 230                                 init1(_cl_ring_element, result[i+j]) (R->_mul(xi,y[j]));
 231                 }
 232                 for (sintL i = xlen-2; i >= 0; i--) {
 233                         var _cl_ring_element xi = x[i];
 234                         for (sintL j = ylen-1; j > 0; j--)
 235                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,y[j]));
 236                         /* j=0 */ init1(_cl_ring_element, result[i]) (R->_mul(xi,y[0]));
 237                 }
 238         } else {
 239                 {
 240                         var sintL j = ylen-1;
 241                         var _cl_ring_element yj = y[j];
 242                         for (sintL i = xlen-1; i >= 0; i--)
 243                                 init1(_cl_ring_element, result[i+j]) (R->_mul(x[i],yj));
 244                 }
 245                 for (sintL j = ylen-2; j >= 0; j--) {
 246                         var _cl_ring_element yj = y[j];
 247                         for (sintL i = xlen-1; i > 0; i--)
 248                                 result[i+j] = R->_plus(result[i+j],R->_mul(x[i],yj));
 249                         /* i=0 */ init1(_cl_ring_element, result[j]) (R->_mul(x[0],yj));
 250                 }
 251         }
 252         // Normalize (not necessary in integral domains).
 253         //gen_normalize(R,result,len);
 254         if (R->_zerop(result[len-1])) cl_abort();
 255         return _cl_UP(UPR, result);
 256 }}
 257
 258 static const _cl_UP gen_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 259 {{
 260         DeclarePoly(cl_SV_ringelt,x);
 261         var cl_heap_ring* R = TheRing(UPR->basering());
 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_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_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_ring_element xi = x[i];
 272                         for (sintL j = i-1; j >= 0; j--)
 273                                 init1(_cl_ring_element, result[i+j]) (R->_mul(xi,x[j]));
 274                 }
 275                 {for (sintL i = xlen-2; i >= 1; i--) {
 276                         var _cl_ring_element xi = x[i];
 277                         for (sintL j = i-1; j >= 1; j--)
 278                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,x[j]));
 279                         /* j=0 */ init1(_cl_ring_element, result[i]) (R->_mul(xi,x[0]));
 280                 }}
 281                 // Double.
 282                 {for (sintL i = len-2; i >= 1; i--)
 283                         result[i] = R->_plus(result[i],result[i]);
 284                 }
 285                 // Add squares.
 286                 init1(_cl_ring_element, result[2*(xlen-1)]) (R->_square(x[xlen-1]));
 287                 for (sintL i = xlen-2; i >= 1; i--)
 288                         result[2*i] = R->_plus(result[2*i],R->_square(x[i]));
 289         }
 290         init1(_cl_ring_element, result[0]) (R->_square(x[0]));
 291         // Normalize (not necessary in integral domains).
 292         //gen_normalize(R,result,len);
 293         if (R->_zerop(result[len-1])) cl_abort();
 294         return _cl_UP(UPR, result);
 295 }}
 296
 297 static const _cl_UP gen_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 gen_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_SV_ringelt,y);
 316         var cl_heap_ring* R = TheRing(UPR->basering());
 317         var sintL ylen = y.length();
 318         if (ylen == 0)
 319                 return _cl_UP(UPR, y);
 320         if (R->zerop(x))
 321                 return _cl_UP(UPR, cl_null_SV_ringelt);
 322         var cl_SV_ringelt result = cl_SV_ringelt(cl_make_heap_SV_ringelt_uninit(ylen));
 323         for (sintL i = ylen-1; i >= 0; i--)
 324                 init1(_cl_ring_element, result[i]) (R->_mul(x,y[i]));
 325         // Normalize (not necessary in integral domains).
 326         //gen_normalize(R,result,ylen);
 327         if (R->_zerop(result[ylen-1])) cl_abort();
 328         return _cl_UP(UPR, result);
 329 }}
 330
 331 static sintL gen_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 332 {
 333         unused UPR;
 334  {      DeclarePoly(cl_SV_ringelt,x);
 335         return (sintL) x.length() - 1;
 336 }}
 337
 338 static const _cl_UP gen_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 339 {
 340         if (!(UPR->basering() == x.ring())) cl_abort();
 341         var cl_heap_ring* R = TheRing(UPR->basering());
 342         if (R->_zerop(x))
 343                 return _cl_UP(UPR, cl_null_SV_ringelt);
 344         else {
 345                 var sintL len = e+1;
 346                 var cl_SV_ringelt result = cl_SV_ringelt(len);
 347                 result[e] = x;
 348                 return _cl_UP(UPR, result);
 349         }
 350 }
 351
 352 static const cl_ring_element gen_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 353 {{
 354         DeclarePoly(cl_SV_ringelt,x);
 355         var cl_heap_ring* R = TheRing(UPR->basering());
 356         if (index < x.length())
 357                 return cl_ring_element(R, x[index]);
 358         else
 359                 return R->zero();
 360 }}
 361
 362 static const _cl_UP gen_create (cl_heap_univpoly_ring* UPR, sintL deg)
 363 {
 364         if (deg < 0)
 365                 return _cl_UP(UPR, cl_null_SV_ringelt);
 366         else {
 367                 var sintL len = deg+1;
 368                 return _cl_UP(UPR, cl_SV_ringelt(len));
 369         }
 370 }
 371
 372 static void gen_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 373 {{
 374         DeclareMutablePoly(cl_SV_ringelt,x);
 375         if (!(UPR->basering() == y.ring())) cl_abort();
 376         if (!(index < x.length())) cl_abort();
 377         x[index] = y;
 378 }}
 379
 380 static void gen_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 381 {{
 382         DeclareMutablePoly(cl_SV_ringelt,x); // NB: x is modified by reference!
 383         var cl_heap_ring* R = TheRing(UPR->basering());
 384         var uintL len = x.length();
 385         if (len > 0)
 386                 gen_normalize(R,x,len);
 387 }}
 388
 389 static const cl_ring_element gen_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_ringelt,x);
 396         var cl_heap_ring* R = TheRing(UPR->basering());
 397         if (!(y.ring() == R)) cl_abort();
 398         var uintL len = x.length();
 399         if (len==0)
 400                 return R->zero();
 401         if (R->_zerop(y))
 402                 return cl_ring_element(R, x[0]);
 403         var sintL i = len-1;
 404         var _cl_ring_element z = x[i];
 405         for ( ; --i >= 0; )
 406                 z = R->_plus(R->_mul(z,y),x[i]);
 407         return cl_ring_element(R, z);
 408 }}
 409
 410 static cl_univpoly_setops gen_setops = {
 411         gen_fprint,
 412         gen_equal
 413 };
 414
 415 static cl_univpoly_addops gen_addops = {
 416         gen_zero,
 417         gen_zerop,
 418         gen_plus,
 419         gen_minus,
 420         gen_uminus
 421 };
 422
 423 static cl_univpoly_mulops gen_mulops = {
 424         gen_one,
 425         gen_canonhom,
 426         gen_mul,
 427         gen_square,
 428         gen_exptpos
 429 };
 430
 431 static cl_univpoly_modulops gen_modulops = {
 432         gen_scalmul
 433 };
 434
 435 static cl_univpoly_polyops gen_polyops = {
 436         gen_degree,
 437         gen_monomial,
 438         gen_coeff,
 439         gen_create,
 440         gen_set_coeff,
 441         gen_finalize,
 442         gen_eval
 443 };
 444
 445 class cl_heap_gen_univpoly_ring : public cl_heap_univpoly_ring {
 446         SUBCLASS_cl_heap_univpoly_ring()
 447 public:
 448         // Constructor.
 449         cl_heap_gen_univpoly_ring (const cl_ring& r)
 450                 : cl_heap_univpoly_ring (r, &gen_setops, &gen_addops, &gen_mulops, &gen_modulops, &gen_polyops) {}
 451 };