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 * This program is free software; you can distribute it and/or modify it 91da177e4SLinus Torvalds * under the terms of the GNU General Public License (Version 2) as 101da177e4SLinus Torvalds * published by the Free Software Foundation. 111da177e4SLinus Torvalds * 121da177e4SLinus Torvalds * This program is distributed in the hope it will be useful, but WITHOUT 131da177e4SLinus Torvalds * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 141da177e4SLinus Torvalds * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 151da177e4SLinus Torvalds * for more details. 161da177e4SLinus Torvalds * 171da177e4SLinus Torvalds * You should have received a copy of the GNU General Public License along 181da177e4SLinus Torvalds * with this program; if not, write to the Free Software Foundation, Inc., 193f7cac41SRalf Baechle * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 201da177e4SLinus Torvalds */ 211da177e4SLinus Torvalds 22cae55066SRalf Baechle #include <linux/compiler.h> 231da177e4SLinus Torvalds 241da177e4SLinus Torvalds #include "ieee754dp.h" 251da177e4SLinus Torvalds 262209bcb1SRalf Baechle int ieee754dp_class(union ieee754dp x) 271da177e4SLinus Torvalds { 281da177e4SLinus Torvalds COMPXDP; 291da177e4SLinus Torvalds EXPLODEXDP; 301da177e4SLinus Torvalds return xc; 311da177e4SLinus Torvalds } 321da177e4SLinus Torvalds 332209bcb1SRalf Baechle int ieee754dp_isnan(union ieee754dp x) 341da177e4SLinus Torvalds { 35c9a10845SMaciej W. Rozycki return ieee754_class_nan(ieee754dp_class(x)); 361da177e4SLinus Torvalds } 371da177e4SLinus Torvalds 38f71baa11SRalf Baechle static inline int ieee754dp_issnan(union ieee754dp x) 391da177e4SLinus Torvalds { 401da177e4SLinus Torvalds assert(ieee754dp_isnan(x)); 41635c9907SRalf Baechle return (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1); 421da177e4SLinus Torvalds } 431da177e4SLinus Torvalds 441da177e4SLinus Torvalds 45d5afa7e9SMaciej W. Rozycki /* 46d5afa7e9SMaciej W. Rozycki * Raise the Invalid Operation IEEE 754 exception 47d5afa7e9SMaciej W. Rozycki * and convert the signaling NaN supplied to a quiet NaN. 48d5afa7e9SMaciej W. Rozycki */ 4990efba36SRalf Baechle union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r) 501da177e4SLinus Torvalds { 51d5afa7e9SMaciej W. Rozycki assert(ieee754dp_issnan(r)); 521da177e4SLinus Torvalds 53d5afa7e9SMaciej W. Rozycki ieee754_setcx(IEEE754_INVALID_OPERATION); 541da177e4SLinus Torvalds return ieee754dp_indef(); 551da177e4SLinus Torvalds } 561da177e4SLinus Torvalds 57de2fc342SRalf Baechle static u64 ieee754dp_get_rounding(int sn, u64 xm) 581da177e4SLinus Torvalds { 591da177e4SLinus Torvalds /* inexact must round of 3 bits 601da177e4SLinus Torvalds */ 611da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 621da177e4SLinus Torvalds switch (ieee754_csr.rm) { 6356a64733SRalf Baechle case FPU_CSR_RZ: 641da177e4SLinus Torvalds break; 6556a64733SRalf Baechle case FPU_CSR_RN: 661da177e4SLinus Torvalds xm += 0x3 + ((xm >> 3) & 1); 671da177e4SLinus Torvalds /* xm += (xm&0x8)?0x4:0x3 */ 681da177e4SLinus Torvalds break; 6956a64733SRalf Baechle case FPU_CSR_RU: /* toward +Infinity */ 701da177e4SLinus Torvalds if (!sn) /* ?? */ 711da177e4SLinus Torvalds xm += 0x8; 721da177e4SLinus Torvalds break; 7356a64733SRalf Baechle case FPU_CSR_RD: /* toward -Infinity */ 741da177e4SLinus Torvalds if (sn) /* ?? */ 751da177e4SLinus Torvalds xm += 0x8; 761da177e4SLinus Torvalds break; 771da177e4SLinus Torvalds } 781da177e4SLinus Torvalds } 791da177e4SLinus Torvalds return xm; 801da177e4SLinus Torvalds } 811da177e4SLinus Torvalds 821da177e4SLinus Torvalds 831da177e4SLinus Torvalds /* generate a normal/denormal number with over,under handling 841da177e4SLinus Torvalds * sn is sign 851da177e4SLinus Torvalds * xe is an unbiased exponent 861da177e4SLinus Torvalds * xm is 3bit extended precision value. 871da177e4SLinus Torvalds */ 882209bcb1SRalf Baechle union ieee754dp ieee754dp_format(int sn, int xe, u64 xm) 891da177e4SLinus Torvalds { 901da177e4SLinus Torvalds assert(xm); /* we don't gen exact zeros (probably should) */ 911da177e4SLinus Torvalds 92ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no execess */ 931da177e4SLinus Torvalds assert(xm & (DP_HIDDEN_BIT << 3)); 941da177e4SLinus Torvalds 951da177e4SLinus Torvalds if (xe < DP_EMIN) { 961da177e4SLinus Torvalds /* strip lower bits */ 971da177e4SLinus Torvalds int es = DP_EMIN - xe; 981da177e4SLinus Torvalds 991da177e4SLinus Torvalds if (ieee754_csr.nod) { 1009e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1019e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1021da177e4SLinus Torvalds 1031da177e4SLinus Torvalds switch(ieee754_csr.rm) { 10456a64733SRalf Baechle case FPU_CSR_RN: 10556a64733SRalf Baechle case FPU_CSR_RZ: 1061da177e4SLinus Torvalds return ieee754dp_zero(sn); 10756a64733SRalf Baechle case FPU_CSR_RU: /* toward +Infinity */ 1081da177e4SLinus Torvalds if (sn == 0) 1091da177e4SLinus Torvalds return ieee754dp_min(0); 1101da177e4SLinus Torvalds else 1111da177e4SLinus Torvalds return ieee754dp_zero(1); 11256a64733SRalf Baechle case FPU_CSR_RD: /* toward -Infinity */ 1131da177e4SLinus Torvalds if (sn == 0) 1141da177e4SLinus Torvalds return ieee754dp_zero(0); 1151da177e4SLinus Torvalds else 1161da177e4SLinus Torvalds return ieee754dp_min(1); 1171da177e4SLinus Torvalds } 1181da177e4SLinus Torvalds } 1191da177e4SLinus Torvalds 120de2fc342SRalf Baechle if (xe == DP_EMIN - 1 && 121de2fc342SRalf Baechle ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3)) 1221da177e4SLinus Torvalds { 1231da177e4SLinus Torvalds /* Not tiny after rounding */ 1249e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 125de2fc342SRalf Baechle xm = ieee754dp_get_rounding(sn, xm); 1261da177e4SLinus Torvalds xm >>= 1; 1271da177e4SLinus Torvalds /* Clear grs bits */ 1281da177e4SLinus Torvalds xm &= ~(DP_MBIT(3) - 1); 1291da177e4SLinus Torvalds xe++; 1301da177e4SLinus Torvalds } 1311da177e4SLinus Torvalds else { 1321da177e4SLinus Torvalds /* sticky right shift es bits 1331da177e4SLinus Torvalds */ 1341da177e4SLinus Torvalds xm = XDPSRS(xm, es); 1351da177e4SLinus Torvalds xe += es; 1361da177e4SLinus Torvalds assert((xm & (DP_HIDDEN_BIT << 3)) == 0); 1371da177e4SLinus Torvalds assert(xe == DP_EMIN); 1381da177e4SLinus Torvalds } 1391da177e4SLinus Torvalds } 1401da177e4SLinus Torvalds if (xm & (DP_MBIT(3) - 1)) { 1419e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1421da177e4SLinus Torvalds if ((xm & (DP_HIDDEN_BIT << 3)) == 0) { 1439e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1441da177e4SLinus Torvalds } 1451da177e4SLinus Torvalds 1461da177e4SLinus Torvalds /* inexact must round of 3 bits 1471da177e4SLinus Torvalds */ 148de2fc342SRalf Baechle xm = ieee754dp_get_rounding(sn, xm); 1491da177e4SLinus Torvalds /* adjust exponent for rounding add overflowing 1501da177e4SLinus Torvalds */ 151ad8fb553SRalf Baechle if (xm >> (DP_FBITS + 3 + 1)) { 1521da177e4SLinus Torvalds /* add causes mantissa overflow */ 1531da177e4SLinus Torvalds xm >>= 1; 1541da177e4SLinus Torvalds xe++; 1551da177e4SLinus Torvalds } 1561da177e4SLinus Torvalds } 1571da177e4SLinus Torvalds /* strip grs bits */ 1581da177e4SLinus Torvalds xm >>= 3; 1591da177e4SLinus Torvalds 160ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 1611da177e4SLinus Torvalds assert(xe >= DP_EMIN); 1621da177e4SLinus Torvalds 1631da177e4SLinus Torvalds if (xe > DP_EMAX) { 1649e8bad1fSRalf Baechle ieee754_setcx(IEEE754_OVERFLOW); 1659e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 1661da177e4SLinus Torvalds /* -O can be table indexed by (rm,sn) */ 1671da177e4SLinus Torvalds switch (ieee754_csr.rm) { 16856a64733SRalf Baechle case FPU_CSR_RN: 1691da177e4SLinus Torvalds return ieee754dp_inf(sn); 17056a64733SRalf Baechle case FPU_CSR_RZ: 1711da177e4SLinus Torvalds return ieee754dp_max(sn); 17256a64733SRalf Baechle case FPU_CSR_RU: /* toward +Infinity */ 1731da177e4SLinus Torvalds if (sn == 0) 1741da177e4SLinus Torvalds return ieee754dp_inf(0); 1751da177e4SLinus Torvalds else 1761da177e4SLinus Torvalds return ieee754dp_max(1); 17756a64733SRalf Baechle case FPU_CSR_RD: /* toward -Infinity */ 1781da177e4SLinus Torvalds if (sn == 0) 1791da177e4SLinus Torvalds return ieee754dp_max(0); 1801da177e4SLinus Torvalds else 1811da177e4SLinus Torvalds return ieee754dp_inf(1); 1821da177e4SLinus Torvalds } 1831da177e4SLinus Torvalds } 1841da177e4SLinus Torvalds /* gen norm/denorm/zero */ 1851da177e4SLinus Torvalds 1861da177e4SLinus Torvalds if ((xm & DP_HIDDEN_BIT) == 0) { 1871da177e4SLinus Torvalds /* we underflow (tiny/zero) */ 1881da177e4SLinus Torvalds assert(xe == DP_EMIN); 1891da177e4SLinus Torvalds if (ieee754_csr.mx & IEEE754_UNDERFLOW) 1909e8bad1fSRalf Baechle ieee754_setcx(IEEE754_UNDERFLOW); 1911da177e4SLinus Torvalds return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm); 1921da177e4SLinus Torvalds } else { 193ad8fb553SRalf Baechle assert((xm >> (DP_FBITS + 1)) == 0); /* no execess */ 1941da177e4SLinus Torvalds assert(xm & DP_HIDDEN_BIT); 1951da177e4SLinus Torvalds 1961da177e4SLinus Torvalds return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); 1971da177e4SLinus Torvalds } 1981da177e4SLinus Torvalds } 199