Re: [R] exactly representable numbers

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Mon 11 Sep 2006 - 15:14:04 GMT

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

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Tue Sep 12 01:21:17 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Mon 11 Sep 2006 - 16:30:31 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.