Convergents of a Continued Fraction

[Thomas Klemm]

Convergents of a Continued Fraction

References

• Continued fraction
• (48SX and 12C): Continued Fractions
• (50G) Fast Continued Fractions
  

Goal

We want to calculate the convergents of a continued fraction.

Formula

For example, here are the convergents for \([0;1,5,2,2]\).

\(
\begin{matrix}
n & −2 & −1 & 0 & 1 & 2 & 3 & 4 \\
a_n & & & 0 & 1 & 5 & 2 & 2 \\
h_n & 0 & 1 & 0 & 1 & 5 & 11 & 27 \\
k_n & 1 & 0 & 1 & 1 & 6 & 13 & 32 \\
\end{matrix}
\)

If successive convergents are found, with numerators \(h_1\), \(h_2\), ... and denominators \(k_1\), \(k_2\), ... then the relevant recursive relation is:

\(
\begin{align}
h_n &= a_n h_{n − 1} + h_{n − 2} \\
k_n &= a_n k_{n − 1} + k_{n − 2} \\
\end{align}
\)

The successive convergents are given by the expression:

\(
\begin{align}
\frac{h_n}{k_n}
\end{align}
\)

Using this formula, we would need to track four values on the stack.
However, I couldn't figure out how to do that without using registers.

But we can divide all the equations by \(k_{n − 1}\) and get:

\(
\begin{align}
\frac{h_n}{k_{n − 1}} &= a_n \frac{h_{n − 1}}{k_{n − 1}} + \frac{h_{n − 2}}{k_{n − 1}} \\
\\
\frac{k_n}{k_{n − 1}} &= a_n + \frac{k_{n − 2}}{k_{n − 1}} \\
\end{align}
\)

Now we introduce three variables \(u\), \(v\) and \(w\):

\(
\begin{align}
u &= \frac{h_{n − 1}}{k_{n − 2}} \\
\\
v &= \frac{k_{n − 1}}{k_{n − 2}} \\
\\
w &= \frac{h_{n − 2}}{k_{n − 2}} \\
\end{align}
\)

With \(u{'}\), \(v{'}\) and \(w{'}\) we denote the successors:

\(
\begin{align}
u{'} &= \frac{h_{n}}{k_{n − 1}} =& a_n w{'} + \frac{w}{v} \\
\\
v{'} &= \frac{k_{n}}{k_{n − 1}} =& a_n + \frac{1}{v} \\
\\
w{'} &= \frac{h_{n − 1}}{k_{n − 1}} =& \frac{u}{v} \\
\end{align}
\)

So we only need to keep track of the three values \(v\), \(w\), and \(w{'}\) on the stack.
The sequence \(w, w{'}, w{''}, \cdots\) are the convergents of the continued fraction.

Programs

This program is for the HP-42S, but it will work on most HP calculators with little or no modification:

Code:

00 { 18-Byte Prgm }
01 X<>Y
02 R↑
03 LASTX
04 ÷
05 LASTX
06 1/X
07 R↑
08▸LBL 00
09 +
10 X<>Y
11 LASTX
12 R↑
13 ×
14 LASTX
15 R↓
16 +
17 X<>Y
18 ÷
19 END

This the same program for the HP-15C:

Code:

   001 {          34 } x↔y
   002 {       43 33 } g R⬆
   003 {       43 36 } g LSTx
   004 {          10 } ÷
   005 {       43 36 } g LSTx
   006 {          15 } 1/x
   007 {       43 33 } g R⬆
   008 {    42 21  0 } f LBL 0
   009 {          40 } +
   010 {          34 } x↔y
   011 {       43 36 } g LSTx
   012 {       43 33 } g R⬆
   013 {          20 } ×
   014 {       43 36 } g LSTx
   015 {          33 } R⬇
   016 {          40 } +
   017 {          34 } x↔y
   018 {          10 } ÷
   019 {       43 32 } g RTN

Stack Diagram

We enter the program at label 00.
At that point you can consider \(\frac{w}{v} = 1\) and \(\frac{1}{v} = 0\).

Code:

| L  | T  | Z   | Y    | X    | 
+----+----+-----+------+------+ 
| v  | w  | w   | w'   | a    |  X<>Y
| v  | w  | w   | a    | w'   |  R↑
| v  | w  | a   | w'   | w    |  LASTX
| v  | a  | w'  | w    | v    |  ÷
| v  | a  | a   | w'   | w/v  |  LASTX
|    | a  | w'  | w/v  | v    |  1/X
|    | a  | w'  | w/v  | 1/v  |  R↑
|    | w' | w/v | 1/v  | a    |  LBL 00
|    | w' | w/v | 1/v  | a    |  +
| a  | w' | w'  | w/v  | v'   |  X<>Y
| a  | w' | w'  | v'   | w/v  |  LASTX
|    | w' | v'  | w/v  | a    |  R↑
|    | v' | w/v | a    | w'   |  ×
| w' | v' | v'  | w/v  | a*w' |  LASTX
|    | v' | w/v | a*w' | w'   |  R↓
|    | w' | v'  | w/v  | a*w' |  +
|    | w' | w'  | v'   | u'   |  X<>Y
|    | w' | w'  | u'   | v'   |  ÷
| v' | w' | w'  | w'   | w"   |

Usage

We start with the first two elements of the continued fraction and place them on the stack.
However, we must enter the values \(1\) and \(0\) in between:

\(a_0\) ENTER
1 ENTER
0 ENTER
\(a_1\) ENTER
XEQ 00

In the following examples we write only briefly:

\(a_0\) 1 0 \(a_1\) XEQ 00

Now we can proceed with:

\(a_2\) R/S
\(a_3\) R/S
\(a_4\) R/S
\(a_5\) R/S
…


Examples

Rational Numbers

This is the finite continued fraction example from above and represents a rational number.

\(
[0; 1, 5, 2, 2]
\)

0 1 0 1 XEQ 00
5 R/S
2 R/S
2 R/S

Code:

    0        0.00000000000                     0                     1
    1        1.00000000000                     1                     1
    5        0.83333333333                     5                     6
    2        0.84615384615                    11                    13
    2        0.84375000000                    27                    32

Golden Ratio

The golden ratio \(\varphi\) has terms equal to \(1\) everywhere—the smallest values possible—which makes \(\varphi\) the most difficult number to approximate rationally.
In this sense, therefore, it is the "most irrational" of all irrational numbers.

\(
\varphi = [\overline{1}]
\)

1 1 0 1 XEQ 00
1 R/S
1 R/S
1 R/S
1 R/S
1 R/S
1 R/S
1 R/S
1 R/S
…

Code:

    1        1.00000000000                     1                     1
    1        2.00000000000                     2                     1
    1        1.50000000000                     3                     2
    1        1.66666666667                     5                     3
    1        1.60000000000                     8                     5
    1        1.62500000000                    13                     8
    1        1.61538461538                    21                    13
    1        1.61904761905                    34                    21
    1        1.61764705882                    55                    34
    1        1.61818181818                    89                    55
    1        1.61797752809                   144                    89
    1        1.61805555556                   233                   144
    1        1.61802575107                   377                   233
    1        1.61803713528                   610                   377
    1        1.61803278689                   987                   610
    1        1.61803444782                  1597                   987
    1        1.61803381340                  2584                  1597
    1        1.61803405573                  4181                  2584
    1        1.61803396317                  6765                  4181
    1        1.61803399852                 10946                  67651

Also it produces the Fibonacci sequence.

Periodic Continued Fractions

The numbers with periodic continued fraction expansion are precisely the irrational solutions of quadratic equations with rational coefficients.

\(
\sqrt{42} = [6; \overline{2, 12}]
\)

6 1 0 2 XEQ 00
12 R/S
2 R/S
12 R/S
2 R/S
12 R/S
2 R/S
12 R/S
2 R/S
12 R/S
2 R/S
…

Code:

    6        6.00000000000                     6                     1
    2        6.50000000000                    13                     2
   12        6.48000000000                   162                    25
    2        6.48076923077                   337                    52
   12        6.48073959938                  4206                   649
    2        6.48074074074                  8749                  1350
   12        6.48074069678                109194                 16849
    2        6.48074069847                227137                 35048
   12        6.48074069841               2834838                437425
    2        6.48074069841               5896813                909898
   12        6.48074069841              73596594              11356201
    2        6.48074069841             153090001              23622300

Other Patterns

While there is no discernible pattern in the simple continued fraction expansion of \(\pi\), there is one for \(e\), the base of the natural logarithm:

\(
e = [2; 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, 1, 1, 10, 1, 1, 12, 1, 1, 14, 1, 1, 16, 1, 1, 18, 1, 1, \cdots]
\)

2 1 0 1 XEQ 00
2 R/S
1 R/S
1 R/S
4 R/S
1 R/S
1 R/S
6 R/S
1 R/S
1 R/S
8 R/S
1 R/S
…

Code:

    2        2.00000000000                     2                     1
    1        3.00000000000                     3                     1
    2        2.66666666667                     8                     3
    1        2.75000000000                    11                     4
    1        2.71428571429                    19                     7
    4        2.71875000000                    87                    32
    1        2.71794871795                   106                    39
    1        2.71830985915                   193                    71
    6        2.71827956989                  1264                   465
    1        2.71828358209                  1457                   536
    1        2.71828171828                  2721                  1001
    8        2.71828183521                 23225                  8544
    1        2.71828182294                 25946                  9545
    1        2.71828182874                 49171                 18089
   10        2.71828182845                517656                190435
    1        2.71828182847                566827                208524
    1        2.71828182846               1084483                398959
   12        2.71828182846              13580623               4996032

\(\pi\)

No pattern has ever been found in this representation.

\(
\pi = [3; 7,15,1,292,1,1,1,2,1,3,1,14,2,1,1,2,2,2,2,1,84,2,1,1,15,3,13, \cdots]
\)

3 1 0 7 XEQ 00
15 R/S
1 R/S
292 R/S
1 R/S
1 R/S
1 R/S
2 R/S
1 R/S
3 R/S
1 R/S
14 R/S
2 R/S
…

Code:

    3        3.00000000000                     3                     1
    7        3.14285714286                    22                     7
   15        3.14150943396                   333                   106
    1        3.14159292035                   355                   113
  292        3.14159265301                103993                 33102
    1        3.14159265392                104348                 33215
    1        3.14159265347                208341                 66317
    1        3.14159265362                312689                 99532
    2        3.14159265358                833719                265381
    1        3.14159265359               1146408                364913
    3        3.14159265359               4272943               1360120
    1        3.14159265359               5419351               1725033
   14        3.14159265359              80143857              25510582

[Albert Chan]

Code work even if continued fraction coefficient is 0. Very nice!

OP HP-42S code:

1 1 0 2 XEQ 00      --> 1 + 1/2 = 1.5
0 R/S                    --> 1 + 1/(2 + 1/0) = 1 + 1/inf = 1
3 R/S                    --> 1 + 1/(2 + 1/(0 + 1/3)) = 1 + 1/(2 + 3) = 1.2

Here is lua code to allow for Generalized Continued Fraction



We simply scale CF numerator to 1, then use OP simple continued fraction code.

x = b0 + a1/(b1 + a2/(b2 + a3/(b3 + ...
   = b0 + 1/(b1/a1 + (a2/a1)/(b2 + a3/(b3 + ...
   = b0 + 1/(b1/a1 + 1/(b2/(a2/a1) + a3/(a2/a1)/(b3 + ...

Code:

function dfc2f(b0)
    local k,u,v,w,w2 = 1,1,0,0
    return function(a, b)
        if not b then a, b = 1, a end
        k = a / k
        b = b / k   -- generalized CF to simple CF
        if not w2 then
            v = b; w2 = 1/b
        else
            u = b*w2 + w/v
            v = b + 1/v
            w, w2 = w2, u/v
        end
        return b0 + w2
    end
end

Lets try tan(pi/4)=1, using generalized CF form, and simple CF form (both are equivalent)

lua> CF1, CF2 = dfc2f(0), dfc2f(0)
lua> x = pi/4
lua> CF1(x,1), CF2(1/x)
0.7853981633974483      0.7853981633974483

lua> for k=1,8 do n=2*k+1; print(CF1(-x*x, n), CF2((-1)^k*n/x)) end
0.988689239934205        0.988689239934205
0.9997876809149683      0.9997876809149683
0.9999978684156949      0.9999978684156949
0.9999999865263551      0.9999999865263551
0.9999999999413254      0.9999999999413254
0.9999999999998131      0.9999999999998131
0.9999999999999997      0.9999999999999997
1.0000000000000002      1.0000000000000002

[Albert Chan]

(03-13-2022 11:41 PM)Thomas Klemm Wrote:  With \(u{'}\), \(v{'}\) and \(w{'}\) we denote the successors:

\(
\begin{align}
u{'} &= \frac{h_{n}}{k_{n − 1}} =& a_n w{'} + \frac{w}{v} \\
\\
v{'} &= \frac{k_{n}}{k_{n − 1}} =& a_n + \frac{1}{v} \\
\\
w{'} &= \frac{h_{n − 1}}{k_{n − 1}} =& \frac{u}{v} \\
\end{align}
\)

So we only need to keep track of the three values \(v\), \(w\), and \(w{'}\) on the stack.
The sequence \(w, w{'}, w{''}, \cdots\) are the convergents of the continued fraction.

Assuming {u,v,w} start out finite. (i.e. a₁ ≠ 0)

This is what happen if we have double-zero CF coeffs, a_n = a_n+1 = 0, n>1

w' = u/v
u' = w/v
v' = 1/v

w'' = u'/v' = (w/v) / (1/v) = w
u'' = w'/v' = (u/v) / (1/v) = u
v'' = 1/v' = 1/(1/v) = v

0 + 1/(0 + 1/CF) = 1/(1/CF) = CF

{u'',v'',w''} = {u,v,w}, matching expectation.

[Albert Chan]

(03-13-2022 11:41 PM)Thomas Klemm Wrote:  For example, here are the convergents for \([0;1,5,2,2]\).

\(
\begin{matrix}
n & −2 & −1 & 0 & 1 & 2 & 3 & 4 \\
a_n & & & 0 & 1 & 5 & 2 & 2 \\
h_n & 0 & 1 & 0 & 1 & 5 & 11 & 27 \\
k_n & 1 & 0 & 1 & 1 & 6 & 13 & 32 \\
\end{matrix}
\)

Just a visual explanation why OP setup only need to track 3 values on the stack.

Scaling the entries does not change convergent ratios.
h(n)/k(n) = w(n)/1, and we don't need to store denominator 1

Code:

a(n)      0 1    5    2     2
h(n) 0  1 0 1    5
k(n) 1  0 1 1    6

            1/6  5/6  11/6   11/13
            ---  ---  ---- = -----
            1/6  1    13/6     1
       
                 5/13 11/13  27/13   27/32
                 ---- -----  ----- = -----
                 6/13   1    32/13     1

But, there is a cost. Convergent integer ratio is lost, only its value remained.
Also, a1≠0 is required, because we can't scale denominator of 0 to 1.
Also, we can't start with a0 alone, because k(-1)=0, unable to scale to 1

[Albert Chan]

(03-13-2022 11:41 PM)Thomas Klemm Wrote:  We enter the program at label 00.
At that point you can consider \(\frac{w}{v} = 1\) and \(\frac{1}{v} = 0\).

We start with the first two elements of the continued fraction and place them on the stack.
However, we must enter the values \(1\) and \(0\) in between:

\(a_0\) ENTER
1 ENTER
0 ENTER
\(a_1\) ENTER
XEQ 00

I was curious if we enter some other value pairs between a0, a1

Code:

            a0     a1          a2
h(n)    z   a0  (a0*a1+z)  (a0*a1+z)*a2 + a0
k(n)    y   1   (a1+y)     (a1+y)*a2 + 1

\(\displaystyle \frac{h_1}{k_1}
= \frac{a_0\;a_1 + z}{a_1 + y}
\;=\; a_0 + \frac{z - a_0\;y}{a_1 + y} \)

\(\displaystyle \frac{h_2}{k_2}
= \frac{(a_0\;a_1+z)×a_2 + a_0}{(a_1+y)×a_2+1}
\;=\; a_0 + \frac{z - a_0\;y}{a_1 + y + \frac{1}{a_2}} \)
...

\(\displaystyle \frac{h_∞}{k_∞}
= a_0 + \frac{z - a_0\;y}{[a_1 + y ;\; a_2 , a_3 ,\; ...]}\)

If {z,y} = {1,0}, this reduced to simple continued fraction \([a_0 ;\; a_1, a_2, a_3, \; ...]\)

Confirm by example:

3 Enter 14 Enter 159 Enter 2654
XEQ 00      → 2.835407038748666903661571276217561
1 R/S         → 2.835465529495380241648898365316275
2 R/S         → 2.835446037199383959246534770761757

lua> a0, z, y, a1 = 3, 14, 159, 2654
lua> f = dfc2f(a0)
lua> f(z-a0*y, a1+y)
2.835407038748667
lua> f(1)
2.83546552949538
lua> f(2)
2.835446037199384

[C.Ret]

(03-13-2022 11:41 PM)Thomas Klemm Wrote:  For example, here are the convergents for \([0;1,5,2,2]\).

\(
\begin{matrix}
n & −2 & −1 & 0 & 1 & 2 & 3 & 4 \\
a_n & & & 0 & 1 & 5 & 2 & 2 \\
h_n & 0 & 1 & 0 & 1 & 5 & 11 & 27 \\
k_n & 1 & 0 & 1 & 1 & 6 & 13 & 32 \\
\end{matrix}
\)

If successive convergents are found, with numerators \(h_1\), \(h_2\), ... and denominators \(k_1\), \(k_2\), ... then the relevant recursive relation is:

\(
\begin{align}
h_n &= a_n h_{n − 1} + h_{n − 2} \\
k_n &= a_n k_{n − 1} + k_{n − 2} \\
\end{align}
\)

The successive convergents are given by the expression:

\(
\begin{align}
\frac{h_n}{k_n}
\end{align}
\)

Using this formula, we would need to track four values on the stack.

Thank you for this very informative post and for sharing such a concise and elegant code.

I also don't know how to manage four elements on the stack. However, since there are only two numerators and two denominators, we only need to use two stack elements for each new \(a_i\). The four variables come from the two numerator and denominator pairs. If we can use just one stack level for each element, we’re safe.

One easy way to achieve this on an HP-15C is to use the complex-stack. I’m not sure if the same trick can be applied to other HP calculators that support complex arithmetic.

Here is a simple program for the HP-15C that uses the complex stack to compute the convergents:

Code:

001 { 42 21 11 } f LBL A    ;  Add next element a(i) to the set [ a₀ ; a₁ , a₂ , …  ] 
002 {       34 }   X↔Y
003 {       36 }   ENTER          ;; make a copy of previous ('h,'k)
004 {       33 }   R↓
005 {       20 }    *
006 {       40 }    +              ;; compute (h,k) = ("h,"k) + ('h,'k)*a(i)
007 {    44  0 }   STO 0
008 {    42 30 } f Re↔Im
009 { 44 10  0 }   STO / 0        ;; store h/k in register R0
010 {    42 30 } f Re↔Im

Code:

011 { 42 21 15 } f LBL E    ; Show Estimated value of convergent h/k
012 {    45  0 }   RCL 0
013 {    42 31 } f PSE
014 {       33 }   R↓
015 {    43 32 } g RTN

Code:

016 { 42 21 13 } f LBL C    ; CLEAR / redo from start
017 { 43  4  8 } g SF 8
018 {        1 }    1
019 {       36 }   ENTER
020 {    42 30 } f Re↔Im
021 {       34 }   X↔Y
022 {    43 32 } g RTN

LBL C: Clears the stack and restarts from the beginning for a new continued fraction.
LBL A or R/S: Adds the next integer \(a_i\) and compute a new convergent (h,k) where h is the real part and k the imaginary part of X: register. Also h/k is store in register R0.
LBL E: Briefly displays an estimate of the current convergent (which is stored in register R0).

At each step, f-(i) can be used to view the denominator; h (Re-part) and k (Im-part). The register R0 contains estimated value of convergent h/k.

However, make sure to leave (h−2,k−2) and (h−1,k−1) in registers Y: and X: respectively before pressing f−A or R/S.

Code:

Example: [ 0 ; 1 , 5 , 2 , 2 ]

ON
f-C                        1.00000  f-(i)   0.00000   
0 f-A    // 0.00000 //     0.00000  f-(i)   1.00000
1 R/S    // 1.00000 //     1.00000  f-(i)   1.00000  
5 R/S    // 0.83333 //     5.00000  f-(i)   6.00000
2 R/S    // 0.84615 //    11.00000  f-(i)  13.00000
2 R/S    // 0.84375 //    27.00000  f-(i)  32.00000
f-E      // 0.84375 //    27.00000
RCL 0                      0.84375

Code:

Example: Euler's Constant γ [ 0; 1, 1, 2, 1, 2, 1, 4, 3, 13, 5, 1, 1, 8, 1, 2, 4, 1, 1, 40, . . . 

f-C                         1.00000  f-(i)      0.00000   
0 R/S    // 0.00000 //      0.00000  f-(i)      1.00000
1 R/S    // 1.00000 //      1.00000  f-(i)      1.00000
1 R/S    // 0.50000 //      1.00000  f-(i)      2.00000
2 R/S    // 0.60000 //      3.00000  f-(i)      5.00000
1 R/S    // 0.57143 //      4.00000  f-(i)      7.00000
2 R/S    // 0.57895 //     11.00000  f-(i)     19.00000
1 R/S    // 0.57692 //     15.00000  f-(i)     26.00000
4 R/S    // 0.57724 //     71.00000  f-(i)    123.00000
3 R/S    // 0.57722 //    228.00000  f-(i)    395.00000
13 R/S   // 0.57722 //  3,035.00000  f-(i)  5,258.00000
5 R/S    // 0.57722 // 15,403.00000  f-(i) 26,685.00000
. . .

I didn’t quite reach my goal, as my code is a bit longer and uses one register (I couldn’t figure out how to compute Re/Im from the stack without a register), whereas your brilliant code is shorter and only uses the stack.

That said, my approach is somewhat simpler for the user, as it only requires entering the continued fraction elements in order and nothing else.

[Thomas Klemm]

(09-20-2024 10:24 PM)C.Ret Wrote:  One easy way to achieve this on an HP-15C is to use the complex-stack. I’m not sure if the same trick can be applied to other HP calculators that support complex arithmetic.

This program is for the HP-42S:

Code:

00 { 27-Byte Prgm }
01▸LBL A
02 RCL× ST Y
03 RCL+ ST Z
04▸LBL E
05 ENTER
06 COMPLEX
07 ÷
08 PSE
09 R↓
10 RTN
11▸LBL C
12 0
13 1
14 COMPLEX
15 1
16 END

To use it set LCLBL and select the CUSTOM menu.