Re: [Rd] Numerical instability in new R Windows development version

From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>
Date: Fri, 27 Jan 2012 17:32:02 +0000

On 27/01/2012 13:26, Duncan Murdoch wrote:
> On 12-01-27 7:23 AM, Hans W Borchers wrote:
>> I have a question concerning the new Windows toolchain for R>= 2.14.2.
>> When trying out my package 'pracma' on the win-builder development
>> version
>> it will stop with the following error message:
>>
>> > f3<- function(x, y) sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1))
>> > dblquad(f3, -1, 1, -1, 1) # 2.094395124 , i.e. 2/3*pi , err = 2e-8
>> Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> Warning in sqrt((1 - (x^2 + y^2)) * (x^2 + y^2<= 1)) : NaNs produced
>> Error in integrate(function(y) f(x, y), ya, yb, subdivisions = subdivs, :
>> non-finite function value
>> Calls: dblquad ...
>> <Anonymous> -> f -> do.call -> mapply -> <Anonymous> -> integrate
>> Execution halted
>> ** running examples for arch 'x64' ... ERROR
>> Running examples in 'pracma-Ex.R' failed
>>
>> This probably means that the following expression got negative for some
>> values x, y:
>>
>> (1 - (x^2 + y^2)) * (x^2 + y^2<= 1)
>
> I think you're right, it's a bug, hopefully easy to fix. Here's a
> simpler version:
>
> x <- 0*(-1)
> sqrt(x)
>
> x is a "negative zero", and the sqrt() function incorrectly produces a
> NaN in the new toolchain.

Well, for some definition of 'incorrectly'. It is clearly what the author of that piece of code intended.

It would be helpful if people would cite definitive references. Someone is going to have to report this on the bugtracker, and at present I don't have enough evidence to do so: the C99/C11 standards do not seem to mandate a particular value (they do say what happens for values less than zero, but C compilers are allowed to have or not have signed zeroes). (Various Unix-alikes say what they do, usually -0, but that's not evidence that other answers are 'incorrect'.)

> Duncan Murdoch
>
>>
>> It appears to be an often used trick in numerical analysis. One
>> advantage is
>> that a function using it is immediately vectorized while an expression
>> such
>> as, e.g., "max(0, 1 - (x^2 + y^2))" is not.
>>
>> The example runs fine on Debian Linux and Mac OS X 32-/64-bit
>> architectures.
>> In my understanding the approach is correct and, as said above, often
>> used in
>> numerical applications.
>>
>> Can someone explain to me why this fails for the Windows 64-bit
>> compiler and
>> what I should use instead. Thanks.
>>
>> Hans Werner Borchers
>> ABB Corporate Research
>>
>> ______________________________________________
>> R-devel_at_r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> R-devel_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

-- 
Brian D. Ripley,                  ripley_at_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_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Received on Fri 27 Jan 2012 - 17:36:02 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

All messages

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 27 Jan 2012 - 19:10:12 GMT.

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

list of date sections of archive