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 * 81da177e4SLinus Torvalds * ######################################################################## 91da177e4SLinus Torvalds * 101da177e4SLinus Torvalds * This program is free software; you can distribute it and/or modify it 111da177e4SLinus Torvalds * under the terms of the GNU General Public License (Version 2) as 121da177e4SLinus Torvalds * published by the Free Software Foundation. 131da177e4SLinus Torvalds * 141da177e4SLinus Torvalds * This program is distributed in the hope it will be useful, but WITHOUT 151da177e4SLinus Torvalds * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 161da177e4SLinus Torvalds * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 171da177e4SLinus Torvalds * for more details. 181da177e4SLinus Torvalds * 191da177e4SLinus Torvalds * You should have received a copy of the GNU General Public License along 201da177e4SLinus Torvalds * with this program; if not, write to the Free Software Foundation, Inc., 211da177e4SLinus Torvalds * 59 Temple Place - Suite 330, Boston MA 02111-1307, USA. 221da177e4SLinus Torvalds * 231da177e4SLinus Torvalds * ######################################################################## 241da177e4SLinus Torvalds */ 251da177e4SLinus Torvalds 261eaf198dSRalf Baechle #include <stdarg.h> 27cae55066SRalf Baechle #include <linux/compiler.h> 281da177e4SLinus Torvalds 291da177e4SLinus Torvalds #include "ieee754dp.h" 301da177e4SLinus Torvalds 312209bcb1SRalf Baechle int ieee754dp_class(union ieee754dp x) 321da177e4SLinus Torvalds { 331da177e4SLinus Torvalds COMPXDP; 341da177e4SLinus Torvalds EXPLODEXDP; 351da177e4SLinus Torvalds return xc; 361da177e4SLinus Torvalds } 371da177e4SLinus Torvalds 382209bcb1SRalf Baechle int ieee754dp_isnan(union ieee754dp x) 391da177e4SLinus Torvalds { 401da177e4SLinus Torvalds return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 411da177e4SLinus Torvalds } 421da177e4SLinus Torvalds 43f71baa11SRalf Baechle static inline int ieee754dp_issnan(union ieee754dp x) 441da177e4SLinus Torvalds { 451da177e4SLinus Torvalds assert(ieee754dp_isnan(x)); 46ad8fb553SRalf Baechle return ((DPMANT(x) & DP_MBIT(DP_FBITS-1)) == DP_MBIT(DP_FBITS-1)); 471da177e4SLinus Torvalds } 481da177e4SLinus Torvalds 491da177e4SLinus Torvalds 502209bcb1SRalf Baechle union ieee754dp __cold ieee754dp_xcpt(union ieee754dp r, const char *op, ...) 511da177e4SLinus Torvalds { 521da177e4SLinus Torvalds struct ieee754xctx ax; 539e8bad1fSRalf Baechle if (!ieee754_tstx()) 541da177e4SLinus Torvalds return r; 551da177e4SLinus Torvalds 561da177e4SLinus Torvalds ax.op = op; 571da177e4SLinus Torvalds ax.rt = IEEE754_RT_DP; 581da177e4SLinus Torvalds ax.rv.dp = r; 591da177e4SLinus Torvalds va_start(ax.ap, op); 601da177e4SLinus Torvalds ieee754_xcpt(&ax); 618142294dSRichard Knutsson va_end(ax.ap); 621da177e4SLinus Torvalds return ax.rv.dp; 631da177e4SLinus Torvalds } 641da177e4SLinus Torvalds 652209bcb1SRalf Baechle union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r, const char *op, ...) 661da177e4SLinus Torvalds { 671da177e4SLinus Torvalds struct ieee754xctx ax; 681da177e4SLinus Torvalds 691da177e4SLinus Torvalds assert(ieee754dp_isnan(r)); 701da177e4SLinus Torvalds 711da177e4SLinus Torvalds if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ 721da177e4SLinus Torvalds return r; 731da177e4SLinus Torvalds 749e8bad1fSRalf Baechle if (!ieee754_setandtestcx(IEEE754_INVALID_OPERATION)) { 751da177e4SLinus Torvalds /* not enabled convert to a quiet NaN */ 76ad8fb553SRalf Baechle DPMANT(r) &= (~DP_MBIT(DP_FBITS-1)); 771da177e4SLinus Torvalds if (ieee754dp_isnan(r)) 781da177e4SLinus Torvalds return r; 791da177e4SLinus Torvalds else 801da177e4SLinus Torvalds return ieee754dp_indef(); 811da177e4SLinus Torvalds } 821da177e4SLinus Torvalds 831da177e4SLinus Torvalds ax.op = op; 841da177e4SLinus Torvalds ax.rt = 0; 851da177e4SLinus Torvalds ax.rv.dp = r; 861da177e4SLinus Torvalds va_start(ax.ap, op); 871da177e4SLinus Torvalds ieee754_xcpt(&ax); 888142294dSRichard Knutsson va_end(ax.ap); 891da177e4SLinus Torvalds return ax.rv.dp; 901da177e4SLinus Torvalds } 911da177e4SLinus Torvalds 921da177e4SLinus Torvalds static u64 get_rounding(int sn, u64 xm) 931da177e4SLinus Torvalds { 941da177e4SLinus Torvalds /* inexact must round of 3 bits 951da177e4SLinus Torvalds */ 961da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 971da177e4SLinus Torvalds switch (ieee754_csr.rm) { 981da177e4SLinus Torvalds case IEEE754_RZ: 991da177e4SLinus Torvalds break; 1001da177e4SLinus Torvalds case IEEE754_RN: 1011da177e4SLinus Torvalds xm += 0x3 + ((xm >> 3) & 1); 1021da177e4SLinus Torvalds /* xm += (xm&0x8)?0x4:0x3 */ 1031da177e4SLinus Torvalds break; 1041da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1051da177e4SLinus Torvalds if (!sn) /* ?? */ 1061da177e4SLinus Torvalds xm += 0x8; 1071da177e4SLinus Torvalds break; 1081da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1091da177e4SLinus Torvalds if (sn) /* ?? */ 1101da177e4SLinus Torvalds xm += 0x8; 1111da177e4SLinus Torvalds break; 1121da177e4SLinus Torvalds } 1131da177e4SLinus Torvalds } 1141da177e4SLinus Torvalds return xm; 1151da177e4SLinus Torvalds } 1161da177e4SLinus Torvalds 1171da177e4SLinus Torvalds 1181da177e4SLinus Torvalds /* generate a normal/denormal number with over,under handling 1191da177e4SLinus Torvalds * sn is sign 1201da177e4SLinus Torvalds * xe is an unbiased exponent 1211da177e4SLinus Torvalds * xm is 3bit extended precision value. 1221da177e4SLinus Torvalds */ 1232209bcb1SRalf Baechle union ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 1241da177e4SLinus Torvalds { 1251da177e4SLinus Torvalds assert(xm); /* we don't gen exact zeros (probably should) */ 1261da177e4SLinus Torvalds 127ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no execess */ 1281da177e4SLinus Torvalds assert(xm & (DP_HIDDEN_BIT << 3)); 1291da177e4SLinus Torvalds 1301da177e4SLinus Torvalds if (xe < DP_EMIN) { 1311da177e4SLinus Torvalds /* strip lower bits */ 1321da177e4SLinus Torvalds int es = DP_EMIN - xe; 1331da177e4SLinus Torvalds 1341da177e4SLinus Torvalds if (ieee754_csr.nod) { 1359e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1369e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1371da177e4SLinus Torvalds 1381da177e4SLinus Torvalds switch(ieee754_csr.rm) { 1391da177e4SLinus Torvalds case IEEE754_RN: 1401da177e4SLinus Torvalds case IEEE754_RZ: 1411da177e4SLinus Torvalds return ieee754dp_zero(sn); 1421da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1431da177e4SLinus Torvalds if (sn == 0) 1441da177e4SLinus Torvalds return ieee754dp_min(0); 1451da177e4SLinus Torvalds else 1461da177e4SLinus Torvalds return ieee754dp_zero(1); 1471da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1481da177e4SLinus Torvalds if (sn == 0) 1491da177e4SLinus Torvalds return ieee754dp_zero(0); 1501da177e4SLinus Torvalds else 1511da177e4SLinus Torvalds return ieee754dp_min(1); 1521da177e4SLinus Torvalds } 1531da177e4SLinus Torvalds } 1541da177e4SLinus Torvalds 1551da177e4SLinus Torvalds if (xe == DP_EMIN - 1 156ad8fb553SRalf Baechle && get_rounding(sn, xm) >> (DP_FBITS + 1 + 3)) 1571da177e4SLinus Torvalds { 1581da177e4SLinus Torvalds /* Not tiny after rounding */ 1599e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1601da177e4SLinus Torvalds xm = get_rounding(sn, xm); 1611da177e4SLinus Torvalds xm >>= 1; 1621da177e4SLinus Torvalds /* Clear grs bits */ 1631da177e4SLinus Torvalds xm &= ~(DP_MBIT(3) - 1); 1641da177e4SLinus Torvalds xe++; 1651da177e4SLinus Torvalds } 1661da177e4SLinus Torvalds else { 1671da177e4SLinus Torvalds /* sticky right shift es bits 1681da177e4SLinus Torvalds */ 1691da177e4SLinus Torvalds xm = XDPSRS(xm, es); 1701da177e4SLinus Torvalds xe += es; 1711da177e4SLinus Torvalds assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 1721da177e4SLinus Torvalds assert(xe == DP_EMIN); 1731da177e4SLinus Torvalds } 1741da177e4SLinus Torvalds } 1751da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 1769e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1771da177e4SLinus Torvalds if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 1789e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1791da177e4SLinus Torvalds } 1801da177e4SLinus Torvalds 1811da177e4SLinus Torvalds /* inexact must round of 3 bits 1821da177e4SLinus Torvalds */ 1831da177e4SLinus Torvalds xm = get_rounding(sn, xm); 1841da177e4SLinus Torvalds /* adjust exponent for rounding add overflowing 1851da177e4SLinus Torvalds */ 186ad8fb553SRalf Baechle if (xm >> (DP_FBITS + 3 + 1)) { 1871da177e4SLinus Torvalds /* add causes mantissa overflow */ 1881da177e4SLinus Torvalds xm >>= 1; 1891da177e4SLinus Torvalds xe++; 1901da177e4SLinus Torvalds } 1911da177e4SLinus Torvalds } 1921da177e4SLinus Torvalds /* strip grs bits */ 1931da177e4SLinus Torvalds xm >>= 3; 1941da177e4SLinus Torvalds 195ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 1961da177e4SLinus Torvalds assert(xe >= DP_EMIN); 1971da177e4SLinus Torvalds 1981da177e4SLinus Torvalds if (xe > DP_EMAX) { 1999e8bad1fSRalf Baechle ieee754_setcx(IEEE754_OVERFLOW); 2009e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 2011da177e4SLinus Torvalds /* -O can be table indexed by (rm,sn) */ 2021da177e4SLinus Torvalds switch (ieee754_csr.rm) { 2031da177e4SLinus Torvalds case IEEE754_RN: 2041da177e4SLinus Torvalds return ieee754dp_inf(sn); 2051da177e4SLinus Torvalds case IEEE754_RZ: 2061da177e4SLinus Torvalds return ieee754dp_max(sn); 2071da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 2081da177e4SLinus Torvalds if (sn == 0) 2091da177e4SLinus Torvalds return ieee754dp_inf(0); 2101da177e4SLinus Torvalds else 2111da177e4SLinus Torvalds return ieee754dp_max(1); 2121da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 2131da177e4SLinus Torvalds if (sn == 0) 2141da177e4SLinus Torvalds return ieee754dp_max(0); 2151da177e4SLinus Torvalds else 2161da177e4SLinus Torvalds return ieee754dp_inf(1); 2171da177e4SLinus Torvalds } 2181da177e4SLinus Torvalds } 2191da177e4SLinus Torvalds /* gen norm/denorm/zero */ 2201da177e4SLinus Torvalds 2211da177e4SLinus Torvalds if ((xm & DP_HIDDEN_BIT) == 0) { 2221da177e4SLinus Torvalds /* we underflow (tiny/zero) */ 2231da177e4SLinus Torvalds assert(xe == DP_EMIN); 2241da177e4SLinus Torvalds if (ieee754_csr.mx & IEEE754_UNDERFLOW) 2259e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 2261da177e4SLinus Torvalds return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 2271da177e4SLinus Torvalds } else { 228ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 2291da177e4SLinus Torvalds assert(xm & DP_HIDDEN_BIT); 2301da177e4SLinus Torvalds 2311da177e4SLinus Torvalds return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 2321da177e4SLinus Torvalds } 2331da177e4SLinus Torvalds } 234