[PATCH] libm -- fix expl kernels

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

[PATCH] libm -- fix expl kernels

Steve Kargl
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
_______________________________________________
[hidden email] mailing list
https://lists.freebsd.org/mailman/listinfo/freebsd-current
To unsubscribe, send any mail to "[hidden email]"