1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak __cexp = cexp 31 32 /* INDENT OFF */ 33 /* 34 * dcomplex cexp(dcomplex z); 35 * 36 * x+iy x 37 * e = e (cos(y)+i*sin(y)) 38 * 39 * Over/underflow issue 40 * -------------------- 41 * exp(x) may be huge but cos(y) or sin(y) may be tiny. So we use 42 * function __k_cexp(x,&n) to return exp(x) = __k_cexp(x,&n)*2**n. 43 * Thus if exp(x+iy) = A + Bi and t = __k_cexp(x,&n), then 44 * A = t*cos(y)*2**n, B = t*sin(y)*2**n 45 * 46 * Purge off all exceptional arguments: 47 * (x,0) --> (exp(x),0) for all x, include inf and NaN 48 * (+inf, y) --> (+inf, NaN) for inf, nan 49 * (-inf, y) --> (+-0, +-0) for y = inf, nan 50 * (x,+-inf/NaN) --> (NaN,NaN) for finite x 51 * For all other cases, return 52 * (x,y) --> exp(x)*cos(y)+i*exp(x)*sin(y)) 53 * 54 * Algorithm for out of range x and finite y 55 * 1. compute exp(x) in factor form (t=__k_cexp(x,&n))*2**n 56 * 2. compute sincos(y,&s,&c) 57 * 3. compute t*s+i*(t*c), then scale back to 2**n and return. 58 */ 59 /* INDENT ON */ 60 61 #include "libm.h" /* exp/scalbn/sincos/__k_cexp */ 62 #include "complex_wrapper.h" 63 64 static const double zero = 0.0; 65 66 dcomplex 67 cexp(dcomplex z) { 68 dcomplex ans; 69 double x, y, t, c, s; 70 int n, ix, iy, hx, hy, lx, ly; 71 72 x = D_RE(z); 73 y = D_IM(z); 74 hx = HI_WORD(x); 75 lx = LO_WORD(x); 76 hy = HI_WORD(y); 77 ly = LO_WORD(y); 78 ix = hx & 0x7fffffff; 79 iy = hy & 0x7fffffff; 80 if ((iy | ly) == 0) { /* y = 0 */ 81 D_RE(ans) = exp(x); 82 D_IM(ans) = y; 83 } else if (ISINF(ix, lx)) { /* x is +-inf */ 84 if (hx < 0) { 85 if (iy >= 0x7ff00000) { 86 D_RE(ans) = zero; 87 D_IM(ans) = zero; 88 } else { 89 sincos(y, &s, &c); 90 D_RE(ans) = zero * c; 91 D_IM(ans) = zero * s; 92 } 93 } else { 94 if (iy >= 0x7ff00000) { 95 D_RE(ans) = x; 96 D_IM(ans) = y - y; 97 } else { 98 (void) sincos(y, &s, &c); 99 D_RE(ans) = x * c; 100 D_IM(ans) = x * s; 101 } 102 } 103 } else { 104 (void) sincos(y, &s, &c); 105 if (ix >= 0x40862E42) { /* |x| > 709.78... ~ log(2**1024) */ 106 t = __k_cexp(x, &n); 107 D_RE(ans) = scalbn(t * c, n); 108 D_IM(ans) = scalbn(t * s, n); 109 } else { 110 t = exp(x); 111 D_RE(ans) = t * c; 112 D_IM(ans) = t * s; 113 } 114 } 115 return (ans); 116 } 117