Index: fusl/src/complex/__cexp.c |
diff --git a/fusl/src/complex/__cexp.c b/fusl/src/complex/__cexp.c |
index 05ac28c75c6271219bc66d3c1b27a8750161b1bd..166ff7c2f991bd2349655beb27ec8e47e0a1d9d4 100644 |
--- a/fusl/src/complex/__cexp.c |
+++ b/fusl/src/complex/__cexp.c |
@@ -27,7 +27,7 @@ |
#include "libm.h" |
-static const uint32_t k = 1799; /* constant for reduction */ |
+static const uint32_t k = 1799; /* constant for reduction */ |
static const double kln2 = 1246.97177782734161156; /* k * ln2 */ |
/* |
@@ -37,22 +37,21 @@ static const double kln2 = 1246.97177782734161156; /* k * ln2 */ |
* 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; |
+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; |
+ /* |
+ * 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; |
} |
/* |
@@ -64,24 +63,24 @@ static double __frexp_exp(double x, int *expt) |
* 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. |
*/ |
-double complex __ldexp_cexp(double complex z, int expt) |
-{ |
- double x, y, exp_x, scale1, scale2; |
- int ex_expt, half_expt; |
+double complex __ldexp_cexp(double complex 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; |
+ 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); |
+ /* |
+ * 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); |
+ return CMPLX(cos(y) * exp_x * scale1 * scale2, |
+ sin(y) * exp_x * scale1 * scale2); |
} |