From murdoch.duncan at gmail.com  Fri Jan 27 19:26:17 2012
From: murdoch.duncan at gmail.com (Duncan Murdoch)
Date: Fri, 27 Jan 2012 13:26:17 -0500
Subject: [Rd] Numerical instability in new R Windows development version
In-Reply-To: <4F22DF92.9000204@stats.ox.ac.uk>
References: <CAML4n3MNRsp_kTpPHarDh7_kLCVwGZJkweOxkoo6S8OMwENfjg@mail.gmail.com>
	<4F22A61A.4020301@gmail.com> <4F22DF92.9000204@stats.ox.ac.uk>
Message-ID: <4F22EC49.80603@gmail.com>

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.

Duncan Murdoch



[Rd] Numerical instability in new R Windows development version
Duncan Murdoch