| home | send comment | send link | add bookmark |

Approximation of irrational numbers

by Stephen R. Schmitt


Definition

A simple continued fraction is an expression like this:
                  1
x = a0 + -------------------
                    1
         a1 + --------------
                       1
              a2 + ---------
                   a3 . . .

The quantities ai are called partial quotients, and the rational number determined by including n terms of the continued fraction is the nth convergent. Every irrational number can be represented as an infinite continued fraction. As n increases, the nth convergent becomes an increasingly accurate rational approximation to an irrational number.

Algorithm

The partial quotients, ai, for a simple continued fraction can be determined using the recursion
r0 = x
rn = 1/(rn-1 - an-1)
where
an = floor( rn )

The convergents are defined as cn = pn/qn. Setting the initial terms pi, qi to

p-2 = 0, q-2 = 1
p-1 = 1, q-1 = 0
p0 = a0, q0 = 1
allows us to compute subsequent terms from the recursion
pn = anpn-1 + pn-2
qn = anqn-1 + qn-2

Zeno source code

Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.

const precision : real := 0.0000000001

program

    var x : real

    put "number to be approximated"...
    get x

    fraction( x )
    
end program

procedure fraction( x : real )
% compute a sequence of partial quotients and convergents

    var a, p, p1, p2, q, q1, q2 : int
    var r : real

    % initialize
    r  := x
    p2 := 0
    q2 := 1
    p1 := 1
    q1 := 0
    
    while true

        % compute partial quotient and convergent
        a := floor( r )
        p := a*p1 + p2
        q := a*q1 + q2

        put (p/q):20:14, " = ", p, "/", q
        exit when abs(p/q - x) < precision

        % advance
        r := 1/(r - a)
        p2 := p1
        q2 := q1
        p1 := p
        q1 := q

    end while

end procedure

function abs( f : real ) : real
% absolute value of real number

    if f < 0.0 then
        f := -f
    end if

    return f

end function

Sample output

number to be approximated? 3.1415926535897932384626433832795
    3.00000000000000 = 3/1
    3.14285714285714 = 22/7
    3.14150943396226 = 333/106
    3.14159292035398 = 355/113
    3.14159265301190 = 103993/33102
    3.14159265392142 = 104348/33215
    3.14159265346744 = 208341/66317
    3.14159265361894 = 312689/99532

AbCd Classics - free on-line books


Copyright © 2005, Stephen R. Schmitt