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

[Eli Friedman]

On Fri, Dec 9, 2011 at 7:58 AM, Edward Meewis
<ed at extraordinarymachine.nl> wrote:
> 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?

I'm not getting quite the same results as you, but your testcase isn't
really testing what you want to test: std::exp is not the same as
::exp.

-Eli