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