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 = 0x40862e42, /* high bits of MAX_EXP * ln2 ~= 710 */ 37219359Sdascexp_ovfl = 0x4096b8e4, /* (MAX_EXP - MIN_DENORM_EXP) * ln2 */ 38219359Sdask = 1799; /* constant for reduction */ 39219359Sdas 40219359Sdasstatic const double 41219359Sdaskln2 = 1246.97177782734161156; /* k * ln2 */ 42219359Sdas 43219359Sdasdouble complex 44219359Sdascexp(double complex z) 45219359Sdas{ 46219359Sdas double x, y, exp_x; 47219359Sdas uint32_t hx, hy, lx, ly; 48219359Sdas int scale; 49219359Sdas 50219359Sdas x = creal(z); 51219359Sdas y = cimag(z); 52219359Sdas 53219359Sdas EXTRACT_WORDS(hy, ly, y); 54219359Sdas hy &= 0x7fffffff; 55219359Sdas 56219359Sdas /* cexp(x + I 0) = exp(x) + I 0 */ 57219359Sdas if ((hy | ly) == 0) 58219359Sdas return (cpack(exp(x), y)); 59219359Sdas if (hy >= 0x7ff00000) { 60219359Sdas EXTRACT_WORDS(hx, lx, x); 61219359Sdas if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { 62219359Sdas /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ 63219359Sdas return (cpack(y - y, y - y)); 64219359Sdas } else if (hx & 0x80000000) { 65219359Sdas /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ 66219359Sdas return (cpack(0.0, 0.0)); 67219359Sdas } else { 68219359Sdas /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ 69219359Sdas return (cpack(x, y - y)); 70219359Sdas } 71219359Sdas } 72219359Sdas 73219359Sdas GET_HIGH_WORD(hx, x); 74219359Sdas if (hx >= exp_ovfl && hx <= cexp_ovfl) { 75219359Sdas /* 76219359Sdas * x is between 709.7 and 1454.3, so we must scale to avoid 77219359Sdas * overflow in exp(x). We use exp(x) = exp(x - kln2) * 2**k, 78219359Sdas * carefully chosen to minimize |exp(kln2) - 2**k|. We also 79219359Sdas * scale the exponent of exp(x) to MANT_DIG to avoid loss of 80219359Sdas * accuracy due to underflow if sin(y) is tiny. 81219359Sdas */ 82219359Sdas exp_x = exp(x - kln2); 83219359Sdas GET_HIGH_WORD(hx, exp_x); 84219359Sdas SET_HIGH_WORD(exp_x, (hx & 0xfffff) | ((0x3ff + 52) << 20)); 85219359Sdas scale = (hx >> 20) - (0x3ff + 52) + k; 86219359Sdas return (cpack(scalbn(cos(y) * exp_x, scale), 87219359Sdas scalbn(sin(y) * exp_x, scale))); 88219359Sdas } else { 89219359Sdas /* 90219359Sdas * Cases covered here: 91219359Sdas * - x < exp_ovfl and exp(x) won't overflow (common case) 92219359Sdas * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 93219359Sdas * - x = +-Inf (generated by exp()) 94219359Sdas * - x = NaN (spurious inexact exception from y) 95219359Sdas */ 96219359Sdas exp_x = exp(x); 97219359Sdas return (cpack(exp_x * cos(y), exp_x * sin(y))); 98219359Sdas } 99219359Sdas} 100