Tonelli-Shanks algorithm: Difference between revisions
Content added Content deleted
m (→{{header|Phix}}: added syntax colouring the hard way) |
(→{{header|C}}: add a new version of tonelli_shanks algorithm) |
||
Line 1,015: | Line 1,015: | ||
</pre> |
</pre> |
||
=={{header|C}}== |
=={{header|C}}== |
||
=== Version 1 === |
|||
{{trans|C#}} |
{{trans|C#}} |
||
<lang c>#include <stdbool.h> |
<lang c>#include <stdbool.h> |
||
Line 1,142: | Line 1,143: | ||
root1 = 30468 |
root1 = 30468 |
||
root2 = 69581</pre> |
root2 = 69581</pre> |
||
=== Version 2 === |
|||
<lang 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. |
|||
} |
|||
</lang> |
|||
A is assumed odd prime, the algorithm requires O(log A + r * r) multiplications modulo A, where r is the power of 2 dividing A − 1. |
|||
=={{header|C sharp|C#}}== |
=={{header|C sharp|C#}}== |