From stephen.b.weston at gmail.com  Fri Jan 27 20:04:01 2012
From: stephen.b.weston at gmail.com (Stephen Weston)
Date: Fri, 27 Jan 2012 14:04:01 -0500
Subject: [Rd] Numerical instability in new R Windows development version
In-Reply-To: <4F22EFB6.9010407@gmail.com>
References: <CAML4n3MNRsp_kTpPHarDh7_kLCVwGZJkweOxkoo6S8OMwENfjg@mail.gmail.com>
	<4F22A61A.4020301@gmail.com> <4F22DF92.9000204@stats.ox.ac.uk>
	<4F22EFB6.9010407@gmail.com>
Message-ID: <CALh21iL37vgi=r2rHSx9okiCihtuc9Z3xccT-XE8=-H7Lh+3ig@mail.gmail.com>

On Fri, Jan 27, 2012 at 1:40 PM, Duncan Murdoch
<murdoch.duncan at gmail.com> 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'.)
>
>
> This page:
>
> http://pubs.opengroup.org/onlinepubs/9699919799/functions/sqrt.html
>
> also says that the correct answer is -0. ?I don't know if that has any
> authority at all...
>
> Duncan Murdoch

Section 5.4.1 "Arithmetic operations" of "IEEE Standard for Floating-
Point Arithmetic", IEEE Std 754-2008 says:

  "The operation squareRoot(x) computes ? x. It has a positive sign for
all operands ? 0, except that squareRoot(?0) shall be ?0."

- Steve



[Rd] Numerical instability in new R Windows development version
Stephen Weston