%   prog4a.m -- Determine Machine Epsilon
%   Written Apr 25, 2005 by Eric Olson for Mathematics 466
%
%   For more information see Step 3 in First Steps in Numerical
%   Analysis, 2nd Edition, by Hosking, Joe, Joyce and Turner

function prog4a
    data=[ 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0, 1000.0 ];
    disp(sprintf('prog4a.m -- Determine Machine Epsilon\n'));
    for i=1:length(data)
        x=data(i);
        delta=maxdelta(x);
        disp(sprintf('x=%g delta=%.20g',x,delta));
    end

function r=maxdelta(x)
    a=0.0;
    if x>0.0
        b=x;
    else
        b=-x;
    end
    deold=a;
    while 1 
        delta=(a+b)/2.0;
        if deold==delta
            break
        end
        r1=x+delta;
        if r1==x
            a=delta;
        else
            b=delta;
        end
        deold=delta;
    end
    r=a;
