1 /* IEEE754 floating point arithmetic 2 * double precision square root 3 */ 4 /* 5 * MIPS floating point support 6 * Copyright (C) 1994-2000 Algorithmics Ltd. 7 * 8 * ######################################################################## 9 * 10 * This program is free software; you can distribute it and/or modify it 11 * under the terms of the GNU General Public License (Version 2) as 12 * published by the Free Software Foundation. 13 * 14 * This program is distributed in the hope it will be useful, but WITHOUT 15 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 16 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 17 * for more details. 18 * 19 * You should have received a copy of the GNU General Public License along 20 * with this program; if not, write to the Free Software Foundation, Inc., 21 * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. 22 * 23 * ######################################################################## 24 */ 25 26 27 #include "ieee754dp.h" 28 29 static const unsigned table[] = { 30 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 31 29598, 36145, 43202, 50740, 58733, 67158, 75992, 32 85215, 83599, 71378, 60428, 50647, 41945, 34246, 33 27478, 21581, 16499, 12183, 8588, 5674, 3403, 34 1742, 661, 130 35 }; 36 37 ieee754dp ieee754dp_sqrt(ieee754dp x) 38 { 39 struct _ieee754_csr oldcsr; 40 ieee754dp y, z, t; 41 unsigned scalx, yh; 42 COMPXDP; 43 44 EXPLODEXDP; 45 CLEARCX; 46 FLUSHXDP; 47 48 /* x == INF or NAN? */ 49 switch (xc) { 50 case IEEE754_CLASS_QNAN: 51 /* sqrt(Nan) = Nan */ 52 return ieee754dp_nanxcpt(x, "sqrt"); 53 case IEEE754_CLASS_SNAN: 54 SETCX(IEEE754_INVALID_OPERATION); 55 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); 56 case IEEE754_CLASS_ZERO: 57 /* sqrt(0) = 0 */ 58 return x; 59 case IEEE754_CLASS_INF: 60 if (xs) { 61 /* sqrt(-Inf) = Nan */ 62 SETCX(IEEE754_INVALID_OPERATION); 63 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); 64 } 65 /* sqrt(+Inf) = Inf */ 66 return x; 67 case IEEE754_CLASS_DNORM: 68 DPDNORMX; 69 /* fall through */ 70 case IEEE754_CLASS_NORM: 71 if (xs) { 72 /* sqrt(-x) = Nan */ 73 SETCX(IEEE754_INVALID_OPERATION); 74 return ieee754dp_nanxcpt(ieee754dp_indef(), "sqrt"); 75 } 76 break; 77 } 78 79 /* save old csr; switch off INX enable & flag; set RN rounding */ 80 oldcsr = ieee754_csr; 81 ieee754_csr.mx &= ~IEEE754_INEXACT; 82 ieee754_csr.sx &= ~IEEE754_INEXACT; 83 ieee754_csr.rm = IEEE754_RN; 84 85 /* adjust exponent to prevent overflow */ 86 scalx = 0; 87 if (xe > 512) { /* x > 2**-512? */ 88 xe -= 512; /* x = x / 2**512 */ 89 scalx += 256; 90 } else if (xe < -512) { /* x < 2**-512? */ 91 xe += 512; /* x = x * 2**512 */ 92 scalx -= 256; 93 } 94 95 y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 96 97 /* magic initial approximation to almost 8 sig. bits */ 98 yh = y.bits >> 32; 99 yh = (yh >> 1) + 0x1ff80000; 100 yh = yh - table[(yh >> 15) & 31]; 101 y.bits = ((u64) yh << 32) | (y.bits & 0xffffffff); 102 103 /* Heron's rule once with correction to improve to ~18 sig. bits */ 104 /* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */ 105 t = ieee754dp_div(x, y); 106 y = ieee754dp_add(y, t); 107 y.bits -= 0x0010000600000000LL; 108 y.bits &= 0xffffffff00000000LL; 109 110 /* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */ 111 /* t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */ 112 z = t = ieee754dp_mul(y, y); 113 t.parts.bexp += 0x001; 114 t = ieee754dp_add(t, z); 115 z = ieee754dp_mul(ieee754dp_sub(x, z), y); 116 117 /* t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; */ 118 t = ieee754dp_div(z, ieee754dp_add(t, x)); 119 t.parts.bexp += 0x001; 120 y = ieee754dp_add(y, t); 121 122 /* twiddle last bit to force y correctly rounded */ 123 124 /* set RZ, clear INEX flag */ 125 ieee754_csr.rm = IEEE754_RZ; 126 ieee754_csr.sx &= ~IEEE754_INEXACT; 127 128 /* t=x/y; ...chopped quotient, possibly inexact */ 129 t = ieee754dp_div(x, y); 130 131 if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) { 132 133 if (!(ieee754_csr.sx & IEEE754_INEXACT)) 134 /* t = t-ulp */ 135 t.bits -= 1; 136 137 /* add inexact to result status */ 138 oldcsr.cx |= IEEE754_INEXACT; 139 oldcsr.sx |= IEEE754_INEXACT; 140 141 switch (oldcsr.rm) { 142 case IEEE754_RP: 143 y.bits += 1; 144 /* drop through */ 145 case IEEE754_RN: 146 t.bits += 1; 147 break; 148 } 149 150 /* y=y+t; ...chopped sum */ 151 y = ieee754dp_add(y, t); 152 153 /* adjust scalx for correctly rounded sqrt(x) */ 154 scalx -= 1; 155 } 156 157 /* py[n0]=py[n0]+scalx; ...scale back y */ 158 y.parts.bexp += scalx; 159 160 /* restore rounding mode, possibly set inexact */ 161 ieee754_csr = oldcsr; 162 163 return y; 164 } 165