k_exp.c revision 284810
144743Smarkm/*- 244743Smarkm * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> 344743Smarkm * All rights reserved. 444743Smarkm * 544743Smarkm * Redistribution and use in source and binary forms, with or without 644743Smarkm * modification, are permitted provided that the following conditions 744743Smarkm * are met: 844743Smarkm * 1. Redistributions of source code must retain the above copyright 944743Smarkm * notice, this list of conditions and the following disclaimer. 1044743Smarkm * 2. Redistributions in binary form must reproduce the above copyright 1144743Smarkm * notice, this list of conditions and the following disclaimer in the 1244743Smarkm * documentation and/or other materials provided with the distribution. 1344743Smarkm * 1444743Smarkm * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 1544743Smarkm * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 1644743Smarkm * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 1744743Smarkm * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 1844743Smarkm * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 1944743Smarkm * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 2044743Smarkm * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 2144743Smarkm * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22311814Sdim * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 2344743Smarkm * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 2444743Smarkm * SUCH DAMAGE. 2544743Smarkm */ 2644743Smarkm 2744743Smarkm#include <sys/cdefs.h> 2844743Smarkm__FBSDID("$FreeBSD: stable/10/lib/msun/src/k_exp.c 284810 2015-06-25 13:01:10Z tijl $"); 2944743Smarkm 3044743Smarkm#include <complex.h> 3144743Smarkm 3244743Smarkm#include "math.h" 3344743Smarkm#include "math_private.h" 3444743Smarkm 3544743Smarkmstatic const uint32_t k = 1799; /* constant for reduction */ 3644743Smarkmstatic const double kln2 = 1246.97177782734161156; /* k * ln2 */ 3744743Smarkm 3844743Smarkm/* 3944743Smarkm * Compute exp(x), scaled to avoid spurious overflow. An exponent is 4044743Smarkm * returned separately in 'expt'. 4144743Smarkm * 4244743Smarkm * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 4344743Smarkm * Output: 2**1023 <= y < 2**1024 4444743Smarkm */ 4544743Smarkmstatic double 4644743Smarkm__frexp_exp(double x, int *expt) 4744743Smarkm{ 4844743Smarkm double exp_x; 4944743Smarkm uint32_t hx; 5044743Smarkm 5144743Smarkm /* 5244743Smarkm * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to 5344743Smarkm * minimize |exp(kln2) - 2**k|. We also scale the exponent of 5444743Smarkm * exp_x to MAX_EXP so that the result can be multiplied by 5544743Smarkm * a tiny number without losing accuracy due to denormalization. 5644743Smarkm */ 5744743Smarkm exp_x = exp(x - kln2); 5844743Smarkm GET_HIGH_WORD(hx, exp_x); 5944743Smarkm *expt = (hx >> 20) - (0x3ff + 1023) + k; 6044743Smarkm SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); 6144743Smarkm return (exp_x); 6244743Smarkm} 6344743Smarkm 6444743Smarkm/* 6544743Smarkm * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. 6644743Smarkm * They are intended for large arguments (real part >= ln(DBL_MAX)) 6744743Smarkm * where care is needed to avoid overflow. 6844743Smarkm * 6944743Smarkm * The present implementation is narrowly tailored for our hyperbolic and 7044743Smarkm * exponential functions. We assume expt is small (0 or -1), and the caller 7144743Smarkm * has filtered out very large x, for which overflow would be inevitable. 7244743Smarkm */ 7344743Smarkm 7444743Smarkmdouble 7544743Smarkm__ldexp_exp(double x, int expt) 7644743Smarkm{ 7744743Smarkm double exp_x, scale; 7844743Smarkm int ex_expt; 7944743Smarkm 8044743Smarkm exp_x = __frexp_exp(x, &ex_expt); 8144743Smarkm expt += ex_expt; 8244743Smarkm INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); 8344743Smarkm return (exp_x * scale); 8444743Smarkm} 8544743Smarkm 8644743Smarkmdouble complex 8744743Smarkm__ldexp_cexp(double complex z, int expt) 88{ 89 double x, y, exp_x, scale1, scale2; 90 int ex_expt, half_expt; 91 92 x = creal(z); 93 y = cimag(z); 94 exp_x = __frexp_exp(x, &ex_expt); 95 expt += ex_expt; 96 97 /* 98 * Arrange so that scale1 * scale2 == 2**expt. We use this to 99 * compensate for scalbn being horrendously slow. 100 */ 101 half_expt = expt / 2; 102 INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); 103 half_expt = expt - half_expt; 104 INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); 105 106 return (CMPLX(cos(y) * exp_x * scale1 * scale2, 107 sin(y) * exp_x * scale1 * scale2)); 108} 109