Re: [Rd] Computing means, variances and sums

From: Duncan Murdoch <murdoch_at_stats.uwo.ca>
Date: Wed 22 Feb 2006 - 12:11:44 GMT

On 2/22/2006 3:52 AM, Prof Brian Ripley wrote:
> 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

I think Windows also preserves the FPU control word across context switches. The problem is that it doesn't do a context switch when you make a call into the system. In particular, the shell services (file dialogs, URL recognition, etc.) involve a large number of DLLs, all running within the same context. MSVC++ normally chooses 53 bit precision as the default, and will set that when a DLL starts up, or whenever a function in it asks to reset the FPU.

Duncan Murdoch
>
> 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
>>
>>
>

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

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