s_cexpf.c revision 219359
1219359Sdas/*- 2219359Sdas * Copyright (c) 2011 David Schultz <das@FreeBSD.ORG> 3219359Sdas * All rights reserved. 4219359Sdas * 5219359Sdas * Redistribution and use in source and binary forms, with or without 6219359Sdas * modification, are permitted provided that the following conditions 7219359Sdas * are met: 8219359Sdas * 1. Redistributions of source code must retain the above copyright 9219359Sdas * notice, this list of conditions and the following disclaimer. 10219359Sdas * 2. Redistributions in binary form must reproduce the above copyright 11219359Sdas * notice, this list of conditions and the following disclaimer in the 12219359Sdas * documentation and/or other materials provided with the distribution. 13219359Sdas * 14219359Sdas * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 15219359Sdas * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 16219359Sdas * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 17219359Sdas * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 18219359Sdas * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 19219359Sdas * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 20219359Sdas * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 21219359Sdas * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22219359Sdas * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23219359Sdas * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24219359Sdas * SUCH DAMAGE. 25219359Sdas */ 26219359Sdas 27219359Sdas#include <sys/cdefs.h> 28219359Sdas__FBSDID("$FreeBSD: head/lib/msun/src/s_cexpf.c 219359 2011-03-07 03:09:24Z das $"); 29219359Sdas 30219359Sdas#include <complex.h> 31219359Sdas#include <math.h> 32219359Sdas 33219359Sdas#include "math_private.h" 34219359Sdas 35219359Sdasstatic const uint32_t 36219359Sdasexp_ovfl = 0x42b17218, /* MAX_EXP * ln2 ~= 88.722839355 */ 37219359Sdascexp_ovfl = 0x43400074, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ 38219359Sdask = 235; /* constant for reduction */ 39219359Sdas 40219359Sdasstatic const float 41219359Sdaskln2 = 162.88958740f; /* k * ln2 */ 42219359Sdas 43219359Sdasfloat complex 44219359Sdascexpf(float complex z) 45219359Sdas{ 46219359Sdas float x, y, exp_x; 47219359Sdas uint32_t hx, hy; 48219359Sdas int scale; 49219359Sdas 50219359Sdas x = crealf(z); 51219359Sdas y = cimagf(z); 52219359Sdas 53219359Sdas GET_FLOAT_WORD(hy, y); 54219359Sdas hy &= 0x7fffffff; 55219359Sdas 56219359Sdas /* cexp(x + I 0) = exp(x) + I 0 */ 57219359Sdas if (hy == 0) 58219359Sdas return (cpackf(expf(x), y)); 59219359Sdas GET_FLOAT_WORD(hx, x); 60219359Sdas if (hy >= 0x7f800000) { 61219359Sdas if ((hx & 0x7fffffff) != 0x7f800000) { 62219359Sdas /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ 63219359Sdas return (cpackf(y - y, y - y)); 64219359Sdas } else if (hx & 0x80000000) { 65219359Sdas /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ 66219359Sdas return (cpackf(0.0, 0.0)); 67219359Sdas } else { 68219359Sdas /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ 69219359Sdas return (cpackf(x, y - y)); 70219359Sdas } 71219359Sdas } 72219359Sdas 73219359Sdas if (hx >= exp_ovfl && hx <= cexp_ovfl) { 74219359Sdas /* 75219359Sdas * x is between 88.7 and 192, so we must scale to avoid 76219359Sdas * overflow in expf(x). We use exp(x) = exp(x - kln2) * 2**k, 77219359Sdas * carefully chosen to minimize |exp(kln2) - 2**k|. We also 78219359Sdas * scale the exponent of exp(x) to MANT_DIG to avoid loss of 79219359Sdas * accuracy due to underflow if sin(y) is tiny. 80219359Sdas */ 81219359Sdas exp_x = expf(x - kln2); 82219359Sdas GET_FLOAT_WORD(hx, exp_x); 83219359Sdas SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 23) << 23)); 84219359Sdas scale = (hx >> 23) - (0x7f + 23) + k; 85219359Sdas return (cpackf(scalbnf(cosf(y) * exp_x, scale), 86219359Sdas scalbnf(sinf(y) * exp_x, scale))); 87219359Sdas } else { 88219359Sdas /* 89219359Sdas * Cases covered here: 90219359Sdas * - x < exp_ovfl and exp(x) won't overflow (common case) 91219359Sdas * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 92219359Sdas * - x = +-Inf (generated by exp()) 93219359Sdas * - x = NaN (spurious inexact exception from y) 94219359Sdas */ 95219359Sdas exp_x = expf(x); 96219359Sdas return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); 97219359Sdas } 98219359Sdas} 99