157 lines
3.1 KiB
D
157 lines
3.1 KiB
D
import std.bigint;
|
|
import std.random;
|
|
import std.stdio;
|
|
|
|
struct PExp {
|
|
BigInt prime;
|
|
int exp;
|
|
}
|
|
|
|
BigInt gcd(BigInt x, BigInt y) {
|
|
if (y == 0) {
|
|
return x;
|
|
}
|
|
return gcd(y, x % y);
|
|
}
|
|
|
|
/// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
|
|
BigInt modPow(BigInt b, BigInt e, BigInt n) {
|
|
if (n == 1) return BigInt(0);
|
|
BigInt result = 1;
|
|
b = b % n;
|
|
while (e > 0) {
|
|
if (e % 2 == 1) {
|
|
result = (result * b) % n;
|
|
}
|
|
e >>= 1;
|
|
b = (b*b) % n;
|
|
}
|
|
return result;
|
|
}
|
|
|
|
BigInt pow(long b, long e) {
|
|
return pow(BigInt(b), BigInt(e));
|
|
}
|
|
BigInt pow(BigInt b, BigInt e) {
|
|
if (e == 0) {
|
|
return BigInt(1);
|
|
}
|
|
|
|
BigInt result = 1;
|
|
while (e > 1) {
|
|
if (e % 2 == 0) {
|
|
b *= b;
|
|
e /= 2;
|
|
} else {
|
|
result *= b;
|
|
b *= b;
|
|
e = (e - 1) / 2;
|
|
}
|
|
}
|
|
|
|
return b * result;
|
|
}
|
|
|
|
BigInt sqrt(BigInt self) {
|
|
BigInt b = self;
|
|
while (true) {
|
|
BigInt a = b;
|
|
b = self / a + a >> 1;
|
|
if (b >= a) return a;
|
|
}
|
|
}
|
|
|
|
long bitLength(BigInt self) {
|
|
BigInt bi = self;
|
|
long length;
|
|
while (bi != 0) {
|
|
length++;
|
|
bi >>= 1;
|
|
}
|
|
return length;
|
|
}
|
|
|
|
PExp[] factor(BigInt n) {
|
|
PExp[] pf;
|
|
BigInt nn = n;
|
|
int b = 0;
|
|
int e = 1;
|
|
while ((nn & e) == 0) {
|
|
e <<= 1;
|
|
b++;
|
|
}
|
|
if (b > 0) {
|
|
nn = nn >> b;
|
|
pf ~= PExp(BigInt(2), b);
|
|
}
|
|
BigInt s = nn.sqrt();
|
|
BigInt d = 3;
|
|
while (nn > 1) {
|
|
if (d > s) d = nn;
|
|
e = 0;
|
|
while (true) {
|
|
BigInt div, rem;
|
|
nn.divMod(d, div, rem);
|
|
if (rem.bitLength > 0) break;
|
|
nn = div;
|
|
e++;
|
|
}
|
|
if (e > 0) {
|
|
pf ~= PExp(d, e);
|
|
s = nn.sqrt();
|
|
}
|
|
d += 2;
|
|
}
|
|
|
|
return pf;
|
|
}
|
|
|
|
BigInt moBachShallit58(BigInt a, BigInt n, PExp[] pf) {
|
|
BigInt n1 = n - 1;
|
|
BigInt mo = 1;
|
|
foreach(pe; pf) {
|
|
BigInt y = n1 / pe.prime.pow(BigInt(pe.exp));
|
|
int o = 0;
|
|
BigInt x = a.modPow(y, n);
|
|
while (x > 1) {
|
|
x = x.modPow(pe.prime, n);
|
|
o++;
|
|
}
|
|
BigInt o1 = pe.prime.pow(BigInt(o));
|
|
o1 = o1 / gcd(mo, o1);
|
|
mo = mo * o1;
|
|
}
|
|
return mo;
|
|
}
|
|
|
|
void moTest(ulong a, ulong n) {
|
|
moTest(BigInt(a), n);
|
|
}
|
|
void moTest(BigInt a, ulong n) {
|
|
// Commented out because the implementations tried all failed for the -2 and -3 tests.
|
|
// if (!n.isProbablePrime()) {
|
|
// writeln("Not computed. Modulus must be prime for this algorithm.");
|
|
// return;
|
|
// }
|
|
if (a.bitLength < 100) {
|
|
write("ord(", a, ")");
|
|
} else {
|
|
write("ord([big])");
|
|
}
|
|
write(" mod ", n, " ");
|
|
BigInt nn = n;
|
|
BigInt mob = moBachShallit58(a, nn, factor(nn - 1));
|
|
writeln("= ", mob);
|
|
}
|
|
|
|
void main() {
|
|
moTest(37, 3343);
|
|
|
|
moTest(pow(10, 100) + 1, 7919);
|
|
moTest(pow(10, 1000) + 1, 15485863);
|
|
moTest(pow(10, 10000) - 1, 22801763489);
|
|
|
|
moTest(1511678068, 7379191741);
|
|
moTest(3047753288, 2257683301);
|
|
}
|