# Chakravala method

The chakravala method (Sanskrit: चक्रवाल विधि) is a cyclic algorithm to solve indeterminate quadratic equations, including Pell's equation. It is commonly attributed to Bhāskara II, (c. 1114 – 1185 CE)[1][2] although some attribute it to Jayadeva (c. 950 ~ 1000 CE).[3] Jayadeva pointed out that Brahmagupta's approach to solving equations of this type could be generalized, and he then described this general method, which was later refined by Bhāskara II in his Bijaganita treatise. He called it the Chakravala method: chakra meaning "wheel" in Sanskrit, a reference to the cyclic nature of the algorithm.[4] E. O. Selenius held that no European performances at the time of Bhāskara, nor much later, exceeded its marvellous height of mathematical complexity.[1][4]

This method is also known as the cyclic method and contains traces of mathematical induction.[5]

## History

Brahmagupta in 628 CE studied indeterminate quadratic equations, including Pell's equation

$\,x^2 = Ny^2 + 1,$

for minimum integers x and y. Brahmagupta could solve it for several N, but not all.

Jayadeva (9th century) and Bhaskara (12th century) offered the first complete solution to the equation, using the chakravala method to find (for the notorious N = 61 case)

$\,x = 1 766 319 049$ and $\,y = 226 153 980.$

This case was first solved in Europe by Brouncker in 1657–58 in response to a challenge by Fermat, and a method first completely described by Lagrange in 1766.[6] Lagrange's method, however, requires the calculation of 21 successive convergents of the continued fraction for the square root of 61, while the chakravala method is much simpler. Selenius, in his assessment of the chakravala method, states

"The method represents a best approximation algorithm of minimal length that, owing to several minimization properties, with minimal effort and avoiding large numbers automatically produces the best solutions to the equation. The chakravala method anticipated the European methods by more than a thousand years. But no European performances in the whole field of algebra at a time much later than Bhaskara's, nay nearly equal up to our times, equalled the marvellous complexity and ingenuity of chakravala."[1][4]

Hermann Hankel calls the chakravala method

"the finest thing achieved in the theory of numbers before Lagrange."[7]

## The method

The chakravala method for solving Pell's equation is based on the observation by Brahmagupta (see Brahmagupta's identity) that

$(x_1^2 - Ny_1^2)(x_2^2 - Ny_2^2) = (x_1x_2 + Ny_1y_2)^2 - N(x_1y_2 + x_2y_1)^2$

This defines a "composition" (samāsa) of two triples $(x_1, y_1, k_1)$ and $(x_2, y_2, k_2)$ that are solutions of $x^2 - Ny^2 = k$, to generate the new triple

$(x_1x_2 + Ny_1y_2 \,,\, x_1y_2 + x_2y_1 \,,\, k_1k_2).$

In the general method, the main idea is that any triple $(a,b,k)$ (that is, one which satisfies $a^2 - Nb^2 = k$) can be composed with the trivial triple $(m, 1, m^2 - N)$ to get the new triple $(am + Nb, a+bm, k(m^2-N))$ for any m. Assuming we started with a triple for which $\gcd(a,b)=1$, this can be scaled down by k (this is Bhaskara's lemma):

$a^2 - Nb^2 = k \implies \left(\frac{am+Nb}{k}\right)^2 - N\left(\frac{a+bm}{k}\right)^2 = \frac{m^2-N}{k}$

Since the signs inside the squares do not matter, the following substitutions are possible:

$a\leftarrow\frac{am+Nb}{|k|}, b\leftarrow\frac{a+bm}{|k|}, k\leftarrow\frac{m^2-N}{k}$

When a positive integer m is chosen so that (a + bm)/k is an integer, so are the other two numbers in the triple. Among such m, the method chooses one that minimizes the absolute value of m2 − N and hence that of (m2 − N)/k. Then the substitution relations are applied for m equal to the chosen value. This results in a new triple (a, b, k). The process is repeated until a triple with $k=1$ is found. This method always terminates with a solution (proved by Lagrange in 1768).[8] Optionally, we can stop when k is ±1, ±2, or ±4, as Brahmagupta's approach gives a solution for those cases. Amendments to the Chakravala Method for Solving Pell's Equation




1) Amendment 1 - Addition of Variable, f, to Vector

The current method uses a vector ( x_current, y_current, k_current) with k_current = x_current^2-Ny_current^2. The first amendment is to add f to the vector ( f_new, x_new, y_new, k_new).




1.1) Definition of New Vector In Terms of the Current Vector

f_new = gcd(x_current, N)
x_new = x_current/gcd(x_current, N)
y_new = y_current

and


k_new = f_new.x_new^2-N.y_new^2/f_new




1.2) How Two New Vectors Combine

The result of combining 2 new vectors (f1, x1, y1, k1) and (f2, x2, y2, k2) is:




(1, f1.x1.f2.x2+N.y1.y2, f1.x1.y2+f2.x2.y1, f1.f2.k1.k2)




1.3) How a New Vector is Reduced

Reduction of a vector is a now 2-step process.




The first step is the current reduction step:




(f , x , y , k ) -> ( f' = f , x' = x /gcd( x , y ), y' = y/gcd(x, y), k' = k /gcd( x , y)^2 )




The second step is:




(f', x', y', k') -> ( f = f'.gcd( x', N/f'), x= x'/gcd( x', N/f'), y = y' , k= k'/gcd( x', N/f' ) )




1.4) Change to Finding a Starter Vector and Finding a Vector Through the BIG M Method

The method for finding a starter vector and for finding a vector through the BIG M method now should search for the optimum vector by iterating over all factors, f, of N where f < sqrt(N).




1.5) Transposing a Vector to Ensure f < N/f

Sometimes after a combination and reduction f > N/f. To keep f < N/f it will be necessary to employ a transpose transformation.
( f, x, y, k) -> ( f_transpose = N/f, x_transpose = y, y_transpose = x, k_transpose = -k)




2) Amendment 2 - Direct Solutions

The current method would require that when k satisfies the condition above that the vector combine with itself (or sometimes - when |k| = 4 - the first vector encountered with |k|=4) until k=1.




The second Amendment is to employ direct solutions.




2.1) Cases That Can Be Solved Directly During the Chakravala Cycle

In the following cases the solution can be obtained directly:




k=1 and f=1 - x_solution = x

                                          y_solution = y


|k| in { 1, 2} and ( f>1 or k<>1 ) - x_solution = (2fx^2-k)/|k|

                                          y_solution = 2xy/|k|


|k| = 4 and f is even - x_solution = (fx^2-2sgn(k))^2/2-1

                                          y_solution = xy(fx^2-2sgn(k))/2


k=4 and f=1 - x_solution = x(x^2-3)/2

                                          y_solution = y(x^2-1)/2


|k| = 4 and f is odd and ( f>1 or k=-4 ) - x_solution = fx^2(fx^2-3sgn(k))^2/2-sgn(k)

                                          y_solution = xy(fx^2-sgn(k))(fx^2-3sgn(k))/2


2.2) Cases That Can Be Solved Without the Need to Even Look for a Starter Vector




If N is of form:




a^2.b^2± a,
a^2.b^2±2a, or
a^2.b^2±4a




x_solution and y_solution can be obtained without looking for a starter by setting f to a, x to b, y to 1, k to -±1, -±2 or -±4 and solving using the direct solutions in 2.1).




Example N=44 is 2^2.3^2+4.2, so f=2, x=3, y=1, k=-4 and

x_solution = (fx^2-2sgn(k))^2/2-1 = 199,
y_solution = xy(fx^2-2sgn(k))/2 = 30.




Also, if N is of form:




a^2.b^2± a^2 and 4b.(2b^2±1) = 0 mod a,
a^2.b^2±2a^2 and 2b.( b^2±1) = 0 mod a,
a^2.b^2±4a^2 and b.( b^2±2) = 0 mod 2a and b is even, or
a^2.b^2±4a^2 and b.( b^2±2).[( b^2±2)^2+1].[( b^2±2)^2+3] = 0 mod 2a and b is odd




x_solution and y_solution can be obtained without looking for a starter by setting (f,x,y,k) to:




( a^2, (2b^2±1)/ a, 2b , 1),
( a^2, ( b^2±1)/ a, b , 1),
( a^2, ( b^2±2)/2a, b/2, 1), or
( a^2, ( b^2±2)/ a, b , 4)
respectively and solving.




Example N=153 is 3^2.4^2+3^2 and 4.4.(2.4^2+1) = 0 mod 3, so f=9, x=11, y=8, k=1 and




x_solution = (2fx^2-k)/|k| = 2177,
y_solution = 2xy/|k| = 176.




2.3 Direct Solutions for General Expressions for N

If N divides a1^2.b1^2 and a2^2.b2^2 in the ratio of:




w1a1:w2a2 and w1+w2 or w1w2 = ±1,±2 or ±4




A general expression for the solution can be obtained.




Example (1) 19 divides 4^2 and 5^2 in the ratio 1:2.




Add c1c2n to 4 and 5
c1c2n+4 c1c2n+5
Multiply the square of each expression by the other's ratio weighting
2(c1c2n+4)^2 (c1c2n+5)^2
Add the expressions together and divide by the sum of the ratio weightings
(3c1^2c2^2n^2+26c1c2n+57)/(1+2)
Set c1 to 3 to ensure coefficient of n is divisible by 3. (this is the purpose of c1)
9c2^2n^2+26c2n+19




This is the general expression for N!




Combine the reduced forms of ( 1, 3c2n+4, 1, -(2c2n+3)) and ( 1, 3c2n+5, 1, 2(2c2n+3))
f for each of them will be 1 because gcd( 3, 4, 9, 26, 19) = 1 and gcd( 3, 5, 9, 26, 19) = 1
( 1, 3c2n+4, 1, -(2c2n+3))
( 1, 3c2n+5, 1, 2(2c2n+3))
( 1, 18c2^2n^2+53c2n+39, 6c2n+9, -2(2c2n+3)^2)
( 1, (2c2n+3)(9c2n+13), 3(2c2n+3), -2(2c2n+3)^2) reduce
( 1, 9c2n+13, 3, -2)




The purpose of c2 is to ensure that after the second reduction step, |k| is 1 or 2, or, failing that, |k| = 4 and x is always odd regardless of n. k=-2 so set c2 to 1 and solve using the appropriate formula in 2.1.




N = 9n^2+26n+19
x_solution = (2fx^2-k)/|k| = (9n+13)^2+1
y_solution = 2xy/|k| = 3(9n+13)




This gives a direct solution for N= 19, 54, 107 etc.




Example (2) 55 divides 5^2 and 7^2 in the ratio 5:-1.




Add c1c2n to 7 and 7
c1c2n+5 c1c2n+7
Multiply the square of each expression by the other's ratio weighting
-(c1c2n+5)^2 5(c1c2n+7)^2
Add the expressions together and divide by the sum of the ratio weightings
(4c1^2c2^2n^2+60c1c2n+220)/(5-1)
Set c1 to 1; coefficient of n is already divisible by 4.
c2^2n^2+15c2n+55




Combine the reduced forms of ( 1, c2n+5, 1, -5(c2n+6)) and ( 1, c2n+7, 1, -(c2n+6))
f for each of them will be 1 because gcd( 1, 5, 1, 15, 55) = 1 and gcd( 1, 7, 1, 15, 55) = 1
( 1, c2n+5, 1, -5(c2n+6))
( 1, c2n+7, 1, - (c2n+6))
( 1, 2c2^2n^2+27c2n+90, 2c2n+12, 5(c2n+6)^2)
( 1, (c2n+6)(2c2n+15), 2(c2n+6), 5(c2n+6)^2) reduce
( 1, 2c2n+15, 2, 5)




The purpose of c2 is to ensure that after the second reduction step, |k| is 1 or 2, or, failing that, |k| = 4 and x is always odd regardless of n. k=5, so set c2 to 5. This will make x and k divisible by 5.
( 1, 10n+15, 2, 5) reduce using second reduction step
( 5, 2n+ 3, 2, 1) |k| is 1, so solve using the appropriate formula in 2.1




N = 25n^2+75n+55
x_solution = (2fx^2-k)/|k| = 10(2n+3)^2-1
y_solution = 2xy/|k| = 4(2n+3)




This gives a direct solution for N= 55, 155, 305 etc.




Example (3) 57 divides 6^2 and 8^2 in the ratio 3:1.




Add c1c2n to 6 and 8
c1c2n+6 c1c2n+8
Multiply the square of each expression by the other's ratio weighting
(c1c2n+6)^2 3(c1c2n+8)^2
Add the expressions together and divide by the sum of the ratio weightings
(4c1^2c2^2n^2+60c1c2n+228)/(3+1)
Set c1 to 1; coefficient of n is already divisible by 4.
c2^2n^2+15c2n+57




Combine the reduced forms of ( 1, c2n+6, 1, -3(c2n+7)) and ( 1, c2n+8, 1, (c2n+7))
f for each of them will be 1 because gcd( 1, 6, 1, 15, 57) = 1 and gcd( 1, 8, 1, 15, 57) = 1
( 1, c2n+6, 1, -3(c2n+7))
( 1, c2n+8, 1, (c2n+7))
( 1, 2c2^2n^2+29c2n+105, 2c2n+14, -3(c2n+7)^2)
( 1, (c2n+7)(2c2n+15), 2(c2n+7), -3(c2n+7)^2) reduce
( 1, 2c2n+15, 2, -3)




The purpose of c2 is to ensure that after the second reduction step, |k| is 1 or 2, or, failing that, |k| = 4 and x is always odd regardless of n. k=-3, so set c2 to 3. This will make x and k divisible by 3.
( 1, 6n+15, 2, -3) reduce using second reduction step
( 3, 2n+ 5, 2, -1) |k| is 1, so solve using the appropriate formula in 2.1




N = 9n^2+45n+57
x_solution = (2fx^2-k)/|k| = 6(2n+5)^2+1
y_solution = 2xy/|k| = 4(2n+5)




This gives a direct solution for N= 57, 111, 183 etc.




Example (4) 61 divides 7^2 and 8^2 in the ratio 4:1.




Add c1c2n to 7 and 8
c1c2n+7 c1c2n+8
Multiply the square of each expression by the other's ratio weighting
(c1c2n+7)^2 4(c1c2n+8)^2
Add the expressions together and divide by the sum of the ratio weightings
(5c1^2c2^2n^2+78c1c2n+305)/(4+1)
Set c1 to 5; to ensure coefficient of n is divisible by 5.
25c2^2n^2+78c2n+61




Combine the reduced forms of ( 1, 5c2n+7, 1, -4(2c2n+3)) and ( 1, 5c2n+8, 1, (2c2n+3))
f for each of them will be 1 because gcd( 5, 7, 25, 78, 61) = 1 and gcd( 5, 8, 25, 78, 61) = 1
( 1, 5c2n+7, 1, -4(2c2n+3))
( 1, 5c2n+8, 1, (2c2n+3))
( 1, 50c2^2n^2+153c2n+117, 10c2n+15, -4(2c2n+3)^2)
( 1, (2c2n+3)(25c2n+39), 5(2c2n+3), -4(2c2n+3)^2) reduce
( 1, 25c2n+39, 5, -4)




The purpose of c2 is to ensure that after the second reduction step, |k| is 1 or 2, or, failing that, |k| = 4 and x is always odd regardless of n. k=-4, so is there a value of c2 that would make gcd( x, |k|) (= gcd( 25c2, 39, 4) )= 2 or 4. No. If c2 were 1, would x always be odd regardless of n (coefficient of n is even and constant is odd)? No. Set c2 to 2 so that x will always be odd.
( 1, 50n+39, 5, -4)




N = 100n^2+156n+61
x_solution = fx^2(fx^2+3)^2/2-sgn(k) = (50n+39)^2[(50n+39)^2+3]^2/2+1
y_solution = xy(fx^2+1)(fx^2+3)/2 = 5(50n+39)[(50n+39)^2+1][(50n+39)^2+3]/2




This gives a direct solution for N= 61, 317, 773 etc.




Example (5) 164 divides 12^2 and 13^2 in the ratio 4:1 (this is the example for 61 "shifted" by c1. 12=7+c1, 13=8+c1)




Add c1c2n to 12 and 13
c1c2n+12 c1c2n+13
Multiply the square of each expression by the other's ratio weighting
(c1c2n+12)^2 4(c1c2n+13)^2
Add the expressions together and divide by the sum of the ratio weightings
(5c1^2c2^2n^2+128c1c2n+820)/(4+1)
Set c1 to 5; to ensure coefficient of n is divisible by 5.
25c2^2n^2+128c2n+164




Combine the reduced forms of ( 1, 5c2n+12, 1, -4(2c2n+5)) and ( 1, 5c2n+13, 1, (2c2n+3))
f for each of them will be 1 because gcd( 5, 12, 25, 128, 164) = 0 and gcd( 5, 13, 25, 128, 164) = 0
( 1, 5c2n+12, 1, -4(2c2n+5))
( 1, 5c2n+13, 1, (2c2n+5))
( 1, 50c2^2n^2+253c2n+320, 10c2n+25, -4(2c2n+5)^2)
( 1, (2c2n+5)(25c2n+64), 5(2c2n+5), -4(2c2n+3)^2) reduce
( 1, 25c2n+64, 5, -4)




The purpose of c2 is to ensure that after the second reduction step, |k| is 1 or 2, or, failing that, |k| = 4 and x is always odd regardless of n. k=-4, so is there a value of c2 that would make gcd( x, |k|) (= gcd( 25c2, 64, 4) )= 2 or 4. Yes, set c2 = 2.




( 1, 50n+64, 5, -4) reduce using second reduction step
( 2, 25n+32, 5, -2) |k| is 2, so solve using the appropriate formula in 2.1




N = 100n^2+256n+164
x_solution = (2fx^2-k)/|k| = 4(25n+32)^2+1
y_solution = 2xy/|k| = 10(25n+32)




This gives a direct solution for N= 164, 520, 1076 etc.




Amendment 3 - Before BIG M, Search Previously Encountered Vectors for a Vector With Which to Combine

With the current method, if |k| is not in { 1, 2, 4}, the only option is to employ the BIG M method to obtain a vector. However, a previously encountered vector would be suitable for combining with the current vector, if:




1) it is not the current vector or its transpose; and




2)|k_prev.k_current| equals




i) 1, 2 or 4 times a prime squared; or
ii) 1, 2 or 4 times a prime cubed; or
iii) 64, 128 or 256 and ( k_prev = 8 or k_curr = 8 ) and

            mod( f_prev.x_prev.y_curr+f_curr.x_curr.y_prev, 4)=0); or


iv) 1, 2 or 4 times prime_1^2.prime_2

           ( and there is a second previously encountered vector where |k| = prime_2 )


In light of this, if no direct solution is possible, one of a pair (or triplet for iii) of candidate vectors that satisfy one of the above conditions is a better starter vector than the one with lowest |k|.




Working Examples

N=28 (example of direct solution with N of form a^2.b^2-4.a)




Current Method

Starter find ( floor(sqrt(N)), 1, floor(sqrt(N))^2-N)

       and  ( ceil(sqrt(N)), 1, ceil(sqrt(N))^2-N)


( 5, 1, -3)
( 6, 1, 8)
neither vector has |k| in { 1, 2, 4}, so choose vector with lower |k|
( 5, 1, -3) |k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-28)
(5M+28, 5+M, -3(M^2-28)) we want mod(y=M+5,|-3|)=0 and |-3(M^2-28)| to be minimised => M=4
( 48, 9, 36) reduce
( 16, 3, 4) |k| = 4, so combine with self
( 16, 3, 4)
( 508, 96, 16) reduce
( 127, 24, 1) k = 1, so solve
x_solution = x = 127
y_solution = y = 24




New Method

N=2^2.3^2-4.2 so
x_solution = (2.3^2-2.1)^2/2-1 = 127
y_solution = 3.1.(2.3^2-2.1)/2 = 24




N=46 (example where a previous candidate satisfies condition ii) in Amendment 3 and BIG M under Current Method would not

    return the vector, but under the New Method it would because of the extra reduction step)


Current Method

Starter find ( floor(sqrt(N)), 1, floor(sqrt(N))^2-N)

       and  ( ceil(sqrt(N)), 1, ceil(sqrt(N))^2-N)


( 6, 1, -10)
( 7, 1, 3)
neither vector has |k| in { 1, 2, 4}, so choose vector with lower |k|
( 7, 1, 3)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-46)
( 7M+46, M+7, 3(M^2-46)) we want mod(y=M+7,|3|)=0 and |3(M^2-46)| to be minimised => M=8
( 102, 15, 54) reduce
( 34, 5, 6)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-46)
( 34M+230, 5M+34, 6(M^2-46)) we want mod(y=5M+34,|6|)=0 and |6(M^2-46)| to be minimised => M=4
( 366, 54, -180)
( 61, 9, -5)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-46)
( 61M+414, 9M+61, -5(M^2-46)) we want mod(y=9M+61,|-5|)=0 and |-5(M^2-46)| to be minimised => M=6
( 780, 115, 50)
( 156, 23, 2) combine with self
( 156, 23, 2)
( 48670, 7176, 4) reduce
( 24335, 3588, 1) k = 1, so solve
x_solution = x = 24335
y_solution = y = 3588




New Method

Starter find ( f, floor(sqrt(N)/f), 1, f.floor(sqrt(N)/f)^2-N/f)

       and  ( f, ceil(sqrt(N)/f), 1, f.ceil(sqrt(N)/f)^2-N/f)


for f = all factors of N < sqrt(N)
Set f=1
( 1, 6, 1, -10) ignored because gcd( x, N/f) > 1
( 1, 7, 1, 3)
Set f=2
( 2, 3, 1, -5)
( 2, 4, 1, 9)
No candidate has |k| in { 1, 2, 4} but ( 1, 7, 1, 3) and ( 2, 4, 1, 9) form a pair that satisfies ii) of Amendment 3
( 1, 7, 1, 3)
( 2, 4, 1, 9)
( 1, 102, 15, 54) reduce
( 2, 17, 5, 3)
( Note equivalent vector under Current Method at this stage had |k| = 6 and would not have allowed BIG M to

 return ( 7, 1, 3), but under New Method if BIG M was used it would return ( 1, 7, 1, 3). )


( 1, 7, 1, 3) |k| is not { 1, 2, 4} but ( 1, 7, 1, 3) satisfies case i) of Amendment 3
( 1, 468, 69, 18)
( 2, 78, 23, 1) |k| is in { 1, 2, 4} so we can stop and use this vector to solve the equation directly
x_solution = (2fx^2-k)/|k| = 24335
y_solution = 2xy/|k| = 3588




N=67 ( example where a previous candidate satisfies condition i) in Amendment 3 that would not have been returned by BIG M,

      because |k| was not a multiple of |k| from current vector)


Current Method

Starter find ( floor(sqrt(N)), 1, floor(sqrt(N))^2-N)

       and  ( ceil(sqrt(N)), 1, ceil(sqrt(N))^2-N)


( 8, 1, -3)
( 9, 1, 14)
neither vector has |k| in { 1, 2, 4}, so choose vector with lower |k|
( 8, 1, -3)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-67)
( 8M+67, M+8, -3(M^2-67)) we want mod(y=M+8,|-3|)=0 and |-3(M^2-67)| to be minimised => M=7
( 123, 15, 54) reduce
( 41, 5, 6)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-67)
( 41M+335, 5M+41, 6(M^2-67)) we want mod(y=5M+41,|6|)=0 and |6(M^2-67)| to be minimised => M=5
( 540, 66, -252) reduce
( 90, 11, -7)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-67)
( 90M+737, 11M+90, -7(M^2-67)) we want mod(y=11M+90,|-7|)=0 and |-7(M^2-67)| to be minimised => M=9
( 1547, 189, -98) reduce
( 221, 27, -2) combine with self
( 221, 27, -2)
( 97684, 11934, 4) reduce
( 48842, 5967, 1) k = 1, so solve
x_solution = x = 48842
y_solution = y = 5967




New Method

Starter find ( f, floor(sqrt(N)/f), 1, f.floor(sqrt(N)/f)^2-N/f)

       and  ( f, ceil(sqrt(N)/f), 1, f.ceil(sqrt(N)/f)^2-N/f)


for f = all factors of N < sqrt(N)
Set f=1
( 1, 8, 1, -3)
( 1, 9, 1, 14)




No candidate has |k| in { 1, 2, 4} and no cases in amendment 3 hold, so starter is vector with lowest |k|
( 1, 8, 1, -3) use BIG M
( F, M, 1, FM^2-67/F)
( 1, 8FM+67, FM+8, -3(F^2M^2-67)) we want gcd(x/F=(8MF+67)/F, N/F=67/F)=1, mod(y=FM+8,|-3|)=0 and

                                         |-3(F.M^2-67/F)| to be minimised => F=1,M=7


( 1, 123, 15, 54) reduce
( 1, 41, 5, 6) |k| is not { 1, 2, 4} but ( 1, 8, 1, -3) satisfies case (i) of Amendment 3
( 1, 8, 1, - 3) ( Note BIG M would not have returned same vector )
( 1, 663, 81, -18)
( 1, 221, 27, - 2) |k| is in { 1, 2, 4} so we can stop and use this vector to solve the equation directly
x_solution = (2fx^2-k)/|k| = 48842
y_solution = 2xy/|k| = 5967




N=97 ( example where a previous candidate satisfies condition iii) in Amendment 3)




Current Method

Starter find ( floor(sqrt(N)), 1, floor(sqrt(N))^2-N)

       and  (  ceil(sqrt(N)), 1,  ceil(sqrt(N))^2-N)


( 9, 1, -16)
(10, 1, 3)
neither vector has |k| in { 1, 2, 4}, so choose vector with lower |k|
(10, 1, 3)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-97)
(10M+97, M+10, 3(M^2-97)) we want mod(y=M+10,|3|)=0 and |3(M^2-97)| to be minimised => M=11
( 207, 21, 72) reduce
( 69, 7, 8)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-97)
(69M+679, 7M+69, 8(M^2-97)) we want mod(y=7M+69,|8|)=0 and |8(M^2-97)| to be minimised => M=5

                           (tie between M=5 and M=13)


( 1024, 104, -576) reduce
( 128, 13, - 9)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-97)
(128M+1261, 13M+128, -9(M^2-97)) we want mod(y=13M+128,|-9|)=0 and |-9(M^2-97)| to be minimised => M=13
( 2925, 297, -648) reduce
( 325, 33, - 8)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-97)
(325M+3201, 33M+325, -8(M^2-97)) we want mod(y=33M+325,|-8|)=0 and |-8(M^2-97)| to be minimised => M=11
(6776, 688, -192) reduce
( 847, 86, - 3)|k| not in { 1, 2, 4}, so use BIG M
( M, 1, M^2-97)
(847M+8342, 86M+847, -3(M^2-97)) we want mod(y=386M+847,|-3|)=0 and |-3(M^2-97)| to be minimised => M=10
(16812, 1707, -9) reduce
( 5604, 569, -1) combine with self
( 5604, 569, -1)
( 62809633, 6377352, 1) k = 1, so solve
x_solution = x = 62809633
y_solution = y = 6377352




New Method

Starter find ( f, floor(sqrt(N)/f), 1, f.floor(sqrt(N)/f)^2-N/f)

       and  ( f, ceil(sqrt(N)/f), 1, f.ceil(sqrt(N)/f)^2-N/f)


for f = all factors of N < sqrt(N)




Set f=1
( 1, 9, 1, -16)
( 1, 10, 1, 3)
No candidate has |k| in { 1, 2, 4} and no cases in amendment 3 hold, so starter is vector with lowest |k|
( 1, 10, 1, 3) use BIG M
( F, M, 1, FM^2-97/F)
( 1, 10FM+97, FM+10, 3(F^2M^2-97)) we want gcd(x/F=(10FM+97)/F, N/F=97/F)=1, mod(y=FM+10,|3|)=0 and

                                         |3(F.M^2-97/F)| to be minimised => F=1,M=11


( 1, 207, 21, 72) reduce
( 1, 69, 7, 8) |k| is not { 1, 2, 4} but ( 1, 9, 1, -16) satisfies case iii) of Amendment 3

                 |8.-16| = 128, mod(1x9x7+1x69x1,4)=0


( 1, 9, 1,-16)
( 1, 1300, 132, -128) reduce
( 1, 325, 33, -8) |k| is not { 1, 2, 4} but ( 1, 69, 7, 8) satisfies case iii) of Amendment 3

                 |-8.8| = 64, mod(1x325x7+1x69x33,4)=0


( 1, 69, 7, 8)
( 1, 44832, 4552, -64) reduce
( 1, 5604, 569, - 1) |k| is in { 1, 2, 4} so we can stop and use this vector to solve the equation directly
x_solution = (2fx^2-k)/|k| = 62809633
y_solution = 2xy/|k| = 6377352

## Examples

### n = 61

The n = 61 case (determining an integer solution satisfying $a^2 - 61b^2 = 1$), issued as a challenge by Fermat many centuries later, was given by Bhaskara as an example.[8]

We start with a solution $a^2 - 61b^2 = k$ for any k found by any means. In this case we can let b be 1, thus, since $8^2 - 61\cdot1^2 = 3$, we have the triple $(a,b,k) = (8, 1, 3)$. Composing it with $(m, 1, m^2-61)$ gives the triple $(8m+61, 8+m, 3(m^2-61))$, which is scaled down (or Bhaskara's lemma is directly used) to get:

$\left( \frac{8m+61}{3}, \frac{8+m}{3}, \frac{m^2-61}{3} \right).$

For 3 to divide $8+m$ and $|m^2-61|$ to be minimal, we choose $m=7$, so that we have the triple $(39, 5, -4)$. Now that k is −4, we can use Brahmagupta's idea: it can be scaled down to the rational solution $(39/2, 5/2, -1)\,$, which composed with itself three times, with $m={7,11,9}$ respectively, when k becomes square and scaling can be applied, this gives $(1523/2, 195/2, 1)\,$. Finally, such procedure can be repeated until the solution is found (requiring 9 additional self-compositions and 4 additional square-scalings): $(1766319049,\, 226153980,\, 1)$. This is the minimal integer solution.

### n = 67

Suppose we are to solve $x^2 - 67y^2 = 1$ for x and y.[9]

We start with a solution $a^2 - 67b^2 = k$ for any k found by any means; in this case we can let b be 1, thus producing $8^2 - 67\cdot1^2 = -3$. At each step, we find an m > 0 such that k divides a + bm, and |m2 − 67| is minimal. We then update a, b, and k to $\frac{am+Nb}{|k|}, \frac{a+bm}{|k|}, \text{ and }\frac{m^2-N}{k}$ respectively.

First iteration

We have $(a,b,k) = (8,1,-3)$. We want a positive integer m such that k divides a + bm, i.e. 3 divides 8 + m, and |m2 − 67| is minimal. The first condition implies that m is of the form 3t + 1 (i.e. 1, 4, 7, 10,… etc.), and among such m, the minimal value is attained for m = 7. Replacing (abk) with $\left(\frac{am+Nb}{|k|}, \frac{a+bm}{|k|}, \frac{m^2-N}{k}\right)$, we get the new values $a = (8\cdot7+67\cdot1)/3 = 41, b = (8 + 1\cdot7)/3 = 5, k = (7^2-67)/(-3) = 6$. That is, we have the new solution:

$41^2 - 67\cdot(5)^2 = 6.$

At this point, one round of the cyclic algorithm is complete.

Second iteration

We now repeat the process. We have $(a,b,k) = (41,5,6)$. We want an m > 0 such that k divides a + bm, i.e. 6 divides 41 + 5m, and |m2 − 67| is minimal. The first condition implies that m is of the form 6t + 5 (i.e. 5, 11, 17,… etc.), and among such m, |m2 − 67| is minimal for m = 5. This leads to the new solution a = (41⋅5 + 67⋅5)/6, etc.:

$90^2 - 67 \cdot 11^2 = -7.$
Third iteration

For 7 to divide 90 + 11m, we must have m = 2 + 7t (i.e. 2, 9, 16,… etc.) and among such m, we pick m = 9.

$221^2 - 67\cdot 27^2 = -2.$
Final solution

At this point, we could continue with the cyclic method (and it would end, after seven iterations), but since the right-hand side is among ±1, ±2, ±4, we can also use Brahmagupta's observation directly. Composing the triple (221, 27, −2) with itself, we get

$\left(\frac{221^2 + 67\cdot27^2}{2}\right)^2 - 67\cdot(221\cdot27)^2 = 1,$

that is, we have the integer solution:

$48842^2 - 67 \cdot 5967^2 = 1.$

This equation approximates $\sqrt{67}$ as $\frac{48842}{5967}$ to within a margin of about $2 \times 10^{-9}$.

## Notes

1. ^ a b c Hoiberg & Ramchandani – Students' Britannica India: Bhaskaracharya II, page 200
2. ^ Kumar, page 23
3. ^ Plofker, page 474
4. ^ a b c Goonatilake, page 127 – 128
5. ^ Cajori (1918), p. 197

"The process of reasoning called "Mathematical Induction" has had several independent origins. It has been traced back to the Swiss Jakob (James) Bernoulli, the Frenchman B. Pascal and P. Fermat, and the Italian F. Maurolycus. [...] By reading a little between the lines one can find traces of mathematical induction still earlier, in the writings of the Hindus and the Greeks, as, for instance, in the "cyclic method" of Bhaskara, and in Euclid's proof that the number of primes is infinite."

6. ^
7. ^ Kaye (1919), p. 337.
8. ^ a b John Stillwell (2002), Mathematics and its history (2 ed.), Springer, pp. 72–76, ISBN 978-0-387-95336-6
9. ^ The example in this section is given (with notation $Q_n$ for k, $P_n$ for m, etc.) in: Michael J. Jacobson; Hugh C. Williams (2009), Solving the Pell equation, Springer, p. 31, ISBN 978-0-387-84922-5

## References

• Florian Cajori (1918), Origin of the Name "Mathematical Induction", The American Mathematical Monthly 25 (5), p. 197-201.
• George Gheverghese Joseph, The Crest of the Peacock: Non-European Roots of Mathematics (1975).
• G. R. Kaye, "Indian Mathematics", Isis 2:2 (1919), p. 326–356.
• C. O. Selenius, "Rationale of the chakravala process of Jayadeva and Bhaskara II", Historia Mathematica 2 (1975), pp. 167-184.
• C. O. Selenius, "Kettenbruch theoretische Erklarung der zyklischen Methode zur Losung der Bhaskara-Pell-Gleichung", Acta Acad. Abo. Math. Phys. 23 (10) (1963).
• Hoiberg, Dale & Ramchandani, Indu (2000). Students' Britannica India. Mumbai: Popular Prakashan. ISBN 0-85229-760-2
• Goonatilake, Susantha (1998). Toward a Global Science: Mining Civilizational Knowledge. Indiana: Indiana University Press. ISBN 0-253-33388-1.
• Kumar, Narendra (2004). Science in Ancient India. Delhi: Anmol Publications Pvt Ltd. ISBN 81-261-2056-8
• Ploker, Kim (2007) "Mathematics in India". The Mathematics of Egypt, Mesopotamia, China, India, and Islam: A Sourcebook New Jersey: Princeton University Press. ISBN 0-691-11485-4
• Edwards, Harold (1977). Fermat's Last Theorem. New York: Springer. ISBN 0-387-90230-9.