will have followup with cexpf and others
-- v3: msvcr120: Add cexp() implementation. libs/musl: Pass math_error callback to exp. msvcr120/tests: Test exp.
From: Daniel Lehman dlehman25@gmail.com
--- dlls/msvcr120/tests/msvcr120.c | 71 ++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+)
diff --git a/dlls/msvcr120/tests/msvcr120.c b/dlls/msvcr120/tests/msvcr120.c index 856e9251e61..7c62a921a7e 100644 --- a/dlls/msvcr120/tests/msvcr120.c +++ b/dlls/msvcr120/tests/msvcr120.c @@ -195,6 +195,23 @@ static _Cancellation_beacon* (__thiscall *p__Cancellation_beacon_ctor)(_Cancella static void (__thiscall *p__Cancellation_beacon_dtor)(_Cancellation_beacon*); static MSVCRT_bool (__thiscall *p__Cancellation_beacon__Confirm_cancel)(_Cancellation_beacon*);
+static inline BOOL compare_double(double f, double g, unsigned int ulps) +{ + ULONGLONG x = *(ULONGLONG *)&f; + ULONGLONG y = *(ULONGLONG *)&g; + + if (f < 0) + x = ~x + 1; + else + x |= ((ULONGLONG)1)<<63; + if (g < 0) + y = ~y + 1; + else + y |= ((ULONGLONG)1)<<63; + + return (x>y ? x-y : y-x) <= ulps; +} + #define SETNOFAIL(x,y) x = (void*)GetProcAddress(module,y) #define SET(x,y) do { SETNOFAIL(x,y); ok(x != NULL, "Export '%s' not found\n", y); } while(0)
@@ -1786,6 +1803,59 @@ static void test__fsopen(void) setlocale(LC_ALL, "C"); }
+static int matherr_called; +static int CDECL matherr_callback(struct _exception *e) +{ + matherr_called = 1; + return 0; +} + +static void test_exp(void) +{ + static const struct { + double x, exp; + errno_t e; + BOOL todo; + } tests[] = { + { NAN, NAN, EDOM, TRUE }, + { -NAN, -NAN, EDOM, TRUE }, + { INFINITY, INFINITY }, + { -INFINITY, 0.0 }, + { 0.0, 1.0 }, + { 1.0, 2.7182818284590451 }, + { 709.7, 1.6549840276802644e+308 }, + { 709.782712893384, 1.7976931348622732e+308 }, + { 709.782712893385, INFINITY, ERANGE }, + }; + errno_t e; + double r; + int i; + + __setusermatherr(matherr_callback); + + for(i=0; i<ARRAY_SIZE(tests); i++) { + errno = 0; + matherr_called = 0; + r = exp(tests[i].x); + e = errno; + if(_isnan(tests[i].exp)) + ok(_isnan(r), "expected NAN, got %0.16e for %d\n", r, i); + else + ok(compare_double(r, tests[i].exp, 16), "expected %0.16e, got %0.16e for %d\n", tests[i].exp, r, i); + ok(signbit(r) == signbit(tests[i].exp), "expected sign %x, got %x for %d\n", + signbit(tests[i].exp), signbit(r), i); + todo_wine_if(tests[i].todo) { + ok(e == tests[i].e, "expected errno %i, but got %i for %d\n", tests[i].e, e, i); + if (tests[i].e) + ok(matherr_called, "matherr wasn't called for %d\n", i); + else + ok(!matherr_called, "matherr was called for %d\n", i); + } + } + + __setusermatherr(NULL); +} + START_TEST(msvcr120) { if (!init()) return; @@ -1813,4 +1883,5 @@ START_TEST(msvcr120) test_strcmp(); test_gmtime64(); test__fsopen(); + test_exp(); }
From: Daniel Lehman dlehman25@gmail.com
exp calls matherr callback. cexp does not cexp calls into exp --- libs/musl/src/internal/libm.h | 1 + libs/musl/src/math/exp.c | 19 ++++++++++++------- 2 files changed, 13 insertions(+), 7 deletions(-)
diff --git a/libs/musl/src/internal/libm.h b/libs/musl/src/internal/libm.h index fa3ec4bd607..e24e10edaa6 100644 --- a/libs/musl/src/internal/libm.h +++ b/libs/musl/src/internal/libm.h @@ -10,6 +10,7 @@ typedef float float_t; typedef double double_t;
+typedef double matherr_t(int, const char *, double, double, double); hidden double math_error(int type, const char *name, double arg1, double arg2, double retval);
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 diff --git a/libs/musl/src/math/exp.c b/libs/musl/src/math/exp.c index db67213b7d3..0c518a31ebd 100644 --- a/libs/musl/src/math/exp.c +++ b/libs/musl/src/math/exp.c @@ -28,7 +28,7 @@ a double. (int32_t)KI is the k used in the argument reduction and exponent adjustment of scale, positive k here means the result may overflow and negative k means the result may underflow. */ -static inline double specialcase(double x, double_t tmp, uint64_t sbits, uint64_t ki) +static inline double specialcase(double x, double_t tmp, uint64_t sbits, uint64_t ki, matherr_t matherr) { double_t scale, y;
@@ -38,7 +38,7 @@ static inline double specialcase(double x, double_t tmp, uint64_t sbits, uint64_ scale = asdouble(sbits); y = 0x1p1009 * (scale + scale * tmp); if (isinf(y)) - return math_error(_OVERFLOW, "exp", x, 0, y); + return matherr(_OVERFLOW, "exp", x, 0, y); return eval_as_double(y); } /* k < 0, need special care in the subnormal range. */ @@ -61,7 +61,7 @@ static inline double specialcase(double x, double_t tmp, uint64_t sbits, uint64_ /* The underflow exception needs to be signaled explicitly. */ fp_force_eval(fp_barrier(0x1p-1022) * 0x1p-1022); y = 0x1p-1022 * y; - return math_error(_UNDERFLOW, "exp", x, 0, y); + return matherr(_UNDERFLOW, "exp", x, 0, y); } y = 0x1p-1022 * y; return eval_as_double(y); @@ -73,7 +73,7 @@ static inline uint32_t top12(double x) return asuint64(x) >> 52; }
-double __cdecl exp(double x) +double __cdecl __exp(double x, matherr_t matherr) { uint32_t abstop; uint64_t ki, idx, top, sbits; @@ -91,9 +91,9 @@ double __cdecl exp(double x) if (abstop >= top12(INFINITY)) return 1.0 + x; if (asuint64(x) >> 63) - return math_error(_UNDERFLOW, "exp", x, 0, fp_barrier(DBL_MIN) * DBL_MIN); + return matherr(_UNDERFLOW, "exp", x, 0, fp_barrier(DBL_MIN) * DBL_MIN); else - return math_error(_OVERFLOW, "exp", x, 0, fp_barrier(DBL_MAX) * DBL_MAX); + return matherr(_OVERFLOW, "exp", x, 0, fp_barrier(DBL_MAX) * DBL_MAX); } /* Large x is special cased below. */ abstop = 0; @@ -118,9 +118,14 @@ double __cdecl exp(double x) /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp. */ tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5); if (predict_false(abstop == 0)) - return specialcase(x, tmp, sbits, ki); + return specialcase(x, tmp, sbits, ki, matherr); scale = asdouble(sbits); /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there is no spurious underflow here even without fma. */ return eval_as_double(scale + scale * tmp); } + +double __cdecl exp(double x) +{ + return __exp(x, math_error); +}
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
This merge request was approved by Piotr Caban.
-
Daniel Lehman -
Daniel Lehman (＠dlehman25)
-
Piotr Caban (＠piotr)