Date: Sat, 19 Sep 2020 16:58:49 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: freebsd-hackers@freebsd.org, freebsd-current@freebsd.org Subject: [PATCH] libm -- fix expl kernels Message-ID: <20200919235849.GA47701@troutmask.apl.washington.edu>
next in thread | raw e-mail | index | archive | help
The following patch fixes the ld80 and ld128 kernels for expl(x).
I can't test ld128, so that patch is untested.
. Micro-optimization: use sincosl(x) instead of a call to cosl(x) and
a call to sinl(x). Argument reduction is done once not twice.
. Use a long double constant instead of an invalid double constant.
. Spell scale2 correctly
--- /usr/src/lib/msun/ld80/k_expl.h 2019-02-23 07:45:04.364785000 -0800
+++ ld80/k_expl.h 2019-12-26 08:15:29.366146000 -0800
@@ -277,7 +277,7 @@
static inline long double complex
__ldexp_cexpl(long double complex z, int expt)
{
- long double exp_x, hi, lo;
+ long double c, exp_x, hi, lo, s;
long double x, y, scale1, scale2;
int half_expt, k;
@@ -285,16 +285,17 @@
y = cimagl(z);
__k_expl(x, &hi, &lo, &k);
- exp_x = (lo + hi) * 0x1p16382;
+ exp_x = (lo + hi) * 0x1p16382L;
expt += k - 16382;
scale1 = 1;
half_expt = expt / 2;
SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
scale2 = 1;
- SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
+ SET_LDBL_EXPSIGN(scale2, BIAS + expt - half_expt);
- return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
- sinl(y) * exp_x * scale1 * scale2));
+ sincosl(y, &s, &c);
+ return (CMPLXL(c * exp_x * scale1 * scale2,
+ s * exp_x * scale1 * scale2));
}
#endif /* _COMPLEX_H */
--- /usr/src/lib/msun/ld128/k_expl.h 2020-03-07 09:48:15.719876000 -0800
+++ ld128/k_expl.h 2020-09-19 16:51:43.577195000 -0700
@@ -300,7 +300,7 @@
static inline long double complex
__ldexp_cexpl(long double complex z, int expt)
{
- long double exp_x, hi, lo;
+ long double c, exp_x, hi, lo, s;
long double x, y, scale1, scale2;
int half_expt, k;
@@ -308,16 +308,17 @@
y = cimagl(z);
__k_expl(x, &hi, &lo, &k);
- exp_x = (lo + hi) * 0x1p16382;
+ exp_x = (lo + hi) * 0x1p16382L;
expt += k - 16382;
scale1 = 1;
half_expt = expt / 2;
SET_LDBL_EXPSIGN(scale1, BIAS + half_expt);
scale2 = 1;
- SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
+ SET_LDBL_EXPSIGN(scale2, BIAS + expt - half_expt);
- return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
- sinl(y) * exp_x * scale1 * scale2));
+ sincosl(y, &s, &c);
+ return (CMPLXL(c * exp_x * scale1 * scale2,
+ s * exp_x * scale1 * scale2));
}
#endif /* _COMPLEX_H */
--
Steve
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20200919235849.GA47701>