125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis /*
2225c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2325c28e83SPiotr Jasiukajtis */
2425c28e83SPiotr Jasiukajtis /*
2525c28e83SPiotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
2625c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2725c28e83SPiotr Jasiukajtis */
2825c28e83SPiotr Jasiukajtis
29*ddc0e0b5SRichard Lowe #pragma weak __cpowf = cpowf
3025c28e83SPiotr Jasiukajtis
3125c28e83SPiotr Jasiukajtis #include "libm.h"
3225c28e83SPiotr Jasiukajtis #include "complex_wrapper.h"
3325c28e83SPiotr Jasiukajtis
3425c28e83SPiotr Jasiukajtis extern void sincospi(double, double *, double *);
3525c28e83SPiotr Jasiukajtis extern void sincospif(float, float *, float *);
3625c28e83SPiotr Jasiukajtis extern double atan2pi(double, double);
3725c28e83SPiotr Jasiukajtis extern float atan2pif(float, float);
3825c28e83SPiotr Jasiukajtis
3925c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
4025c28e83SPiotr Jasiukajtis extern int __swapRP(int);
4125c28e83SPiotr Jasiukajtis #endif
4225c28e83SPiotr Jasiukajtis
4325c28e83SPiotr Jasiukajtis static const double
4425c28e83SPiotr Jasiukajtis dpi = 3.1415926535897931160E0, /* Hex 2^ 1 * 1.921FB54442D18 */
4525c28e83SPiotr Jasiukajtis dhalf = 0.5,
4625c28e83SPiotr Jasiukajtis dsqrt2 = 1.41421356237309514547, /* 3FF6A09E 667F3BCD */
4725c28e83SPiotr Jasiukajtis dinvpi = 0.3183098861837906715377675;
4825c28e83SPiotr Jasiukajtis
4925c28e83SPiotr Jasiukajtis static const float one = 1.0F, zero = 0.0F;
5025c28e83SPiotr Jasiukajtis
5125c28e83SPiotr Jasiukajtis #define hiinf 0x7f800000
5225c28e83SPiotr Jasiukajtis
5325c28e83SPiotr Jasiukajtis fcomplex
cpowf(fcomplex z,fcomplex w)5425c28e83SPiotr Jasiukajtis cpowf(fcomplex z, fcomplex w) {
5525c28e83SPiotr Jasiukajtis fcomplex ans;
5625c28e83SPiotr Jasiukajtis float x, y, u, v, t, c, s;
5725c28e83SPiotr Jasiukajtis double dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
5825c28e83SPiotr Jasiukajtis int ix, iy, hx, hy, hv, hu, iu, iv, j;
5925c28e83SPiotr Jasiukajtis
6025c28e83SPiotr Jasiukajtis x = F_RE(z);
6125c28e83SPiotr Jasiukajtis y = F_IM(z);
6225c28e83SPiotr Jasiukajtis u = F_RE(w);
6325c28e83SPiotr Jasiukajtis v = F_IM(w);
6425c28e83SPiotr Jasiukajtis hx = THE_WORD(x);
6525c28e83SPiotr Jasiukajtis hy = THE_WORD(y);
6625c28e83SPiotr Jasiukajtis hu = THE_WORD(u);
6725c28e83SPiotr Jasiukajtis hv = THE_WORD(v);
6825c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff;
6925c28e83SPiotr Jasiukajtis iy = hy & 0x7fffffff;
7025c28e83SPiotr Jasiukajtis iu = hu & 0x7fffffff;
7125c28e83SPiotr Jasiukajtis iv = hv & 0x7fffffff;
7225c28e83SPiotr Jasiukajtis
7325c28e83SPiotr Jasiukajtis j = 0;
7425c28e83SPiotr Jasiukajtis if (iv == 0) { /* z**(real) */
7525c28e83SPiotr Jasiukajtis if (hu == 0x3f800000) { /* (anything) ** 1 is itself */
7625c28e83SPiotr Jasiukajtis F_RE(ans) = x;
7725c28e83SPiotr Jasiukajtis F_IM(ans) = y;
7825c28e83SPiotr Jasiukajtis } else if (iu == 0) { /* (anything) ** 0 is 1 */
7925c28e83SPiotr Jasiukajtis F_RE(ans) = one;
8025c28e83SPiotr Jasiukajtis F_IM(ans) = zero;
8125c28e83SPiotr Jasiukajtis } else if (iy == 0) { /* (real)**(real) */
8225c28e83SPiotr Jasiukajtis F_IM(ans) = zero;
8325c28e83SPiotr Jasiukajtis if (hx < 0 && ix < hiinf && iu < hiinf) {
8425c28e83SPiotr Jasiukajtis /* -x ** u is exp(i*pi*u)*pow(x,u) */
8525c28e83SPiotr Jasiukajtis t = powf(-x, u);
8625c28e83SPiotr Jasiukajtis sincospif(u, &s, &c);
8725c28e83SPiotr Jasiukajtis F_RE(ans) = (c == zero)? c: c * t;
8825c28e83SPiotr Jasiukajtis F_IM(ans) = (s == zero)? s: s * t;
8925c28e83SPiotr Jasiukajtis } else {
9025c28e83SPiotr Jasiukajtis F_RE(ans) = powf(x, u);
9125c28e83SPiotr Jasiukajtis }
9225c28e83SPiotr Jasiukajtis } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
9325c28e83SPiotr Jasiukajtis if (ix > hiinf || iy > hiinf || iu > hiinf) {
9425c28e83SPiotr Jasiukajtis F_RE(ans) = F_IM(ans) = x + y + u;
9525c28e83SPiotr Jasiukajtis } else {
9625c28e83SPiotr Jasiukajtis v = fabsf(y);
9725c28e83SPiotr Jasiukajtis if (ix != 0)
9825c28e83SPiotr Jasiukajtis v += fabsf(x);
9925c28e83SPiotr Jasiukajtis t = atan2pif(y, x);
10025c28e83SPiotr Jasiukajtis sincospif(t * u, &s, &c);
10125c28e83SPiotr Jasiukajtis F_RE(ans) = (c == zero)? c: c * v;
10225c28e83SPiotr Jasiukajtis F_IM(ans) = (s == zero)? s: s * v;
10325c28e83SPiotr Jasiukajtis }
10425c28e83SPiotr Jasiukajtis } else if (ix == iy) { /* if |x| == |y| */
10525c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
10625c28e83SPiotr Jasiukajtis int rp = __swapRP(fp_extended);
10725c28e83SPiotr Jasiukajtis #endif
10825c28e83SPiotr Jasiukajtis dx = (double)x;
10925c28e83SPiotr Jasiukajtis du = (double)u;
11025c28e83SPiotr Jasiukajtis dt = (hx >= 0)? 0.25 : 0.75;
11125c28e83SPiotr Jasiukajtis if (hy < 0)
11225c28e83SPiotr Jasiukajtis dt = -dt;
11325c28e83SPiotr Jasiukajtis dr = pow(dsqrt2 * dx, du);
11425c28e83SPiotr Jasiukajtis sincospi(dt * du, &ds, &dc);
11525c28e83SPiotr Jasiukajtis F_RE(ans) = (float)(dr * dc);
11625c28e83SPiotr Jasiukajtis F_IM(ans) = (float)(dr * ds);
11725c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
11825c28e83SPiotr Jasiukajtis if (rp != fp_extended)
11925c28e83SPiotr Jasiukajtis (void) __swapRP(rp);
12025c28e83SPiotr Jasiukajtis #endif
12125c28e83SPiotr Jasiukajtis } else {
12225c28e83SPiotr Jasiukajtis j = 1;
12325c28e83SPiotr Jasiukajtis }
12425c28e83SPiotr Jasiukajtis if (j == 0)
12525c28e83SPiotr Jasiukajtis return (ans);
12625c28e83SPiotr Jasiukajtis }
12725c28e83SPiotr Jasiukajtis if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
12825c28e83SPiotr Jasiukajtis /*
12925c28e83SPiotr Jasiukajtis * non-zero imaginery part(s) with inf component(s) yields NaN
13025c28e83SPiotr Jasiukajtis */
13125c28e83SPiotr Jasiukajtis t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
13225c28e83SPiotr Jasiukajtis F_RE(ans) = F_IM(ans) = t - t;
13325c28e83SPiotr Jasiukajtis } else {
13425c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
13525c28e83SPiotr Jasiukajtis int rp = __swapRP(fp_extended);
13625c28e83SPiotr Jasiukajtis #endif
13725c28e83SPiotr Jasiukajtis /* INDENT OFF */
13825c28e83SPiotr Jasiukajtis /*
13925c28e83SPiotr Jasiukajtis * r = u*log(hypot(x,y))-v*atan2(y,x),
14025c28e83SPiotr Jasiukajtis * q = u*atan2(y,x)+v*log(hypot(x,y))
14125c28e83SPiotr Jasiukajtis * or
14225c28e83SPiotr Jasiukajtis * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
14325c28e83SPiotr Jasiukajtis * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
14425c28e83SPiotr Jasiukajtis * ans = exp(r)*(cospi(q/pi) + i sinpi(q/pi))
14525c28e83SPiotr Jasiukajtis */
14625c28e83SPiotr Jasiukajtis /* INDENT ON */
14725c28e83SPiotr Jasiukajtis dx = (double)x;
14825c28e83SPiotr Jasiukajtis dy = (double)y;
14925c28e83SPiotr Jasiukajtis du = (double)u;
15025c28e83SPiotr Jasiukajtis dv = (double)v;
15125c28e83SPiotr Jasiukajtis if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */
15225c28e83SPiotr Jasiukajtis dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
15325c28e83SPiotr Jasiukajtis else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
15425c28e83SPiotr Jasiukajtis dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);
15525c28e83SPiotr Jasiukajtis else
15625c28e83SPiotr Jasiukajtis dt = dhalf * log(dx * dx + dy * dy);
15725c28e83SPiotr Jasiukajtis dp = atan2pi(dy, dx);
15825c28e83SPiotr Jasiukajtis if (iv == 0) { /* dv = 0 */
15925c28e83SPiotr Jasiukajtis dr = exp(du * dt);
16025c28e83SPiotr Jasiukajtis dq = du * dp;
16125c28e83SPiotr Jasiukajtis } else {
16225c28e83SPiotr Jasiukajtis dr = exp(du * dt - dv * dp * dpi);
16325c28e83SPiotr Jasiukajtis dq = du * dp + dv * dt * dinvpi;
16425c28e83SPiotr Jasiukajtis }
16525c28e83SPiotr Jasiukajtis sincospi(dq, &ds, &dc);
16625c28e83SPiotr Jasiukajtis F_RE(ans) = (float)(dr * dc);
16725c28e83SPiotr Jasiukajtis F_IM(ans) = (float)(dr * ds);
16825c28e83SPiotr Jasiukajtis #if defined(__i386) && !defined(__amd64)
16925c28e83SPiotr Jasiukajtis if (rp != fp_extended)
17025c28e83SPiotr Jasiukajtis (void) __swapRP(rp);
17125c28e83SPiotr Jasiukajtis #endif
17225c28e83SPiotr Jasiukajtis }
17325c28e83SPiotr Jasiukajtis return (ans);
17425c28e83SPiotr Jasiukajtis }
175