From: Daniel Lehman dlehman25@gmail.com
--- dlls/msvcr120/msvcr120.spec | 2 +- dlls/msvcr120/tests/msvcr120.c | 99 ++++++++++++++++++++++ dlls/msvcr120_app/msvcr120_app.spec | 2 +- dlls/ucrtbase/ucrtbase.spec | 2 +- include/msvcrt/complex.h | 2 + libs/musl/Makefile.in | 5 +- libs/musl/src/complex/__cexp.c | 87 +++++++++++++++++++ libs/musl/src/complex/cexp.c | 117 ++++++++++++++++++++++++++ libs/musl/src/internal/complex_impl.h | 18 ++++ 9 files changed, 330 insertions(+), 4 deletions(-) create mode 100644 libs/musl/src/complex/__cexp.c create mode 100644 libs/musl/src/complex/cexp.c create mode 100644 libs/musl/src/internal/complex_impl.h
diff --git a/dlls/msvcr120/msvcr120.spec b/dlls/msvcr120/msvcr120.spec index 5c588bd5889..56b26852ce1 100644 --- a/dlls/msvcr120/msvcr120.spec +++ b/dlls/msvcr120/msvcr120.spec @@ -2073,7 +2073,7 @@ @ stub ccosl @ cdecl ceil(double) @ cdecl -arch=!i386 ceilf(float) -@ stub cexp +@ cdecl -norelay cexp(int128) @ stub cexpf @ stub cexpl @ cdecl cimag(int128) diff --git a/dlls/msvcr120/tests/msvcr120.c b/dlls/msvcr120/tests/msvcr120.c index 7c62a921a7e..5a200d5e705 100644 --- a/dlls/msvcr120/tests/msvcr120.c +++ b/dlls/msvcr120/tests/msvcr120.c @@ -1852,6 +1852,104 @@ static void test_exp(void) ok(!matherr_called, "matherr was called for %d\n", i); } } +} + +static void test_cexp(void) +{ + static const struct { + double r, i; + double rexp, iexp; + errno_t e; + } tests[] = { + { INFINITY, 0.0, INFINITY, 0.0 }, + { INFINITY, -0.0, INFINITY, -0.0 }, + { -INFINITY, 0.0, 0.0, 0.0 }, + { -INFINITY, -0.0, 0.0, -0.0 }, + { 0.0, INFINITY, NAN, NAN, EDOM }, + { -0.0, INFINITY, NAN, NAN, EDOM }, + { 0.0, -INFINITY, NAN, NAN, EDOM }, + { -0.0, -INFINITY, NAN, NAN, EDOM }, + { 100.0, INFINITY, NAN, NAN, EDOM }, + { -100.0, INFINITY, NAN, NAN, EDOM }, + { 100.0, -INFINITY, NAN, NAN, EDOM }, + { -100.0, -INFINITY, NAN, NAN, EDOM }, + { -INFINITY, 2.0, -0.0, 0.0 }, + { -INFINITY, 4.0, -0.0, -0.0 }, + { INFINITY, 2.0, -INFINITY, INFINITY }, + { INFINITY, 4.0, -INFINITY, -INFINITY }, + { INFINITY, INFINITY, INFINITY, NAN, EDOM }, + { INFINITY, -INFINITY, INFINITY, NAN, EDOM }, + { -INFINITY, INFINITY, 0.0, 0.0 }, + { -INFINITY, -INFINITY, 0.0, -0.0 }, + { -INFINITY, NAN, 0.0, 0.0 }, + { INFINITY, NAN, INFINITY, NAN }, + { NAN, 0.0, NAN, 0.0 }, + { NAN, -0.0, NAN, -0.0 }, + { NAN, 1.0, NAN, NAN }, + { NAN, INFINITY, NAN, NAN }, + { 0.0, NAN, NAN, NAN }, + { 1.0, NAN, NAN, NAN }, + { NAN, NAN, NAN, NAN }, + }; + static const struct { + double r, i; + double rexp, iexp; + errno_t e; + } tests2[] = { + { 0.0, M_PI, -1.0, 1.2246467991473532e-016 }, + { 709.7, 0.0, 1.6549840276802644e+308, 0.0 }, + { 709.78271289338397, 0.0, 1.7976931348622734e+308, 0.0 }, + { 709.8, 0.0, INFINITY, 0.0, ERANGE }, + { 746.0, 0.0, INFINITY, 0.0, ERANGE }, + { 747.0, 0.0, INFINITY, 0.0, ERANGE }, + { 1454.3, 0.0, INFINITY, 0.0, ERANGE }, + { 709.7, M_PI, -1.6549840276802644e+308, 2.0267708921386307e+292 }, + { 709.78271289338397, M_PI, -1.7976931348622734e+308, 2.2015391434582545e+292 }, + { 709.8, M_PI, -INFINITY, 2.2399282475936458e+292, ERANGE }, + { 746.0, M_PI, -INFINITY, 1.1794902399837951e+308, ERANGE }, + { 747.0, M_PI, -INFINITY, INFINITY, ERANGE }, + { 1454.3, M_PI, -INFINITY, INFINITY, ERANGE }, + }; + _Dcomplex c, r; + errno_t e; + int i; + + __setusermatherr(matherr_callback); + + for(i=0; i<ARRAY_SIZE(tests); i++) { + c = _Cbuild(tests[i].r, tests[i].i); + errno = 0; + matherr_called = 0; + r = cexp(c); + e = errno; + if(_isnan(tests[i].rexp)) + ok(_isnan(r._Val[0]), "expected NAN, got %0.16e for %d\n", r._Val[0], i); + else + ok(r._Val[0] == tests[i].rexp, "expected %0.16e, got %0.16e for %d\n", tests[i].rexp, r._Val[0], i); + if(_isnan(tests[i].iexp)) + ok(_isnan(r._Val[1]), "expected NAN, got %0.16e for %d\n", r._Val[1], i); + else + ok(r._Val[1] == tests[i].iexp, "expected %0.16e, got %0.16e for %d\n", tests[i].iexp, r._Val[1], i); + ok(signbit(r._Val[0]) == signbit(tests[i].rexp), "expected sign %x, got %x for %d\n", + signbit(tests[i].rexp), signbit(r._Val[0]), i); + ok((signbit(r._Val[1]) == signbit(tests[i].iexp)) || + broken(tests[i].r == -INFINITY && tests[i].i == -0.0) /* older win10 */, + "expected sign %x, got %x for %d\n", signbit(tests[i].iexp), signbit(r._Val[1]), i); + ok(e == tests[i].e, "expected errno %i, but got %i for %d\n", tests[i].e, e, i); + ok(!matherr_called, "matherr was called for %d\n", i); + } + + for(i=0; i<ARRAY_SIZE(tests2); i++) { + errno = 0; + matherr_called = 0; + c = _Cbuild(tests2[i].r, tests2[i].i); + r = cexp(c); + e = errno; + ok(compare_double(r._Val[0], tests2[i].rexp, 16), "expected %0.16e, got %0.16e for real %d\n", tests2[i].rexp, r._Val[0], i); + ok(compare_double(r._Val[1], tests2[i].iexp, 16), "expected %0.16e, got %0.16e for imag %d\n", tests2[i].iexp, r._Val[1], i); + ok(e == tests2[i].e, "expected errno %i, but got %i for %d\n", tests2[i].e, e, i); + ok(!matherr_called, "matherr was called for %d\n", i); + }
__setusermatherr(NULL); } @@ -1884,4 +1982,5 @@ START_TEST(msvcr120) test_gmtime64(); test__fsopen(); test_exp(); + test_cexp(); } diff --git a/dlls/msvcr120_app/msvcr120_app.spec b/dlls/msvcr120_app/msvcr120_app.spec index fc2ac655c40..c7cfe4e071e 100644 --- a/dlls/msvcr120_app/msvcr120_app.spec +++ b/dlls/msvcr120_app/msvcr120_app.spec @@ -1740,7 +1740,7 @@ @ stub ccosl @ cdecl ceil(double) msvcr120.ceil @ cdecl -arch=!i386 ceilf(float) msvcr120.ceilf -@ stub cexp +@ cdecl -norelay cexp(int128) msvcr120.cexp @ stub cexpf @ stub cexpl @ cdecl cimag(int128) msvcr120.cimag diff --git a/dlls/ucrtbase/ucrtbase.spec b/dlls/ucrtbase/ucrtbase.spec index e03e0b51d60..230f5ceb1d7 100644 --- a/dlls/ucrtbase/ucrtbase.spec +++ b/dlls/ucrtbase/ucrtbase.spec @@ -2218,7 +2218,7 @@ @ stub ccosl @ cdecl ceil(double) @ cdecl -arch=!i386 ceilf(float) -@ stub cexp +@ cdecl -norelay cexp(int128) @ stub cexpf @ stub cexpl @ cdecl cimag(int128) diff --git a/include/msvcrt/complex.h b/include/msvcrt/complex.h index f8c6fdd47bb..78ddaf5088c 100644 --- a/include/msvcrt/complex.h +++ b/include/msvcrt/complex.h @@ -25,6 +25,8 @@ typedef _C_double_complex _Dcomplex; typedef _C_float_complex _Fcomplex;
_ACRTIMP _Dcomplex __cdecl _Cbuild(double, double); +_ACRTIMP _Dcomplex __cdecl cexp(_Dcomplex); + _ACRTIMP double __cdecl cimag(_Dcomplex); _ACRTIMP double __cdecl creal(_Dcomplex);
diff --git a/libs/musl/Makefile.in b/libs/musl/Makefile.in index bcb14ffef78..aed423e29b8 100644 --- a/libs/musl/Makefile.in +++ b/libs/musl/Makefile.in @@ -1,8 +1,11 @@ EXTLIB = libmusl.a EXTRAINCL = -I$(srcdir)/src/internal -I$(srcdir)/arch/generic -EXTRADEFS = -D_ACRTIMP= -D_NO_CRT_MATH_INLINE +EXTRADEFS = -D_ACRTIMP= -D_NO_CRT_MATH_INLINE \ + -fno-builtin
SOURCES = \ + src/complex/__cexp.c \ + src/complex/cexp.c \ src/math/__cos.c \ src/math/__cosdf.c \ src/math/__expo2.c \ diff --git a/libs/musl/src/complex/__cexp.c b/libs/musl/src/complex/__cexp.c new file mode 100644 index 00000000000..eb9f9a82790 --- /dev/null +++ b/libs/musl/src/complex/__cexp.c @@ -0,0 +1,87 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/k_exp.c */ +/*- + * Copyright (c) 2011 David Schultz das@FreeBSD.ORG + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include "complex_impl.h" + +static const uint32_t k = 1799; /* constant for reduction */ +static const double kln2 = 1246.97177782734161156; /* k * ln2 */ + +/* + * Compute exp(x), scaled to avoid spurious overflow. An exponent is + * returned separately in 'expt'. + * + * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 + * Output: 2**1023 <= y < 2**1024 + */ +static double __frexp_exp(double x, int *expt) +{ + double exp_x; + uint32_t hx; + + /* + * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to + * minimize |exp(kln2) - 2**k|. We also scale the exponent of + * exp_x to MAX_EXP so that the result can be multiplied by + * a tiny number without losing accuracy due to denormalization. + */ + exp_x = exp(x - kln2); + GET_HIGH_WORD(hx, exp_x); + *expt = (hx >> 20) - (0x3ff + 1023) + k; + SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); + return exp_x; +} + +/* + * __ldexp_cexp(x, expt) compute exp(x) * 2**expt. + * It is intended for large arguments (real part >= ln(DBL_MAX)) + * where care is needed to avoid overflow. + * + * The present implementation is narrowly tailored for our hyperbolic and + * exponential functions. We assume expt is small (0 or -1), and the caller + * has filtered out very large x, for which overflow would be inevitable. + */ +_Dcomplex __ldexp_cexp(_Dcomplex z, int expt) +{ + double x, y, exp_x, scale1, scale2; + int ex_expt, half_expt; + + x = creal(z); + y = cimag(z); + exp_x = __frexp_exp(x, &ex_expt); + expt += ex_expt; + + /* + * Arrange so that scale1 * scale2 == 2**expt. We use this to + * compensate for scalbn being horrendously slow. + */ + half_expt = expt / 2; + INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); + half_expt = expt - half_expt; + INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); + + return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2); +} diff --git a/libs/musl/src/complex/cexp.c b/libs/musl/src/complex/cexp.c new file mode 100644 index 00000000000..9330898a2cb --- /dev/null +++ b/libs/musl/src/complex/cexp.c @@ -0,0 +1,117 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_cexp.c */ +/*- + * Copyright (c) 2011 David Schultz das@FreeBSD.ORG + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include "complex_impl.h" + +static const uint32_t +exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ +cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ + +extern double __cdecl __exp(double x, matherr_t matherr); + +static double cexp_callback( int type, const char *name, double arg1, double arg2, double retval ) +{ + + switch (type) { + case _DOMAIN: + errno = EDOM; + break; + case _SING: + case _OVERFLOW: + errno = ERANGE; + break; + } + return retval; +} + +_Dcomplex __cdecl cexp(_Dcomplex z) +{ + _Dcomplex r; + double x, y, exp_x; + uint32_t hx, hy, lx, ly; + + x = creal(z); + y = cimag(z); + + EXTRACT_WORDS(hy, ly, y); + hy &= 0x7fffffff; + + /* cexp(x + I 0) = exp(x) + I 0 */ + if ((hy | ly) == 0) { + if (isnan(x)) return CMPLX(x, y); + return CMPLX(__exp(x, cexp_callback), y); + } + EXTRACT_WORDS(hx, lx, x); + /* cexp(0 + I y) = cos(y) + I sin(y) */ + if (((hx & 0x7fffffff) | lx) == 0) { + if (isinf(y)) { + errno = EDOM; + return CMPLX(NAN, NAN); + } + return CMPLX(cos(y), sin(y)); + } + + if (hy >= 0x7ff00000) { + if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { + /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ + if (isinf(y)) { + if (!isnan(x)) errno = EDOM; + return CMPLX(NAN, NAN); + } + return CMPLX(y - y, y - y); + } else if (hx & 0x80000000) { + /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ + return CMPLX(0.0, copysign(0.0, y)); + } else { + /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ + if (!isnan(y)) errno = EDOM; + return CMPLX(x, NAN); + } + } + + if (hx >= exp_ovfl && hx <= cexp_ovfl) { + /* + * x is between 709.7 and 1454.3, so we must scale to avoid + * overflow in exp(x). + */ + r = __ldexp_cexp(z, 0); + if (isinf(r._Val[0])) errno = ERANGE; + return r; + } else { + /* + * Cases covered here: + * - x < exp_ovfl and exp(x) won't overflow (common case) + * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 + * - x = +-Inf (generated by exp()) + * - x = NaN (spurious inexact exception from y) + */ + if (isnan(x)) + return CMPLX(NAN, NAN); + exp_x = __exp(x, cexp_callback); + return CMPLX(exp_x * cos(y), exp_x * sin(y)); + } +} diff --git a/libs/musl/src/internal/complex_impl.h b/libs/musl/src/internal/complex_impl.h new file mode 100644 index 00000000000..c1ae69f45f2 --- /dev/null +++ b/libs/musl/src/internal/complex_impl.h @@ -0,0 +1,18 @@ +#ifndef _COMPLEX_IMPL_H +#define _COMPLEX_IMPL_H + +#include <complex.h> +#include "libm.h" + +#undef CMPLX +#undef CMPLXF +#undef CMPLXL + +#define CMPLX(x, y) _Cbuild(x, y) +#define CMPLXF(x, y) _FCbuild(x, y) +#define CMPLXL(x, y) _LCbuild(x, y) + +hidden _Dcomplex __ldexp_cexp(_Dcomplex,int); +hidden _Fcomplex __ldexp_cexpf(_Fcomplex,int); + +#endif