1*7c478bd9Sstevel@tonic-gate /* 2*7c478bd9Sstevel@tonic-gate * CDDL HEADER START 3*7c478bd9Sstevel@tonic-gate * 4*7c478bd9Sstevel@tonic-gate * The contents of this file are subject to the terms of the 5*7c478bd9Sstevel@tonic-gate * Common Development and Distribution License, Version 1.0 only 6*7c478bd9Sstevel@tonic-gate * (the "License"). You may not use this file except in compliance 7*7c478bd9Sstevel@tonic-gate * with the License. 8*7c478bd9Sstevel@tonic-gate * 9*7c478bd9Sstevel@tonic-gate * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 10*7c478bd9Sstevel@tonic-gate * or http://www.opensolaris.org/os/licensing. 11*7c478bd9Sstevel@tonic-gate * See the License for the specific language governing permissions 12*7c478bd9Sstevel@tonic-gate * and limitations under the License. 13*7c478bd9Sstevel@tonic-gate * 14*7c478bd9Sstevel@tonic-gate * When distributing Covered Code, include this CDDL HEADER in each 15*7c478bd9Sstevel@tonic-gate * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 16*7c478bd9Sstevel@tonic-gate * If applicable, add the following below this CDDL HEADER, with the 17*7c478bd9Sstevel@tonic-gate * fields enclosed by brackets "[]" replaced with your own identifying 18*7c478bd9Sstevel@tonic-gate * information: Portions Copyright [yyyy] [name of copyright owner] 19*7c478bd9Sstevel@tonic-gate * 20*7c478bd9Sstevel@tonic-gate * CDDL HEADER END 21*7c478bd9Sstevel@tonic-gate */ 22*7c478bd9Sstevel@tonic-gate /* 23*7c478bd9Sstevel@tonic-gate * Copyright 2003 Sun Microsystems, Inc. All rights reserved. 24*7c478bd9Sstevel@tonic-gate * Use is subject to license terms. 25*7c478bd9Sstevel@tonic-gate */ 26*7c478bd9Sstevel@tonic-gate 27*7c478bd9Sstevel@tonic-gate #pragma ident "%Z%%M% %I% %E% SMI" 28*7c478bd9Sstevel@tonic-gate 29*7c478bd9Sstevel@tonic-gate #include "quad.h" 30*7c478bd9Sstevel@tonic-gate 31*7c478bd9Sstevel@tonic-gate static const double C[] = { 32*7c478bd9Sstevel@tonic-gate 0.0, 33*7c478bd9Sstevel@tonic-gate 0.5, 34*7c478bd9Sstevel@tonic-gate 1.0, 35*7c478bd9Sstevel@tonic-gate 68719476736.0, 36*7c478bd9Sstevel@tonic-gate 536870912.0, 37*7c478bd9Sstevel@tonic-gate 48.0, 38*7c478bd9Sstevel@tonic-gate 16.0, 39*7c478bd9Sstevel@tonic-gate 1.52587890625000000000e-05, 40*7c478bd9Sstevel@tonic-gate 2.86102294921875000000e-06, 41*7c478bd9Sstevel@tonic-gate 5.96046447753906250000e-08, 42*7c478bd9Sstevel@tonic-gate 3.72529029846191406250e-09, 43*7c478bd9Sstevel@tonic-gate 1.70530256582424044609e-13, 44*7c478bd9Sstevel@tonic-gate 7.10542735760100185871e-15, 45*7c478bd9Sstevel@tonic-gate 8.67361737988403547206e-19, 46*7c478bd9Sstevel@tonic-gate 2.16840434497100886801e-19, 47*7c478bd9Sstevel@tonic-gate 1.27054942088145050860e-21, 48*7c478bd9Sstevel@tonic-gate 1.21169035041947413311e-27, 49*7c478bd9Sstevel@tonic-gate 9.62964972193617926528e-35, 50*7c478bd9Sstevel@tonic-gate 4.70197740328915003187e-38 51*7c478bd9Sstevel@tonic-gate }; 52*7c478bd9Sstevel@tonic-gate 53*7c478bd9Sstevel@tonic-gate #define zero C[0] 54*7c478bd9Sstevel@tonic-gate #define half C[1] 55*7c478bd9Sstevel@tonic-gate #define one C[2] 56*7c478bd9Sstevel@tonic-gate #define two36 C[3] 57*7c478bd9Sstevel@tonic-gate #define two29 C[4] 58*7c478bd9Sstevel@tonic-gate #define three2p4 C[5] 59*7c478bd9Sstevel@tonic-gate #define two4 C[6] 60*7c478bd9Sstevel@tonic-gate #define twom16 C[7] 61*7c478bd9Sstevel@tonic-gate #define three2m20 C[8] 62*7c478bd9Sstevel@tonic-gate #define twom24 C[9] 63*7c478bd9Sstevel@tonic-gate #define twom28 C[10] 64*7c478bd9Sstevel@tonic-gate #define three2m44 C[11] 65*7c478bd9Sstevel@tonic-gate #define twom47 C[12] 66*7c478bd9Sstevel@tonic-gate #define twom60 C[13] 67*7c478bd9Sstevel@tonic-gate #define twom62 C[14] 68*7c478bd9Sstevel@tonic-gate #define three2m71 C[15] 69*7c478bd9Sstevel@tonic-gate #define three2m91 C[16] 70*7c478bd9Sstevel@tonic-gate #define twom113 C[17] 71*7c478bd9Sstevel@tonic-gate #define twom124 C[18] 72*7c478bd9Sstevel@tonic-gate 73*7c478bd9Sstevel@tonic-gate static const unsigned 74*7c478bd9Sstevel@tonic-gate fsr_re = 0x00000000u, 75*7c478bd9Sstevel@tonic-gate fsr_rn = 0xc0000000u; 76*7c478bd9Sstevel@tonic-gate 77*7c478bd9Sstevel@tonic-gate #ifdef __sparcv9 78*7c478bd9Sstevel@tonic-gate 79*7c478bd9Sstevel@tonic-gate /* 80*7c478bd9Sstevel@tonic-gate * _Qp_sqrt(pz, x) sets *pz = sqrt(*x). 81*7c478bd9Sstevel@tonic-gate */ 82*7c478bd9Sstevel@tonic-gate void 83*7c478bd9Sstevel@tonic-gate _Qp_sqrt(union longdouble *pz, const union longdouble *x) 84*7c478bd9Sstevel@tonic-gate 85*7c478bd9Sstevel@tonic-gate #else 86*7c478bd9Sstevel@tonic-gate 87*7c478bd9Sstevel@tonic-gate /* 88*7c478bd9Sstevel@tonic-gate * _Q_sqrt(x) returns sqrt(*x). 89*7c478bd9Sstevel@tonic-gate */ 90*7c478bd9Sstevel@tonic-gate union longdouble 91*7c478bd9Sstevel@tonic-gate _Q_sqrt(const union longdouble *x) 92*7c478bd9Sstevel@tonic-gate 93*7c478bd9Sstevel@tonic-gate #endif /* __sparcv9 */ 94*7c478bd9Sstevel@tonic-gate 95*7c478bd9Sstevel@tonic-gate { 96*7c478bd9Sstevel@tonic-gate union longdouble z; 97*7c478bd9Sstevel@tonic-gate union xdouble u; 98*7c478bd9Sstevel@tonic-gate double c, d, rr, r[2], tt[3], xx[4], zz[5]; 99*7c478bd9Sstevel@tonic-gate unsigned int xm, fsr, lx, wx[3]; 100*7c478bd9Sstevel@tonic-gate unsigned int msw, frac2, frac3, frac4, rm; 101*7c478bd9Sstevel@tonic-gate int ex, ez; 102*7c478bd9Sstevel@tonic-gate 103*7c478bd9Sstevel@tonic-gate if (QUAD_ISZERO(*x)) { 104*7c478bd9Sstevel@tonic-gate Z = *x; 105*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 106*7c478bd9Sstevel@tonic-gate } 107*7c478bd9Sstevel@tonic-gate 108*7c478bd9Sstevel@tonic-gate xm = x->l.msw; 109*7c478bd9Sstevel@tonic-gate 110*7c478bd9Sstevel@tonic-gate __quad_getfsrp(&fsr); 111*7c478bd9Sstevel@tonic-gate 112*7c478bd9Sstevel@tonic-gate /* handle nan and inf cases */ 113*7c478bd9Sstevel@tonic-gate if ((xm & 0x7fffffff) >= 0x7fff0000) { 114*7c478bd9Sstevel@tonic-gate if ((x->l.msw & 0xffff) | x->l.frac2 | x->l.frac3 | 115*7c478bd9Sstevel@tonic-gate x->l.frac4) { 116*7c478bd9Sstevel@tonic-gate if (!(x->l.msw & 0x8000)) { 117*7c478bd9Sstevel@tonic-gate /* snan, signal invalid */ 118*7c478bd9Sstevel@tonic-gate if (fsr & FSR_NVM) { 119*7c478bd9Sstevel@tonic-gate __quad_fsqrtq(x, &Z); 120*7c478bd9Sstevel@tonic-gate } else { 121*7c478bd9Sstevel@tonic-gate Z = *x; 122*7c478bd9Sstevel@tonic-gate Z.l.msw |= 0x8000; 123*7c478bd9Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA | 124*7c478bd9Sstevel@tonic-gate FSR_NVC; 125*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr); 126*7c478bd9Sstevel@tonic-gate } 127*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 128*7c478bd9Sstevel@tonic-gate } 129*7c478bd9Sstevel@tonic-gate Z = *x; 130*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 131*7c478bd9Sstevel@tonic-gate } 132*7c478bd9Sstevel@tonic-gate if (x->l.msw & 0x80000000) { 133*7c478bd9Sstevel@tonic-gate /* sqrt(-inf), signal invalid */ 134*7c478bd9Sstevel@tonic-gate if (fsr & FSR_NVM) { 135*7c478bd9Sstevel@tonic-gate __quad_fsqrtq(x, &Z); 136*7c478bd9Sstevel@tonic-gate } else { 137*7c478bd9Sstevel@tonic-gate Z.l.msw = 0x7fffffff; 138*7c478bd9Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff; 139*7c478bd9Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC; 140*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr); 141*7c478bd9Sstevel@tonic-gate } 142*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 143*7c478bd9Sstevel@tonic-gate } 144*7c478bd9Sstevel@tonic-gate /* sqrt(inf), return inf */ 145*7c478bd9Sstevel@tonic-gate Z = *x; 146*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 147*7c478bd9Sstevel@tonic-gate } 148*7c478bd9Sstevel@tonic-gate 149*7c478bd9Sstevel@tonic-gate /* handle negative numbers */ 150*7c478bd9Sstevel@tonic-gate if (xm & 0x80000000) { 151*7c478bd9Sstevel@tonic-gate if (fsr & FSR_NVM) { 152*7c478bd9Sstevel@tonic-gate __quad_fsqrtq(x, &Z); 153*7c478bd9Sstevel@tonic-gate } else { 154*7c478bd9Sstevel@tonic-gate Z.l.msw = 0x7fffffff; 155*7c478bd9Sstevel@tonic-gate Z.l.frac2 = Z.l.frac3 = Z.l.frac4 = 0xffffffff; 156*7c478bd9Sstevel@tonic-gate fsr = (fsr & ~FSR_CEXC) | FSR_NVA | FSR_NVC; 157*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr); 158*7c478bd9Sstevel@tonic-gate } 159*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 160*7c478bd9Sstevel@tonic-gate } 161*7c478bd9Sstevel@tonic-gate 162*7c478bd9Sstevel@tonic-gate /* now x is finite, positive */ 163*7c478bd9Sstevel@tonic-gate __quad_setfsrp((unsigned *)&fsr_re); 164*7c478bd9Sstevel@tonic-gate 165*7c478bd9Sstevel@tonic-gate /* get the normalized significand and exponent */ 166*7c478bd9Sstevel@tonic-gate ex = (int)(xm >> 16); 167*7c478bd9Sstevel@tonic-gate lx = xm & 0xffff; 168*7c478bd9Sstevel@tonic-gate if (ex) { 169*7c478bd9Sstevel@tonic-gate lx |= 0x10000; 170*7c478bd9Sstevel@tonic-gate wx[0] = x->l.frac2; 171*7c478bd9Sstevel@tonic-gate wx[1] = x->l.frac3; 172*7c478bd9Sstevel@tonic-gate wx[2] = x->l.frac4; 173*7c478bd9Sstevel@tonic-gate } else { 174*7c478bd9Sstevel@tonic-gate if (lx | (x->l.frac2 & 0xfffe0000)) { 175*7c478bd9Sstevel@tonic-gate wx[0] = x->l.frac2; 176*7c478bd9Sstevel@tonic-gate wx[1] = x->l.frac3; 177*7c478bd9Sstevel@tonic-gate wx[2] = x->l.frac4; 178*7c478bd9Sstevel@tonic-gate ex = 1; 179*7c478bd9Sstevel@tonic-gate } else if (x->l.frac2 | (x->l.frac3 & 0xfffe0000)) { 180*7c478bd9Sstevel@tonic-gate lx = x->l.frac2; 181*7c478bd9Sstevel@tonic-gate wx[0] = x->l.frac3; 182*7c478bd9Sstevel@tonic-gate wx[1] = x->l.frac4; 183*7c478bd9Sstevel@tonic-gate wx[2] = 0; 184*7c478bd9Sstevel@tonic-gate ex = -31; 185*7c478bd9Sstevel@tonic-gate } else if (x->l.frac3 | (x->l.frac4 & 0xfffe0000)) { 186*7c478bd9Sstevel@tonic-gate lx = x->l.frac3; 187*7c478bd9Sstevel@tonic-gate wx[0] = x->l.frac4; 188*7c478bd9Sstevel@tonic-gate wx[1] = wx[2] = 0; 189*7c478bd9Sstevel@tonic-gate ex = -63; 190*7c478bd9Sstevel@tonic-gate } else { 191*7c478bd9Sstevel@tonic-gate lx = x->l.frac4; 192*7c478bd9Sstevel@tonic-gate wx[0] = wx[1] = wx[2] = 0; 193*7c478bd9Sstevel@tonic-gate ex = -95; 194*7c478bd9Sstevel@tonic-gate } 195*7c478bd9Sstevel@tonic-gate while ((lx & 0x10000) == 0) { 196*7c478bd9Sstevel@tonic-gate lx = (lx << 1) | (wx[0] >> 31); 197*7c478bd9Sstevel@tonic-gate wx[0] = (wx[0] << 1) | (wx[1] >> 31); 198*7c478bd9Sstevel@tonic-gate wx[1] = (wx[1] << 1) | (wx[2] >> 31); 199*7c478bd9Sstevel@tonic-gate wx[2] <<= 1; 200*7c478bd9Sstevel@tonic-gate ex--; 201*7c478bd9Sstevel@tonic-gate } 202*7c478bd9Sstevel@tonic-gate } 203*7c478bd9Sstevel@tonic-gate ez = ex - 0x3fff; 204*7c478bd9Sstevel@tonic-gate if (ez & 1) { 205*7c478bd9Sstevel@tonic-gate /* make exponent even */ 206*7c478bd9Sstevel@tonic-gate lx = (lx << 1) | (wx[0] >> 31); 207*7c478bd9Sstevel@tonic-gate wx[0] = (wx[0] << 1) | (wx[1] >> 31); 208*7c478bd9Sstevel@tonic-gate wx[1] = (wx[1] << 1) | (wx[2] >> 31); 209*7c478bd9Sstevel@tonic-gate wx[2] <<= 1; 210*7c478bd9Sstevel@tonic-gate ez--; 211*7c478bd9Sstevel@tonic-gate } 212*7c478bd9Sstevel@tonic-gate 213*7c478bd9Sstevel@tonic-gate /* extract the significands into doubles */ 214*7c478bd9Sstevel@tonic-gate c = twom16; 215*7c478bd9Sstevel@tonic-gate xx[0] = (double)((int)lx) * c; 216*7c478bd9Sstevel@tonic-gate 217*7c478bd9Sstevel@tonic-gate c *= twom24; 218*7c478bd9Sstevel@tonic-gate xx[0] += (double)((int)(wx[0] >> 8)) * c; 219*7c478bd9Sstevel@tonic-gate 220*7c478bd9Sstevel@tonic-gate c *= twom24; 221*7c478bd9Sstevel@tonic-gate xx[1] = (double)((int)(((wx[0] << 16) | (wx[1] >> 16)) & 222*7c478bd9Sstevel@tonic-gate 0xffffff)) * c; 223*7c478bd9Sstevel@tonic-gate 224*7c478bd9Sstevel@tonic-gate c *= twom24; 225*7c478bd9Sstevel@tonic-gate xx[2] = (double)((int)(((wx[1] << 8) | (wx[2] >> 24)) & 226*7c478bd9Sstevel@tonic-gate 0xffffff)) * c; 227*7c478bd9Sstevel@tonic-gate 228*7c478bd9Sstevel@tonic-gate c *= twom24; 229*7c478bd9Sstevel@tonic-gate xx[3] = (double)((int)(wx[2] & 0xffffff)) * c; 230*7c478bd9Sstevel@tonic-gate 231*7c478bd9Sstevel@tonic-gate /* approximate the divisor for the Newton iteration */ 232*7c478bd9Sstevel@tonic-gate c = xx[0] + xx[1]; 233*7c478bd9Sstevel@tonic-gate c = __quad_dp_sqrt(&c); 234*7c478bd9Sstevel@tonic-gate rr = half / c; 235*7c478bd9Sstevel@tonic-gate 236*7c478bd9Sstevel@tonic-gate /* compute the first five "digits" of the square root */ 237*7c478bd9Sstevel@tonic-gate zz[0] = (c + two29) - two29; 238*7c478bd9Sstevel@tonic-gate tt[0] = zz[0] + zz[0]; 239*7c478bd9Sstevel@tonic-gate r[0] = (xx[0] - zz[0] * zz[0]) + xx[1]; 240*7c478bd9Sstevel@tonic-gate 241*7c478bd9Sstevel@tonic-gate zz[1] = (rr * (r[0] + xx[2]) + three2p4) - three2p4; 242*7c478bd9Sstevel@tonic-gate tt[1] = zz[1] + zz[1]; 243*7c478bd9Sstevel@tonic-gate r[0] -= tt[0] * zz[1]; 244*7c478bd9Sstevel@tonic-gate r[1] = xx[2] - zz[1] * zz[1]; 245*7c478bd9Sstevel@tonic-gate c = (r[1] + three2m20) - three2m20; 246*7c478bd9Sstevel@tonic-gate r[0] += c; 247*7c478bd9Sstevel@tonic-gate r[1] = (r[1] - c) + xx[3]; 248*7c478bd9Sstevel@tonic-gate 249*7c478bd9Sstevel@tonic-gate zz[2] = (rr * (r[0] + r[1]) + three2m20) - three2m20; 250*7c478bd9Sstevel@tonic-gate tt[2] = zz[2] + zz[2]; 251*7c478bd9Sstevel@tonic-gate r[0] -= tt[0] * zz[2]; 252*7c478bd9Sstevel@tonic-gate r[1] -= tt[1] * zz[2]; 253*7c478bd9Sstevel@tonic-gate c = (r[1] + three2m44) - three2m44; 254*7c478bd9Sstevel@tonic-gate r[0] += c; 255*7c478bd9Sstevel@tonic-gate r[1] = (r[1] - c) - zz[2] * zz[2]; 256*7c478bd9Sstevel@tonic-gate 257*7c478bd9Sstevel@tonic-gate zz[3] = (rr * (r[0] + r[1]) + three2m44) - three2m44; 258*7c478bd9Sstevel@tonic-gate r[0] = ((r[0] - tt[0] * zz[3]) + r[1]) - tt[1] * zz[3]; 259*7c478bd9Sstevel@tonic-gate r[1] = -tt[2] * zz[3]; 260*7c478bd9Sstevel@tonic-gate c = (r[1] + three2m91) - three2m91; 261*7c478bd9Sstevel@tonic-gate r[0] += c; 262*7c478bd9Sstevel@tonic-gate r[1] = (r[1] - c) - zz[3] * zz[3]; 263*7c478bd9Sstevel@tonic-gate 264*7c478bd9Sstevel@tonic-gate zz[4] = (rr * (r[0] + r[1]) + three2m71) - three2m71; 265*7c478bd9Sstevel@tonic-gate 266*7c478bd9Sstevel@tonic-gate /* reduce to three doubles, making sure zz[1] is positive */ 267*7c478bd9Sstevel@tonic-gate zz[0] += zz[1] - twom47; 268*7c478bd9Sstevel@tonic-gate zz[1] = twom47 + zz[2] + zz[3]; 269*7c478bd9Sstevel@tonic-gate zz[2] = zz[4]; 270*7c478bd9Sstevel@tonic-gate 271*7c478bd9Sstevel@tonic-gate /* if the third term might lie on a rounding boundary, perturb it */ 272*7c478bd9Sstevel@tonic-gate if (zz[2] == (twom62 + zz[2]) - twom62) { 273*7c478bd9Sstevel@tonic-gate /* here we just need to get the sign of the remainder */ 274*7c478bd9Sstevel@tonic-gate c = (((((r[0] - tt[0] * zz[4]) - tt[1] * zz[4]) + r[1]) 275*7c478bd9Sstevel@tonic-gate - tt[2] * zz[4]) - (zz[3] + zz[3]) * zz[4]) - zz[4] * zz[4]; 276*7c478bd9Sstevel@tonic-gate if (c < zero) 277*7c478bd9Sstevel@tonic-gate zz[2] -= twom124; 278*7c478bd9Sstevel@tonic-gate else if (c > zero) 279*7c478bd9Sstevel@tonic-gate zz[2] += twom124; 280*7c478bd9Sstevel@tonic-gate } 281*7c478bd9Sstevel@tonic-gate 282*7c478bd9Sstevel@tonic-gate /* 283*7c478bd9Sstevel@tonic-gate * propagate carries/borrows, using round-to-negative-infinity mode 284*7c478bd9Sstevel@tonic-gate * to make all terms nonnegative (note that we can't encounter a 285*7c478bd9Sstevel@tonic-gate * borrow so large that the roundoff is unrepresentable because 286*7c478bd9Sstevel@tonic-gate * we took care to make zz[1] positive above) 287*7c478bd9Sstevel@tonic-gate */ 288*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr_rn); 289*7c478bd9Sstevel@tonic-gate c = zz[1] + zz[2]; 290*7c478bd9Sstevel@tonic-gate zz[2] += (zz[1] - c); 291*7c478bd9Sstevel@tonic-gate zz[1] = c; 292*7c478bd9Sstevel@tonic-gate c = zz[0] + zz[1]; 293*7c478bd9Sstevel@tonic-gate zz[1] += (zz[0] - c); 294*7c478bd9Sstevel@tonic-gate zz[0] = c; 295*7c478bd9Sstevel@tonic-gate 296*7c478bd9Sstevel@tonic-gate /* adjust exponent and strip off integer bit */ 297*7c478bd9Sstevel@tonic-gate ez = (ez >> 1) + 0x3fff; 298*7c478bd9Sstevel@tonic-gate zz[0] -= one; 299*7c478bd9Sstevel@tonic-gate 300*7c478bd9Sstevel@tonic-gate /* the first 48 bits of fraction come from zz[0] */ 301*7c478bd9Sstevel@tonic-gate u.d = d = two36 + zz[0]; 302*7c478bd9Sstevel@tonic-gate msw = u.l.lo; 303*7c478bd9Sstevel@tonic-gate zz[0] -= (d - two36); 304*7c478bd9Sstevel@tonic-gate 305*7c478bd9Sstevel@tonic-gate u.d = d = two4 + zz[0]; 306*7c478bd9Sstevel@tonic-gate frac2 = u.l.lo; 307*7c478bd9Sstevel@tonic-gate zz[0] -= (d - two4); 308*7c478bd9Sstevel@tonic-gate 309*7c478bd9Sstevel@tonic-gate /* the next 32 come from zz[0] and zz[1] */ 310*7c478bd9Sstevel@tonic-gate u.d = d = twom28 + (zz[0] + zz[1]); 311*7c478bd9Sstevel@tonic-gate frac3 = u.l.lo; 312*7c478bd9Sstevel@tonic-gate zz[0] -= (d - twom28); 313*7c478bd9Sstevel@tonic-gate 314*7c478bd9Sstevel@tonic-gate /* condense the remaining fraction; errors here won't matter */ 315*7c478bd9Sstevel@tonic-gate c = zz[0] + zz[1]; 316*7c478bd9Sstevel@tonic-gate zz[1] = ((zz[0] - c) + zz[1]) + zz[2]; 317*7c478bd9Sstevel@tonic-gate zz[0] = c; 318*7c478bd9Sstevel@tonic-gate 319*7c478bd9Sstevel@tonic-gate /* get the last word of fraction */ 320*7c478bd9Sstevel@tonic-gate u.d = d = twom60 + (zz[0] + zz[1]); 321*7c478bd9Sstevel@tonic-gate frac4 = u.l.lo; 322*7c478bd9Sstevel@tonic-gate zz[0] -= (d - twom60); 323*7c478bd9Sstevel@tonic-gate 324*7c478bd9Sstevel@tonic-gate /* keep track of what's left for rounding; note that the error */ 325*7c478bd9Sstevel@tonic-gate /* in computing c will be non-negative due to rounding mode */ 326*7c478bd9Sstevel@tonic-gate c = zz[0] + zz[1]; 327*7c478bd9Sstevel@tonic-gate 328*7c478bd9Sstevel@tonic-gate /* get the rounding mode */ 329*7c478bd9Sstevel@tonic-gate rm = fsr >> 30; 330*7c478bd9Sstevel@tonic-gate 331*7c478bd9Sstevel@tonic-gate /* round and raise exceptions */ 332*7c478bd9Sstevel@tonic-gate fsr &= ~FSR_CEXC; 333*7c478bd9Sstevel@tonic-gate if (c != zero) { 334*7c478bd9Sstevel@tonic-gate fsr |= FSR_NXC; 335*7c478bd9Sstevel@tonic-gate 336*7c478bd9Sstevel@tonic-gate /* decide whether to round the fraction up */ 337*7c478bd9Sstevel@tonic-gate if (rm == FSR_RP || (rm == FSR_RN && (c > twom113 || 338*7c478bd9Sstevel@tonic-gate (c == twom113 && ((frac4 & 1) || (c - zz[0] != zz[1])))))) { 339*7c478bd9Sstevel@tonic-gate /* round up and renormalize if necessary */ 340*7c478bd9Sstevel@tonic-gate if (++frac4 == 0) 341*7c478bd9Sstevel@tonic-gate if (++frac3 == 0) 342*7c478bd9Sstevel@tonic-gate if (++frac2 == 0) 343*7c478bd9Sstevel@tonic-gate if (++msw == 0x10000) { 344*7c478bd9Sstevel@tonic-gate msw = 0; 345*7c478bd9Sstevel@tonic-gate ez++; 346*7c478bd9Sstevel@tonic-gate } 347*7c478bd9Sstevel@tonic-gate } 348*7c478bd9Sstevel@tonic-gate } 349*7c478bd9Sstevel@tonic-gate 350*7c478bd9Sstevel@tonic-gate /* stow the result */ 351*7c478bd9Sstevel@tonic-gate z.l.msw = (ez << 16) | msw; 352*7c478bd9Sstevel@tonic-gate z.l.frac2 = frac2; 353*7c478bd9Sstevel@tonic-gate z.l.frac3 = frac3; 354*7c478bd9Sstevel@tonic-gate z.l.frac4 = frac4; 355*7c478bd9Sstevel@tonic-gate 356*7c478bd9Sstevel@tonic-gate if ((fsr & FSR_CEXC) & (fsr >> 23)) { 357*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr); 358*7c478bd9Sstevel@tonic-gate __quad_fsqrtq(x, &Z); 359*7c478bd9Sstevel@tonic-gate } else { 360*7c478bd9Sstevel@tonic-gate Z = z; 361*7c478bd9Sstevel@tonic-gate fsr |= (fsr & 0x1f) << 5; 362*7c478bd9Sstevel@tonic-gate __quad_setfsrp(&fsr); 363*7c478bd9Sstevel@tonic-gate } 364*7c478bd9Sstevel@tonic-gate QUAD_RETURN(Z); 365*7c478bd9Sstevel@tonic-gate } 366