[Rd] Friday question: negative zero

[Petr Savicky]

On Fri, Aug 31, 2007 at 08:39:02PM -0400, Duncan Murdoch wrote:
[snip]
> The other day I came across one in complex numbers, and it took me a 
> while to figure out that negative zero was what was happening:
> 
>   > x <- complex(real = -1)
>   > x
> [1] -1+0i
>   > 1/x
> [1] -1+0i
>   > x^(1/3)
> [1] 0.5+0.8660254i
>   > (1/x)^(1/3)
> [1] 0.5-0.8660254i
> 
> (The imaginary part of 1/x is negative zero.)
> 
> As a Friday question:  are there other ways to create and detect 
> negative zero in R?
> 
> And another somewhat more serious question:  is the behaviour of 
> negative zero consistent across platforms?  (The calculations above were 
> done in Windows in R-devel.)

Let me contribute to the statistics over platforms, which already contains
powerpc-apple-darwin8.9.1 (Steven McKinney) and two different results
on AMD 64 (deepayan.sarkar, Jeffrey Horner).

I tried  1/Im(1/complex(real = -1)) on three different platforms, all
of which are Linux, but with different CPU and different gcc version.
In all cases, R-devel_2007-08-31 is used.

The results differ. 

1. [...]$ grep name /proc/cpuinfo
   model name      : Pentium III (Coppermine)
   gcc version 4.0.2 20050901 (prerelease) (SUSE Linux)

   1/Im(1/complex(real = -1)) # Inf

2. [...]$ grep name /proc/cpuinfo
   model name      : Celeron (Coppermine)
   gcc version 4.1.0 (SUSE Linux)

   1/Im(1/complex(real = -1)) # -Inf

3. [...]$ grep name /proc/cpuinfo # two processor machine
   model name      : Intel(R) Xeon(TM) CPU 2.80GHz
   model name      : Intel(R) Xeon(TM) CPU 2.80GHz
   gcc version 3.3.3 20040412 (Red Hat Linux 3.3.3-7)

   1/Im(1/complex(real = -1)) # -Inf

The result does not depend on option -ffloat-store. On Xeon CPU,
it does not depend on, whether SSE (IEEE 754) arithmetic is
used (-march=pentium4 -mfpmath=sse) or not.

Let us try something more.

  z <- complex(real = -1)
  1/Im(z)        # [1] Inf
  1/Im(Conj(z))  # [1] -Inf

  z <- complex(real = -1, imaginary = -0)
  1/Im(z)        # [1] -Inf
  1/Im(Conj(z))  # [1] Inf

These results are consistent across platforms.
Let us look at 1/Im(1/z).

Platform 1 (Pentium III):

  z <- complex(real = -1)
  1/Im(1/z)      # [1] Inf

  z <- complex(real = -1, imaginary = -0)
  1/Im(1/z)      # [1] Inf

Platform 2 (Celeron):

  z <- complex(real = -1)
  1/Im(1/z)      # [1] -Inf

  z <- complex(real = -1, imaginary = -0)
  1/Im(1/z)      # [1] -Inf   (*)

Platform 3 (Xeon):

  z <- complex(real = -1)
  1/Im(1/z)      # [1] -Inf

  z <- complex(real = -1, imaginary = -0)
  1/Im(1/z)      # [1] Inf   (*)

So, besides the first difference between 1 and 2+3,
there is also a difference between 2 and 3 (marked (*)).

Petr Savicky.
===