1226597Sdas/*- 2226597Sdas * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> 3226597Sdas * All rights reserved. 4226597Sdas * 5226597Sdas * Redistribution and use in source and binary forms, with or without 6226597Sdas * modification, are permitted provided that the following conditions 7226597Sdas * are met: 8226597Sdas * 1. Redistributions of source code must retain the above copyright 9226597Sdas * notice, this list of conditions and the following disclaimer. 10226597Sdas * 2. Redistributions in binary form must reproduce the above copyright 11226597Sdas * notice, this list of conditions and the following disclaimer in the 12226597Sdas * documentation and/or other materials provided with the distribution. 13226597Sdas * 14226597Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15226597Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16226597Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17226597Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18226597Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19226597Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20226597Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21226597Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22226597Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23226597Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24226597Sdas * SUCH DAMAGE. 25226597Sdas */ 26226597Sdas 27226597Sdas#include <sys/cdefs.h> 28226597Sdas__FBSDID("$FreeBSD$"); 29226597Sdas 30226597Sdas#include <complex.h> 31226597Sdas 32226597Sdas#include "math.h" 33226597Sdas#include "math_private.h" 34226597Sdas 35226597Sdasstatic const uint32_t k = 1799; /* constant for reduction */ 36226597Sdasstatic const double kln2 = 1246.97177782734161156; /* k * ln2 */ 37226597Sdas 38226597Sdas/* 39226597Sdas * Compute exp(x), scaled to avoid spurious overflow. An exponent is 40226597Sdas * returned separately in 'expt'. 41226597Sdas * 42226597Sdas * Input: ln(DBL_MAX) <= x < ln(2 * DBL_MAX / DBL_MIN_DENORM) ~= 1454.91 43226597Sdas * Output: 2**1023 <= y < 2**1024 44226597Sdas */ 45226597Sdasstatic double 46226597Sdas__frexp_exp(double x, int *expt) 47226597Sdas{ 48226597Sdas double exp_x; 49226597Sdas uint32_t hx; 50226597Sdas 51226597Sdas /* 52226597Sdas * We use exp(x) = exp(x - kln2) * 2**k, carefully chosen to 53226597Sdas * minimize |exp(kln2) - 2**k|. We also scale the exponent of 54226597Sdas * exp_x to MAX_EXP so that the result can be multiplied by 55226597Sdas * a tiny number without losing accuracy due to denormalization. 56226597Sdas */ 57226597Sdas exp_x = exp(x - kln2); 58226597Sdas GET_HIGH_WORD(hx, exp_x); 59226597Sdas *expt = (hx >> 20) - (0x3ff + 1023) + k; 60226597Sdas SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 1023) << 20)); 61226597Sdas return (exp_x); 62226597Sdas} 63226597Sdas 64226597Sdas/* 65226597Sdas * __ldexp_exp(x, expt) and __ldexp_cexp(x, expt) compute exp(x) * 2**expt. 66226597Sdas * They are intended for large arguments (real part >= ln(DBL_MAX)) 67226597Sdas * where care is needed to avoid overflow. 68226597Sdas * 69226597Sdas * The present implementation is narrowly tailored for our hyperbolic and 70226597Sdas * exponential functions. We assume expt is small (0 or -1), and the caller 71226597Sdas * has filtered out very large x, for which overflow would be inevitable. 72226597Sdas */ 73226597Sdas 74226597Sdasdouble 75226597Sdas__ldexp_exp(double x, int expt) 76226597Sdas{ 77226597Sdas double exp_x, scale; 78226597Sdas int ex_expt; 79226597Sdas 80226597Sdas exp_x = __frexp_exp(x, &ex_expt); 81226597Sdas expt += ex_expt; 82226597Sdas INSERT_WORDS(scale, (0x3ff + expt) << 20, 0); 83226597Sdas return (exp_x * scale); 84226597Sdas} 85226597Sdas 86226597Sdasdouble complex 87226597Sdas__ldexp_cexp(double complex z, int expt) 88226597Sdas{ 89226597Sdas double x, y, exp_x, scale1, scale2; 90226597Sdas int ex_expt, half_expt; 91226597Sdas 92226597Sdas x = creal(z); 93226597Sdas y = cimag(z); 94226597Sdas exp_x = __frexp_exp(x, &ex_expt); 95226597Sdas expt += ex_expt; 96226597Sdas 97226597Sdas /* 98226597Sdas * Arrange so that scale1 * scale2 == 2**expt. We use this to 99226597Sdas * compensate for scalbn being horrendously slow. 100226597Sdas */ 101226597Sdas half_expt = expt / 2; 102226597Sdas INSERT_WORDS(scale1, (0x3ff + half_expt) << 20, 0); 103226597Sdas half_expt = expt - half_expt; 104226597Sdas INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); 105226597Sdas 106226597Sdas return (cpack(cos(y) * exp_x * scale1 * scale2, 107226597Sdas sin(y) * exp_x * scale1 * scale2)); 108226597Sdas} 109