1/*
2 * psqrt - calculate square roots modulo a prime
3 *
4 * Copyright (C) 1999  David I. Bell
5 *
6 * Calc is open software; you can redistribute it and/or modify it under
7 * the terms of the version 2.1 of the GNU Lesser General Public License
8 * as published by the Free Software Foundation.
9 *
10 * Calc is distributed in the hope that it will be useful, but WITHOUT
11 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
13 * Public License for more details.
14 *
15 * A copy of version 2.1 of the GNU Lesser General Public License is
16 * distributed with calc under the filename COPYING-LGPL.  You should have
17 * received a copy with calc; if not, write to Free Software Foundation, Inc.
18 * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
19 *
20 * Under source code control:	1990/02/15 01:50:35
21 * File existed as early as:	before 1990
22 *
23 * Share and enjoy!  :-)	http://www.isthe.com/chongo/tech/comp/calc/
24 */
25
26/*
27 * Returns null if number is not prime or if there is no square root.
28 * The smaller square root is always returned.
29 */
30
31
32define psqrt(u, p)
33{
34	local	p1, q, n, y, r, v, w, t, k;
35
36	p1 = p - 1;
37	r = lowbit(p1);
38	q = p >> r;
39	t = 1 << (r - 1);
40	for (n = 2; ; n++) {
41		if (ptest(n, 1) == 0)
42			continue;
43		y = pmod(n, q, p);
44		k = pmod(y, t, p);
45		if (k == 1)
46			continue;
47		if (k != p1)
48			return;
49		break;
50	}
51	t = pmod(u, (q - 1) / 2, p);
52	v = (t * u) % p;
53	w = (t^2 * u) % p;
54	while (w != 1) {
55		k = 0;
56		t = w;
57		do {
58			k++;
59			t = t^2 % p;
60		} while (t != 1);
61		if (k == r)
62			return;
63		t = pmod(y, 1 << (r - k - 1), p);
64		y = t^2 % p;
65		v = (v * t) % p;
66		w = (w * y) % p;
67		r = k;
68	}
69	return min(v, p - v);
70}
71