# macheps.r                         J F Monahan, June 2001
 #                                   revised August 2007, 2009
 # compute machine epsilon
 #   *** these results on R version 2.9.0 on Linux       ***
 #   *** Vista version 2.9.1 gives old 2**(-53)+2**(-63) ***
 #
 # first try some values to get in ballpark
  eps <- 2**(-56)
  "2**(-56)"
  eps
  "Is 1+eps bigger than 1 for eps = j*2**(-56), j = 1 to 16"
  (1!=(1+eps*(1:16)))      # F if different, T if not 
 #
 # look at representations with sprintf and "%a"
  sprintf("%a",(1+eps*(1:16)))
 #
 # 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 than 1 for eps = j*2**(-59), j = 60 to 70"
  eps <- 2**(-59)
  (1!=(1+eps*(60:70)))
 #
 # look at representations 
  sprintf("%a",(1+eps*(60:70)))
 #
 # another try
  eps <- 2**(-53)          # first guess + smaller
  (1!=(1+eps*(1 + 2**(-51:-55))))
 #
 # look at representations
  sprintf("%a",(1+eps*(1 + 2**(-51:-55))))
 #
 # this is different -- the real machine epsilon
  "Is 1+eps bigger than 1 for eps = 2**(-53) + 2**(-105)"
  eps <- 2**(-53) + 2**(-105)
  eps
 # this should do it
  sprintf("%a",eps)
  (1!=(1+eps))
 #
 # this is the same -- just a hair too small
  "Is 1+eps bigger than 1 for eps = 2**(-53) + 2**(-106)"
  eps <- 2**(-53) + 2**(-106)
  eps
 # this should be too much 
  sprintf("%a",eps)
 # test it
  (1!=(1+eps))
 #
  rm(list=ls())
  q()