% prog1a.m -- Determine Machine Epsilon
% Written Sep 2, 2007 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 prog1a
    data=[ 1 2 3 4 5 10 100 1000 ];
    disp(sprintf('prog1a.m -- Determine Machine Epsilon\n'));
    for i=1:length(data)
        x=data(i);
        delta=maxdelta(x);
        disp(sprintf('x=%g delta=%g',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;