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 2005 Sun Microsystems, Inc. All rights reserved.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * Use is subject to license terms.
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern void sincospi(double, double *, double *);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern void sincospif(float, float *, float *);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern double atan2pi(double, double);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern float atan2pif(float, float);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisextern int __swapRP(int);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const double
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dpi = 3.1415926535897931160E0, /* Hex 2^ 1 * 1.921FB54442D18 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dsqrt2 = 1.41421356237309514547, /* 3FF6A09E 667F3BCD */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtisstatic const float one = 1.0F, zero = 0.0F;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis float x, y, u, v, t, c, s;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis double dx, dy, du, dv, dt, dc, ds, dp, dq, dr;
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (hu == 0x3f800000) { /* (anything) ** 1 is itself */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iu == 0) { /* (anything) ** 0 is 1 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (iy == 0) { /* (real)**(real) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* -x ** u is exp(i*pi*u)*pow(x,u) */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix == 0 || ix >= hiinf || iy >= hiinf) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix > hiinf || iy > hiinf || iu > hiinf) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis sincospif(t * u, &s, &c);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis } else if (ix == iy) { /* if |x| == |y| */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (iu >= hiinf || iv >= hiinf || ix >= hiinf || iy >= hiinf) {
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * non-zero imaginery part(s) with inf component(s) yields NaN
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis t = fabsf(x) + fabsf(y) + fabsf(u) + fabsf(v);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT OFF */
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 * r = u*log(hypot(x,y))-v*pi*atan2pi(y,x),
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * q/pi = u*atan2pi(y,x)+v*log(hypot(x,y))/pi
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis * ans = exp(r)*(cospi(q/pi) + i sinpi(q/pi))
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis /* INDENT ON */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis if (ix > 0x3f000000 && ix < 0x40000000) /* .5 < |x| < 2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dt = dhalf * log1p((dx - 1.0) * (dx + 1.0) + dy * dy);
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis else if (iy > 0x3f000000 && iy < 0x40000000) /* .5 < |y| < 2 */
25c28e83beb90e7c80452a7c818c5e6f73a07dc8Piotr Jasiukajtis dt = dhalf * log1p((dy - 1.0) * (dy + 1.0) + dx * dx);