HP Articles Forum
[Return to the Index ]
[ Previous | Next ]
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.