1/*
2 * pell - solve Pell's equation
3 *
4 * Copyright (C) 1999,2021  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:34
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 * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
28 * Type the solution to Pell's equation for a particular D.
29 */
30
31
32define pell(D)
33{
34	local X, Y;
35
36	X = pellx(D);
37	if (isnull(X)) {
38		print "D=":D:" is square";
39		return;
40	}
41	Y = isqrt((X^2 - 1) / D);
42	print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
43}
44
45
46/*
47 * Function to solve Pell's equation
48 * Returns the solution X to:
49 *	X^2 - D * Y^2 = 1
50 */
51define pellx(D)
52{
53	local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n;
54	local mat ans[2,2];
55	local mat tmp[2,2];
56
57	R = isqrt(D);
58	Vp = D - R^2;
59	if (Vp == 0)
60		return;
61	Rp = R + R;
62	U = Rp;
63	Up = U;
64	V = 1;
65	A = 0;
66	n = 0;
67	ans[0,0] = 1;
68	ans[1,1] = 1;
69	tmp[0,1] = 1;
70	tmp[1,0] = 1;
71	do {
72		T = V;
73		V = A * (Up - U) + Vp;
74		Vp = T;
75		A = U // V;
76		Up = U;
77		U = Rp - U % V;
78		tmp[0,0] = A;
79		ans *= tmp;
80		n++;
81	} while (A != Rp);
82	Q2 = ans[[1]];
83	Q1 = isqrt(Q2^2 * D + 1);
84	if (isodd(n)) {
85		T = Q1^2 + D * Q2^2;
86		Q2 = Q1 * Q2 * 2;
87		Q1 = T;
88	}
89	return Q1;
90}
91