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 26cae55066SRalf Baechle #include <linux/compiler.h> 271da177e4SLinus Torvalds 281da177e4SLinus Torvalds #include "ieee754dp.h" 291da177e4SLinus Torvalds 302209bcb1SRalf Baechle int ieee754dp_class(union ieee754dp x) 311da177e4SLinus Torvalds { 321da177e4SLinus Torvalds COMPXDP; 331da177e4SLinus Torvalds EXPLODEXDP; 341da177e4SLinus Torvalds return xc; 351da177e4SLinus Torvalds } 361da177e4SLinus Torvalds 372209bcb1SRalf Baechle int ieee754dp_isnan(union ieee754dp x) 381da177e4SLinus Torvalds { 391da177e4SLinus Torvalds return ieee754dp_class(x) >= IEEE754_CLASS_SNAN; 401da177e4SLinus Torvalds } 411da177e4SLinus Torvalds 42f71baa11SRalf Baechle static inline int ieee754dp_issnan(union ieee754dp x) 431da177e4SLinus Torvalds { 441da177e4SLinus Torvalds assert(ieee754dp_isnan(x)); 45ad8fb553SRalf Baechle return ((DPMANT(x) & DP_MBIT(DP_FBITS-1)) == DP_MBIT(DP_FBITS-1)); 461da177e4SLinus Torvalds } 471da177e4SLinus Torvalds 481da177e4SLinus Torvalds 4990efba36SRalf Baechle union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r) 501da177e4SLinus Torvalds { 511da177e4SLinus Torvalds assert(ieee754dp_isnan(r)); 521da177e4SLinus Torvalds 531da177e4SLinus Torvalds if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */ 541da177e4SLinus Torvalds return r; 551da177e4SLinus Torvalds 569e8bad1fSRalf Baechle if (!ieee754_setandtestcx(IEEE754_INVALID_OPERATION)) { 571da177e4SLinus Torvalds /* not enabled convert to a quiet NaN */ 58ad8fb553SRalf Baechle DPMANT(r) &= (~DP_MBIT(DP_FBITS-1)); 591da177e4SLinus Torvalds if (ieee754dp_isnan(r)) 601da177e4SLinus Torvalds return r; 611da177e4SLinus Torvalds else 621da177e4SLinus Torvalds return ieee754dp_indef(); 631da177e4SLinus Torvalds } 641da177e4SLinus Torvalds 6590efba36SRalf Baechle return r; 661da177e4SLinus Torvalds } 671da177e4SLinus Torvalds 68de2fc342SRalf Baechle static u64 ieee754dp_get_rounding(int sn, u64 xm) 691da177e4SLinus Torvalds { 701da177e4SLinus Torvalds /* inexact must round of 3 bits 711da177e4SLinus Torvalds */ 721da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 731da177e4SLinus Torvalds switch (ieee754_csr.rm) { 741da177e4SLinus Torvalds case IEEE754_RZ: 751da177e4SLinus Torvalds break; 761da177e4SLinus Torvalds case IEEE754_RN: 771da177e4SLinus Torvalds xm += 0x3 + ((xm >> 3) & 1); 781da177e4SLinus Torvalds /* xm += (xm&0x8)?0x4:0x3 */ 791da177e4SLinus Torvalds break; 801da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 811da177e4SLinus Torvalds if (!sn) /* ?? */ 821da177e4SLinus Torvalds xm += 0x8; 831da177e4SLinus Torvalds break; 841da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 851da177e4SLinus Torvalds if (sn) /* ?? */ 861da177e4SLinus Torvalds xm += 0x8; 871da177e4SLinus Torvalds break; 881da177e4SLinus Torvalds } 891da177e4SLinus Torvalds } 901da177e4SLinus Torvalds return xm; 911da177e4SLinus Torvalds } 921da177e4SLinus Torvalds 931da177e4SLinus Torvalds 941da177e4SLinus Torvalds /* generate a normal/denormal number with over,under handling 951da177e4SLinus Torvalds * sn is sign 961da177e4SLinus Torvalds * xe is an unbiased exponent 971da177e4SLinus Torvalds * xm is 3bit extended precision value. 981da177e4SLinus Torvalds */ 992209bcb1SRalf Baechle union ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 1001da177e4SLinus Torvalds { 1011da177e4SLinus Torvalds assert(xm); /* we don't gen exact zeros (probably should) */ 1021da177e4SLinus Torvalds 103ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no execess */ 1041da177e4SLinus Torvalds assert(xm & (DP_HIDDEN_BIT << 3)); 1051da177e4SLinus Torvalds 1061da177e4SLinus Torvalds if (xe < DP_EMIN) { 1071da177e4SLinus Torvalds /* strip lower bits */ 1081da177e4SLinus Torvalds int es = DP_EMIN - xe; 1091da177e4SLinus Torvalds 1101da177e4SLinus Torvalds if (ieee754_csr.nod) { 1119e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1129e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1131da177e4SLinus Torvalds 1141da177e4SLinus Torvalds switch(ieee754_csr.rm) { 1151da177e4SLinus Torvalds case IEEE754_RN: 1161da177e4SLinus Torvalds case IEEE754_RZ: 1171da177e4SLinus Torvalds return ieee754dp_zero(sn); 1181da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1191da177e4SLinus Torvalds if (sn == 0) 1201da177e4SLinus Torvalds return ieee754dp_min(0); 1211da177e4SLinus Torvalds else 1221da177e4SLinus Torvalds return ieee754dp_zero(1); 1231da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1241da177e4SLinus Torvalds if (sn == 0) 1251da177e4SLinus Torvalds return ieee754dp_zero(0); 1261da177e4SLinus Torvalds else 1271da177e4SLinus Torvalds return ieee754dp_min(1); 1281da177e4SLinus Torvalds } 1291da177e4SLinus Torvalds } 1301da177e4SLinus Torvalds 131de2fc342SRalf Baechle if (xe == DP_EMIN - 1 && 132de2fc342SRalf Baechle ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3)) 1331da177e4SLinus Torvalds { 1341da177e4SLinus Torvalds /* Not tiny after rounding */ 1359e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 136de2fc342SRalf Baechle xm = ieee754dp_get_rounding(sn, xm); 1371da177e4SLinus Torvalds xm >>= 1; 1381da177e4SLinus Torvalds /* Clear grs bits */ 1391da177e4SLinus Torvalds xm &= ~(DP_MBIT(3) - 1); 1401da177e4SLinus Torvalds xe++; 1411da177e4SLinus Torvalds } 1421da177e4SLinus Torvalds else { 1431da177e4SLinus Torvalds /* sticky right shift es bits 1441da177e4SLinus Torvalds */ 1451da177e4SLinus Torvalds xm = XDPSRS(xm, es); 1461da177e4SLinus Torvalds xe += es; 1471da177e4SLinus Torvalds assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 1481da177e4SLinus Torvalds assert(xe == DP_EMIN); 1491da177e4SLinus Torvalds } 1501da177e4SLinus Torvalds } 1511da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 1529e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1531da177e4SLinus Torvalds if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 1549e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1551da177e4SLinus Torvalds } 1561da177e4SLinus Torvalds 1571da177e4SLinus Torvalds /* inexact must round of 3 bits 1581da177e4SLinus Torvalds */ 159de2fc342SRalf Baechle xm = ieee754dp_get_rounding(sn, xm); 1601da177e4SLinus Torvalds /* adjust exponent for rounding add overflowing 1611da177e4SLinus Torvalds */ 162ad8fb553SRalf Baechle if (xm >> (DP_FBITS + 3 + 1)) { 1631da177e4SLinus Torvalds /* add causes mantissa overflow */ 1641da177e4SLinus Torvalds xm >>= 1; 1651da177e4SLinus Torvalds xe++; 1661da177e4SLinus Torvalds } 1671da177e4SLinus Torvalds } 1681da177e4SLinus Torvalds /* strip grs bits */ 1691da177e4SLinus Torvalds xm >>= 3; 1701da177e4SLinus Torvalds 171ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 1721da177e4SLinus Torvalds assert(xe >= DP_EMIN); 1731da177e4SLinus Torvalds 1741da177e4SLinus Torvalds if (xe > DP_EMAX) { 1759e8bad1fSRalf Baechle ieee754_setcx(IEEE754_OVERFLOW); 1769e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1771da177e4SLinus Torvalds /* -O can be table indexed by (rm,sn) */ 1781da177e4SLinus Torvalds switch (ieee754_csr.rm) { 1791da177e4SLinus Torvalds case IEEE754_RN: 1801da177e4SLinus Torvalds return ieee754dp_inf(sn); 1811da177e4SLinus Torvalds case IEEE754_RZ: 1821da177e4SLinus Torvalds return ieee754dp_max(sn); 1831da177e4SLinus Torvalds case IEEE754_RU: /* toward +Infinity */ 1841da177e4SLinus Torvalds if (sn == 0) 1851da177e4SLinus Torvalds return ieee754dp_inf(0); 1861da177e4SLinus Torvalds else 1871da177e4SLinus Torvalds return ieee754dp_max(1); 1881da177e4SLinus Torvalds case IEEE754_RD: /* toward -Infinity */ 1891da177e4SLinus Torvalds if (sn == 0) 1901da177e4SLinus Torvalds return ieee754dp_max(0); 1911da177e4SLinus Torvalds else 1921da177e4SLinus Torvalds return ieee754dp_inf(1); 1931da177e4SLinus Torvalds } 1941da177e4SLinus Torvalds } 1951da177e4SLinus Torvalds /* gen norm/denorm/zero */ 1961da177e4SLinus Torvalds 1971da177e4SLinus Torvalds if ((xm & DP_HIDDEN_BIT) == 0) { 1981da177e4SLinus Torvalds /* we underflow (tiny/zero) */ 1991da177e4SLinus Torvalds assert(xe == DP_EMIN); 2001da177e4SLinus Torvalds if (ieee754_csr.mx & IEEE754_UNDERFLOW) 2019e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 2021da177e4SLinus Torvalds return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 2031da177e4SLinus Torvalds } else { 204ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 2051da177e4SLinus Torvalds assert(xm & DP_HIDDEN_BIT); 2061da177e4SLinus Torvalds 2071da177e4SLinus Torvalds return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 2081da177e4SLinus Torvalds } 2091da177e4SLinus Torvalds } 210