// gcc main.c -std=c99 -lm <-m32> #include <stdio.h> #include <math.h> #include <complex.h> #include <fenv.h> #include <stdlib.h> // abs function for int argument #define ireal(a) (((int*)&(a))[0]) #define iimag(a) (((int*)&(a))[1]) #define real(a) (((float*)&(a))[0]) #define imag(a) (((float*)&(a))[1]) int main() { float _Complex a,b,r1,r2; ireal(a)=0xbf0e38e4; iimag(a)=0xbde38e3e; //argument ireal(b)=0x3f0e38e3; iimag(b)=0xbde38e3e; //argument fesetround(FE_TONEAREST); r1 = cpowf(a,b); //actual result r2 = cpow(a,b); //actual result 2 printf("inputs: a.real = (0x%08x)\t%e\n", ireal(a), real(a)); printf(" a.imag = (0x%08x)\t%e\n", iimag(a), imag(a)); printf("\n"); printf("inputs: b.real = (0x%08x)\t%e\n", ireal(b), real(b)); printf(" b.imag = (0x%08x)\t%e\n", iimag(b), imag(b)); printf("\n"); printf("actual.real = (0x%08x)\t%.13e\n", ireal(r1), real(r1)); printf("expected.real = (0x%08x)\t%.13e\n", ireal(r2), real(r2)); printf("error.real = %d ulp\n", abs(ireal(r1) - ireal(r2))); printf("\n"); printf("actual.imag = (0x%08x)\t%.13e\n", iimag(r1), imag(r1)); printf("expected.imag = (0x%08x)\t%.13e\n", iimag(r2), imag(r2)); printf("error.imag = %d ulp\n", abs(iimag(r1) - iimag(r2))); return 0; } Result x86_32 (core2): inputs: a.real = (0xbf0e38e4) -5.555556e-01 a.imag = (0xbde38e3e) -1.111111e-01 inputs: b.real = (0x3f0e38e3) 5.555555e-01 b.imag = (0xbde38e3e) -1.111111e-01 actual.real = (0xba6f8923) -9.1375614283606e-04 expected.real = (0xba6f8d77) -9.1382063692436e-04 error.real = 1108 ulp actual.imag = (0xbf069c6d) -5.2582436800003e-01 expected.imag = (0xbf069c6d) -5.2582436800003e-01 error.imag = 0 ulp Result x86_64 (core2): inputs: a.real = (0xbf0e38e4) -5.555556e-01 a.imag = (0xbde38e3e) -1.111111e-01 inputs: b.real = (0x3f0e38e3) 5.555555e-01 b.imag = (0xbde38e3e) -1.111111e-01 actual.real = (0xba6f8923) -9.1375614283606e-04 expected.real = (0xba6f8d77) -9.1382063692436e-04 error.real = 1108 ulp actual.imag = (0xbf069c6d) -5.2582436800003e-01 expected.imag = (0xbf069c6d) -5.2582436800003e-01 error.imag = 0 ulp
Again, please give testcases that use hex floats rather than using type punning, so that the test inputs can more readily be extracted for use elsewhere. cpowf, cpow, cpowl are generically inaccurate for all architectures and floating-point formats - there's nothing x86-specific. They are also by far the hardest complex functions to make accurate; the errors in this example are tiny compared to those you can get with properly chosen examples. In general for long double you need to compute log and atan2 results to over 16000 places to get good results for cpowl. Consider cpowl (LDBL_MAX, I * LDBL_MAX) for example - you need to compute LDBL_MAX * log (LDBL_MAX), reduced mod 2pi. Many cases will overflow or underflow but you still need to determine the signs of zero or infinity in those cases.
Reformatted test case: #include <stdio.h> #include <math.h> #include <complex.h> #include <fenv.h> #define real(a) (((float*)&(a))[0]) #define imag(a) (((float*)&(a))[1]) int main() { float _Complex a,b,r1; double _Complex r2; a=-0x1.1c71c8p-1 + I * -0x1.c71c7cp-4; //argument b= 0x1.1c71c6p-1 + I * -0x1.c71c7cp-4; //argument fesetround(FE_TONEAREST); r1 = cpowf(a,b); //actual result r2 = cpow(a,b); //actual result 2 printf("inputs: a = %15.6a + I * %15.6a\n", creal(a), cimag(a)); printf(" b = %15.6a + I * %15.6a\n", creal(b), cimag(b)); printf("actual = %15.6a + I * %15.6a\n", creal(r1),cimag(r1)); printf("expected = %15.6a + I * %15.6a\n", creal(r2),cimag(r2)); printf("error = %15f + I * %15f ulp\n", fabs(creal(r1) - creal(r2)) * exp2(23.0-logbf(creal(r1))), fabs(cimag(r1) - cimag(r2)) * exp2(23.0-logbf(cimag(r1))) ); return 0; } Results (same x86 and x86_64) inputs: a = -0x1.1c71c8p-1 + I * -0x1.c71c7cp-4 b = 0x1.1c71c6p-1 + I * -0x1.c71c7cp-4 actual = -0x1.df1246p-11 + I * -0x1.0d38dap-1 expected = -0x1.df1aefp-11 + I * -0x1.0d38d9p-1 error = 1108.410103 + I * 0.495932 ulp
Hello, as an additional data point, using cpow to compute i^2 yields a number with a non-zero imaginary part, as was noticed in later messages of the conversation started here: http://gcc.gnu.org/ml/libstdc++/2013-01/msg00058.html
Related LSB test suite failure "/olver/math/cexp/tests/cpow_spec" https://lsbbugs.linuxfoundation.org/show_bug.cgi?id=2912
There are at least two classes of accuracy issues with cpow(z,z2): (1) Insufficient precision in intermediate results, affecting primarily "large" arguments (2) Bad handling of |z1| == 1. Sample program below. This is actually easy to solve: use atan2pi, sinpi, and cospi instead of atan2, sin, and cos. #include <complex.h> #include <stdio.h> #include <math.h> int dont_optimize = 0; int main (int argc, char **argv) { double complex z1, z2, z; /* * z here should be real. I see an imaginary part of * 2.54571e-17, except when gcc is allowed to optimize * the whole thing at compile time. Exact result should * be -exp(-pi/2) =~ -0.20788 * * Surprise: Mathematica gets it wrong too. I see the * same value as gcc for N[I^(I+2)]. */ z1 = dont_optimize + I; z2 = 2 + I; z = cpow (z1, z2); printf ("%g + %g i\n", creal (z), cimag (z)); /* * Really scary version. Exact result is exp(-pi/2) * =~ 0.20788. I see "0.177058 + -0.108924 i". */ z2 = ldexp (1,53) + I; z = cpow (z1, z2); printf ("%g + %g i\n", creal (z), cimag (z)); return 0; }
Ok, make that three classes of inaccuracy: (3) Bad handling of real arguments. I see different values in this case: z1 = dont_optimize + 607/128; z2 = 1234567 / 8192; z = cpow (z1, z2); printf ("%.20g + %g i\n", creal (z), cimag (z)); printf ("%.20g\n", pow (creal (z1), creal (z2))); I get: 2.0370359763344673153e+90 + 0 i 2.0370359763344860863e+90 For real arguments cpow should defer to pow. (Or at least use a method that is no less accurate than pow; I think the current evidence is that pow is better.)