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.
Links:
[1] http://gnu.org/s/guile
[2] http://www.open-std.org/jtc1/sc22/wg14/www/docs/n1570.pdf