1/* 2 * sumsq - find unique two positive integers whose squares sum to a given 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:37 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 * Determine the unique two positive integers whose squares sum to the 28 * specified prime. This is always possible for all primes of the form 29 * 4N+1, and always impossible for primes of the form 4N-1. 30 */ 31 32 33define ss(p) 34{ 35 local a, b, i, p4; 36 37 if (p == 2) { 38 print "1^2 + 1^2 = 2"; 39 return; 40 } 41 if ((p % 4) != 1) { 42 print p, "is not of the form 4N+1"; 43 return; 44 } 45 if (!ptest(p, min(p-2, 10))) { 46 print p, "is not a prime"; 47 return; 48 } 49 p4 = (p - 1) / 4; 50 i = 2; 51 do { 52 a = pmod(i++, p4, p); 53 } while ((a^2 % p) == 1); 54 b = p; 55 while (b^2 > p) { 56 i = b % a; 57 b = a; 58 a = i; 59 } 60 print a : "^2 +" , b : "^2 =" , a^2 + b^2; 61} 62