[cfe-dev] [libcxx/test] Puzzling floating point behaviour

[Edward Meewis]

Hi all,

When running the libcxx I got asserts in complex.transcendentals. After 
a lot of digging I have the following:

-- pow_complex_test.cpp:
#include <complex>
#include <cmath>
#include <cstdio>

int main()
{
     const std::complex<long double> a = std::complex<long double>(2, 3);
     const std::complex<long double> b = std::complex<long double>(2, 0);
     std::complex<long double> x = std::complex<long double>(-5, 12);

     std::complex<long double> c = pow(a, b);

     printf("c = pow(a,b) = %.20Le + i*%.20Le\n", real(c), imag(c));
     printf("x = %.20Le + i*%.20Le\n\n", real(x), imag(x));
     printf("err = c - x = %.20Le + i*%.20Le\n\n", real(c-x), imag(c-x));

     long double r_a = sqrt(real(a)*real(a) + imag(a)*imag(a));
     long double th_a = atan(imag(a)/real(a));

     /* bw = b*log(a) = 2*log(r_a) + i*2*th_a */

     std::complex<long double> bw = b*log(a);

     /* d = exp(x) = exp(xr) * (cos(xi) + i*sin(xi) */

     std::complex<long double> d = exp(b*log(a));
     std::complex<long double> dx0 = exp(std::complex<long 
double>(bw.real(), bw.imag()));
     std::complex<long double> dx1 =
         std::complex<long double>(exp(bw.real()) * cos(bw.imag()),
                                   exp(bw.real()) * sin(bw.imag()));
     std::complex<long double> dx2 =
         (long double)exp(bw.real()) *
         std::complex<long double>(cos(bw.imag()), sin(bw.imag()));

     printf("d = exp(b*log(a)) = %.20Le + i*%.20Le\n", real(d), imag(d));
     printf("dx0 = %.20Le + i*%.20Le\n", real(dx0), imag(dx0));
     printf("dx1 = %.20Le + i*%.20Le\n", real(dx1), imag(dx1));
     printf("dx2 = %.20Le + i*%.20Le\n\n", real(dx2), imag(dx2));
     printf("err = d - x = %.20Le + i*%.20Le\n", real(d-x), imag(d-x));
     printf("err = dx0 - x = %.20Le + i*%.20Le\n", real(dx0-x), 
imag(dx0-x));
     printf("err = dx1 - x = %.20Le + i*%.20Le\n", real(dx1-x), 
imag(dx1-x));
     printf("err = dx2 - x = %.20Le + i*%.20Le\n", real(dx2-x), 
imag(dx2-x));

}
--
clang -std=c++0x -nostdinc++ -I/home/emeewis/contrib/libcxx/include 
-nodefaultlibs -L/home/emeewis/contrib/libcxx/lib -lc++ 
-L/usr/home/emeewis/contrib/FreeBSD-HEAD/lib/libcxxrt -lcxxrt -lm 
pow_complex_test.cpp && ./a.out

c = pow(a,b) = -4.99999999981561327189e+00 + i*1.19999999999030784181e+01
x = -5.00000000000000000000e+00 + i*1.20000000000000000000e+01

err = c - x = 1.84386728108165698359e-10 + i*-9.69215818713564658538e-11

d = exp(b*log(a)) = -4.99999999981561327189e+00 + 
i*1.19999999999030784181e+01
dx0 = -4.99999999981561327189e+00 + i*1.19999999999030784181e+01
dx1 = -5.00000000000000000000e+00 + i*1.20000000000000000000e+01
dx2 = -5.00000000000000000000e+00 + i*1.20000000000000000000e+01

err = d - x = 1.84386728108165698359e-10 + i*-9.69215818713564658538e-11
err = dx0 - x = 1.84386728108165698359e-10 + i*-9.69215818713564658538e-11
err = dx1 - x = 0.00000000000000000000e+00 + i*0.00000000000000000000e+00
err = dx2 - x = 0.00000000000000000000e+00 + i*0.00000000000000000000e+00

So, calculating c = pow(a,b) ever more directly eventually getting the 
correct result.

Compiling against libstdc++ gets correct results.

I'm guessing long doubles get demoted somewhere. Does anyone have an 
idea what the problem may be?

Regards, Edward