[R] exactly representable numbers
Duncan Murdoch
murdoch at stats.uwo.ca
Mon Sep 11 15:40:28 CEST 2006
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
More information about the R-help
mailing list