diff options
Diffstat (limited to 'usr/src/libm/src/complex/cpowf.c')
| -rw-r--r-- | usr/src/libm/src/complex/cpowf.c | 168 |
1 files changed, 168 insertions, 0 deletions
diff --git a/usr/src/libm/src/complex/cpowf.c b/usr/src/libm/src/complex/cpowf.c new file mode 100644 index 0000000..1063be0 --- /dev/null +++ b/usr/src/libm/src/complex/cpowf.c @@ -0,0 +1,168 @@ +/* + * CDDL HEADER START + * + * The contents of this file are subject to the terms of the + * Common Development and Distribution License (the "License"). + * You may not use this file except in compliance with the License. + * + * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE + * or http://www.opensolaris.org/os/licensing. + * See the License for the specific language governing permissions + * and limitations under the License. + * + * When distributing Covered Code, include this CDDL HEADER in each + * file and include the License file at usr/src/OPENSOLARIS.LICENSE. + * If applicable, add the following below this CDDL HEADER, with the + * fields enclosed by brackets "[]" replaced with your own identifying + * information: Portions Copyright [yyyy] [name of copyright owner] + * + * CDDL HEADER END + */ +/* + * Copyright 2005 Sun Microsystems, Inc. All rights reserved. + * Use is subject to license terms. + */ + +#pragma ident "@(#)cpowf.c 1.6 06/01/23 SMI" + +#pragma weak cpowf = __cpowf + +#include "libm.h" +#include "complex_wrapper.h" + +#if defined(__i386) && !defined(__amd64) +extern int __swapRP(int); +#endif + +static const double + dpi = 3.1415926535897931160E0, /* Hex 2^ 1 * 1.921FB54442D18 */ + dhalf = 0.5, + dsqrt2 = 1.41421356237309514547, /* 3FF6A09E 667F3BCD */ + dinvpi = 0.3183098861837906715377675; + +static const float one = 1.0F, zero = 0.0F; + +#define hiinf 0x7f800000 + +fcomplex +cpowf(fcomplex z, fcomplex w) { + fcomplex ans; + float x, y, u, v, t, c, s; + double dx, dy, du, dv, dt, dc, ds, dp, dq, dr; + int ix, iy, hx, hy, hv, hu, iu, iv, j; + + x = F_RE(z); + y = F_IM(z); + u = F_RE(w); + v = F_IM(w); + hx = THE_WORD(x); + hy = THE_WORD(y); + hu = THE_WORD(u); + hv = THE_WORD(v); + ix = hx & 0x7fffffff; + iy = hy & 0x7fffffff; + iu = hu & 0x7fffffff; + iv = hv & 0x7fffffff; + + j = 0; + if (iv == 0) { /* z**(real) */ + if (hu == 0x3f800000) { /* (anything) ** 1 is itself */ + F_RE(ans) = x; + F_IM(ans) = y; + } else if (iu == 0) { /* (anything) ** 0 is 1 */ + F_RE(ans) = one; + F_IM(ans) = zero; + } else if (iy == 0) { /* (real)**(real) */ + F_IM(ans) = zero; + if (hx < 0 && ix < hiinf && iu < hiinf) { + /* -x ** u is exp(i*pi*u)*pow(x,u) */ + t = powf(-x, u); + sincospif(u, &s, &c); + F_RE(ans) = (c == zero)? c: c * t; + F_IM(ans) = (s == zero)? s: s * t; + } else { + F_RE(ans) = powf(x, u); + } + } else if (ix == 0 || ix >= hiinf || iy >= hiinf) { + if (ix > hiinf || iy > hiinf || iu > hiinf) { + F_RE(ans) = F_IM(ans) = x + y + u; + } else { + v = fabsf(y); + if (ix != 0) + v += fabsf(x); + t = atan2pif(y, x); + sincospif(t * u, &s, &c); + F_RE(ans) = (c == zero)? c: c * v; + F_IM(ans) = (s == zero)? s: s * v; + } + } else if (ix == iy) { /* if |x| == |y| */ +#if defined(__i386) && !defined(__amd64) + int rp = __swapRP(fp_extended); +#endif + dx = (double)x; + du = (double)u; + dt = (hx >= 0)? 0.25 : 0.75; + if (hy < 0) + dt = -dt; + dr = pow(dsqrt2 * dx, du); + sincospi(dt * du, &ds, &dc); + F_RE(ans) = (float)(dr * dc); + F_IM(ans) = (float)(dr * ds); +#if defined(__i386) && !defined(__amd64) + if (rp != fp_extended) + (void) __swapRP(rp); +#endif + } else { + j = 1; + } + if (j == 0) + return (ans); + } + if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) { + /* + * non-zero imaginery part(s) with inf component(s) yields NaN + */ + t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v); + F_RE(ans) = F_IM(ans) = t - t; + } else { +#if defined(__i386) && !defined(__amd64) + int rp = __swapRP(fp_extended); +#endif + /* INDENT OFF */ + /* + * r = u*log(hypot(x,y))-v*atan2(y,x), + * q = u*atan2(y,x)+v*log(hypot(x,y)) + * or + * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x), + * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi + * ans = exp(r)*(cospi(q/pi) + i sinpi(q/pi)) + */ + /* INDENT ON */ + dx = (double)x; + dy = (double)y; + du = (double)u; + dv = (double)v; + if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */ + dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy); + else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */ + dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx); + else + dt = dhalf * log(dx * dx + dy * dy); + dp = atan2pi(dy, dx); + if (iv == 0) { /* dv = 0 */ + dr = exp(du * dt); + dq = du * dp; + } else { + dr = exp(du * dt - dv * dp * dpi); + dq = du * dp + dv * dt * dinvpi; + } + sincospi(dq, &ds, &dc); + F_RE(ans) = (float)(dr * dc); + F_IM(ans) = (float)(dr * ds); +#if defined(__i386) && !defined(__amd64) + if (rp != fp_extended) + (void) __swapRP(rp); +#endif + } + return (ans); +} |