The Four Meanings of "Accurate to 3 Places"

[Joe Horn]

(12-29-2017 10:06 PM)acapde Wrote:  Additionally, I would like to play more with the continued fraction method. It seems that the implementation is not difficult but (being lazy) does anyone happen to have the algorithm already in (Sys) RPL / RPN or any other simple format?

In 1994, Branimir Buljan asked me via email how the "Continued Fraction by Fast Recursion Algorithm" works. Everything below is an excerpt from my reply to him. Please note that this posting is ONLY intended to elucidate the continued-fraction algorithm itself. This is NOT the best way to code it. The code below introduces round-off errors due to executing 1/FP(real), and the algorithm itself misses ALL intermediate convergents. Avoiding these problems takes more code; cf. my PDQ Algorithm posting for details.

-----

Key the following two little programs into your HP48:

SETUP: << (0,1) (1,0) >>

ITER8: << ROT IP LASTARG FP DUP { INV } IFT 4 ROLLD OVER * ROT + >>

To use, type any decimal number, press SETUP once, then press ITER8 (pronounced "iterate") repeatedly, carefully watching the stack. You'll see the "fast recursion formula" in action.

For example, Press PI (with flag -2 set, of course) and then SETUP:

3: 3.14159265359
2: (0,1)
1: (1,0)

Level 1 contains the current approximation in (numerator, denominator) format. As you can see, it's 1/0 at setup (represents infinity).

Press ITER8:

3: 7.06251330592
2: (1,0)
1: (3,1)

This means that 3/1 is the first approximation of pi. The first approximation is always less than the input decimal. The second approximation is always greater than the input decimal:

Press ITER8:

3: 15.9965944095
2: (3,1)
1: (22,7)

So 22/7 is the next approximation of pi, and is greater than pi, as expected, but closer than the previous approximation. Each time you press ITER8, the previous approximation goes up to level 2. This is necessary for the algorithm, but it's also nice to visually compare levels 1 & 2.

The algorithm is easy: take the integer part (IP) of level 3, multiply it by level 1, and add level 2. That's the new approximation. The old approximation goes to level 2, and level 3 gets replaced with 1/FP(x) where x is the former contents of level 3.

So the next iteration, by this definition, would be equal to:

3: 1/0.9965944095
2: (22,7)
1: 15*(22,7)+(3,1)

Press ITER8 and see it:

3: 1.00341722818
2: (22,7)
1: (333,106)

So 333/106 is the next approximation of pi, and is less than pi, as expected, and even closer than the previous approximation. Each press of ITER8 will produce the next approximation, alternating less/greater than pi, and always getting closer to it. Eventually (after 9 iterations in this case) the fraction will be equal to the input as far as the HP 48 is concerned (12 digits max).

The HP 48 ->Q function uses a similar technique, iterating until the difference between the input and the approximation equals zero (rounded to the current display format).

If level 3 ever becomes equal to 0, then stop iterating; the exact fraction has been found. (Running ITER8 further will just perform a SWAP).