Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Get distance to next largest floating point number in R [duplicate]

Is there any implementation of functionality in R, such that it is possible to get the next representable floating point number from a given floating point number. This would be similar to the nextafter function in the C standard library. Schemes such as number + .Machine$double.eps don't work in general.

like image 624
Rónán Daly Avatar asked Nov 20 '22 16:11

Rónán Daly


1 Answers

No, but there are two ways you can make it:

Using C

If you want the exact functionality of the nextafter() function, you can write a C function that works as an interface to the function such that the following two constraints are met:

  • The function does not return a value. All work is accomplished as a "side effect" (changing the values of arguments).
  • All the arguments are pointers. Even scalars are vectors (of length one) in R.

That function should then be compiled as a shared library:

R CMD SHLIB foo.c

for UNIX-like OSs. The shared library can be called using dyn.load("foo.so"). You can then call the function from inside R using the .C() function

.C("foo", ...)

A more in depth treatment of calling C from R is here.

Using R

number + .Machine$double.eps is the way to go but you have to consider edge cases, such as if x - y < .Machine$double.eps or if x == y. I would write the function like this:

nextafter <- function(x, y){
  # Appropriate type checking and bounds checking goes here
  delta = y - x
  if(x > 0){
    factor = 2^floor(log2(x)) + ifelse(x >= 4, 1, 0)
      } else if (x < 0) {
    factor = 65
  }
  if (delta > .Machine$double.eps){
    return(x + factor * .Machine$double.eps)
  } else if (delta < .Machine$double.eps){
    return(x - factor * .Machine$double.eps)
  } else {
    return(x)
  }
}

Now, unlike C, if you want to check integers, you can do so in the same function but you need to change the increment based on the type.

UPDATE The previous code did not perform as expected for numbers larger than 2. There is a factor that needs to be multiplied by the .Machine$double.eps to make it large enough to cause the numbers to be different. It is related to the nearest power of 2 plus one. You can get an idea of how this works with the below code:

n <- -100
factor <- vector('numeric', 100)
for(i in 1:n){
  j = 0
  while(TRUE){
    j = j + 1
    if(i - j * .Machine$double.eps != i) break()
  }
  factor[i] = j
}  
like image 162
Christopher Louden Avatar answered Nov 22 '22 06:11

Christopher Louden