95 lines
3.1 KiB
C
95 lines
3.1 KiB
C
// return (a * b) % mod, avoiding overflow errors while doing modular multiplication.
|
|
static unsigned multiplication_modulo(unsigned a, unsigned b, const unsigned mod) {
|
|
unsigned res = 0, tmp;
|
|
for (b %= mod; a; a & 1 ? b >= mod - res ? res -= mod : 0, res += b : 0, a >>= 1, (tmp = b) >= mod - b ? tmp -= mod : 0, b += tmp);
|
|
return res % mod;
|
|
}
|
|
|
|
// return (n ^ exp) % mod
|
|
static unsigned mod_pow(unsigned n, unsigned exp, const unsigned mod) {
|
|
unsigned res = 1;
|
|
for (n %= mod; exp; exp & 1 ? res = multiplication_modulo(res, n, mod) : 0, n = multiplication_modulo(n, n, mod), exp >>= 1);
|
|
return res;
|
|
}
|
|
|
|
static unsigned tonelli_shanks_1(const unsigned n, const unsigned mod) {
|
|
// return root such that (root * root) % mod congruent to n % mod.
|
|
// return 0 if no solution to the congruence exists.
|
|
// mod is assumed odd prime.
|
|
const unsigned a = n % mod;
|
|
unsigned res, b, c, d, e, f, g, h;
|
|
if (mod_pow(a, (mod - 1) >> 1, mod) != 1)
|
|
res = 0;
|
|
else
|
|
switch (mod & 7) {
|
|
case 3 : case 7 :
|
|
res = mod_pow(a, (mod + 1) >> 2, mod);
|
|
break;
|
|
case 5 :
|
|
res = mod_pow(a, (mod + 3) >> 3, mod);
|
|
if (multiplication_modulo(res, res, mod) != a){
|
|
b = mod_pow(2, (mod - 1) >> 2, mod);
|
|
res = multiplication_modulo(res, b, mod);
|
|
}
|
|
break;
|
|
default :
|
|
if (a == 1)
|
|
res = 1;
|
|
else {
|
|
for (c = mod - 1, d = 2; d < mod && mod_pow(d, c >> 1, mod) != c; ++d);
|
|
for (e = 0; !(c & 1); ++e, c >>= 1);
|
|
f = mod_pow(a, c, mod);
|
|
b = mod_pow(d, c, mod);
|
|
for (h = 0, g = 0; h < e; h++) {
|
|
d = mod_pow(b, g, mod);
|
|
d = multiplication_modulo(d, f, mod);
|
|
d = mod_pow(d, 1 << (e - 1 - h), mod);
|
|
if (d == mod - 1)
|
|
g += 1 << h;
|
|
}
|
|
f = mod_pow(a, (c + 1) >> 1, mod);
|
|
b = mod_pow(b, g >> 1, mod);
|
|
res = multiplication_modulo(f, b, mod);
|
|
}
|
|
}
|
|
return res;
|
|
}
|
|
|
|
// return root such that (root * root) % mod congruent to n % mod.
|
|
// return 0 (the default value of a) if no solution to the congruence exists.
|
|
static unsigned tonelli_shanks_2(unsigned n, const unsigned mod) {
|
|
unsigned a = 0, b = mod - 1, c, d = b, e = 0, f = 2, g;
|
|
if (mod_pow(n, b >> 1, mod) == 1) {
|
|
for (; !(d & 1); ++e, d >>= 1);
|
|
if (e == 1)
|
|
a = mod_pow(n, (mod + 1) >> 2, mod);
|
|
else {
|
|
for (; b != mod_pow(f, b >> 1, mod); ++f);
|
|
for (b = mod_pow(f, d, mod), a = mod_pow(n, (d + 1) >> 1, mod), c = mod_pow(n, d, mod), g = e; c != 1; g = d) {
|
|
for (d = 0, e = c, --g; e != 1 && d < g; ++d)
|
|
e = multiplication_modulo(e, e, mod);
|
|
for (f = b, n = g - d; n--;)
|
|
f = multiplication_modulo(f, f, mod);
|
|
a = multiplication_modulo(a, f, mod);
|
|
b = multiplication_modulo(f, f, mod);
|
|
c = multiplication_modulo(c, b, mod);
|
|
}
|
|
}
|
|
}
|
|
return a;
|
|
}
|
|
|
|
#include <assert.h>
|
|
int main() {
|
|
unsigned n, mod, root ; /* root_2 = mod - root */
|
|
|
|
n = 27875, mod = 26371, root = tonelli_shanks_1(n, mod);
|
|
assert(root == 14320); // 14320 * 14320 mod 26371 = 1504 and 1504 = 27875 mod 26371
|
|
|
|
n = 1111111111, mod = 1111111121, root = tonelli_shanks_1(n, mod);
|
|
assert(root == 88664850);
|
|
|
|
n = 5258, mod = 3851, root = tonelli_shanks_1(n, mod);
|
|
assert(root == 0); // no solution to the congruence exists.
|
|
}
|