cpow.c revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8
/*
* 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
* 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 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/* INDENT OFF */
/*
* dcomplex cpow(dcomplex z);
*
* z**w analytically equivalent to
*
* cpow(z,w) = cexp(w clog(z))
*
* Let z = x+iy, w = u+iv.
* Since
* _________
* / 2 2 -1 y
* log(x+iy) = log(\/ x + y ) + i tan (---)
* x
*
* 1 2 2 -1 y
* = --- log(x + y ) + i tan (---)
* 2 x
* u 2 2 -1 y
* (u+iv)* log(x+iy) = --- log(x + y ) - v tan (---) + (1)
* 2 x
*
* v 2 2 -1 y
* i * [ --- log(x + y ) + u tan (---) ] (2)
* 2 x
*
* = r + i q
*
* Therefore,
* w r+iq r
* z = e = e (cos(q)+i*sin(q))
* _______
* / 2 2
* r \/ x + y -v*atan2(y,x)
* Here e can be expressed as: u * e
*
* Special cases (in the order of appearance):
* 1. (anything) ** 0 is 1
* 2. (anything) ** 1 is itself
* 3. When v = 0, y = 0:
* If x is finite and negative, and u is finite, then
* x ** u = exp(u*pi i) * pow(|x|, u);
* otherwise,
* x ** u = pow(x, u);
* 4. When v = 0, x = 0 or |x| = |y| or x is inf or y is inf:
* (x + y i) ** u = r * exp(q i)
* where
* r = hypot(x,y) ** u
* q = u * atan2pi(y, x)
*
* 5. otherwise, z**w is NAN if any x, y, u, v is a Nan or inf
*
* Note: many results of special cases are obtained in terms of
* polar coordinate. In the conversion from polar to rectangle:
* r exp(q i) = r * cos(q) + r * sin(q) i,
* we regard r * 0 is 0 except when r is a NaN.
*/
/* INDENT ON */
#include "complex_wrapper.h"
extern void sincospi(double, double *, double *);
static const double
huge = 1e300,
tiny = 1e-300,
invln2 = 1.44269504088896338700e+00,
one = 1.0,
zero = 0.0;
static const int hiinf = 0x7ff00000;
extern double atan2pi(double, double);
/*
* Assuming |t[0]| > |t[1]| and |t[2]| > |t[3]|, sum4fp subroutine
* compute t[0] + t[1] + t[2] + t[3] into two double fp numbers.
*/
static double
/*
* Rearrange ti so that |t1| >= |t2| >= |t3| >= |t4|
*/
} else {
}
} else
t3 = t;
}
/* summing r = t1 + t2 + t3 + t4 to w1 + w2 */
return (w1);
}
int i, j, k;
x = D_RE(z);
y = D_IM(z);
u = D_RE(w);
v = D_IM(w);
j = 0;
/* -x ** u is exp(i*pi*u)*pow(x,u) */
r = pow(-x, u);
sincospi(u, &s, &c);
} else
else {
r = fabs(y);
else
t = atan2pi(y, x);
sincospi(t * u, &s, &c);
}
if (hx >= 0) {
sincospi(t * u, &s, &c);
} else if ((lu & 3) == 0) {
sincospi(t * u, &s, &c);
} else {
r = (hy >= 0)? u : -u;
t = -0.25 * r;
w1 = r + t;
}
/* |x| >= 1 or |u| < 1023 */
else /* special treatment */
j = 2;
if (j == 0) {
}
} else
j = 1;
if (j == 0)
return (ans);
}
/*
* non-zero imag part(s) with inf component(s) yields NaN
*/
} else {
k = 0; /* no scaling */
u *= .0009765625; /* scale 2**-10 to avoid overflow */
v *= .0009765625;
k = 1; /* scale by 2**-10 */
}
/*
* Use similated higher precision arithmetic to compute:
* r = u * log(hypot(x, y)) - v * atan2(y, x)
* q = u * atan2(y, x) + v * log(hypot(x, y))
*/
u1 = u;
v1 = v;
/* compute q = u * atan2(y, x) + v * log(hypot(x, y)) */
if (j != 2) {
if (j == 1) { /* v = 0 */
w1 = b[0] + b[1];
} else {
}
if (k == 1)
/*
* square (cos(q) + i sin(q)) k times to get
* (cos(2^k * q + i sin(2^k * q)
*/
for (i = 0; i < 10; i++) {
t1 = s * c;
c = (c + s) * (c - s);
}
}
/* compute r = u * (t1, t2) - v * (t3, t4) */
if (j == 1) { /* v = 0 */
w1 = b[0] + b[1];
} else {
}
}
/* compute exp(w1 + w2) */
r = one;
else { /* compute exp(w1 + w2) */
t1 = (double) k;
}
if (c != zero) c *= r;
if (s != zero) s *= r;
if (k != 0) {
c = scalbn(c, k);
s = scalbn(s, k);
}
}
return (ans);
}