[PATCH] libm -- fix expl kernels

From: Steve Kargl <sgk_at_troutmask.apl.washington.edu>
Date: Sat, 19 Sep 2020 16:58:49 -0700
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
_at__at_ -277,7 +277,7 _at__at_
 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;
 
_at__at_ -285,16 +285,17 _at__at_
 	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
_at__at_ -300,7 +300,7 _at__at_
 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;
 
_at__at_ -308,16 +308,17 _at__at_
 	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
Received on Sat Sep 19 2020 - 21:58:59 UTC

This archive was generated by hypermail 2.4.0 : Wed May 19 2021 - 11:41:25 UTC