s_cexpf.c revision 226413
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 226413 2011-10-16 05:37:01Z 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); 60226413Sdas /* cexp(0 + I y) = cos(y) + I sin(y) */ 61226413Sdas if ((hx & 0x7fffffff) == 0) 62226413Sdas return (cpackf(cosf(y), sinf(y))); 63226413Sdas 64219359Sdas if (hy >= 0x7f800000) { 65219359Sdas if ((hx & 0x7fffffff) != 0x7f800000) { 66219359Sdas /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ 67219359Sdas return (cpackf(y - y, y - y)); 68219359Sdas } else if (hx & 0x80000000) { 69219359Sdas /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ 70219359Sdas return (cpackf(0.0, 0.0)); 71219359Sdas } else { 72219359Sdas /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ 73219359Sdas return (cpackf(x, y - y)); 74219359Sdas } 75219359Sdas } 76219359Sdas 77219359Sdas if (hx >= exp_ovfl && hx <= cexp_ovfl) { 78219359Sdas /* 79219359Sdas * x is between 88.7 and 192, so we must scale to avoid 80219359Sdas * overflow in expf(x). We use exp(x) = exp(x - kln2) * 2**k, 81219359Sdas * carefully chosen to minimize |exp(kln2) - 2**k|. We also 82219359Sdas * scale the exponent of exp(x) to MANT_DIG to avoid loss of 83219359Sdas * accuracy due to underflow if sin(y) is tiny. 84219359Sdas */ 85219359Sdas exp_x = expf(x - kln2); 86219359Sdas GET_FLOAT_WORD(hx, exp_x); 87219359Sdas SET_FLOAT_WORD(exp_x, (hx & 0x7fffff) | ((0x7f + 23) << 23)); 88219359Sdas scale = (hx >> 23) - (0x7f + 23) + k; 89219359Sdas return (cpackf(scalbnf(cosf(y) * exp_x, scale), 90219359Sdas scalbnf(sinf(y) * exp_x, scale))); 91219359Sdas } else { 92219359Sdas /* 93219359Sdas * Cases covered here: 94219359Sdas * - x < exp_ovfl and exp(x) won't overflow (common case) 95219359Sdas * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 96219359Sdas * - x = +-Inf (generated by exp()) 97219359Sdas * - x = NaN (spurious inexact exception from y) 98219359Sdas */ 99219359Sdas exp_x = expf(x); 100219359Sdas return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); 101219359Sdas } 102219359Sdas} 103