[Rd] Numerical instability in new R Windows development version

From: Hans W Borchers <hwborchers_at_googlemail.com>
Date: Fri, 27 Jan 2012 13:23:45 +0100


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

This probably means that the following expression got negative for some values x, y:

  (1 - (x^2 + y^2)) * (x^2 + y^2 <= 1)

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 Received on Fri 27 Jan 2012 - 12:40:49 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 Mon 30 Jan 2012 - 22:40:13 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