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$"); 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 */ 37226597Sdascexp_ovfl = 0x43400074; /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ 38219359Sdas 39219359Sdasfloat complex 40219359Sdascexpf(float complex z) 41219359Sdas{ 42219359Sdas float x, y, exp_x; 43219359Sdas uint32_t hx, hy; 44219359Sdas 45219359Sdas x = crealf(z); 46219359Sdas y = cimagf(z); 47219359Sdas 48219359Sdas GET_FLOAT_WORD(hy, y); 49219359Sdas hy &= 0x7fffffff; 50219359Sdas 51219359Sdas /* cexp(x + I 0) = exp(x) + I 0 */ 52219359Sdas if (hy == 0) 53219359Sdas return (cpackf(expf(x), y)); 54219359Sdas GET_FLOAT_WORD(hx, x); 55226413Sdas /* cexp(0 + I y) = cos(y) + I sin(y) */ 56226413Sdas if ((hx & 0x7fffffff) == 0) 57226413Sdas return (cpackf(cosf(y), sinf(y))); 58226413Sdas 59219359Sdas if (hy >= 0x7f800000) { 60219359Sdas if ((hx & 0x7fffffff) != 0x7f800000) { 61219359Sdas /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ 62219359Sdas return (cpackf(y - y, y - y)); 63219359Sdas } else if (hx & 0x80000000) { 64219359Sdas /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ 65219359Sdas return (cpackf(0.0, 0.0)); 66219359Sdas } else { 67219359Sdas /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ 68219359Sdas return (cpackf(x, y - y)); 69219359Sdas } 70219359Sdas } 71219359Sdas 72219359Sdas if (hx >= exp_ovfl && hx <= cexp_ovfl) { 73219359Sdas /* 74219359Sdas * x is between 88.7 and 192, so we must scale to avoid 75226597Sdas * overflow in expf(x). 76219359Sdas */ 77226597Sdas return (__ldexp_cexpf(z, 0)); 78219359Sdas } else { 79219359Sdas /* 80219359Sdas * Cases covered here: 81219359Sdas * - x < exp_ovfl and exp(x) won't overflow (common case) 82219359Sdas * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 83219359Sdas * - x = +-Inf (generated by exp()) 84219359Sdas * - x = NaN (spurious inexact exception from y) 85219359Sdas */ 86219359Sdas exp_x = expf(x); 87219359Sdas return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); 88219359Sdas } 89219359Sdas} 90