src/modinteger/cl_MI_std.h

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