1 // m > 1, standard representation, no tricks
5 static void std_fprint (cl_heap_modint_ring* R, std::ostream& stream, const _cl_MI &x)
7 fprint(stream,R->_retract(x));
8 fprint(stream," mod ");
9 fprint(stream,R->modulus);
12 static const cl_I std_reduce_modulo (cl_heap_modint_ring* R, const cl_I& x)
14 return mod(x,R->modulus);
17 static const _cl_MI std_canonhom (cl_heap_modint_ring* R, const cl_I& x)
19 return _cl_MI(R, mod(x,R->modulus));
22 static const cl_I std_retract (cl_heap_modint_ring* R, const _cl_MI& x)
28 static const _cl_MI std_random (cl_heap_modint_ring* R, random_state& randomstate)
30 return _cl_MI(R, random_I(randomstate,R->modulus));
33 static const _cl_MI std_zero (cl_heap_modint_ring* R)
38 static bool std_zerop (cl_heap_modint_ring* R, const _cl_MI& x)
44 static const _cl_MI std_plus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
46 var cl_I zr = x.rep + y.rep;
47 return _cl_MI(R, (zr >= R->modulus ? zr - R->modulus : zr));
50 static const _cl_MI std_minus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
52 var cl_I zr = x.rep - y.rep;
53 return _cl_MI(R, (minusp(zr) ? zr + R->modulus : zr));
56 static const _cl_MI std_uminus (cl_heap_modint_ring* R, const _cl_MI& x)
58 return _cl_MI(R, (zerop(x.rep) ? (cl_I)0 : R->modulus - x.rep));
61 static const _cl_MI std_one (cl_heap_modint_ring* R)
66 static const _cl_MI std_mul (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
68 return _cl_MI(R, mod(x.rep * y.rep, R->modulus));
71 static const _cl_MI std_square (cl_heap_modint_ring* R, const _cl_MI& x)
73 return _cl_MI(R, mod(square(x.rep), R->modulus));
76 static const cl_MI_x std_recip (cl_heap_modint_ring* R, const _cl_MI& x)
78 var const cl_I& xr = x.rep;
80 var cl_I g = xgcd(xr,R->modulus,&u,&v);
81 // g = gcd(x,m) = x*u+m*v
83 return cl_MI(R, (minusp(u) ? u + R->modulus : u));
85 throw division_by_0_exception();
86 return cl_notify_composite(R,xr);
89 static const cl_MI_x std_div (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
91 var const cl_I& yr = y.rep;
93 var cl_I g = xgcd(yr,R->modulus,&u,&v);
94 // g = gcd(y,m) = y*u+m*v
96 return cl_MI(R, mod(x.rep * (minusp(u) ? u + R->modulus : u), R->modulus));
98 throw division_by_0_exception();
99 return cl_notify_composite(R,yr);
102 static uint8 const ord2_table[256] = // maps i -> ord2(i) for i>0
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
122 static uint8 const odd_table[256] = // maps i -> i/2^ord2(i) for i>0
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
142 static const _cl_MI std_expt_pos (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
146 // Right-Left Binary, [Cohen, Algorithm 1.2.1.]
148 // Solange b gerade, setze a:=a*a, b:=b/2. (Invariante: a^b = x^y.)
150 // Solange b:=floor(b/2) >0 ist,
151 // setze a:=a*a, und falls b ungerade, setze c:=a*c.
155 while (!oddp(b)) { a = R->_square(a); b = b >> 1; } // a^b = x^y
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:
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)
174 var uintC nn = integer_length(n);
177 // k = 1, normal Left-Right Binary algorithm.
178 var uintL _n = FN_to_UV(n);
180 for (var int i = nn-2; i >= 0; i--) {
187 // General Left-Right base 2^k algorithm.
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 uintC 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.
209 var const uintD* n_LSDptr;
210 var const uintD* n_MSDptr;
211 I_to_NDS_nocopy(n, n_MSDptr=,,n_LSDptr=,false,);
212 var const uintL k_mask = bit(k)-1;
214 var unsigned int carrybits = 0;
215 for (var uintC i = 0; i < nnk; i++) {
216 if (carrybits >= k) {
217 n_digits[i] = carry & k_mask;
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);
229 // Compute maximum odd base-2^k digit.
230 var uintL maxodd = 1;
232 for (var uintC i = 0; i < nnk; i++) {
233 var uintL d = n_digits[i];
236 if (d > maxodd) maxodd = d;
240 for (var uintC i = 0; i < nnk; i++) {
241 var uintL d = n_digits[i];
243 var uintL d2; ord2_32(d,d2=);
245 if (d > maxodd) maxodd = d;
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());
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));
258 // Compute a = x^n_digits[nnk-1].
260 var uintL d = n_digits[nnk-1];
261 if (d == 0) throw runtime_exception();
267 d = d>>d2; // d := d/2^ord2(d)
269 a = x_oddpow[d]; // x^d
271 if (d==0 && maxodd > 1 && d2>0) {
274 if (!(d2 < k)) throw runtime_exception();
278 for (var sintC i = nnk-2; i >= 0; i--) {
279 // Compute a := a^(2^k) * x^n_digits[i].
280 var uintL d = n_digits[i];
287 d = d>>d2; // d/2^ord2(d)
289 for (var sintL j = k-d2; j>0; j--)
291 a = R->_mul(a,x_oddpow[d]);
295 if (!(d2 <= k)) throw runtime_exception();
300 for (var uintL i = 0; i < maxodd; i++)
301 x_oddpow[i].~_cl_MI();
308 static const cl_MI_x std_expt (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
314 return cl_MI(R,R->_expt_pos(x,y));
316 return R->_recip(R->_expt_pos(x,-y));
319 static cl_modint_setops std_setops = {
324 static cl_modint_addops std_addops = {
331 static cl_modint_mulops std_mulops = {
344 class cl_heap_modint_ring_std : public cl_heap_modint_ring {
345 SUBCLASS_cl_heap_modint_ring()
348 cl_heap_modint_ring_std (const cl_I& m);
349 // Virtual destructor.
350 ~cl_heap_modint_ring_std () {}
353 static void cl_heap_modint_ring_std_destructor (cl_heap* pointer)
355 (*(cl_heap_modint_ring_std*)pointer).~cl_heap_modint_ring_std();
358 cl_class cl_class_modint_ring_std = {
359 cl_heap_modint_ring_std_destructor,
360 cl_class_flags_modint_ring
364 inline cl_heap_modint_ring_std::cl_heap_modint_ring_std (const cl_I& m)
365 : cl_heap_modint_ring (m, &std_setops, &std_addops, &std_mulops)
367 type = &cl_class_modint_ring_std;