cexp[f] patch

Steve Kargl sgk at troutmask.apl.washington.edu
Sat Feb 2 00:25:01 UTC 2019


The following patch utilizes sincos[f] in the computation
of cexp[f].  For 20 million random z=x+Iy drawn in the 
box defined by x,y in [0,MAX_EXP*LN2], this amounts to
a 8.7% and 11.4% speed improvement over computing sin[f]
and cos[f] individually.

Index: msun/src/s_cexp.c
===================================================================
--- msun/src/s_cexp.c	(revision 2135)
+++ msun/src/s_cexp.c	(working copy)
@@ -41,7 +41,7 @@
 double complex
 cexp(double complex z)
 {
-	double x, y, exp_x;
+	double c, exp_x, s, x, y;
 	uint32_t hx, hy, lx, ly;
 
 	x = creal(z);
@@ -55,8 +55,10 @@
 		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 (((hx & 0x7fffffff) | lx) == 0) {
+		sincos(y, &s, &c);
+		return (CMPLX(c, s));
+	}
 
 	if (hy >= 0x7ff00000) {
 		if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
@@ -86,6 +88,7 @@
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = exp(x);
-		return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
+		sincos(y, &s, &c);
+		return (CMPLX(exp_x * c, exp_x * s));
 	}
 }
Index: msun/src/s_cexpf.c
===================================================================
--- msun/src/s_cexpf.c	(revision 2135)
+++ msun/src/s_cexpf.c	(working copy)
@@ -41,7 +41,7 @@
 float complex
 cexpf(float complex z)
 {
-	float x, y, exp_x;
+	float c, exp_x, s, x, y;
 	uint32_t hx, hy;
 
 	x = crealf(z);
@@ -55,8 +55,10 @@
 		return (CMPLXF(expf(x), y));
 	GET_FLOAT_WORD(hx, x);
 	/* cexp(0 + I y) = cos(y) + I sin(y) */
-	if ((hx & 0x7fffffff) == 0)
-		return (CMPLXF(cosf(y), sinf(y)));
+	if ((hx & 0x7fffffff) == 0) {
+		sincosf(y, &s, &c);
+		return (CMPLXF(c, s));
+	}
 
 	if (hy >= 0x7f800000) {
 		if ((hx & 0x7fffffff) != 0x7f800000) {
@@ -86,6 +88,7 @@
 		 *  -  x = NaN (spurious inexact exception from y)
 		 */
 		exp_x = expf(x);
-		return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
+		sincosf(y, &s, &c);
+		return (CMPLXF(exp_x * c, exp_x * s));
 	}
 }

-- 
Steve


More information about the freebsd-numerics mailing list