# macheps.s                         J F Monahan, June 2001
 # compute machine epsilon
 #
# first try some values to get in ballpark
eps <- 2**(-56)
 "2**(-56)"
eps
 "Is 1+eps bigger 1 for eps = j*2**(-56), j = 1 to 16"
d <- (1!=(1+eps*(1:16)))      # F if different, T if not 
d
# at this point it's around 8*2**(-56) or 2**(-53)
#  or ... 64*2**(-59), so try 60 to 70
 "Is 1+eps bigger 1 for eps = j*2**(-59), j = 60 to 70"
eps <- 2**(-59)
d <- (1!=(1+eps*(60:70)))
d
# this is different -- the real machine epsilon
 "Is 1+eps bigger 1 for eps = 2**(-53) + 2**(-105)"
eps <- 2**(-53) + 2**(-105)
eps
d <- (1!=(1+eps))
d
# this is the same -- just a hair too big
 "Is 1+eps bigger 1 for eps = 2**(-53) + 2**(-106)"
eps <- 2**(-53) + 2**(-106)
d <- (1!=(1+eps))
d
 rm(list=ls())