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