 The Museum of HP Calculators

HP Articles Forum

The continued fractions algorithm (was: Re: Converting a decimal # to a fraction #)

Posted by rupert9000 on 22 Apr 2000, 6:09 p.m.

Well, for my current practical purposes I developed the simple algorithm previously posted, but now I have seized the opportunity to re-implement a more efficient algorithm based on the continued fractions. It is well explained, for practical purposes, on the HP67 Mat I Pac manual (afterall this is the HP Museum...) :

for any given number, let's say decnum, the continued fractions algorithm can generate a converging serie of fractional approximations of it. Each fractional approximation is computed as a numerator n(i) and a denominator d(i) as follows:

n(i) = a(i) * n(i-1) + n(i-2) d(i) = a(i) * d(i-1) + d(i-2)

The following values are assigned by definition: n(-1) = 0, d(-1) = 1, n(0) = 1, d(0) = 0.

The values for a(i) coefficients are computed as follows:

a(i+1) = int( x(i) / y(i) ), where x(i+1) = y(i) and y(i+1) = x(i) - a(i+1) * y(i).

Again, the following values are assigned by definition: a(1) = int( decnum ), x(1) = 1, y(1) = frac( decnum ).

For a reference on the web look at http://archives.math.utk.edu/articles/atuyl/confrac/.

Here is a listing (on this site it might not be properly displayed), written in plain BASIC because it allows to easily test an algorithm, is self-explaining and easily adaptable to any kind of calculator/computer.

---------------------------------------

REM Continued fractions algorithm

DIM a1, a2, x1, y1, x2, y2 AS DOUBLE

DIM dm1, d0, d1, nm1, n0, n1 AS DOUBLE

DIM decnum, delta, f AS DOUBLE

DIM i AS LONG

CLS

PRINT ""

INPUT "Number "; decnum

nm1 = 0

dm1 = 1

n0 = 1

d0 = 0

a1 = int(decnum)

x1 = 1

y1 = decnum - a1 : REM y1 = frac(decnum)

i = 1

DO

a2 = INT(x1 / y1)

x2 = y1

y2 = x1 - a2 * y1

n1 = a1 * n0 + nm1

d1 = a1 * d0 + dm1

f = n1 / d1

delta = ABS(decnum - f)

PRINT "i="; i; " N(i)="; n1; " D(i)="; d1; " N(i)/D(i)="; f; " Delta="; delta

PRINT " press any key..."

DO

LOOP WHILE INKEY\$ = ""

a1 = a2

x1 = x2

y1 = y2

nm1 = n0

dm1 = d0

n0 = n1

d0 = d1

i = i + 1

LOOP WHILE y1 > 0

PRINT ""

PRINT "END PRG"

SYSTEM

---------------------------------------

As I previously wrote, this algorithm is more efficient than those based on trial and error. Moreover, it doesn't generate results like 0.3695=3,695/10,000 as some other algorithms do.

In order to make a practical example I can describe the typical use I do of such kind of calculation while drawing and dimensioning mechanical parts with a CAD. Let's say a threading diameter of 87.3 mm = 3.437007874 inches must be expressed in inches plus a fraction of an inch. Entering 0.437007874 as decnum the program generates the following serie of fractional approximations:

i= 1 N(i)= 0 D(i)= 1 N(i)/D(i)= 0 Delta= .4370079 i= 2 N(i)= 1 D(i)= 2 N(i)/D(i)= .5 Delta= 6.299213E-02 i= 3 N(i)= 3 D(i)= 7 N(i)/D(i)= .4285714285714285 Delta= 8.436446E-03 i= 4 N(i)= 7 D(i)= 16 N(i)/D(i)= .4375 Delta= 4.921257E-04 i= 5 N(i)= 52 D(i)= 119 N(i)/D(i)= .4369747899159664 Delta= 3.308434E-05 i= 6 N(i)= 59 D(i)= 135 N(i)/D(i)= .4370370370370371 Delta= 2.916279E-05 i= 7 N(i)= 111 D(i)= 254 N(i)/D(i)= .437007874015748 Delta= 2.34664E-10 i= 8 N(i)= 7331720 D(i)= 1.677709E+07 N(i)/D(i)= .4370078742504138 Delta= 1.776357E-15

The threading diameter will then be dimensioned as 3" 7/16.

Someone might prefer to enter a tolerance and display just the single fractional approximation that satisfies the condition abs( n(i)/d(i) - decnum ) 0? 78 GTO`ITERAT` 79`END PRG` 80 AVIEW 81 BEEP 82 RTN

Reg. Use

00 a1 01 a2 02 x1 03 y1 04 i 05 y2 06 dm1 07 d0 08 d1 09 nm1 10 n0 11 n1 12 decnum

Well, I hope my English is sufficiently intelligible and that what I tried to explain above can be found useful.