/* * Copyright (c) 1993 David I. Bell * Permission is granted to use, distribute, or modify this source, * provided that this copyright notice remains intact. * * Calculate square roots modulo a prime. * * Returns null if number is not prime or if there is no square root. * The smaller square root is always returned. */ define psqrt(u, p) { local p1, q, n, y, r, v, w, t, k; p1 = p - 1; r = lowbit(p1); q = p >> r; t = 1 << (r - 1); for (n = 2; ; n++) { if (ptest(n, 1) == 0) continue; y = pmod(n, q, p); k = pmod(y, t, p); if (k == 1) continue; if (k != p1) return; break; } t = pmod(u, (q - 1) / 2, p); v = (t * u) % p; w = (t^2 * u) % p; while (w != 1) { k = 0; t = w; do { k++; t = t^2 % p; } while (t != 1); if (k == r) return; t = pmod(y, 1 << (r - k - 1), p); y = t^2 % p; v = (v * t) % p; w = (w * y) % p; r = k; } return min(v, p - v); } global lib_debug; if (!isnum(lib_debug) || lib_debug>0) print "psqrt(u, p) defined";