[Rd] Numerical instability in new R Windows development version

[Duncan Murdoch]

On 27/01/2012 1:26 PM, Duncan Murdoch wrote:
> On 27/01/2012 12:32 PM, Prof Brian Ripley wrote:
> >  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'.)
>
> Section 6.3 of IEEE 754-2008 says
>
> Except that squareRoot(−0) shall be −0, every numeric squareRoot result
> shall have a positive sign.

I believe the corresponding ISO standard is ISO/IEC/IEEE 60559:2011, but I don't have a copy, and I don't think my library does.

Duncan Murdoch