src/modinteger/cl_MI_std.h

   1 // m > 1, standard representation, no tricks
   2
   3 static void std_fprint (cl_heap_modint_ring* R, cl_ostream stream, const _cl_MI &x)
   4 {
   5         fprint(stream,R->_retract(x));
   6         fprint(stream," mod ");
   7         fprint(stream,R->modulus);
   8 }
   9
  10 static const cl_I std_reduce_modulo (cl_heap_modint_ring* R, const cl_I& x)
  11 {
  12         return mod(x,R->modulus);
  13 }
  14
  15 static const _cl_MI std_canonhom (cl_heap_modint_ring* R, const cl_I& x)
  16 {
  17         return _cl_MI(R, mod(x,R->modulus));
  18 }
  19
  20 static const cl_I std_retract (cl_heap_modint_ring* R, const _cl_MI& x)
  21 {
  22         unused R;
  23         return x.rep;
  24 }
  25
  26 static const _cl_MI std_random (cl_heap_modint_ring* R, cl_random_state& randomstate)
  27 {
  28         return _cl_MI(R, random_I(randomstate,R->modulus));
  29 }
  30
  31 static const _cl_MI std_zero (cl_heap_modint_ring* R)
  32 {
  33         return _cl_MI(R, 0);
  34 }
  35
  36 static cl_boolean std_zerop (cl_heap_modint_ring* R, const _cl_MI& x)
  37 {
  38         unused R;
  39         return zerop(x.rep);
  40 }
  41
  42 static const _cl_MI std_plus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
  43 {
  44         var cl_I zr = x.rep + y.rep;
  45         return _cl_MI(R, (zr >= R->modulus ? zr - R->modulus : zr));
  46 }
  47
  48 static const _cl_MI std_minus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
  49 {
  50         var cl_I zr = x.rep - y.rep;
  51         return _cl_MI(R, (minusp(zr) ? zr + R->modulus : zr));
  52 }
  53
  54 static const _cl_MI std_uminus (cl_heap_modint_ring* R, const _cl_MI& x)
  55 {
  56         return _cl_MI(R, (zerop(x.rep) ? (cl_I)0 : R->modulus - x.rep));
  57 }
  58
  59 static const _cl_MI std_one (cl_heap_modint_ring* R)
  60 {
  61         return _cl_MI(R, 1);
  62 }
  63
  64 static const _cl_MI std_mul (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
  65 {
  66         return _cl_MI(R, mod(x.rep * y.rep, R->modulus));
  67 }
  68
  69 static const _cl_MI std_square (cl_heap_modint_ring* R, const _cl_MI& x)
  70 {
  71         return _cl_MI(R, mod(square(x.rep), R->modulus));
  72 }
  73
  74 static const cl_MI_x std_recip (cl_heap_modint_ring* R, const _cl_MI& x)
  75 {
  76         var const cl_I& xr = x.rep;
  77         var cl_I u, v;
  78         var cl_I g = xgcd(xr,R->modulus,&u,&v);
  79         // g = gcd(x,m) = x*u+m*v
  80         if (eq(g,1))
  81                 return cl_MI(R, (minusp(u) ? u + R->modulus : u));
  82         if (zerop(xr))
  83                 cl_error_division_by_0();
  84         return cl_notify_composite(R,xr);
  85 }
  86
  87 static const cl_MI_x std_div (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
  88 {
  89         var const cl_I& yr = y.rep;
  90         var cl_I u, v;
  91         var cl_I g = xgcd(yr,R->modulus,&u,&v);
  92         // g = gcd(y,m) = y*u+m*v
  93         if (eq(g,1))
  94                 return cl_MI(R, mod(x.rep * (minusp(u) ? u + R->modulus : u), R->modulus));
  95         if (zerop(yr))
  96                 cl_error_division_by_0();
  97         return cl_notify_composite(R,yr);
  98 }
  99
 100 static uint8 const ord2_table[256] = // maps i -> ord2(i) for i>0
 101 {
 102  63,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 103   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 104   5,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 105   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 106   6,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 107   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 108   5,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 109   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 110   7,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 111   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 112   5,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 113   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 114   6,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 115   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 116   5,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0,
 117   4,  0,  1,  0,  2,  0,  1,  0,  3,  0,  1,  0,  2,  0,  1,  0
 118 };
 119
 120 static uint8 const odd_table[256] = // maps i -> i/2^ord2(i) for i>0
 121 {
 122   0,   1,   1,   3,   1,   5,   3,   7,   1,   9,   5,  11,   3,  13,   7,  15,
 123   1,  17,   9,  19,   5,  21,  11,  23,   3,  25,  13,  27,   7,  29,  15,  31,
 124   1,  33,  17,  35,   9,  37,  19,  39,   5,  41,  21,  43,  11,  45,  23,  47,
 125   3,  49,  25,  51,  13,  53,  27,  55,   7,  57,  29,  59,  15,  61,  31,  63,
 126   1,  65,  33,  67,  17,  69,  35,  71,   9,  73,  37,  75,  19,  77,  39,  79,
 127   5,  81,  41,  83,  21,  85,  43,  87,  11,  89,  45,  91,  23,  93,  47,  95,
 128   3,  97,  49,  99,  25, 101,  51, 103,  13, 105,  53, 107,  27, 109,  55, 111,
 129   7, 113,  57, 115,  29, 117,  59, 119,  15, 121,  61, 123,  31, 125,  63, 127,
 130   1, 129,  65, 131,  33, 133,  67, 135,  17, 137,  69, 139,  35, 141,  71, 143,
 131   9, 145,  73, 147,  37, 149,  75, 151,  19, 153,  77, 155,  39, 157,  79, 159,
 132   5, 161,  81, 163,  41, 165,  83, 167,  21, 169,  85, 171,  43, 173,  87, 175,
 133  11, 177,  89, 179,  45, 181,  91, 183,  23, 185,  93, 187,  47, 189,  95, 191,
 134   3, 193,  97, 195,  49, 197,  99, 199,  25, 201, 101, 203,  51, 205, 103, 207,
 135  13, 209, 105, 211,  53, 213, 107, 215,  27, 217, 109, 219,  55, 221, 111, 223,
 136   7, 225, 113, 227,  57, 229, 115, 231,  29, 233, 117, 235,  59, 237, 119, 239,
 137  15, 241, 121, 243,  61, 245, 123, 247,  31, 249, 125, 251,  63, 253, 127, 255
 138 };
 139
 140 static const _cl_MI std_expt_pos (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
 141 {
 142 #if 0
 143         // Methode:
 144         // Right-Left Binary, [Cohen, Algorithm 1.2.1.]
 145         //   a:=x, b:=y.
 146         //   Solange b gerade, setze a:=a*a, b:=b/2. (Invariante: a^b = x^y.)
 147         //   c:=a.
 148         //   Solange b:=floor(b/2) >0 ist,
 149         //     setze a:=a*a, und falls b ungerade, setze c:=a*c.
 150         //   Liefere c.
 151         var _cl_MI a = x;
 152         var cl_I b = y;
 153         while (!oddp(b)) { a = R->_square(a); b = b >> 1; } // a^b = x^y
 154         var _cl_MI c = a;
 155         until (eq(b,1)) {
 156                 b = b >> 1;
 157                 a = R->_square(a);
 158                 // a^b*c = x^y
 159                 if (oddp(b))
 160                         c = R->_mul(a,c);
 161         }
 162         return c;
 163 #else
 164         // Methode:
 165         // Left-Right base 2^k, [Cohen, Algorithm 1.2.4.]
 166         // Good values of k, depending on the size nn of the exponent n:
 167         //   k = 1 for nn <= 8
 168         //   k = 2 for nn <= 24
 169         //   k = 3 for nn <= 69.8
 170         //   ... k for nn <= k*(k+1)*2^(2k)/(2^(k+1)-k-2)
 171         var cl_I n = y;
 172         var uintL nn = integer_length(n);
 173         // n has nn bits.
 174         if (nn <= 8) {
 175                 // k = 1, normal Left-Right Binary algorithm.
 176                 var uintL _n = FN_to_UL(n);
 177                 var _cl_MI a = x;
 178                 for (var int i = nn-2; i >= 0; i--) {
 179                         a = R->_square(a);
 180                         if (_n & bit(i))
 181                                 a = R->_mul(a,x);
 182                 }
 183                 return cl_MI(R,a);
 184         } else {
 185                 // General Left-Right base 2^k algorithm.
 186                 CL_ALLOCA_STACK;
 187                 var uintL k;
 188                      if (nn <= 24)      k = 2;
 189                 else if (nn <= 69)      k = 3;
 190                 else if (nn <= 196)     k = 4;
 191                 else if (nn <= 538)     k = 5;
 192                 else if (nn <= 1433)    k = 6;
 193                 else if (nn <= 3714)    k = 7;
 194                 else if (nn <= 9399)    k = 8;
 195                 else if (nn <= 23290)   k = 9;
 196                 else if (nn <= 56651)   k = 10;
 197                 else if (nn <= 135598)  k = 11;
 198                 else if (nn <= 320034)  k = 12;
 199                 else if (nn <= 746155)  k = 13;
 200                 else if (nn <= 1721160) k = 14;
 201                 else if (nn <= 3933180) k = 15;
 202                 else /* if (nn <= 8914120) */ k = 16;
 203                 var uintL nnk = ceiling(nn,k); // number of base-2^k digits in n
 204                 var uint16* n_digits = cl_alloc_array(uint16,nnk);
 205                 // Split n into base-2^k digits.
 206                 {
 207                         var const uintD* n_LSDptr;
 208                         var const uintD* n_MSDptr;
 209                         I_to_NDS_nocopy(n, n_MSDptr=,,n_LSDptr=,cl_false,);
 210                         var const uintL k_mask = bit(k)-1;
 211                         var uintD carry = 0;
 212                         var unsigned int carrybits = 0;
 213                         for (var uintL i = 0; i < nnk; i++) {
 214                                 if (carrybits >= k) {
 215                                         n_digits[i] = carry & k_mask;
 216                                         carry = carry >> k;
 217                                         carrybits -= k;
 218                                 } else {
 219                                         var uintD next =
 220                                           (n_LSDptr==n_MSDptr ? 0 : lsprefnext(n_LSDptr));
 221                                         n_digits[i] = (carry | (next << carrybits)) & k_mask;
 222                                         carry = next >> (k-carrybits);
 223                                         carrybits = intDsize - (k-carrybits);
 224                                 }
 225                         }
 226                 }
 227                 // Compute maximum odd base-2^k digit.
 228                 var uintL maxodd = 1;
 229                 if (k <= 8) {
 230                         for (var uintL i = 0; i < nnk; i++) {
 231                                 var uintL d = n_digits[i];
 232                                 if (d > 0) {
 233                                         d = odd_table[d];
 234                                         if (d > maxodd) maxodd = d;
 235                                 }
 236                         }
 237                 } else {
 238                         for (var uintL i = 0; i < nnk; i++) {
 239                                 var uintL d = n_digits[i];
 240                                 if (d > 0) {
 241                                         var uintL d2; ord2_32(d,d2=);
 242                                         d = d>>d2;
 243                                         if (d > maxodd) maxodd = d;
 244                                 }
 245                         }
 246                 }
 247                 maxodd = (maxodd+1)/2; // number of odd powers we need
 248                 var _cl_MI* x_oddpow = cl_alloc_array(_cl_MI,maxodd);
 249                 var _cl_MI x2 = (maxodd > 1 ? R->_square(x) : R->_zero());
 250                 {
 251                         init1(_cl_MI, x_oddpow[0]) (x);
 252                         for (var uintL i = 1; i < maxodd; i++)
 253                                 init1(_cl_MI, x_oddpow[i]) (R->_mul(x_oddpow[i-1],x2));
 254                 }
 255                 var _cl_MI a;
 256                 // Compute a = x^n_digits[nnk-1].
 257                 {
 258                         var uintL d = n_digits[nnk-1];
 259                         if (d == 0) cl_abort();
 260                         var uintL d2;
 261                         if (k <= 8)
 262                                 d2 = ord2_table[d];
 263                         else
 264                                 ord2_32(d,d2=);
 265                         d = d>>d2; // d := d/2^ord2(d)
 266                         d = floor(d,2);
 267                         a = x_oddpow[d]; // x^d
 268                         // Square d2 times.
 269                         if (d==0 && maxodd > 1 && d2>0) {
 270                                 a = x2; d2--;
 271                         }
 272                         if (!(d2 < k)) cl_abort();
 273                         for ( ; d2>0; d2--)
 274                                 a = R->_square(a);
 275                 }
 276                 for (var sintL i = nnk-2; i >= 0; i--) {
 277                         // Compute a := a^(2^k) * x^n_digits[i].
 278                         var uintL d = n_digits[i];
 279                         var uintL d2;
 280                         if (d > 0) {
 281                                 if (k <= 8)
 282                                         d2 = ord2_table[d];
 283                                 else
 284                                         ord2_32(d,d2=);
 285                                 d = d>>d2; // d/2^ord2(d)
 286                                 d = floor(d,2);
 287                                 for (var sintL j = k-d2; j>0; j--)
 288                                         a = R->_square(a);
 289                                 a = R->_mul(a,x_oddpow[d]);
 290                         } else
 291                                 d2 = k;
 292                         // Square d2 times.
 293                         if (!(d2 <= k)) cl_abort();
 294                         for ( ; d2>0; d2--)
 295                                 a = R->_square(a);
 296                 }
 297                 {
 298                         for (var uintL i = 0; i < maxodd; i++)
 299                                 x_oddpow[i].~_cl_MI();
 300                 }
 301                 return cl_MI(R,a);
 302         }
 303 #endif
 304 }
 305
 306 static const cl_MI_x std_expt (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
 307 {
 308         if (!minusp(y)) {
 309                 if (zerop(y))
 310                         return R->one();
 311                 else
 312                         return cl_MI(R,R->_expt_pos(x,y));
 313         } else
 314                 return R->_recip(R->_expt_pos(x,-y));
 315 }
 316
 317 static cl_modint_setops std_setops = {
 318         std_fprint,
 319         modint_equal,
 320         std_random
 321 };
 322 static cl_modint_addops std_addops = {
 323         std_zero,
 324         std_zerop,
 325         std_plus,
 326         std_minus,
 327         std_uminus
 328 };
 329 static cl_modint_mulops std_mulops = {
 330         std_one,
 331         std_canonhom,
 332         std_mul,
 333         std_square,
 334         std_expt_pos,
 335         std_recip,
 336         std_div,
 337         std_expt,
 338         std_reduce_modulo,
 339         std_retract
 340 };
 341
 342 class cl_heap_modint_ring_std : public cl_heap_modint_ring {
 343         SUBCLASS_cl_heap_modint_ring()
 344 public:
 345         // Constructor.
 346         cl_heap_modint_ring_std (const cl_I& m)
 347                 : cl_heap_modint_ring (m, &std_setops, &std_addops, &std_mulops) {}
 348         // Virtual destructor.
 349         ~cl_heap_modint_ring_std () {}
 350 };