diff options
Diffstat (limited to 'lib')
| -rw-r--r-- | lib/msun/Makefile | 2 | ||||
| -rw-r--r-- | lib/msun/src/k_exp.c | 108 | ||||
| -rw-r--r-- | lib/msun/src/k_expf.c | 87 | ||||
| -rw-r--r-- | lib/msun/src/math_private.h | 8 | ||||
| -rw-r--r-- | lib/msun/src/s_cexp.c | 19 | ||||
| -rw-r--r-- | lib/msun/src/s_cexpf.c | 19 |
6 files changed, 210 insertions, 33 deletions
diff --git a/lib/msun/Makefile b/lib/msun/Makefile index 7542f63..f1385f1 100644 --- a/lib/msun/Makefile +++ b/lib/msun/Makefile @@ -49,7 +49,7 @@ COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \ e_pow.c e_powf.c e_rem_pio2.c \ e_rem_pio2f.c e_remainder.c e_remainderf.c e_scalb.c e_scalbf.c \ e_sinh.c e_sinhf.c e_sqrt.c e_sqrtf.c fenv.c \ - k_cos.c k_cosf.c k_rem_pio2.c k_sin.c k_sinf.c \ + k_cos.c k_cosf.c k_exp.c k_expf.c k_rem_pio2.c k_sin.c k_sinf.c \ k_tan.c k_tanf.c \ s_asinh.c s_asinhf.c s_atan.c s_atanf.c s_carg.c s_cargf.c s_cargl.c \ s_cbrt.c s_cbrtf.c s_ceil.c s_ceilf.c \ diff --git a/lib/msun/src/k_exp.c b/lib/msun/src/k_exp.c new file mode 100644 index 0000000..f592f69 --- /dev/null +++ b/lib/msun/src/k_exp.c @@ -0,0 +1,108 @@ +/*- + * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <complex.h> + +#include "math.h" +#include "math_private.h" + +static const uint32_t k = 1799; /* constant for reduction */ +static const double kln2 = 1246.97177782734161156; /* k * ln2 */ + +/* + * Compute exp(x), scaled to avoid spurious overflow. An exponent is + * returned separately in 'expt'. + * + * 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; + + /* + * 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); +} + +/* + * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. + * They are intended for large arguments (real part >= ln(DBL_MAX)) + * where care is needed to avoid overflow. + * + * The present implementation is narrowly tailored for our hyperbolic and + * 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 +__ldexp_exp(double x, int expt) +{ + double exp_x, scale; + int ex_expt; + + exp_x = __frexp_exp(x, &ex_expt); + expt += ex_expt; + INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); + return (exp_x * scale); +} + +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; + + /* + * 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 (cpack(cos(y) * exp_x * scale1 * scale2, + sin(y) * exp_x * scale1 * scale2)); +} diff --git a/lib/msun/src/k_expf.c b/lib/msun/src/k_expf.c new file mode 100644 index 0000000..a860b9f --- /dev/null +++ b/lib/msun/src/k_expf.c @@ -0,0 +1,87 @@ +/*- + * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +#include <complex.h> + +#include "math.h" +#include "math_private.h" + +static const uint32_t k = 235; /* constant for reduction */ +static const float kln2 = 162.88958740F; /* k * ln2 */ + +/* + * See k_exp.c for details. + * + * Input: ln(FLT_MAX) <= x < ln(2 * FLT_MAX / FLT_MIN_DENORM) ~= 192.7 + * Output: 2**127 <= y < 2**128 + */ +static float +__frexp_expf(float x, int *expt) +{ + double exp_x; + uint32_t hx; + + exp_x = expf(x - kln2); + GET_FLOAT_WORD(hx, exp_x); + *expt = (hx >> 23) - (0x7f + 127) + k; + SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23)); + return (exp_x); +} + +float +__ldexp_expf(float x, int expt) +{ + float exp_x, scale; + int ex_expt; + + exp_x = __frexp_expf(x, &ex_expt); + expt += ex_expt; + SET_FLOAT_WORD(scale, (0x7f + expt) << 23); + return (exp_x * scale); +} + +float complex +__ldexp_cexpf(float complex z, int expt) +{ + float x, y, exp_x, scale1, scale2; + int ex_expt, half_expt; + + x = crealf(z); + y = cimagf(z); + exp_x = __frexp_expf(x, &ex_expt); + expt += ex_expt; + + half_expt = expt / 2; + SET_FLOAT_WORD(scale1, (0x7f + half_expt) << 23); + half_expt = expt - half_expt; + SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); + + return (cpackf(cosf(y) * exp_x * scale1 * scale2, + sinf(y) * exp_x * scale1 * scale2)); +} diff --git a/lib/msun/src/math_private.h b/lib/msun/src/math_private.h index 2f8a9db..79280e3 100644 --- a/lib/msun/src/math_private.h +++ b/lib/msun/src/math_private.h @@ -397,6 +397,10 @@ int __ieee754_rem_pio2(double,double*); double __kernel_sin(double,double,int); double __kernel_cos(double,double); double __kernel_tan(double,double,int); +double __ldexp_exp(double,int); +#ifdef _COMPLEX_H +double complex __ldexp_cexp(double complex,int); +#endif /* float precision kernel functions */ #ifdef INLINE_REM_PIO2F @@ -415,6 +419,10 @@ float __kernel_cosdf(double); __inline #endif float __kernel_tandf(double,int); +float __ldexp_expf(float,int); +#ifdef _COMPLEX_H +float complex __ldexp_cexpf(float complex,int); +#endif /* long double precision kernel functions */ long double __kernel_sinl(long double, long double, int); diff --git a/lib/msun/src/s_cexp.c b/lib/msun/src/s_cexp.c index b907e1d..abe178f 100644 --- a/lib/msun/src/s_cexp.c +++ b/lib/msun/src/s_cexp.c @@ -34,18 +34,13 @@ __FBSDID("$FreeBSD$"); static const uint32_t exp_ovfl = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ -cexp_ovfl = 0x4096b8e4, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ -k = 1799; /* constant for reduction */ - -static const double -kln2 = 1246.97177782734161156; /* k * ln2 */ +cexp_ovfl = 0x4096b8e4; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ double complex cexp(double complex z) { double x, y, exp_x; uint32_t hx, hy, lx, ly; - int scale; x = creal(z); y = cimag(z); @@ -77,17 +72,9 @@ cexp(double complex z) if (hx >= exp_ovfl && hx <= cexp_ovfl) { /* * x is between 709.7 and 1454.3, so we must scale to avoid - * overflow in 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 MANT_DIG to avoid loss of - * accuracy due to underflow if sin(y) is tiny. + * overflow in exp(x). */ - exp_x = exp(x - kln2); - GET_HIGH_WORD(hx, exp_x); - SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 52) << 20)); - scale = (hx >> 20) - (0x3ff + 52) + k; - return (cpack(scalbn(cos(y) * exp_x, scale), - scalbn(sin(y) * exp_x, scale))); + return (__ldexp_cexp(z, 0)); } else { /* * Cases covered here: diff --git a/lib/msun/src/s_cexpf.c b/lib/msun/src/s_cexpf.c index 08ec545..0e30d08 100644 --- a/lib/msun/src/s_cexpf.c +++ b/lib/msun/src/s_cexpf.c @@ -34,18 +34,13 @@ __FBSDID("$FreeBSD$"); static const uint32_t exp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */ -cexp_ovfl = 0x43400074, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ -k = 235; /* constant for reduction */ - -static const float -kln2 = 162.88958740f; /* k * ln2 */ +cexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ float complex cexpf(float complex z) { float x, y, exp_x; uint32_t hx, hy; - int scale; x = crealf(z); y = cimagf(z); @@ -77,17 +72,9 @@ cexpf(float complex z) if (hx >= exp_ovfl && hx <= cexp_ovfl) { /* * x is between 88.7 and 192, so we must scale to avoid - * overflow in expf(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 MANT_DIG to avoid loss of - * accuracy due to underflow if sin(y) is tiny. + * overflow in expf(x). */ - exp_x = expf(x - kln2); - GET_FLOAT_WORD(hx, exp_x); - SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 23) << 23)); - scale = (hx >> 23) - (0x7f + 23) + k; - return (cpackf(scalbnf(cosf(y) * exp_x, scale), - scalbnf(sinf(y) * exp_x, scale))); + return (__ldexp_cexpf(z, 0)); } else { /* * Cases covered here: |
