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