src/polynomial/elem/cl_UP_MI.h

   1 // Univariate Polynomials over a ring of modular integers.
   2
   3 #include "cln/GV_modinteger.h"
   4 #include "cln/modinteger.h"
   5 #include "cln/abort.h"
   6
   7 namespace cln {
   8
   9 // Assume a ring is a modint ring.
  10   inline cl_heap_modint_ring* TheModintRing (const cl_ring& R)
  11   { return (cl_heap_modint_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 modint_normalize (cl_heap_modint_ring* R, cl_GV_MI& 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_GV_MI newresult = cl_GV_MI(len,R);
  25                 #if 0
  26                 for (var sintL i = len-1; i >= 0; i--)
  27                         newresult[i] = result[i];
  28                 #else
  29                 cl_GV_MI::copy_elements(result,0,newresult,0,len);
  30                 #endif
  31                 result = newresult;
  32         }
  33 }
  34
  35 static void modint_fprint (cl_heap_univpoly_ring* UPR, cl_ostream stream, const _cl_UP& x)
  36 {{
  37         DeclarePoly(cl_GV_MI,x);
  38         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  39         var sintL xlen = x.length();
  40         if (xlen == 0)
  41                 fprint(stream, "0");
  42         else {
  43                 var cl_string varname = get_varname(UPR);
  44                 for (var sintL i = xlen-1; i >= 0; i--)
  45                         if (!R->_zerop(x[i])) {
  46                                 if (i < xlen-1)
  47                                         fprint(stream, " + ");
  48                                 fprint(stream, "(");
  49                                 R->_fprint(stream, x[i]);
  50                                 fprint(stream, ")");
  51                                 if (i > 0) {
  52                                         fprint(stream, "*");
  53                                         fprint(stream, varname);
  54                                         if (i != 1) {
  55                                                 fprint(stream, "^");
  56                                                 fprintdecimal(stream, i);
  57                                         }
  58                                 }
  59                         }
  60         }
  61 }}
  62
  63 static cl_boolean modint_equal (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  64 {{
  65         DeclarePoly(cl_GV_MI,x);
  66         DeclarePoly(cl_GV_MI,y);
  67         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  68         var sintL xlen = x.length();
  69         var sintL ylen = y.length();
  70         if (!(xlen == ylen))
  71                 return cl_false;
  72         for (var sintL i = xlen-1; i >= 0; i--)
  73                 if (!R->_equal(x[i],y[i]))
  74                         return cl_false;
  75         return cl_true;
  76 }}
  77
  78 static const _cl_UP modint_zero (cl_heap_univpoly_ring* UPR)
  79 {
  80         return _cl_UP(UPR, cl_null_GV_I);
  81 }
  82
  83 static cl_boolean modint_zerop (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
  84 {
  85         unused UPR;
  86  {      DeclarePoly(cl_GV_MI,x);
  87         var sintL xlen = x.length();
  88         if (xlen == 0)
  89                 return cl_true;
  90         else
  91                 return cl_false;
  92 }}
  93
  94 static const _cl_UP modint_plus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
  95 {{
  96         DeclarePoly(cl_GV_MI,x);
  97         DeclarePoly(cl_GV_MI,y);
  98         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
  99         var sintL xlen = x.length();
 100         var sintL ylen = y.length();
 101         if (xlen == 0)
 102                 return _cl_UP(UPR, y);
 103         if (ylen == 0)
 104                 return _cl_UP(UPR, x);
 105         // Now xlen > 0, ylen > 0.
 106         if (xlen > ylen) {
 107                 var cl_GV_MI result = cl_GV_MI(xlen,R);
 108                 var sintL i;
 109                 #if 0
 110                 for (i = xlen-1; i >= ylen; i--)
 111                         result[i] = x[i];
 112                 #else
 113                 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
 114                 #endif
 115                 for (i = ylen-1; i >= 0; i--)
 116                         result[i] = R->_plus(x[i],y[i]);
 117                 return _cl_UP(UPR, result);
 118         }
 119         if (xlen < ylen) {
 120                 var cl_GV_MI result = cl_GV_MI(ylen,R);
 121                 var sintL i;
 122                 #if 0
 123                 for (i = ylen-1; i >= xlen; i--)
 124                         result[i] = y[i];
 125                 #else
 126                 cl_GV_MI::copy_elements(y,xlen,result,xlen,ylen-xlen);
 127                 #endif
 128                 for (i = xlen-1; i >= 0; i--)
 129                         result[i] = R->_plus(x[i],y[i]);
 130                 return _cl_UP(UPR, result);
 131         }
 132         // Now xlen = ylen > 0. Add and normalize simultaneously.
 133         for (var sintL i = xlen-1; i >= 0; i--) {
 134                 var _cl_MI hicoeff = R->_plus(x[i],y[i]);
 135                 if (!R->_zerop(hicoeff)) {
 136                         var cl_GV_MI result = cl_GV_MI(i+1,R);
 137                         result[i] = hicoeff;
 138                         for (i-- ; i >= 0; i--)
 139                                 result[i] = R->_plus(x[i],y[i]);
 140                         return _cl_UP(UPR, result);
 141                 }
 142         }
 143         return _cl_UP(UPR, cl_null_GV_I);
 144 }}
 145
 146 static const _cl_UP modint_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 147 {{
 148         DeclarePoly(cl_GV_MI,x);
 149         DeclarePoly(cl_GV_MI,y);
 150         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 151         var sintL xlen = x.length();
 152         var sintL ylen = y.length();
 153         if (xlen == 0)
 154                 return _cl_UP(UPR, y);
 155         if (ylen == 0)
 156                 return _cl_UP(UPR, x);
 157         // Now xlen > 0, ylen > 0.
 158         if (xlen > ylen) {
 159                 var cl_GV_MI result = cl_GV_MI(xlen,R);
 160                 var sintL i;
 161                 #if 0
 162                 for (i = xlen-1; i >= ylen; i--)
 163                         result[i] = x[i];
 164                 #else
 165                 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
 166                 #endif
 167                 for (i = ylen-1; i >= 0; i--)
 168                         result[i] = R->_minus(x[i],y[i]);
 169                 return _cl_UP(UPR, result);
 170         }
 171         if (xlen < ylen) {
 172                 var cl_GV_MI result = cl_GV_MI(ylen,R);
 173                 var sintL i;
 174                 for (i = ylen-1; i >= xlen; i--)
 175                         result[i] = R->_uminus(y[i]);
 176                 for (i = xlen-1; i >= 0; i--)
 177                         result[i] = R->_minus(x[i],y[i]);
 178                 return _cl_UP(UPR, result);
 179         }
 180         // Now xlen = ylen > 0. Add and normalize simultaneously.
 181         for (var sintL i = xlen-1; i >= 0; i--) {
 182                 var _cl_MI hicoeff = R->_minus(x[i],y[i]);
 183                 if (!R->_zerop(hicoeff)) {
 184                         var cl_GV_MI result = cl_GV_MI(i+1,R);
 185                         result[i] = hicoeff;
 186                         for (i-- ; i >= 0; i--)
 187                                 result[i] = R->_minus(x[i],y[i]);
 188                         return _cl_UP(UPR, result);
 189                 }
 190         }
 191         return _cl_UP(UPR, cl_null_GV_I);
 192 }}
 193
 194 static const _cl_UP modint_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 195 {{
 196         DeclarePoly(cl_GV_MI,x);
 197         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 198         var sintL xlen = x.length();
 199         if (xlen == 0)
 200                 return _cl_UP(UPR, x);
 201         // Now xlen > 0.
 202         // Negate. No normalization necessary, since the degree doesn't change.
 203         var sintL i = xlen-1;
 204         var _cl_MI hicoeff = R->_uminus(x[i]);
 205         if (R->_zerop(hicoeff)) cl_abort();
 206         var cl_GV_MI result = cl_GV_MI(xlen,R);
 207         result[i] = hicoeff;
 208         for (i-- ; i >= 0; i--)
 209                 result[i] = R->_uminus(x[i]);
 210         return _cl_UP(UPR, result);
 211 }}
 212
 213 static const _cl_UP modint_one (cl_heap_univpoly_ring* UPR)
 214 {
 215         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 216         var cl_GV_MI result = cl_GV_MI(1,R);
 217         result[0] = R->_one();
 218         return _cl_UP(UPR, result);
 219 }
 220
 221 static const _cl_UP modint_canonhom (cl_heap_univpoly_ring* UPR, const cl_I& x)
 222 {
 223         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 224         var cl_GV_MI result = cl_GV_MI(1,R);
 225         result[0] = R->_canonhom(x);
 226         return _cl_UP(UPR, result);
 227 }
 228
 229 static const _cl_UP modint_mul (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 230 {{
 231         DeclarePoly(cl_GV_MI,x);
 232         DeclarePoly(cl_GV_MI,y);
 233         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 234         var sintL xlen = x.length();
 235         var sintL ylen = y.length();
 236         if (xlen == 0)
 237                 return _cl_UP(UPR, x);
 238         if (ylen == 0)
 239                 return _cl_UP(UPR, y);
 240         // Multiply.
 241         var sintL len = xlen + ylen - 1;
 242         var cl_GV_MI result = cl_GV_MI(len,R);
 243         if (xlen < ylen) {
 244                 {
 245                         var sintL i = xlen-1;
 246                         var _cl_MI xi = x[i];
 247                         for (sintL j = ylen-1; j >= 0; j--)
 248                                 result[i+j] = R->_mul(xi,y[j]);
 249                 }
 250                 for (sintL i = xlen-2; i >= 0; i--) {
 251                         var _cl_MI xi = x[i];
 252                         for (sintL j = ylen-1; j > 0; j--)
 253                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,y[j]));
 254                         /* j=0 */ result[i] = R->_mul(xi,y[0]);
 255                 }
 256         } else {
 257                 {
 258                         var sintL j = ylen-1;
 259                         var _cl_MI yj = y[j];
 260                         for (sintL i = xlen-1; i >= 0; i--)
 261                                 result[i+j] = R->_mul(x[i],yj);
 262                 }
 263                 for (sintL j = ylen-2; j >= 0; j--) {
 264                         var _cl_MI yj = y[j];
 265                         for (sintL i = xlen-1; i > 0; i--)
 266                                 result[i+j] = R->_plus(result[i+j],R->_mul(x[i],yj));
 267                         /* i=0 */ result[j] = R->_mul(x[0],yj);
 268                 }
 269         }
 270         // Normalize (not necessary in integral domains).
 271         //modint_normalize(R,result,len);
 272         if (R->_zerop(result[len-1])) cl_abort();
 273         return _cl_UP(UPR, result);
 274 }}
 275
 276 static const _cl_UP modint_square (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 277 {{
 278         DeclarePoly(cl_GV_MI,x);
 279         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 280         var sintL xlen = x.length();
 281         if (xlen == 0)
 282                 return cl_UP(UPR, x);
 283         var sintL len = 2*xlen-1;
 284         var cl_GV_MI result = cl_GV_MI(len,R);
 285         if (xlen > 1) {
 286                 // Loop through all 0 <= j < i <= xlen-1.
 287                 {
 288                         var sintL i = xlen-1;
 289                         var _cl_MI xi = x[i];
 290                         for (sintL j = i-1; j >= 0; j--)
 291                                 result[i+j] = R->_mul(xi,x[j]);
 292                 }
 293                 {for (sintL i = xlen-2; i >= 1; i--) {
 294                         var _cl_MI xi = x[i];
 295                         for (sintL j = i-1; j >= 1; j--)
 296                                 result[i+j] = R->_plus(result[i+j],R->_mul(xi,x[j]));
 297                         /* j=0 */ result[i] = R->_mul(xi,x[0]);
 298                 }}
 299                 // Double.
 300                 {for (sintL i = len-2; i >= 1; i--)
 301                         result[i] = R->_plus(result[i],result[i]);
 302                 }
 303                 // Add squares.
 304                 result[2*(xlen-1)] = R->_square(x[xlen-1]);
 305                 for (sintL i = xlen-2; i >= 1; i--)
 306                         result[2*i] = R->_plus(result[2*i],R->_square(x[i]));
 307         }
 308         result[0] = R->_square(x[0]);
 309         // Normalize (not necessary in integral domains).
 310         //modint_normalize(R,result,len);
 311         if (R->_zerop(result[len-1])) cl_abort();
 312         return _cl_UP(UPR, result);
 313 }}
 314
 315 static const _cl_UP modint_exptpos (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_I& y)
 316 {
 317         var _cl_UP a = x;
 318         var cl_I b = y;
 319         while (!oddp(b)) { a = UPR->_square(a); b = b >> 1; }
 320         var _cl_UP c = a;
 321         until (b == 1)
 322           { b = b >> 1;
 323             a = UPR->_square(a);
 324             if (oddp(b)) { c = UPR->_mul(a,c); }
 325           }
 326         return c;
 327 }
 328
 329 static const _cl_UP modint_scalmul (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, const _cl_UP& y)
 330 {
 331         if (!(UPR->basering() == x.ring())) cl_abort();
 332  {
 333         DeclarePoly(_cl_MI,x);
 334         DeclarePoly(cl_GV_MI,y);
 335         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 336         var sintL ylen = y.length();
 337         if (ylen == 0)
 338                 return _cl_UP(UPR, y);
 339         if (R->_zerop(x))
 340                 return _cl_UP(UPR, cl_null_GV_I);
 341         // Now ylen > 0.
 342         // No normalization necessary, since the degree doesn't change.
 343         var cl_GV_MI result = cl_GV_MI(ylen,R);
 344         for (sintL i = ylen-1; i >= 0; i--)
 345                 result[i] = R->_mul(x,y[i]);
 346         return _cl_UP(UPR, result);
 347 }}
 348
 349 static sintL modint_degree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 350 {
 351         unused UPR;
 352  {      DeclarePoly(cl_GV_MI,x);
 353         return (sintL) x.length() - 1;
 354 }}
 355
 356 static const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 357 {
 358         if (!(UPR->basering() == x.ring())) cl_abort();
 359  {      DeclarePoly(_cl_MI,x);
 360         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 361         if (R->_zerop(x))
 362                 return _cl_UP(UPR, cl_null_GV_I);
 363         else {
 364                 var sintL len = e+1;
 365                 var cl_GV_MI result = cl_GV_MI(len,R);
 366                 result[e] = x;
 367                 return _cl_UP(UPR, result);
 368         }
 369 }}
 370
 371 static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 372 {{
 373         DeclarePoly(cl_GV_MI,x);
 374         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 375         if (index < x.length())
 376                 return cl_MI(R, x[index]);
 377         else
 378                 return R->zero();
 379 }}
 380
 381 static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
 382 {
 383         if (deg < 0)
 384                 return _cl_UP(UPR, cl_null_GV_I);
 385         else {
 386                 var sintL len = deg+1;
 387                 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 388                 return _cl_UP(UPR, cl_GV_MI(len,R));
 389         }
 390 }
 391
 392 static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 393 {{
 394         DeclareMutablePoly(cl_GV_MI,x);
 395         if (!(UPR->basering() == y.ring())) cl_abort();
 396   {     DeclarePoly(_cl_MI,y);
 397         if (!(index < x.length())) cl_abort();
 398         x[index] = y;
 399 }}}
 400
 401 static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 402 {{
 403         DeclareMutablePoly(cl_GV_MI,x); // NB: x is modified by reference!
 404         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 405         var uintL len = x.length();
 406         if (len > 0)
 407                 modint_normalize(R,x,len);
 408 }}
 409
 410 static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 411 {{
 412         // Method:
 413         // If x = 0, return 0.
 414         // If y = 0, return x[0].
 415         // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
 416         DeclarePoly(cl_GV_MI,x);
 417         if (!(UPR->basering() == y.ring())) cl_abort();
 418   {     DeclarePoly(_cl_MI,y);
 419         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 420         var uintL len = x.length();
 421         if (len==0)
 422                 return R->zero();
 423         if (R->_zerop(y))
 424                 return cl_MI(R, x[0]);
 425         var sintL i = len-1;
 426         var _cl_MI z = x[i];
 427         for ( ; --i >= 0; )
 428                 z = R->_plus(R->_mul(z,y),x[i]);
 429         return cl_MI(R, z);
 430 }}}
 431
 432 static cl_univpoly_setops modint_setops = {
 433         modint_fprint,
 434         modint_equal
 435 };
 436
 437 static cl_univpoly_addops modint_addops = {
 438         modint_zero,
 439         modint_zerop,
 440         modint_plus,
 441         modint_minus,
 442         modint_uminus
 443 };
 444
 445 static cl_univpoly_mulops modint_mulops = {
 446         modint_one,
 447         modint_canonhom,
 448         modint_mul,
 449         modint_square,
 450         modint_exptpos
 451 };
 452
 453 static cl_univpoly_modulops modint_modulops = {
 454         modint_scalmul
 455 };
 456
 457 static cl_univpoly_polyops modint_polyops = {
 458         modint_degree,
 459         modint_monomial,
 460         modint_coeff,
 461         modint_create,
 462         modint_set_coeff,
 463         modint_finalize,
 464         modint_eval
 465 };
 466
 467 class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
 468         SUBCLASS_cl_heap_univpoly_ring()
 469 public:
 470         // Constructor.
 471         cl_heap_modint_univpoly_ring (const cl_ring& r)
 472                 : cl_heap_univpoly_ring (r, &modint_setops, &modint_addops, &modint_mulops, &modint_modulops, &modint_polyops) {}
 473 };
 474
 475 }  // namespace cln