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, std::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_uminus (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 147 {{
 148         DeclarePoly(cl_GV_MI,x);
 149         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 150         var sintL xlen = x.length();
 151         if (xlen == 0)
 152                 return _cl_UP(UPR, x);
 153         // Now xlen > 0.
 154         // Negate. No normalization necessary, since the degree doesn't change.
 155         var sintL i = xlen-1;
 156         var _cl_MI hicoeff = R->_uminus(x[i]);
 157         if (R->_zerop(hicoeff)) cl_abort();
 158         var cl_GV_MI result = cl_GV_MI(xlen,R);
 159         result[i] = hicoeff;
 160         for (i-- ; i >= 0; i--)
 161                 result[i] = R->_uminus(x[i]);
 162         return _cl_UP(UPR, result);
 163 }}
 164
 165 static const _cl_UP modint_minus (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const _cl_UP& y)
 166 {{
 167         DeclarePoly(cl_GV_MI,x);
 168         DeclarePoly(cl_GV_MI,y);
 169         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 170         var sintL xlen = x.length();
 171         var sintL ylen = y.length();
 172         if (ylen == 0)
 173                 return _cl_UP(UPR, x);
 174         if (xlen == 0)
 175                 return modint_uminus(UPR, _cl_UP(UPR, y));
 176         // Now xlen > 0, ylen > 0.
 177         if (xlen > ylen) {
 178                 var cl_GV_MI result = cl_GV_MI(xlen,R);
 179                 var sintL i;
 180                 #if 0
 181                 for (i = xlen-1; i >= ylen; i--)
 182                         result[i] = x[i];
 183                 #else
 184                 cl_GV_MI::copy_elements(x,ylen,result,ylen,xlen-ylen);
 185                 #endif
 186                 for (i = ylen-1; i >= 0; i--)
 187                         result[i] = R->_minus(x[i],y[i]);
 188                 return _cl_UP(UPR, result);
 189         }
 190         if (xlen < ylen) {
 191                 var cl_GV_MI result = cl_GV_MI(ylen,R);
 192                 var sintL i;
 193                 for (i = ylen-1; i >= xlen; i--)
 194                         result[i] = R->_uminus(y[i]);
 195                 for (i = xlen-1; i >= 0; i--)
 196                         result[i] = R->_minus(x[i],y[i]);
 197                 return _cl_UP(UPR, result);
 198         }
 199         // Now xlen = ylen > 0. Add and normalize simultaneously.
 200         for (var sintL i = xlen-1; i >= 0; i--) {
 201                 var _cl_MI hicoeff = R->_minus(x[i],y[i]);
 202                 if (!R->_zerop(hicoeff)) {
 203                         var cl_GV_MI result = cl_GV_MI(i+1,R);
 204                         result[i] = hicoeff;
 205                         for (i-- ; i >= 0; i--)
 206                                 result[i] = R->_minus(x[i],y[i]);
 207                         return _cl_UP(UPR, result);
 208                 }
 209         }
 210         return _cl_UP(UPR, cl_null_GV_I);
 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 sintL modint_ldegree (cl_heap_univpoly_ring* UPR, const _cl_UP& x)
 357 {{
 358         DeclarePoly(cl_GV_MI,x);
 359         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 360         var sintL xlen = x.length();
 361         for (sintL i = 0; i < xlen; i++) {
 362                 if (!R->_zerop(x[i]))
 363                         return i;
 364         }
 365         return -1;
 366 }}
 367
 368 static const _cl_UP modint_monomial (cl_heap_univpoly_ring* UPR, const cl_ring_element& x, uintL e)
 369 {
 370         if (!(UPR->basering() == x.ring())) cl_abort();
 371  {      DeclarePoly(_cl_MI,x);
 372         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 373         if (R->_zerop(x))
 374                 return _cl_UP(UPR, cl_null_GV_I);
 375         else {
 376                 var sintL len = e+1;
 377                 var cl_GV_MI result = cl_GV_MI(len,R);
 378                 result[e] = x;
 379                 return _cl_UP(UPR, result);
 380         }
 381 }}
 382
 383 static const cl_ring_element modint_coeff (cl_heap_univpoly_ring* UPR, const _cl_UP& x, uintL index)
 384 {{
 385         DeclarePoly(cl_GV_MI,x);
 386         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 387         if (index < x.length())
 388                 return cl_MI(R, x[index]);
 389         else
 390                 return R->zero();
 391 }}
 392
 393 static const _cl_UP modint_create (cl_heap_univpoly_ring* UPR, sintL deg)
 394 {
 395         if (deg < 0)
 396                 return _cl_UP(UPR, cl_null_GV_I);
 397         else {
 398                 var sintL len = deg+1;
 399                 var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 400                 return _cl_UP(UPR, cl_GV_MI(len,R));
 401         }
 402 }
 403
 404 static void modint_set_coeff (cl_heap_univpoly_ring* UPR, _cl_UP& x, uintL index, const cl_ring_element& y)
 405 {{
 406         DeclareMutablePoly(cl_GV_MI,x);
 407         if (!(UPR->basering() == y.ring())) cl_abort();
 408   {     DeclarePoly(_cl_MI,y);
 409         if (!(index < x.length())) cl_abort();
 410         x[index] = y;
 411 }}}
 412
 413 static void modint_finalize (cl_heap_univpoly_ring* UPR, _cl_UP& x)
 414 {{
 415         DeclareMutablePoly(cl_GV_MI,x); // NB: x is modified by reference!
 416         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 417         var uintL len = x.length();
 418         if (len > 0)
 419                 modint_normalize(R,x,len);
 420 }}
 421
 422 static const cl_ring_element modint_eval (cl_heap_univpoly_ring* UPR, const _cl_UP& x, const cl_ring_element& y)
 423 {{
 424         // Method:
 425         // If x = 0, return 0.
 426         // If y = 0, return x[0].
 427         // Else compute (...(x[len-1]*y+x[len-2])*y ...)*y + x[0].
 428         DeclarePoly(cl_GV_MI,x);
 429         if (!(UPR->basering() == y.ring())) cl_abort();
 430   {     DeclarePoly(_cl_MI,y);
 431         var cl_heap_modint_ring* R = TheModintRing(UPR->basering());
 432         var uintL len = x.length();
 433         if (len==0)
 434                 return R->zero();
 435         if (R->_zerop(y))
 436                 return cl_MI(R, x[0]);
 437         var sintL i = len-1;
 438         var _cl_MI z = x[i];
 439         for ( ; --i >= 0; )
 440                 z = R->_plus(R->_mul(z,y),x[i]);
 441         return cl_MI(R, z);
 442 }}}
 443
 444 static cl_univpoly_setops modint_setops = {
 445         modint_fprint,
 446         modint_equal
 447 };
 448
 449 static cl_univpoly_addops modint_addops = {
 450         modint_zero,
 451         modint_zerop,
 452         modint_plus,
 453         modint_minus,
 454         modint_uminus
 455 };
 456
 457 static cl_univpoly_mulops modint_mulops = {
 458         modint_one,
 459         modint_canonhom,
 460         modint_mul,
 461         modint_square,
 462         modint_exptpos
 463 };
 464
 465 static cl_univpoly_modulops modint_modulops = {
 466         modint_scalmul
 467 };
 468
 469 static cl_univpoly_polyops modint_polyops = {
 470         modint_degree,
 471         modint_ldegree,
 472         modint_monomial,
 473         modint_coeff,
 474         modint_create,
 475         modint_set_coeff,
 476         modint_finalize,
 477         modint_eval
 478 };
 479
 480 class cl_heap_modint_univpoly_ring : public cl_heap_univpoly_ring {
 481         SUBCLASS_cl_heap_univpoly_ring()
 482 public:
 483         // Constructor.
 484         cl_heap_modint_univpoly_ring (const cl_ring& r);
 485         // Destructor.
 486         ~cl_heap_modint_univpoly_ring () {}
 487 };
 488
 489 static void cl_heap_modint_univpoly_ring_destructor (cl_heap* pointer)
 490 {
 491         (*(cl_heap_modint_univpoly_ring*)pointer).~cl_heap_modint_univpoly_ring();
 492 }
 493
 494 cl_class cl_class_modint_univpoly_ring = {
 495         cl_heap_modint_univpoly_ring_destructor,
 496         0
 497 };
 498
 499 // Constructor.
 500 inline cl_heap_modint_univpoly_ring::cl_heap_modint_univpoly_ring (const cl_ring& r)
 501         : cl_heap_univpoly_ring (r, &modint_setops, &modint_addops, &modint_mulops, &modint_modulops, &modint_polyops)
 502 {
 503         type = &cl_class_modint_univpoly_ring;
 504 }
 505
 506 }  // namespace cln