Convergents of a Continued Fraction

[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.