/*
 * pell - solve Pell's equation
 *
 * Copyright (C) 1999  David I. Bell
 *
 * Calc is open software; you can redistribute it and/or modify it under
 * the terms of the version 2.1 of the GNU Lesser General Public License
 * as published by the Free Software Foundation.
 *
 * Calc is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 * or FITNESS FOR A PARTICULAR PURPOSE.	 See the GNU Lesser General
 * Public License for more details.
 *
 * A copy of version 2.1 of the GNU Lesser General Public License is
 * distributed with calc under the filename COPYING-LGPL.  You should have
 * received a copy with calc; if not, write to Free Software Foundation, Inc.
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 *
 * @(#) $Revision: 29.1 $
 * @(#) $Id: pell.cal,v 29.1 1999/12/14 09:15:31 chongo Exp $
 * @(#) $Source: /usr/local/src/cmd/calc/cal/RCS/pell.cal,v $
 *
 * Under source code control:	1990/02/15 01:50:34
 * File existed as early as:	before 1990
 *
 * Share and enjoy!  :-)	http://reality.sgi.com/chongo/tech/comp/calc/
 */

/*
 * Solve Pell's equation; Returns the solution X to: X^2 - D * Y^2 = 1.
 * Type the solution to pells equation for a particular D.
 */


define pell(D)
{
	local X, Y;

	X = pellx(D);
	if (isnull(X)) {
		print "D=":D:" is square";
		return;
	}
	Y = isqrt((X^2 - 1) / D);
	print X : "^2 - " : D : "*" : Y : "^2 = " : X^2 - D*Y^2;
}


/*
 * Function to solve Pell's equation
 * Returns the solution X to:
 *	X^2 - D * Y^2 = 1
 */
define pellx(D)
{
	local R, Rp, U, Up, V, Vp, A, T, Q1, Q2, n;
	local mat ans[2,2];
	local mat tmp[2,2];

	R = isqrt(D);
	Vp = D - R^2;
	if (Vp == 0)
		return;
	Rp = R + R;
	U = Rp;
	Up = U;
	V = 1;
	A = 0;
	n = 0;
	ans[0,0] = 1;
	ans[1,1] = 1;
	tmp[0,1] = 1;
	tmp[1,0] = 1;
	do {
		T = V;
		V = A * (Up - U) + Vp;
		Vp = T;
		A = U // V;
		Up = U;
		U = Rp - U % V;
		tmp[0,0] = A;
		ans *= tmp;
		n++;
	} while (A != Rp);
	Q2 = ans[[1]];
	Q1 = isqrt(Q2^2 * D + 1);
	if (isodd(n)) {
		T = Q1^2 + D * Q2^2;
		Q2 = Q1 * Q2 * 2;
		Q1 = T;
	}
	return Q1;
}