Anonymous user
Tonelli-Shanks algorithm: Difference between revisions
→{{header|C}}: add a new version of tonelli_shanks algorithm
m (→{{header|Phix}}: added syntax colouring the hard way) |
(→{{header|C}}: add a new version of tonelli_shanks algorithm) |
||
Line 1,015:
</pre>
=={{header|C}}==
=== Version 1 ===
{{trans|C#}}
<lang c>#include <stdbool.h>
Line 1,142 ⟶ 1,143:
root1 = 30468
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#}}==
| |||