25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER START
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * The contents of this file are subject to the terms of the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Common Development and Distribution License (the "License").
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You may not use this file except in compliance with the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * See the License for the specific language governing permissions
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * and limitations under the License.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * CDDL HEADER END
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT OFF */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * dcomplex cpow(dcomplex z);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z**w analytically equivalent to
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * cpow(z,w) = cexp(w clog(z))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Let z = x+iy, w = u+iv.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * log(x+iy) = log(\/ x + y ) + i tan (---)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * = --- log(x + y ) + i tan (---)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (u+iv)* log(x+iy) = --- log(x + y ) - v tan (---) + (1)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * i * [ --- log(x + y ) + u tan (---) ] (2)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * z = e = e (cos(q)+i*sin(q))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * r \/ x + y -v*atan2(y,x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Here e can be expressed as: u * e
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Special cases (in the order of appearance):
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 1. (anything) ** 0 is 1
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 2. (anything) ** 1 is itself
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 3. When v = 0, y = 0:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * If x is finite and negative, and u is finite, then
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x ** u = exp(u*pi i) * pow(|x|, u);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * x ** u = pow(x, u);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 4. When v = 0, x = 0 or |x| = |y| or x is inf or y is inf:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (x + y i) ** u = r * exp(q i)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * r = hypot(x,y) ** u
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * q = u * atan2pi(y, x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * 5. otherwise, z**w is NAN if any x, y, u, v is a Nan or inf
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Note: many results of special cases are obtained in terms of
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * polar coordinate. In the conversion from polar to rectangle:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * r exp(q i) = r * cos(q) + r * sin(q) i,
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * we regard r * 0 is 0 except when r is a NaN.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis/* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis#include "libm.h" /* atan2/exp/fabs/hypot/log/pow/scalbn */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* atan2pi/exp2/sincos/sincospi/__k_clog_r/__k_atan2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern void sincospi(double, double *, double *);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern double atan2pi(double, double);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Assuming |t[0]| > |t[1]| and |t[2]| > |t[3]|, sum4fp subroutine
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * compute t[0] + t[1] + t[2] + t[3] into two double fp numbers.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t1 = ta[0]; t2 = ta[1]; t3 = ta[2]; t4 = ta[3];
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Rearrange ti so that |t1| >= |t2| >= |t3| >= |t4|
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* summing r = t1 + t2 + t3 + t4 to w1 + w2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double b[4], t1, t2, t3, t4, w1, w2, u1, v1, x1, y1;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis int ix, iy, hx, lx, hy, ly, hv, hu, iu, iv, lu, lv;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (((hu - 0x3ff00000) | lu) == 0) { /* z ** 1 = z */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if ((iu | lu) == 0) { /* z ** 0 = 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if ((iy | ly) == 0) { /* (real)**(real) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* -x ** u is exp(i*pi*u)*pow(x,u) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (((ix | lx) == 0) || ix >= hiinf || iy >= hiinf) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sincospi(t * u, &s, &c);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (((ix - iy) | (lx - ly)) == 0) { /* |x| = |y| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sincospi(t * u, &s, &c);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sincospi(t * u, &s, &c);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis r = (hy >= 0)? u : -u;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = -0.25 * r;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (ix >= 0x3ff00000 || iu < 0x408ff800)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* |x| >= 1 or |u| < 1023 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else /* special treatment */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * non-zero imag part(s) with inf component(s) yields NaN
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = fabs(x) + fabs(y) + fabs(u) + fabs(v);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis k = 0; /* no scaling */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis u *= .0009765625; /* scale 2**-10 to avoid overflow */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis v *= .0009765625;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use similated higher precision arithmetic to compute:
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * r = u * log(hypot(x, y)) - v * atan2(y, x)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * q = u * atan2(y, x) + v * log(hypot(x, y))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis x2 = t2 - (x1 - t1); /* log(hypot(x,y)) = x1 + x2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis y2 = t4 - (y1 - t3); /* atan2(y,x) = y1 + y2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute q = u * atan2(y, x) + v * log(hypot(x, y)) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * square (cos(q) + i sin(q)) k times to get
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * (cos(2^k * q + i sin(2^k * q)
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis for (i = 0; i < 10; i++) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis c = (c + s) * (c - s);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute r = u * (t1, t2) - v * (t3, t4) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* check over/underflow for exp(w1 + w2) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* compute exp(w1 + w2) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix < 0x3c900000) /* exp(tiny < 2**-54) = 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (ix >= 0x40880000) /* overflow/underflow */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else { /* compute exp(w1 + w2) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis k = (int) (invln2 * w1 + ((hx >= 0)? 0.5 : -0.5));
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t1 = (double) k;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (c != zero) c *= r;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (s != zero) s *= r;