Re: [R] exactly representable numbers

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Mon 11 Sep 2006 - 13:40:28 GMT

On 9/11/2006 9:01 AM, Prof Brian Ripley wrote:

> On Mon, 11 Sep 2006, Duncan Murdoch wrote:
> 

>> On 9/11/2006 8:20 AM, Prof Brian Ripley wrote:
>> > On Mon, 11 Sep 2006, Duncan Murdoch wrote:
>> >
>> > > On 9/11/2006 6:15 AM, Robin Hankin wrote:
>> > > > Hi Sundar
>> > > >
>> > > >
>> > > > thanks for this. But I didn't make it clear that I'm interested in
>> > > > extreme numbers
>> > > > such as 1e300 and 1e-300.

>> >
>> > That's not relevant, unless you are interested in extremely small numbers.
>> >
>> > > > Then
>> > > >
>> > > > > f(1e300)
>> > > > [1] 7.911257e+283
>> >
>> > (That is inaccurate)
>> >
>> > > > is different from
>> > > >
>> > > > 1e300*.Machine$double.eps
>> >
>> > Yes, since 1e300 is not a power of two. However, Sundar is right in the
>> > sense that this is an upper bound for normalized numbers.
>> >
>> > f(1) is .Machine$double.neg.eps, but so it is for all 1 <= x < 2.
>> > This gives you the answer: .Machine$double.neg.eps * 2^floor(log2(x))
>>
>> I'm not sure what is going wrong, but that is too small (on my machine, at
>> least):
>>
>> > f1 <- function(x) .Machine$double.neg.eps * 2^floor(log2(x))
>> > f1(1e300)
>> [1] 7.435085e+283
>> > 1e300 + f1(1e300) == 1e300
>> [1] TRUE
>>
>> Notice the difference in the 3rd decimal place from the empirical answer from
>> my bisection search below.
> 
> I wasn't going into that much detail: just what accuracy are we looking 
> for here?  .Machine$double.neg.eps isn't that accurate (as its help page 
> says).
> 

>> > Similarly for going below (but carefully as you get an extra halving on the
>> > powers of two).
>> >
>> > These results hold for all but denormalized numbers (those below 1e-308).
>> >
>> >
>> > > > [I'm interested in the gap between successive different exactly
>> > > > representable
>> > > > numbers right across the IEEE range]
>> > >
>> > > I'm not sure the result you're looking for is well defined, because on at
>> > > least the Windows platform, R makes use of 80 bit temporaries as well as
>> > > 64 bit double precision reals. I don't know any, but would guess there
>> > > exist examples of apparently equivalent formulations of your question that
>> > > give different answers because one uses the temporaries and the other
>> > > doesn't.

>> >
>> > Not at R level. For something to get stored in a real vector, it will be a
>> > standard 64-bit double.
>>
>> I don't think that's a proof, since R level code can call C functions, and
>> there are an awful lot of callable functions in R, but I don't have a
>> counter-example.
> 
> Oh, it is proof: the storage is declared as C double, and extended 
> precision double only exists on the chip.  Since R has to look up 'x' 
> whenever it sees the symbol, it looks at the stored value, not on a 
> register.  A compiler could do better, but the current code cannot.

It's a proof of something, but "apparently equivalent formulations" is a vague requirement. I would guess that if I did come up with a counter-example I would have to push the bounds a bit, and it would be questionable whether the formulations were equivalent.

For example, it would clearly be outside the bounds for me to write a function which did the entire computation in C (which would likely give a different answer than R gives). But what if some package already contains that function, and I just call it from R?

Duncan Murdoch



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 Mon Sep 11 23:47:43 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 - 14:30:04 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.