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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23 */ 24 /* 25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26 * Use is subject to license terms. 27 */ 28 29 #pragma weak __cpowf = cpowf 30 31 #include "libm.h" 32 #include "complex_wrapper.h" 33 34 extern void sincospi(double, double *, double *); 35 extern void sincospif(float, float *, float *); 36 extern double atan2pi(double, double); 37 extern float atan2pif(float, float); 38 39 #if defined(__i386) && !defined(__amd64) 40 extern int __swapRP(int); 41 #endif 42 43 static const double 44 dpi = 3.1415926535897931160E0, /* Hex 2^ 1 * 1.921FB54442D18 */ 45 dhalf = 0.5, 46 dsqrt2 = 1.41421356237309514547, /* 3FF6A09E 667F3BCD */ 47 dinvpi = 0.3183098861837906715377675; 48 49 static const float one = 1.0F, zero = 0.0F; 50 51 #define hiinf 0x7f800000 52 53 fcomplex 54 cpowf(fcomplex z, fcomplex w) { 55 fcomplex ans; 56 float x, y, u, v, t, c, s; 57 double dx, dy, du, dv, dt, dc, ds, dp, dq, dr; 58 int ix, iy, hx, hy, hv, hu, iu, iv, j; 59 60 x = F_RE(z); 61 y = F_IM(z); 62 u = F_RE(w); 63 v = F_IM(w); 64 hx = THE_WORD(x); 65 hy = THE_WORD(y); 66 hu = THE_WORD(u); 67 hv = THE_WORD(v); 68 ix = hx & 0x7fffffff; 69 iy = hy & 0x7fffffff; 70 iu = hu & 0x7fffffff; 71 iv = hv & 0x7fffffff; 72 73 j = 0; 74 if (iv == 0) { /* z**(real) */ 75 if (hu == 0x3f800000) { /* (anything) ** 1 is itself */ 76 F_RE(ans) = x; 77 F_IM(ans) = y; 78 } else if (iu == 0) { /* (anything) ** 0 is 1 */ 79 F_RE(ans) = one; 80 F_IM(ans) = zero; 81 } else if (iy == 0) { /* (real)**(real) */ 82 F_IM(ans) = zero; 83 if (hx < 0 && ix < hiinf && iu < hiinf) { 84 /* -x ** u is exp(i*pi*u)*pow(x,u) */ 85 t = powf(-x, u); 86 sincospif(u, &s, &c); 87 F_RE(ans) = (c == zero)? c: c * t; 88 F_IM(ans) = (s == zero)? s: s * t; 89 } else { 90 F_RE(ans) = powf(x, u); 91 } 92 } else if (ix == 0 || ix >= hiinf || iy >= hiinf) { 93 if (ix > hiinf || iy > hiinf || iu > hiinf) { 94 F_RE(ans) = F_IM(ans) = x + y + u; 95 } else { 96 v = fabsf(y); 97 if (ix != 0) 98 v += fabsf(x); 99 t = atan2pif(y, x); 100 sincospif(t * u, &s, &c); 101 F_RE(ans) = (c == zero)? c: c * v; 102 F_IM(ans) = (s == zero)? s: s * v; 103 } 104 } else if (ix == iy) { /* if |x| == |y| */ 105 #if defined(__i386) && !defined(__amd64) 106 int rp = __swapRP(fp_extended); 107 #endif 108 dx = (double)x; 109 du = (double)u; 110 dt = (hx >= 0)? 0.25 : 0.75; 111 if (hy < 0) 112 dt = -dt; 113 dr = pow(dsqrt2 * dx, du); 114 sincospi(dt * du, &ds, &dc); 115 F_RE(ans) = (float)(dr * dc); 116 F_IM(ans) = (float)(dr * ds); 117 #if defined(__i386) && !defined(__amd64) 118 if (rp != fp_extended) 119 (void) __swapRP(rp); 120 #endif 121 } else { 122 j = 1; 123 } 124 if (j == 0) 125 return (ans); 126 } 127 if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) { 128 /* 129 * non-zero imaginery part(s) with inf component(s) yields NaN 130 */ 131 t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v); 132 F_RE(ans) = F_IM(ans) = t - t; 133 } else { 134 #if defined(__i386) && !defined(__amd64) 135 int rp = __swapRP(fp_extended); 136 #endif 137 /* INDENT OFF */ 138 /* 139 * r = u*log(hypot(x,y))-v*atan2(y,x), 140 * q = u*atan2(y,x)+v*log(hypot(x,y)) 141 * or 142 * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x), 143 * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi 144 * ans = exp(r)*(cospi(q/pi) + i sinpi(q/pi)) 145 */ 146 /* INDENT ON */ 147 dx = (double)x; 148 dy = (double)y; 149 du = (double)u; 150 dv = (double)v; 151 if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */ 152 dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy); 153 else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */ 154 dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx); 155 else 156 dt = dhalf * log(dx * dx + dy * dy); 157 dp = atan2pi(dy, dx); 158 if (iv == 0) { /* dv = 0 */ 159 dr = exp(du * dt); 160 dq = du * dp; 161 } else { 162 dr = exp(du * dt - dv * dp * dpi); 163 dq = du * dp + dv * dt * dinvpi; 164 } 165 sincospi(dq, &ds, &dc); 166 F_RE(ans) = (float)(dr * dc); 167 F_IM(ans) = (float)(dr * ds); 168 #if defined(__i386) && !defined(__amd64) 169 if (rp != fp_extended) 170 (void) __swapRP(rp); 171 #endif 172 } 173 return (ans); 174 } 175