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