[Rd] Numerical instability in new R Windows development version

Duncan Murdoch murdoch.duncan at gmail.com
Fri Jan 27 14:26:50 CET 2012


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.

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



More information about the R-devel mailing list