complex number compiler and libc bugs (cexp+conj) on OSX and with the intel compiler (icc)

Today a bug in complex number handling surfaced in guile [1] which only appeared on OSX.

This is a short note just to make sure that the bug is reported somewhere.

Test-code (written mostly by Mark Weaver who also analyzed the bug - I only ran the code on a few platforms I happened to have access to):

// test.c
// compile with gcc -O0 -o test test.c -lm
// or with icc -O0 -o test test.c -lm
#include <complex.h>
#include <stdio.h>

int
main (int argc, char **argv)
{
  double complex z = conj (1.0);
  double complex result;

  if (argc == 1)
    z = conj (0.0);

  result = cexp (z);

  printf ("cexp (%f + %f i) => %f + %f i\n",
          creal (z), cimag (z), creal (result), cimag (result));
  result = conj(result);
  printf ("conj(cexp (%f + %f i)) => %f + %f i\n",
          creal (z), cimag (z), creal (result), cimag (result));

  return 0;
}

As by the C-11 standard [2] (pages 561 and 216) this should return:

cexp (0.000000 + -0.000000 i) => 1.000000 + -0.000000 i

conj(cexp (0.000000 + -0.000000 i)) => 1.000000 + 0.000000 i

Page 561:

— cexp(conj(z)) = conj(cexp(z)).

Page 216:

The conj functions compute the complex conjugate of z, by reversing the sign of its imaginary part.

On OSX it returns (compiled with GCC):

TODO: Check the second line!

cexp (0.000000 + -0.000000 i) => 1.000000 + 0.000000 i

With the intel compiler it returns:

cexp (0.000000 + 0.000000 i) => 1.000000 + 0.000000 i

conj(cexp (0.000000 + 0.000000 i)) => 1.000000 + 0.000000 i

In short: On OSX cexp seems broken. With the intel compiler conj seems broken.

icc --version
# => icc (ICC) 13.1.3 20130607
# => Copyright (C) 1985-2013 Intel Corporation.  All rights reserved.

The OSX compiler is GCC 4.8.2 from MacPorts.

[taylanub] ArneBab: You might want to add that compiler optimizations can result in cexp() calls where there are none (which is how this bug surfaced in our case).

[mark_weaver] cexp(z) = e^z = e^(a+bi) = e^a * e^(bi) = e^a * (cos(b) + i*sin(b))

[mark_weaver] for real 'b', e^(bi) is a point on the unit circle on the complex plane.

[mark_weaver] so cexp(bi) can be used to compute cos(b) and sin(b) simultaneously, and probably faster than calling 'sin' and 'cos' separately.