[PATCH 0/2] MR7748: ucrtbase: Add cexp() implementation.
will have followup with cexpf and others -- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748
From: Daniel Lehman <dlehman25(a)gmail.com> --- include/msvcrt/complex.h | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/include/msvcrt/complex.h b/include/msvcrt/complex.h index 2d7009c470c..53a72b6a000 100644 --- a/include/msvcrt/complex.h +++ b/include/msvcrt/complex.h @@ -24,4 +24,10 @@ typedef struct _C_float_complex typedef _C_double_complex _Dcomplex; typedef _C_float_complex _Fcomplex; +double __cdecl cimag(_Dcomplex); +double __cdecl creal(_Dcomplex); + +float __cdecl cimagf(_Fcomplex); +float __cdecl crealf(_Fcomplex); + #endif /* _COMPLEX_H_DEFINED */ -- GitLab https://gitlab.winehq.org/wine/wine/-/merge_requests/7748
From: Daniel Lehman <dlehman25(a)gmail.com> with following changes: /double complex/_Dcomplex/ --- dlls/msvcr120/msvcr120.spec | 2 +- dlls/msvcr120/tests/msvcr120.c | 24 ++++++++ dlls/msvcr120_app/msvcr120_app.spec | 2 +- dlls/ucrtbase/ucrtbase.spec | 2 +- include/msvcrt/complex.h | 5 +- libs/musl/Makefile.in | 4 +- libs/musl/src/complex/__cexp.c | 87 +++++++++++++++++++++++++++ libs/musl/src/complex/cexp.c | 83 +++++++++++++++++++++++++ libs/musl/src/internal/complex_impl.h | 22 +++++++ 9 files changed, 225 insertions(+), 6 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..19d396f7d90 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 cexp(int128) @ stub cexpf @ stub cexpl @ cdecl cimag(int128) diff --git a/dlls/msvcr120/tests/msvcr120.c b/dlls/msvcr120/tests/msvcr120.c index 21b1c00b4c3..2db3d4b0d2f 100644 --- a/dlls/msvcr120/tests/msvcr120.c +++ b/dlls/msvcr120/tests/msvcr120.c @@ -233,6 +233,7 @@ static double (__cdecl *p_creal)(_Dcomplex); static float (__cdecl *p_crealf)(_Fcomplex); static double (__cdecl *p_cimag)(_Dcomplex); static float (__cdecl *p_cimagf)(_Fcomplex); +static _Dcomplex (__cdecl *p_cexp)(_Dcomplex); static double (__cdecl *p_nexttoward)(double, double); static float (__cdecl *p_nexttowardf)(float, double); static double (__cdecl *p_nexttowardl)(double, double); @@ -299,6 +300,12 @@ static _Fcomplex __cdecl i386_FCbuild(float r, float i) } #endif +static inline BOOL almost_equal(double d1, double d2) { + if(d1-d2>-1e-15 && d1-d2<1e-15) + return TRUE; + return FALSE; +} + #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) @@ -357,6 +364,7 @@ static BOOL init(void) SET(p__Cbuild, "_Cbuild"); SET(p_creal, "creal"); SET(p_cimag, "cimag"); + SET(p_cexp, "cexp"); SET(p__FCbuild, "_FCbuild"); SET(p_crealf, "crealf"); SET(p_cimagf, "cimagf"); @@ -1956,6 +1964,21 @@ static void test__fsopen(void) p_setlocale(LC_ALL, "C"); } +static void test_cexp(void) +{ + _Dcomplex c, r; + + c = p__Cbuild(0.0, M_PI); + r = p_cexp(c); + ok(r.r == -1.0, "r.r = %lf\n", r.r); + ok(almost_equal(r.i, 0.0), "got %e\n", r.i); + + c = p__Cbuild(0.0, INFINITY); + r = p_cexp(c); + ok(isnan(r.r), "r.r = %lf\n", r.r); + ok(isnan(r.i), "r.i = %lf\n", r.i); +} + START_TEST(msvcr120) { if (!init()) return; @@ -1983,4 +2006,5 @@ START_TEST(msvcr120) test_strcmp(); test_gmtime64(); test__fsopen(); + test_cexp(); } diff --git a/dlls/msvcr120_app/msvcr120_app.spec b/dlls/msvcr120_app/msvcr120_app.spec index fc2ac655c40..22d79412696 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 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..4d0fb42ad32 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 cexp(int128) @ stub cexpf @ stub cexpl @ cdecl cimag(int128) diff --git a/include/msvcrt/complex.h b/include/msvcrt/complex.h index 53a72b6a000..f2680a258fa 100644 --- a/include/msvcrt/complex.h +++ b/include/msvcrt/complex.h @@ -24,8 +24,9 @@ typedef struct _C_float_complex typedef _C_double_complex _Dcomplex; typedef _C_float_complex _Fcomplex; -double __cdecl cimag(_Dcomplex); -double __cdecl creal(_Dcomplex); +_Dcomplex __cdecl cexp(_Dcomplex); +double __cdecl cimag(_Dcomplex); +double __cdecl creal(_Dcomplex); float __cdecl cimagf(_Fcomplex); float __cdecl crealf(_Fcomplex); diff --git a/libs/musl/Makefile.in b/libs/musl/Makefile.in index bcb14ffef78..2cfe9685334 100644 --- a/libs/musl/Makefile.in +++ b/libs/musl/Makefile.in @@ -1,8 +1,10 @@ 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(a)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..4c5d8e1b068 --- /dev/null +++ b/libs/musl/src/complex/cexp.c @@ -0,0 +1,83 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/s_cexp.c */ +/*- + * Copyright (c) 2011 David Schultz <das(a)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 */ + +_Dcomplex cexp(_Dcomplex z) +{ + 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) + return CMPLX(exp(x), y); + EXTRACT_WORDS(hx, lx, x); + /* cexp(0 + I y) = cos(y) + I sin(y) */ + if (((hx & 0x7fffffff) | lx) == 0) + return CMPLX(cos(y), sin(y)); + + if (hy >= 0x7ff00000) { + if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { + /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ + return CMPLX(y - y, y - y); + } else if (hx & 0x80000000) { + /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ + return CMPLX(0.0, 0.0); + } else { + /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ + return CMPLX(x, y - y); + } + } + + 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). + */ + return __ldexp_cexp(z, 0); + } 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) + */ + exp_x = exp(x); + 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..47382c0e73c --- /dev/null +++ b/libs/musl/src/internal/complex_impl.h @@ -0,0 +1,22 @@ +#ifndef _COMPLEX_IMPL_H +#define _COMPLEX_IMPL_H + +#include <complex.h> +#include "libm.h" + +#undef __CMPLX +#undef CMPLX +#undef CMPLXF +#undef CMPLXL + +#define __CMPLX(x, y, t) \ + ((union { t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z) + +#define CMPLX(x, y) __CMPLX(x, y, _Dcomplex) +#define CMPLXF(x, y) __CMPLX(x, y, _Fcomplex) +#define CMPLXL(x, y) __CMPLX(x, y, long double) + +hidden _Dcomplex __ldexp_cexp(_Dcomplex,int); +hidden _Fcomplex __ldexp_cexpf(_Fcomplex,int); + +#endif -- GitLab https://gitlab.winehq.org/wine/wine/-/merge_requests/7748
Piotr Caban (@piotr) commented about include/msvcrt/complex.h:
typedef _C_double_complex _Dcomplex; typedef _C_float_complex _Fcomplex;
+double __cdecl cimag(_Dcomplex);
_ACRTIMP double __cdecl cimag(_Dcomplex);
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100113
Piotr Caban (@piotr) commented about dlls/msvcr120/tests/msvcr120.c:
} #endif
+static inline BOOL almost_equal(double d1, double d2) { Please use something like `msvcrt/tests/string.c:compare_double()`instead.
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100114
Piotr Caban (@piotr) commented about dlls/msvcr120/tests/msvcr120.c:
p_setlocale(LC_ALL, "C"); }
+static void test_cexp(void) +{ + _Dcomplex c, r; + + c = p__Cbuild(0.0, M_PI); + r = p_cexp(c); + ok(r.r == -1.0, "r.r = %lf\n", r.r); + ok(almost_equal(r.i, 0.0), "got %e\n", r.i); + + c = p__Cbuild(0.0, INFINITY);
This call should not invoke matherr callback. It's usually not that important but NAN sign bit should also be different. -- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100115
Piotr Caban (@piotr) commented about dlls/msvcr120/msvcr120.spec:
@ stub ccosl @ cdecl ceil(double) @ cdecl -arch=!i386 ceilf(float) -@ stub cexp +@ cdecl cexp(int128)
@ cdecl -norelay cexp(int128)
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100116
Piotr Caban (@piotr) commented about libs/musl/Makefile.in:
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 Why do you need this change?
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100117
Piotr Caban (@piotr) commented about libs/musl/src/internal/complex_impl.h:
+#ifndef _COMPLEX_IMPL_H +#define _COMPLEX_IMPL_H + +#include <complex.h> +#include "libm.h" + +#undef __CMPLX +#undef CMPLX +#undef CMPLXF +#undef CMPLXL + +#define __CMPLX(x, y, t) \ + ((union { t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z) While this code does what you want it looks non portable / over complicated. At least __xy type should be changed to something like:
#define __CMPLX(x, y, t, t2) \
((union { t __z; t2 __xy[2]; }){.__xy = {(x),(y)}}.__z)
#define CMPLX(x, y) __CMPLX(x, y, _Dcomplex, double)
It could be simplified even further to: ```c #define __CMPLX(x, y, t) ((t){(x), (y)}) ``` I'm afraid that none of this code is portable (at least it will not compile with MSVC). -- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100118
On Mon Apr 7 10:02:30 2025 +0000, Piotr Caban wrote:
Why do you need this change? these warnings:
In file included from /home/me/stuff/wine/libs/musl/src/complex/__cexp.c:28:
In file included from /home/me/stuff/wine/libs/musl/src/internal/complex_impl.h:4:
/home/me/stuff/wine/include/msvcrt/complex.h:27:28: warning: incompatible redeclaration of library function 'cexp' [-Wincompatible-library-redeclaration]
27 | _ACRTIMP _Dcomplex __cdecl cexp(_Dcomplex);
| ^
/home/me/stuff/wine/include/msvcrt/complex.h:27:28: note: 'cexp' is a builtin with type '_Complex double (_Complex double)'
/home/me/stuff/wine/include/msvcrt/complex.h:28:25: warning: incompatible redeclaration of library function 'cimag' [-Wincompatible-library-redeclaration]
28 | _ACRTIMP double __cdecl cimag(_Dcomplex);
| ^
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100521
On Thu Apr 10 02:25:17 2025 +0000, Daniel Lehman wrote:
these warnings: ``` In file included from /home/me/stuff/wine/libs/musl/src/complex/__cexp.c:28: In file included from /home/me/stuff/wine/libs/musl/src/internal/complex_impl.h:4: /home/me/stuff/wine/include/msvcrt/complex.h:27:28: warning: incompatible redeclaration of library function 'cexp' [-Wincompatible-library-redeclaration] 27 | _ACRTIMP _Dcomplex __cdecl cexp(_Dcomplex); | ^ /home/me/stuff/wine/include/msvcrt/complex.h:27:28: note: 'cexp' is a builtin with type '_Complex double (_Complex double)' /home/me/stuff/wine/include/msvcrt/complex.h:28:25: warning: incompatible redeclaration of library function 'cimag' [-Wincompatible-library-redeclaration] 28 | _ACRTIMP double __cdecl cimag(_Dcomplex); | ^ ``` Try using -fno-builtin-cexp -fno-builtin-cimag instead. -fno-builtin is a quite big tool, affecting all kinds of stuff from isdigit to memcpy.
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100543
On Mon Apr 7 10:02:30 2025 +0000, Piotr Caban wrote:
While this code does what you want it looks non portable / over complicated. At least __xy type should be changed to something like: ```c #define __CMPLX(x, y, t, t2) \ ((union { t __z; t2 __xy[2]; }){.__xy = {(x),(y)}}.__z) #define CMPLX(x, y) __CMPLX(x, y, _Dcomplex, double) ``` It could be simplified even further to: ```c #define __CMPLX(x, y, t) ((t){(x), (y)}) ``` I'm afraid that none of this code is portable (at least it will not compile with MSVC). It's probably best to define the macros using `_Cbuild` and `_FCbuild` functions:
#define CMPLX(x, y) _Cbuild(x, y)
-- https://gitlab.winehq.org/wine/wine/-/merge_requests/7748#note_100545
participants (4)
-
Alfred Agrell (＠Alcaro) -
Daniel Lehman
-
Daniel Lehman (＠dlehman25)
-
Piotr Caban (＠piotr)