1.4 Machine Epsilon

The minimum value (or one close to it) which can be differentiated from one plus it.

1 + 1e-100 == 1 # TRUE

## [1] TRUE

1 + .Machine$double.eps == 1

## [1] FALSE

1 + .Machine$double.eps / 2 == 1

## [1] TRUE

To get machine epsilon.

MachEps

## function () 
## {
##     eps = 1
##     while (1 + eps > 1) eps = eps/2
##     return(2 * eps)
## }
## <bytecode: 0x000001e2ab941b48>
## <environment: namespace:math>

MachEps()

## [1] 2.22e-16

.Machine$double.eps

## [1] 2.22e-16

A little more generalized epsilon

MachEps2

## function (x = 1, Start = 1) 
## {
##     eps = Start
##     while (x + eps > x) eps = eps/2
##     return(2 * eps)
## }
## <bytecode: 0x000001e2abdc7358>
## <environment: namespace:math>

MachEps2(100, 1) # Different from .Machine$doubple.eps. Larger x, larger epsilon.

## [1] 1.421e-14

(eps2 = MachEps2(1, 0.8))

## [1] 1.776e-16

1 + eps2 == 1 # FALSE and less than .Machine$double.eps

## [1] FALSE

Get really the smallest epsilon by bisectional search

MachEps3

## function (MaxIter = 1000) 
## {
##     lb = 0
##     ub = 1
##     NextEps = (lb + ub)/2
##     i = 0
##     while (NextEps != lb & NextEps != ub & i < MaxIter) {
##         i = i + 1
##         CurEps = NextEps
##         if (1 + CurEps > 1) {
##             ub = CurEps
##         }
##         else {
##             lb = CurEps
##         }
##         NextEps = (lb + ub)/2
##     }
##     attr(CurEps, "iteration") = i
##     return(CurEps)
## }
## <bytecode: 0x000001e2ac25b610>
## <environment: namespace:math>

(eps3 = MachEps3()) # smaller than .Machine$double.eps

## [1] 1.11e-16
## attr(,"iteration")
## [1] 105

1 + eps3 == 1 # FALSE

## [1] FALSE

format(eps3, digits=22)

## [1] "1.110223024625156786971e-16"

format(.Machine$double.eps, digits=22)

## [1] "2.220446049250313080847e-16"

How about the smallest epsilon for 1 - eps < 1

MachEps4

## function (MaxIter = 1000) 
## {
##     lb = -1
##     ub = 0
##     NextEps = (lb + ub)/2
##     i = 0
##     while (NextEps != lb & NextEps != ub & i < MaxIter) {
##         i = i + 1
##         CurEps = NextEps
##         if (1 + CurEps < 1) {
##             lb = CurEps
##         }
##         else {
##             ub = CurEps
##         }
##         NextEps = (lb + ub)/2
##     }
##     attr(CurEps, "iteration") = i
##     return(CurEps)
## }
## <bytecode: 0x000001e2ac6bf7a0>
## <environment: namespace:math>

(eps4 = MachEps4()) # different not only in sign, but also in magnitude!

## [1] -5.551e-17
## attr(,"iteration")
## [1] 106

1 + eps4 == 1

## [1] FALSE

How about the smallest positive value?

NextEps = 1
while(NextEps > 0) {
  eps = NextEps
  NextEps = NextEps/2
}
format(eps, digits=22) # much less than .Machine$double.xmin !!!

## [1] "4.940656458412465440715e-324"

format(.Machine$double.xmin, digits=22) # R specification of smaleest positive value

## [1] "2.22507385850720138280e-308"

format(Bin2Dec(c(rep(0, 63), 1)), digits=22) # IEEE 754 smallest positive value

## [1] "1.112536929253600691400e-308"