// m > 1, standard representation, no tricks
namespace cln {
static void std_fprint (cl_heap_modint_ring* R, std::ostream& stream, const _cl_MI &x)
{
fprint(stream,R->_retract(x));
fprint(stream," mod ");
fprint(stream,R->modulus);
}
static const cl_I std_reduce_modulo (cl_heap_modint_ring* R, const cl_I& x)
{
return mod(x,R->modulus);
}
static const _cl_MI std_canonhom (cl_heap_modint_ring* R, const cl_I& x)
{
return _cl_MI(R, mod(x,R->modulus));
}
static const cl_I std_retract (cl_heap_modint_ring* R, const _cl_MI& x)
{
unused R;
return x.rep;
}
static const _cl_MI std_random (cl_heap_modint_ring* R, random_state& randomstate)
{
return _cl_MI(R, random_I(randomstate,R->modulus));
}
static const _cl_MI std_zero (cl_heap_modint_ring* R)
{
return _cl_MI(R, 0);
}
static cl_boolean std_zerop (cl_heap_modint_ring* R, const _cl_MI& x)
{
unused R;
return zerop(x.rep);
}
static const _cl_MI std_plus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
{
var cl_I zr = x.rep + y.rep;
return _cl_MI(R, (zr >= R->modulus ? zr - R->modulus : zr));
}
static const _cl_MI std_minus (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
{
var cl_I zr = x.rep - y.rep;
return _cl_MI(R, (minusp(zr) ? zr + R->modulus : zr));
}
static const _cl_MI std_uminus (cl_heap_modint_ring* R, const _cl_MI& x)
{
return _cl_MI(R, (zerop(x.rep) ? (cl_I)0 : R->modulus - x.rep));
}
static const _cl_MI std_one (cl_heap_modint_ring* R)
{
return _cl_MI(R, 1);
}
static const _cl_MI std_mul (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
{
return _cl_MI(R, mod(x.rep * y.rep, R->modulus));
}
static const _cl_MI std_square (cl_heap_modint_ring* R, const _cl_MI& x)
{
return _cl_MI(R, mod(square(x.rep), R->modulus));
}
static const cl_MI_x std_recip (cl_heap_modint_ring* R, const _cl_MI& x)
{
var const cl_I& xr = x.rep;
var cl_I u, v;
var cl_I g = xgcd(xr,R->modulus,&u,&v);
// g = gcd(x,m) = x*u+m*v
if (eq(g,1))
return cl_MI(R, (minusp(u) ? u + R->modulus : u));
if (zerop(xr))
cl_error_division_by_0();
return cl_notify_composite(R,xr);
}
static const cl_MI_x std_div (cl_heap_modint_ring* R, const _cl_MI& x, const _cl_MI& y)
{
var const cl_I& yr = y.rep;
var cl_I u, v;
var cl_I g = xgcd(yr,R->modulus,&u,&v);
// g = gcd(y,m) = y*u+m*v
if (eq(g,1))
return cl_MI(R, mod(x.rep * (minusp(u) ? u + R->modulus : u), R->modulus));
if (zerop(yr))
cl_error_division_by_0();
return cl_notify_composite(R,yr);
}
static uint8 const ord2_table[256] = // maps i -> ord2(i) for i>0
{
63, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
7, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0,
4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
};
static uint8 const odd_table[256] = // maps i -> i/2^ord2(i) for i>0
{
0, 1, 1, 3, 1, 5, 3, 7, 1, 9, 5, 11, 3, 13, 7, 15,
1, 17, 9, 19, 5, 21, 11, 23, 3, 25, 13, 27, 7, 29, 15, 31,
1, 33, 17, 35, 9, 37, 19, 39, 5, 41, 21, 43, 11, 45, 23, 47,
3, 49, 25, 51, 13, 53, 27, 55, 7, 57, 29, 59, 15, 61, 31, 63,
1, 65, 33, 67, 17, 69, 35, 71, 9, 73, 37, 75, 19, 77, 39, 79,
5, 81, 41, 83, 21, 85, 43, 87, 11, 89, 45, 91, 23, 93, 47, 95,
3, 97, 49, 99, 25, 101, 51, 103, 13, 105, 53, 107, 27, 109, 55, 111,
7, 113, 57, 115, 29, 117, 59, 119, 15, 121, 61, 123, 31, 125, 63, 127,
1, 129, 65, 131, 33, 133, 67, 135, 17, 137, 69, 139, 35, 141, 71, 143,
9, 145, 73, 147, 37, 149, 75, 151, 19, 153, 77, 155, 39, 157, 79, 159,
5, 161, 81, 163, 41, 165, 83, 167, 21, 169, 85, 171, 43, 173, 87, 175,
11, 177, 89, 179, 45, 181, 91, 183, 23, 185, 93, 187, 47, 189, 95, 191,
3, 193, 97, 195, 49, 197, 99, 199, 25, 201, 101, 203, 51, 205, 103, 207,
13, 209, 105, 211, 53, 213, 107, 215, 27, 217, 109, 219, 55, 221, 111, 223,
7, 225, 113, 227, 57, 229, 115, 231, 29, 233, 117, 235, 59, 237, 119, 239,
15, 241, 121, 243, 61, 245, 123, 247, 31, 249, 125, 251, 63, 253, 127, 255
};
static const _cl_MI std_expt_pos (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
{
#if 0
// Methode:
// Right-Left Binary, [Cohen, Algorithm 1.2.1.]
// a:=x, b:=y.
// Solange b gerade, setze a:=a*a, b:=b/2. (Invariante: a^b = x^y.)
// c:=a.
// Solange b:=floor(b/2) >0 ist,
// setze a:=a*a, und falls b ungerade, setze c:=a*c.
// Liefere c.
var _cl_MI a = x;
var cl_I b = y;
while (!oddp(b)) { a = R->_square(a); b = b >> 1; } // a^b = x^y
var _cl_MI c = a;
until (eq(b,1)) {
b = b >> 1;
a = R->_square(a);
// a^b*c = x^y
if (oddp(b))
c = R->_mul(a,c);
}
return c;
#else
// Methode:
// Left-Right base 2^k, [Cohen, Algorithm 1.2.4.]
// Good values of k, depending on the size nn of the exponent n:
// k = 1 for nn <= 8
// k = 2 for nn <= 24
// k = 3 for nn <= 69.8
// ... k for nn <= k*(k+1)*2^(2k)/(2^(k+1)-k-2)
var cl_I n = y;
var uintL nn = integer_length(n);
// n has nn bits.
if (nn <= 8) {
// k = 1, normal Left-Right Binary algorithm.
var uintL _n = FN_to_UL(n);
var _cl_MI a = x;
for (var int i = nn-2; i >= 0; i--) {
a = R->_square(a);
if (_n & bit(i))
a = R->_mul(a,x);
}
return cl_MI(R,a);
} else {
// General Left-Right base 2^k algorithm.
CL_ALLOCA_STACK;
var uintL k;
if (nn <= 24) k = 2;
else if (nn <= 69) k = 3;
else if (nn <= 196) k = 4;
else if (nn <= 538) k = 5;
else if (nn <= 1433) k = 6;
else if (nn <= 3714) k = 7;
else if (nn <= 9399) k = 8;
else if (nn <= 23290) k = 9;
else if (nn <= 56651) k = 10;
else if (nn <= 135598) k = 11;
else if (nn <= 320034) k = 12;
else if (nn <= 746155) k = 13;
else if (nn <= 1721160) k = 14;
else if (nn <= 3933180) k = 15;
else /* if (nn <= 8914120) */ k = 16;
var uintL nnk = ceiling(nn,k); // number of base-2^k digits in n
var uint16* n_digits = cl_alloc_array(uint16,nnk);
// Split n into base-2^k digits.
{
var const uintD* n_LSDptr;
var const uintD* n_MSDptr;
I_to_NDS_nocopy(n, n_MSDptr=,,n_LSDptr=,cl_false,);
var const uintL k_mask = bit(k)-1;
var uintD carry = 0;
var unsigned int carrybits = 0;
for (var uintL i = 0; i < nnk; i++) {
if (carrybits >= k) {
n_digits[i] = carry & k_mask;
carry = carry >> k;
carrybits -= k;
} else {
var uintD next =
(n_LSDptr==n_MSDptr ? 0 : lsprefnext(n_LSDptr));
n_digits[i] = (carry | (next << carrybits)) & k_mask;
carry = next >> (k-carrybits);
carrybits = intDsize - (k-carrybits);
}
}
}
// Compute maximum odd base-2^k digit.
var uintL maxodd = 1;
if (k <= 8) {
for (var uintL i = 0; i < nnk; i++) {
var uintL d = n_digits[i];
if (d > 0) {
d = odd_table[d];
if (d > maxodd) maxodd = d;
}
}
} else {
for (var uintL i = 0; i < nnk; i++) {
var uintL d = n_digits[i];
if (d > 0) {
var uintL d2; ord2_32(d,d2=);
d = d>>d2;
if (d > maxodd) maxodd = d;
}
}
}
maxodd = (maxodd+1)/2; // number of odd powers we need
var _cl_MI* x_oddpow = cl_alloc_array(_cl_MI,maxodd);
var _cl_MI x2 = (maxodd > 1 ? R->_square(x) : R->_zero());
{
init1(_cl_MI, x_oddpow[0]) (x);
for (var uintL i = 1; i < maxodd; i++)
init1(_cl_MI, x_oddpow[i]) (R->_mul(x_oddpow[i-1],x2));
}
var _cl_MI a;
// Compute a = x^n_digits[nnk-1].
{
var uintL d = n_digits[nnk-1];
if (d == 0) cl_abort();
var uintL d2;
if (k <= 8)
d2 = ord2_table[d];
else
ord2_32(d,d2=);
d = d>>d2; // d := d/2^ord2(d)
d = floor(d,2);
a = x_oddpow[d]; // x^d
// Square d2 times.
if (d==0 && maxodd > 1 && d2>0) {
a = x2; d2--;
}
if (!(d2 < k)) cl_abort();
for ( ; d2>0; d2--)
a = R->_square(a);
}
for (var sintL i = nnk-2; i >= 0; i--) {
// Compute a := a^(2^k) * x^n_digits[i].
var uintL d = n_digits[i];
var uintL d2;
if (d > 0) {
if (k <= 8)
d2 = ord2_table[d];
else
ord2_32(d,d2=);
d = d>>d2; // d/2^ord2(d)
d = floor(d,2);
for (var sintL j = k-d2; j>0; j--)
a = R->_square(a);
a = R->_mul(a,x_oddpow[d]);
} else
d2 = k;
// Square d2 times.
if (!(d2 <= k)) cl_abort();
for ( ; d2>0; d2--)
a = R->_square(a);
}
{
for (var uintL i = 0; i < maxodd; i++)
x_oddpow[i].~_cl_MI();
}
return cl_MI(R,a);
}
#endif
}
static const cl_MI_x std_expt (cl_heap_modint_ring* R, const _cl_MI& x, const cl_I& y)
{
if (!minusp(y)) {
if (zerop(y))
return R->one();
else
return cl_MI(R,R->_expt_pos(x,y));
} else
return R->_recip(R->_expt_pos(x,-y));
}
static cl_modint_setops std_setops = {
std_fprint,
modint_equal,
std_random
};
static cl_modint_addops std_addops = {
std_zero,
std_zerop,
std_plus,
std_minus,
std_uminus
};
static cl_modint_mulops std_mulops = {
std_one,
std_canonhom,
std_mul,
std_square,
std_expt_pos,
std_recip,
std_div,
std_expt,
std_reduce_modulo,
std_retract
};
class cl_heap_modint_ring_std : public cl_heap_modint_ring {
SUBCLASS_cl_heap_modint_ring()
public:
// Constructor.
cl_heap_modint_ring_std (const cl_I& m);
// Virtual destructor.
~cl_heap_modint_ring_std () {}
};
static void cl_heap_modint_ring_std_destructor (cl_heap* pointer)
{
(*(cl_heap_modint_ring_std*)pointer).~cl_heap_modint_ring_std();
}
cl_class cl_class_modint_ring_std = {
cl_heap_modint_ring_std_destructor,
cl_class_flags_modint_ring
};
// Constructor.
inline cl_heap_modint_ring_std::cl_heap_modint_ring_std (const cl_I& m)
: cl_heap_modint_ring (m, &std_setops, &std_addops, &std_mulops)
{
type = &cl_class_modint_ring_std;
}
} // namespace cln
syntax highlighted by Code2HTML, v. 0.9.1