2N/A * The contents of this file are subject to the terms of the 2N/A * Common Development and Distribution License, Version 1.0 only 2N/A * (the "License"). You may not use this file except in compliance 2N/A * See the License for the specific language governing permissions 2N/A * and limitations under the License. 2N/A * When distributing Covered Code, include this CDDL HEADER in each 2N/A * If applicable, add the following below this CDDL HEADER, with the 2N/A * fields enclosed by brackets "[]" replaced with your own identifying 2N/A * information: Portions Copyright [yyyy] [name of copyright owner] 2N/A * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 2N/A * Use is subject to license terms. 2N/A#
pragma ident "%Z%%M% %I% %E% SMI" 2N/Astatic const double C[] = {
2N/A 1.52587890625000000000e-05,
2N/A 2.86102294921875000000e-06,
2N/A 5.96046447753906250000e-08,
2N/A 3.72529029846191406250e-09,
2N/A 1.70530256582424044609e-13,
2N/A 7.10542735760100185871e-15,
2N/A 8.67361737988403547206e-19,
2N/A 2.16840434497100886801e-19,
2N/A 1.27054942088145050860e-21,
2N/A 1.21169035041947413311e-27,
2N/A 9.62964972193617926528e-35,
2N/A 4.70197740328915003187e-38 2N/Astatic const unsigned 2N/A * _Qp_sqrt(pz, x) sets *pz = sqrt(*x). 2N/A * _Q_sqrt(x) returns sqrt(*x). 2N/A#
endif /* __sparcv9 */ 2N/A /* handle nan and inf cases */ 2N/A if ((
xm &
0x7fffffff) >=
0x7fff0000) {
2N/A /* snan, signal invalid */ 2N/A /* sqrt(-inf), signal invalid */ 2N/A /* sqrt(inf), return inf */ 2N/A /* handle negative numbers */ 2N/A /* now x is finite, positive */ 2N/A /* get the normalized significand and exponent */ 2N/A while ((
lx &
0x10000) == 0) {
2N/A /* make exponent even */ 2N/A /* extract the significands into doubles */ 2N/A xx[0] += (
double)((
int)(
wx[0] >>
8)) * c;
2N/A xx[
1] = (
double)((
int)(((
wx[0] <<
16) | (
wx[
1] >>
16)) &
2N/A xx[
2] = (
double)((
int)(((
wx[
1] <<
8) | (
wx[
2] >>
24)) &
2N/A xx[
3] = (
double)((
int)(
wx[
2] &
0xffffff)) * c;
2N/A /* approximate the divisor for the Newton iteration */ 2N/A /* compute the first five "digits" of the square root */ 2N/A r[
1] = (r[
1] - c) +
xx[
3];
2N/A r[0] = ((r[0] -
tt[0] *
zz[
3]) + r[
1]) -
tt[
1] *
zz[
3];
2N/A /* reduce to three doubles, making sure zz[1] is positive */ 2N/A /* if the third term might lie on a rounding boundary, perturb it */ 2N/A /* here we just need to get the sign of the remainder */ 2N/A c = (((((r[0] -
tt[0] *
zz[
4]) -
tt[
1] *
zz[
4]) + r[
1])
2N/A * to make all terms nonnegative (note that we can't encounter a 2N/A * borrow so large that the roundoff is unrepresentable because 2N/A * we took care to make zz[1] positive above) 2N/A /* adjust exponent and strip off integer bit */ 2N/A /* the first 48 bits of fraction come from zz[0] */ 2N/A /* the next 32 come from zz[0] and zz[1] */ 2N/A /* condense the remaining fraction; errors here won't matter */ 2N/A /* get the last word of fraction */ 2N/A /* keep track of what's left for rounding; note that the error */ 2N/A /* in computing c will be non-negative due to rounding mode */ 2N/A /* get the rounding mode */ 2N/A /* round and raise exceptions */ 2N/A /* decide whether to round the fraction up */ 2N/A /* round up and renormalize if necessary */ 2N/A /* stow the result */