Re: [Rd] Computing means, variances and sums

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Wed 22 Feb 2006 - 08:52:13 GMT

I've managed to track this down. The setting of the FPU control word on a ix86 machine changes the precision of (gcc) long double calculations in the FPU, as well as those of double. So if it gets changed to PC_53, long doubles lose accuracy even though 10 bytes are carried around.

For some discussion of extended-precision issues, see Priest's annex to Goldberg's paper at http://www.validlab.com/goldberg/paper.pdf.

Many other OSes other than Windows on ix86 do consider the FPU control word when doing context-switching. There is (often heated) debate as to what the correct default should be, see e.g. the thread begining

http://gcc.gnu.org/ml/gcc/1998-12/msg00473.html

Whereas Linux has selected extended precision, apparently FreeBSD selected 53-bit mantissa, and Solaris x86 used the equivalent of -ffloat-store to ensure stricter compliant to the IEEE754 'double' model (and cynics say, to make the Sparc Solaris performance look good).

It will be interesting to see what MacIntel has selected (I suspect the FreeBSD solution).

This does mean that using long doubles for accumulators is not necessarily always a win although the potential loss is likely to be small.

On Sun, 19 Feb 2006, Duncan Murdoch wrote:

> On 2/19/2006 3:18 PM, Prof Brian Ripley wrote:
>> On Sun, 19 Feb 2006, hadley wickham wrote:
>>
>>>> p.s. If my computations are correct, 0.2 = 0*/2 + 0/4 + 1/8 + 1/16 +
>>>> 0/32 + 0/64 + 1/128 + 1/256 + 0/512 + 0/1024 + 1/2048 + 1/4096 + ... =
>>>> 0.3333333333333h. Perhaps someone can extend this to an FAQ to help
>>>> explain finite precision arithmetic and rounding issues.
>>> This is drifting a bit off topic, but the other day I discovered this
>>> rather nice illustration of the perils of finite precision arithmetic
>>> while creating a contrast matrix:
>>>
>>>> n <- 13
>>>> a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
>>>> rowSums(a)
>>> [1] 2.775558e-16 2.775558e-16 5.551115e-17 5.551115e-17 5.551115e-17
>>> [6] 5.551115e-17 0.000000e+00 -5.551115e-17 0.000000e+00 5.551115e-17
>>> [11] 1.110223e-16 1.665335e-16 2.220446e-16
>>>
>>> Not only do most of the rows not sum to 0, they do not even sum to the
>>> same number! It is hard to remember the familiar rules of arithmetic
>>> do not always apply.
>>
>> I think you will find this example does give all 0's in R-devel, even on
>> platforms like Sparc.
>
> Only until the fpu precision gets changed:
>
>> n <- 13
>> a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
>> rowSums(a)
> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0
>> RSiteSearch('junk')
> A search query has been submitted to http://search.r-project.org
> The results page should open in your browser shortly
>
>> n <- 13
>> a <- matrix(-1/n, ncol=n, nrow=n) + diag(n)
>> rowSums(a)
> [1] 2.775558e-16 2.775558e-16 5.551115e-17 5.551115e-17 5.551115e-17
> [6] 5.551115e-17 0.000000e+00 -5.551115e-17 0.000000e+00 5.551115e-17
> [11] 1.110223e-16 1.665335e-16 2.220446e-16
>
> We still need to protect against these changes. I'll put something together,
> unless you're already working on it.
>
> The approach I'm thinking of is to define a macro to be called in risky
> situations. On platforms where this isn't an issue, the macro would be null;
> on Windows, it would reset the fpu to full precision.
>
> For example, RSiteSearch causes damage in the ShellExecute call in
> do_shellexec called from browseURL, so I'd add protection there. I think we
> should also add detection code somewhere in the evaluation loop to help in
> diagnosing these problems.
>
>> But users do need to remember that computer arithmetic is inexact except in
>> rather narrowly delimited cases.
>
> Yes, that too.
>
> Duncan Murdoch
>
>

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Wed Feb 22 19:54:53 2006

This archive was generated by hypermail 2.1.8 : Wed 22 Feb 2006 - 19:30:50 GMT