[R] exactly representable numbers

[Duncan Murdoch]

On 9/11/2006 11:01 AM, Robin Hankin wrote:
> Hi Duncan
> 
> 
> [snip]
> 
> On 11 Sep 2006, at 12:12, Duncan Murdoch wrote:
> 
>> Here's my version:  not tested too much.
>>
>> f <- function(x) {
>>    u <- x
>>    l <- 0
>>    mid <- u/2
>>    while (l < mid && mid < u) {
>>      if (x < x + mid) u <- mid
>>      else l <- mid
>>      mid <- (l + u)/2
>>    }
>>    u
>> }
>>
> 
> 
> 
> thanks for this.  Wouldn't it be a good idea to have some function
> that returns "the smallest exactly representable number strictly  
> greater than x"?
> 
> Or does this exist already?

I don't know if it exists.  I wouldn't have much use for it, but if you 
would, maybe it's worth writing.

If you try to do it based on my code above, here are some bugs I've 
noticed since sending that:

  - it doesn't work for negative x
  - it doesn't work for denormal x (e.g. 1.e-310)

I'm not sure what goes wrong in the latter case, but it reports numbers 
which are unnecessarily large:

 > f(1.e-310)
[1] 4.940656e-324
 > 1.e-310 == 1.e-310 + 4e-324
[1] FALSE
 > 1.e-310 == 1.e-310 + 3e-324
[1] TRUE

So the right answer is somewhere between 3e-324 and 4e-324, but my 
function says something bigger.

I suspect the best way to do this accurately is to look at the bit 
patterns of the stored numbers, and add a 1 to the least significant 
bit, but that's a bit too much work for me.

Duncan Murdoch


> 
> 
> 
> Robin Hankin
> Uncertainty Analyst
> National Oceanography Centre, Southampton
> European Way, Southampton SO14 3ZH, UK
>   tel  023-8059-7743