% Script: numbers
% find machine epsilon
x = 1; p = 0; y = 1; z = x + y;
while x ~= z, y = y/2; p = p+1; z = x+y; end
disp(sprintf('Machine epsilon = %22.15e and exponent = %1.0f',...
   2*y,p))
% find underflow exponent
x = 1; q = 0;
while x > 0, x = x/2; q = q+1; z = x+y; end
disp(sprintf('underflow expt: base 2 = %1.0f,base 10 = %1.0f',...
   -q/2,-q/2/log(10)))
x = 1; r = 0;
while x ~= inf, x = x*2; r = r+1; z = x+y; end
disp(sprintf('overflow expt: base 2 = %1.0f, base 10 = %1.0f',...
   r/2,r/2/log(10)))  