11da177e4SLinus Torvalds /* IEEE754 floating point arithmetic 21da177e4SLinus Torvalds * double precision: common utilities 31da177e4SLinus Torvalds */ 41da177e4SLinus Torvalds /* 51da177e4SLinus Torvalds * MIPS floating point support 61da177e4SLinus Torvalds * Copyright (C) 1994-2000 Algorithmics Ltd. 71da177e4SLinus Torvalds * http://www.algor.co.uk 81da177e4SLinus Torvalds * 91da177e4SLinus Torvalds * ######################################################################## 101da177e4SLinus Torvalds * 111da177e4SLinus Torvalds * This program is free software; you can distribute it and/or modify it 121da177e4SLinus Torvalds * under the terms of the GNU General Public License (Version 2) as 131da177e4SLinus Torvalds * published by the Free Software Foundation. 141da177e4SLinus Torvalds * 151da177e4SLinus Torvalds * This program is distributed in the hope it will be useful, but WITHOUT 161da177e4SLinus Torvalds * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 171da177e4SLinus Torvalds * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 181da177e4SLinus Torvalds * for more details. 191da177e4SLinus Torvalds * 201da177e4SLinus Torvalds * You should have received a copy of the GNU General Public License along 211da177e4SLinus Torvalds * with this program; if not, write to the Free Software Foundation, Inc., 221da177e4SLinus Torvalds * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. 231da177e4SLinus Torvalds * 241da177e4SLinus Torvalds * ######################################################################## 251da177e4SLinus Torvalds */ 261da177e4SLinus Torvalds 271da177e4SLinus Torvalds 281da177e4SLinus Torvalds #include "ieee754dp.h" 291da177e4SLinus Torvalds 301da177e4SLinus Torvalds int ieee754dp_class(ieee754dp x) 311da177e4SLinus Torvalds { 321da177e4SLinus Torvalds COMPXDP; 331da177e4SLinus Torvalds EXPLODEXDP; 341da177e4SLinus Torvalds return xc; 351da177e4SLinus Torvalds } 361da177e4SLinus Torvalds 371da177e4SLinus Torvalds int ieee754dp_isnan(ieee754dp x) 381da177e4SLinus Torvalds { 391da177e4SLinus Torvalds return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 401da177e4SLinus Torvalds } 411da177e4SLinus Torvalds 421da177e4SLinus Torvalds int ieee754dp_issnan(ieee754dp x) 431da177e4SLinus Torvalds { 441da177e4SLinus Torvalds assert(ieee754dp_isnan(x)); 451da177e4SLinus Torvalds return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1)); 461da177e4SLinus Torvalds } 471da177e4SLinus Torvalds 481da177e4SLinus Torvalds 491da177e4SLinus Torvalds ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...) 501da177e4SLinus Torvalds { 511da177e4SLinus Torvalds struct ieee754xctx ax; 521da177e4SLinus Torvalds if (!TSTX()) 531da177e4SLinus Torvalds return r; 541da177e4SLinus Torvalds 551da177e4SLinus Torvalds ax.op = op; 561da177e4SLinus Torvalds ax.rt = IEEE754_RT_DP; 571da177e4SLinus Torvalds ax.rv.dp = r; 581da177e4SLinus Torvalds va_start(ax.ap, op); 591da177e4SLinus Torvalds ieee754_xcpt(&ax); 601da177e4SLinus Torvalds return ax.rv.dp; 611da177e4SLinus Torvalds } 621da177e4SLinus Torvalds 631da177e4SLinus Torvalds ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...) 641da177e4SLinus Torvalds { 651da177e4SLinus Torvalds struct ieee754xctx ax; 661da177e4SLinus Torvalds 671da177e4SLinus Torvalds assert(ieee754dp_isnan(r)); 681da177e4SLinus Torvalds 691da177e4SLinus Torvalds if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ 701da177e4SLinus Torvalds return r; 711da177e4SLinus Torvalds 721da177e4SLinus Torvalds if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) { 731da177e4SLinus Torvalds /* not enabled convert to a quiet NaN */ 741da177e4SLinus Torvalds DPMANT(r) &= (~DP_MBIT(DP_MBITS-1)); 751da177e4SLinus Torvalds if (ieee754dp_isnan(r)) 761da177e4SLinus Torvalds return r; 771da177e4SLinus Torvalds else 781da177e4SLinus Torvalds return ieee754dp_indef(); 791da177e4SLinus Torvalds } 801da177e4SLinus Torvalds 811da177e4SLinus Torvalds ax.op = op; 821da177e4SLinus Torvalds ax.rt = 0; 831da177e4SLinus Torvalds ax.rv.dp = r; 841da177e4SLinus Torvalds va_start(ax.ap, op); 851da177e4SLinus Torvalds ieee754_xcpt(&ax); 861da177e4SLinus Torvalds return ax.rv.dp; 871da177e4SLinus Torvalds } 881da177e4SLinus Torvalds 891da177e4SLinus Torvalds ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y) 901da177e4SLinus Torvalds { 911da177e4SLinus Torvalds assert(ieee754dp_isnan(x)); 921da177e4SLinus Torvalds assert(ieee754dp_isnan(y)); 931da177e4SLinus Torvalds 941da177e4SLinus Torvalds if (DPMANT(x) > DPMANT(y)) 951da177e4SLinus Torvalds return x; 961da177e4SLinus Torvalds else 971da177e4SLinus Torvalds return y; 981da177e4SLinus Torvalds } 991da177e4SLinus Torvalds 1001da177e4SLinus Torvalds 1011da177e4SLinus Torvalds static u64 get_rounding(int sn, u64 xm) 1021da177e4SLinus Torvalds { 1031da177e4SLinus Torvalds /* inexact must round of 3 bits 1041da177e4SLinus Torvalds */ 1051da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 1061da177e4SLinus Torvalds switch (ieee754_csr.rm) { 1071da177e4SLinus Torvalds case IEEE754_RZ: 1081da177e4SLinus Torvalds break; 1091da177e4SLinus Torvalds case IEEE754_RN: 1101da177e4SLinus Torvalds xm += 0x3 + ((xm >> 3) & 1); 1111da177e4SLinus Torvalds /* xm += (xm&0x8)?0x4:0x3 */ 1121da177e4SLinus Torvalds break; 1131da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1141da177e4SLinus Torvalds if (!sn) /* ?? */ 1151da177e4SLinus Torvalds xm += 0x8; 1161da177e4SLinus Torvalds break; 1171da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1181da177e4SLinus Torvalds if (sn) /* ?? */ 1191da177e4SLinus Torvalds xm += 0x8; 1201da177e4SLinus Torvalds break; 1211da177e4SLinus Torvalds } 1221da177e4SLinus Torvalds } 1231da177e4SLinus Torvalds return xm; 1241da177e4SLinus Torvalds } 1251da177e4SLinus Torvalds 1261da177e4SLinus Torvalds 1271da177e4SLinus Torvalds /* generate a normal/denormal number with over,under handling 1281da177e4SLinus Torvalds * sn is sign 1291da177e4SLinus Torvalds * xe is an unbiased exponent 1301da177e4SLinus Torvalds * xm is 3bit extended precision value. 1311da177e4SLinus Torvalds */ 1321da177e4SLinus Torvalds ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 1331da177e4SLinus Torvalds { 1341da177e4SLinus Torvalds assert(xm); /* we don't gen exact zeros (probably should) */ 1351da177e4SLinus Torvalds 1361da177e4SLinus Torvalds assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */ 1371da177e4SLinus Torvalds assert(xm & (DP_HIDDEN_BIT << 3)); 1381da177e4SLinus Torvalds 1391da177e4SLinus Torvalds if (xe < DP_EMIN) { 1401da177e4SLinus Torvalds /* strip lower bits */ 1411da177e4SLinus Torvalds int es = DP_EMIN - xe; 1421da177e4SLinus Torvalds 1431da177e4SLinus Torvalds if (ieee754_csr.nod) { 1441da177e4SLinus Torvalds SETCX(IEEE754_UNDERFLOW); 1451da177e4SLinus Torvalds SETCX(IEEE754_INEXACT); 1461da177e4SLinus Torvalds 1471da177e4SLinus Torvalds switch(ieee754_csr.rm) { 1481da177e4SLinus Torvalds case IEEE754_RN: 1491da177e4SLinus Torvalds return ieee754dp_zero(sn); 1501da177e4SLinus Torvalds case IEEE754_RZ: 1511da177e4SLinus Torvalds return ieee754dp_zero(sn); 1521da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1531da177e4SLinus Torvalds if(sn == 0) 1541da177e4SLinus Torvalds return ieee754dp_min(0); 1551da177e4SLinus Torvalds else 1561da177e4SLinus Torvalds return ieee754dp_zero(1); 1571da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1581da177e4SLinus Torvalds if(sn == 0) 1591da177e4SLinus Torvalds return ieee754dp_zero(0); 1601da177e4SLinus Torvalds else 1611da177e4SLinus Torvalds return ieee754dp_min(1); 1621da177e4SLinus Torvalds } 1631da177e4SLinus Torvalds } 1641da177e4SLinus Torvalds 1651da177e4SLinus Torvalds if (xe == DP_EMIN - 1 1661da177e4SLinus Torvalds && get_rounding(sn, xm) >> (DP_MBITS + 1 + 3)) 1671da177e4SLinus Torvalds { 1681da177e4SLinus Torvalds /* Not tiny after rounding */ 1691da177e4SLinus Torvalds SETCX(IEEE754_INEXACT); 1701da177e4SLinus Torvalds xm = get_rounding(sn, xm); 1711da177e4SLinus Torvalds xm >>= 1; 1721da177e4SLinus Torvalds /* Clear grs bits */ 1731da177e4SLinus Torvalds xm &= ~(DP_MBIT(3) - 1); 1741da177e4SLinus Torvalds xe++; 1751da177e4SLinus Torvalds } 1761da177e4SLinus Torvalds else { 1771da177e4SLinus Torvalds /* sticky right shift es bits 1781da177e4SLinus Torvalds */ 1791da177e4SLinus Torvalds xm = XDPSRS(xm, es); 1801da177e4SLinus Torvalds xe += es; 1811da177e4SLinus Torvalds assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 1821da177e4SLinus Torvalds assert(xe == DP_EMIN); 1831da177e4SLinus Torvalds } 1841da177e4SLinus Torvalds } 1851da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 1861da177e4SLinus Torvalds SETCX(IEEE754_INEXACT); 1871da177e4SLinus Torvalds if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 1881da177e4SLinus Torvalds SETCX(IEEE754_UNDERFLOW); 1891da177e4SLinus Torvalds } 1901da177e4SLinus Torvalds 1911da177e4SLinus Torvalds /* inexact must round of 3 bits 1921da177e4SLinus Torvalds */ 1931da177e4SLinus Torvalds xm = get_rounding(sn, xm); 1941da177e4SLinus Torvalds /* adjust exponent for rounding add overflowing 1951da177e4SLinus Torvalds */ 1961da177e4SLinus Torvalds if (xm >> (DP_MBITS + 3 + 1)) { 1971da177e4SLinus Torvalds /* add causes mantissa overflow */ 1981da177e4SLinus Torvalds xm >>= 1; 1991da177e4SLinus Torvalds xe++; 2001da177e4SLinus Torvalds } 2011da177e4SLinus Torvalds } 2021da177e4SLinus Torvalds /* strip grs bits */ 2031da177e4SLinus Torvalds xm >>= 3; 2041da177e4SLinus Torvalds 2051da177e4SLinus Torvalds assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 2061da177e4SLinus Torvalds assert(xe >= DP_EMIN); 2071da177e4SLinus Torvalds 2081da177e4SLinus Torvalds if (xe > DP_EMAX) { 2091da177e4SLinus Torvalds SETCX(IEEE754_OVERFLOW); 2101da177e4SLinus Torvalds SETCX(IEEE754_INEXACT); 2111da177e4SLinus Torvalds /* -O can be table indexed by (rm,sn) */ 2121da177e4SLinus Torvalds switch (ieee754_csr.rm) { 2131da177e4SLinus Torvalds case IEEE754_RN: 2141da177e4SLinus Torvalds return ieee754dp_inf(sn); 2151da177e4SLinus Torvalds case IEEE754_RZ: 2161da177e4SLinus Torvalds return ieee754dp_max(sn); 2171da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 2181da177e4SLinus Torvalds if (sn == 0) 2191da177e4SLinus Torvalds return ieee754dp_inf(0); 2201da177e4SLinus Torvalds else 2211da177e4SLinus Torvalds return ieee754dp_max(1); 2221da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 2231da177e4SLinus Torvalds if (sn == 0) 2241da177e4SLinus Torvalds return ieee754dp_max(0); 2251da177e4SLinus Torvalds else 2261da177e4SLinus Torvalds return ieee754dp_inf(1); 2271da177e4SLinus Torvalds } 2281da177e4SLinus Torvalds } 2291da177e4SLinus Torvalds /* gen norm/denorm/zero */ 2301da177e4SLinus Torvalds 2311da177e4SLinus Torvalds if ((xm & DP_HIDDEN_BIT) == 0) { 2321da177e4SLinus Torvalds /* we underflow (tiny/zero) */ 2331da177e4SLinus Torvalds assert(xe == DP_EMIN); 2341da177e4SLinus Torvalds if (ieee754_csr.mx & IEEE754_UNDERFLOW) 2351da177e4SLinus Torvalds SETCX(IEEE754_UNDERFLOW); 2361da177e4SLinus Torvalds return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 2371da177e4SLinus Torvalds } else { 2381da177e4SLinus Torvalds assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */ 2391da177e4SLinus Torvalds assert(xm & DP_HIDDEN_BIT); 2401da177e4SLinus Torvalds 2411da177e4SLinus Torvalds return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 2421da177e4SLinus Torvalds } 2431da177e4SLinus Torvalds } 244