## 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 an be the number of continued fraction expansion,
pn/qn be the partial sum up to n-th expansion,
and introduce the auxiliary parameter Pn, Qn.
Computing n-th terms of the following recurrence formula.

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

When Qn = 1, then pn-1, qn-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

previous index next

E-mail : kc2h-msm@asahi-net.or.jp
Hisanori Mishima