Solves the Brahmagupta (Pell) equation using the chakravala method

chakravala.gp

/*
Finds the minimal solution to a Pell's equation, a^2 - nb^2 = 1 for given n

Inputs:
	n: Non-negative integer
	
Output:
	a,b: Minimal solution to Pell's equation
	
Example:
chakravala(13)

Written by: John Peach 14-Mar-2021 
Wild Peaches               

*/

chakravala(n) = 
{
	\\ Check that n is not square
	if(issquare(n),error("n = ", n, " is square."));

	\\ Initial values for a,b,k
	sqrt_n = sqrt(n);
	a = ceil(sqrt_n);
	b = 1;
	k = a^2 - n;
	
	steps = 1;
		
	\\ Repeat until k = 1
	while(k!=1,
		
		abs_k = abs(k);
		
		\\ Find m s.t. k | (a+bm), select m that minimizes |m^2-n|		
		if(abs_k == 1, m = round(sqrt_n),			
			r = -a % abs_k;    \\ Negative of a mod(k)
			s = (1/b) % abs_k; \\ Inverse of b mod(k) (PARI calculates extended Euclidean internally)
			t = (r*s) % abs_k; \\ Remainder of (a + b x m) when divided by k
			m = t + floor((sqrt_n-t)/k)*k; \\ Approximate value for m, may need slight adjustment
			\\ Adjust m until m < sqrt(n) < m+k
			while(m > sqrt_n, m -= abs_k);
			while(m + abs_k < sqrt_n, m += abs_k);
			\\ Select closest value to n
			if(abs(m^2 - n) > abs((m+abs_k)^2 - n),m = m + abs_k) ); 
			
		\\ Print current triple
		print([a,b,k]);
	
		\\ Update a,b,k
		alpha = a;
		a = (m*a + n*b)/abs_k;
		b = (alpha + m*b)/abs_k;
		k = (m^2 - n)/k;
		
		\\ Increment step counter
		steps += 1;
		
	);
	
	\\ Print final result
	print(a "^2 - " n " x " b "^2 = 1 in " steps " calculations.");
	
	\\ Return vector of [a,b,k]
	return([a,b,k]);	
}