2. Solution by continued fraction expansion

Continued fraction

We cannot solve Pell equation for large n by this primitive way.
For larger n, using the continued fraction expansion.

Following is the example of continued fraction expansion.
For sqrt(2), we can transform as follows;

sqrt(2) = 1 + sqrt(2) -1

                     1
        = 1 + ---------------
                     1
                -----------
                 sqrt(2)-1

                     1
        = 1 + ---------------
                 sqrt(2)+1
                -----------
                    2-1

                   1
        = 1 + -----------
               1+sqrt(2)

replacing "sqrt(2)" by itself and get the following fraction;

                      1
sqrt(2) = 1 + ----------------
                        1
             2 + -------------
                         1
                 2 + ---------

                     2 + ....

This is called as Continued fraction expansion.
In this chapter, representing this in abbreviated way;

sqrt(2) = 1 [2]
( [ ] denotes the repetition)

Another example : sqrt(3)

sqrt(3) = 1 + sqrt(3) -1

                     1
        = 1 + ---------------
                     1
                -----------
                 sqrt(3)-1

                     1
        = 1 + ---------------
                 sqrt(3)+1
                -----------
                    3-1

                       1
        = 1 + ------------------
                    sqrt(3)-1
               1 + -----------
                        2

                        1
        = 1 + ---------------------
                         3-1
               1 + --------------
                    2(sqrt(3)+1)

                       1
        = 1 + ------------------
                        1
               1 + -----------
                    1+sqrt(3)

Replacing sqrt(3) by itself

sqrt(3) = 1 [1,2]

Algorithm for continued fraction expansion and solution of Pell equation

Let a_n be the number of continued fraction expansion,
p_n/q_n be the partial sum up to n-th expansion,
and introduce the auxiliary parameter P_n, Q_n.
Computing n-th terms of the following recurrence formula.

an	a1 = int(sqrt(n)) ...... an = int((a1+Pn)/Q1)
Pn	P1 = 0 P2 = a1 ...... Pn = an-1 * Qn-1-Pn-1
Qn	Q1 = 1 Q2 = n-a12 ...... Qn = (n-Pn2) / Qn-1
pn	p1 = a1 p2 = a1 * a2+1 ...... pn = an * pn-1+pn-2
qn	q1 = 1 q2 = a2 ...... qn = an * qn-1+qn-2

When Q_n = 1, then p_n-1, q_n-1is the solution of Pell equation.

Program

 10   ' pell equation's solution by continue fraction method
 20   input N
 30   ' n=1
 40   Lp1=0:Lq1=1:A1=isqrt(N):Sp1=A1:Sq1=1
 50   ' n=2
 60   Lp2=A1:Lq2=N-A1^2:A2=int((A1+Lp2)/Lq2):Sp2=A1*A2+1:Sq2=A2
 70   if Lq2=1 then F=1:goto 150
 80   ' n
 90   Lp3=A2*Lq2-Lp2:Lq3=(N-Lp3^2)\Lq2:A3=int((isqrt(N)+Lp3)/Lq3)
100   Sp3=A3*Sp2+Sp1:Sq3=A3*Sq2+Sq1
110   if Lq3=1 then F=2:goto 150
120   Lp1=Lp2:Lp2=Lp3:Lq1=Lq2:Lq2=Lq3:Sp1=Sp2:Sp2=Sp3:Sq1=Sq2:Sq2=Sq3:A2=A3
130   goto 90
140   '
150   if F=1 then print Sp1,Sq1,Sp1^2-N*Sq1^2:end
160   print Sp2,Sq2,Sp2^2-N*Sq2^2:end

where

Lp := P
Lq := Q
Sp := p
Sq := q