xref: /openbmc/linux/arch/mips/math-emu/ieee754dp.c (revision 9d5a6349)
19d5a6349SThomas Gleixner // SPDX-License-Identifier: GPL-2.0-only
21da177e4SLinus Torvalds /* IEEE754 floating point arithmetic
31da177e4SLinus Torvalds  * double precision: common utilities
41da177e4SLinus Torvalds  */
51da177e4SLinus Torvalds /*
61da177e4SLinus Torvalds  * MIPS floating point support
71da177e4SLinus Torvalds  * Copyright (C) 1994-2000 Algorithmics Ltd.
81da177e4SLinus Torvalds  */
91da177e4SLinus Torvalds 
10cae55066SRalf Baechle #include <linux/compiler.h>
111da177e4SLinus Torvalds 
121da177e4SLinus Torvalds #include "ieee754dp.h"
131da177e4SLinus Torvalds 
ieee754dp_class(union ieee754dp x)142209bcb1SRalf Baechle int ieee754dp_class(union ieee754dp x)
151da177e4SLinus Torvalds {
161da177e4SLinus Torvalds 	COMPXDP;
171da177e4SLinus Torvalds 	EXPLODEXDP;
181da177e4SLinus Torvalds 	return xc;
191da177e4SLinus Torvalds }
201da177e4SLinus Torvalds 
ieee754dp_isnan(union ieee754dp x)21e06b530bSMaciej W. Rozycki static inline int ieee754dp_isnan(union ieee754dp x)
221da177e4SLinus Torvalds {
23c9a10845SMaciej W. Rozycki 	return ieee754_class_nan(ieee754dp_class(x));
241da177e4SLinus Torvalds }
251da177e4SLinus Torvalds 
ieee754dp_issnan(union ieee754dp x)26f71baa11SRalf Baechle static inline int ieee754dp_issnan(union ieee754dp x)
271da177e4SLinus Torvalds {
2890d53a91SMaciej W. Rozycki 	int qbit;
2990d53a91SMaciej W. Rozycki 
301da177e4SLinus Torvalds 	assert(ieee754dp_isnan(x));
3190d53a91SMaciej W. Rozycki 	qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1);
3290d53a91SMaciej W. Rozycki 	return ieee754_csr.nan2008 ^ qbit;
331da177e4SLinus Torvalds }
341da177e4SLinus Torvalds 
351da177e4SLinus Torvalds 
36d5afa7e9SMaciej W. Rozycki /*
37d5afa7e9SMaciej W. Rozycki  * Raise the Invalid Operation IEEE 754 exception
38d5afa7e9SMaciej W. Rozycki  * and convert the signaling NaN supplied to a quiet NaN.
39d5afa7e9SMaciej W. Rozycki  */
ieee754dp_nanxcpt(union ieee754dp r)4090efba36SRalf Baechle union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
411da177e4SLinus Torvalds {
42d5afa7e9SMaciej W. Rozycki 	assert(ieee754dp_issnan(r));
431da177e4SLinus Torvalds 
44d5afa7e9SMaciej W. Rozycki 	ieee754_setcx(IEEE754_INVALID_OPERATION);
45acd9e20cSMaciej W. Rozycki 	if (ieee754_csr.nan2008) {
4690d53a91SMaciej W. Rozycki 		DPMANT(r) |= DP_MBIT(DP_FBITS - 1);
47acd9e20cSMaciej W. Rozycki 	} else {
48acd9e20cSMaciej W. Rozycki 		DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1);
49acd9e20cSMaciej W. Rozycki 		if (!ieee754dp_isnan(r))
50acd9e20cSMaciej W. Rozycki 			DPMANT(r) |= DP_MBIT(DP_FBITS - 2);
51acd9e20cSMaciej W. Rozycki 	}
5290d53a91SMaciej W. Rozycki 
5390d53a91SMaciej W. Rozycki 	return r;
541da177e4SLinus Torvalds }
551da177e4SLinus Torvalds 
ieee754dp_get_rounding(int sn,u64 xm)56de2fc342SRalf Baechle static u64 ieee754dp_get_rounding(int sn, u64 xm)
571da177e4SLinus Torvalds {
581da177e4SLinus Torvalds 	/* inexact must round of 3 bits
591da177e4SLinus Torvalds 	 */
601da177e4SLinus Torvalds 	if (xm & (DP_MBIT(3) - 1)) {
611da177e4SLinus Torvalds 		switch (ieee754_csr.rm) {
6256a64733SRalf Baechle 		case FPU_CSR_RZ:
631da177e4SLinus Torvalds 			break;
6456a64733SRalf Baechle 		case FPU_CSR_RN:
651da177e4SLinus Torvalds 			xm += 0x3 + ((xm >> 3) & 1);
661da177e4SLinus Torvalds 			/* xm += (xm&0x8)?0x4:0x3 */
671da177e4SLinus Torvalds 			break;
6856a64733SRalf Baechle 		case FPU_CSR_RU:	/* toward +Infinity */
691da177e4SLinus Torvalds 			if (!sn)	/* ?? */
701da177e4SLinus Torvalds 				xm += 0x8;
711da177e4SLinus Torvalds 			break;
7256a64733SRalf Baechle 		case FPU_CSR_RD:	/* toward -Infinity */
731da177e4SLinus Torvalds 			if (sn) /* ?? */
741da177e4SLinus Torvalds 				xm += 0x8;
751da177e4SLinus Torvalds 			break;
761da177e4SLinus Torvalds 		}
771da177e4SLinus Torvalds 	}
781da177e4SLinus Torvalds 	return xm;
791da177e4SLinus Torvalds }
801da177e4SLinus Torvalds 
811da177e4SLinus Torvalds 
821da177e4SLinus Torvalds /* generate a normal/denormal number with over,under handling
831da177e4SLinus Torvalds  * sn is sign
841da177e4SLinus Torvalds  * xe is an unbiased exponent
851da177e4SLinus Torvalds  * xm is 3bit extended precision value.
861da177e4SLinus Torvalds  */
ieee754dp_format(int sn,int xe,u64 xm)872209bcb1SRalf Baechle union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
881da177e4SLinus Torvalds {
891da177e4SLinus Torvalds 	assert(xm);		/* we don't gen exact zeros (probably should) */
901da177e4SLinus Torvalds 
9192a76f6dSAdam Buchbinder 	assert((xm >> (DP_FBITS + 1 + 3)) == 0);	/* no excess */
921da177e4SLinus Torvalds 	assert(xm & (DP_HIDDEN_BIT << 3));
931da177e4SLinus Torvalds 
941da177e4SLinus Torvalds 	if (xe < DP_EMIN) {
951da177e4SLinus Torvalds 		/* strip lower bits */
961da177e4SLinus Torvalds 		int es = DP_EMIN - xe;
971da177e4SLinus Torvalds 
981da177e4SLinus Torvalds 		if (ieee754_csr.nod) {
999e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_UNDERFLOW);
1009e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_INEXACT);
1011da177e4SLinus Torvalds 
1021da177e4SLinus Torvalds 			switch(ieee754_csr.rm) {
10356a64733SRalf Baechle 			case FPU_CSR_RN:
10456a64733SRalf Baechle 			case FPU_CSR_RZ:
1051da177e4SLinus Torvalds 				return ieee754dp_zero(sn);
10656a64733SRalf Baechle 			case FPU_CSR_RU:    /* toward +Infinity */
1071da177e4SLinus Torvalds 				if (sn == 0)
1081da177e4SLinus Torvalds 					return ieee754dp_min(0);
1091da177e4SLinus Torvalds 				else
1101da177e4SLinus Torvalds 					return ieee754dp_zero(1);
11156a64733SRalf Baechle 			case FPU_CSR_RD:    /* toward -Infinity */
1121da177e4SLinus Torvalds 				if (sn == 0)
1131da177e4SLinus Torvalds 					return ieee754dp_zero(0);
1141da177e4SLinus Torvalds 				else
1151da177e4SLinus Torvalds 					return ieee754dp_min(1);
1161da177e4SLinus Torvalds 			}
1171da177e4SLinus Torvalds 		}
1181da177e4SLinus Torvalds 
119de2fc342SRalf Baechle 		if (xe == DP_EMIN - 1 &&
120de2fc342SRalf Baechle 		    ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
1211da177e4SLinus Torvalds 		{
1221da177e4SLinus Torvalds 			/* Not tiny after rounding */
1239e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_INEXACT);
124de2fc342SRalf Baechle 			xm = ieee754dp_get_rounding(sn, xm);
1251da177e4SLinus Torvalds 			xm >>= 1;
1261da177e4SLinus Torvalds 			/* Clear grs bits */
1271da177e4SLinus Torvalds 			xm &= ~(DP_MBIT(3) - 1);
1281da177e4SLinus Torvalds 			xe++;
1291da177e4SLinus Torvalds 		}
1301da177e4SLinus Torvalds 		else {
1311da177e4SLinus Torvalds 			/* sticky right shift es bits
1321da177e4SLinus Torvalds 			 */
1331da177e4SLinus Torvalds 			xm = XDPSRS(xm, es);
1341da177e4SLinus Torvalds 			xe += es;
1351da177e4SLinus Torvalds 			assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
1361da177e4SLinus Torvalds 			assert(xe == DP_EMIN);
1371da177e4SLinus Torvalds 		}
1381da177e4SLinus Torvalds 	}
1391da177e4SLinus Torvalds 	if (xm & (DP_MBIT(3) - 1)) {
1409e8bad1fSRalf Baechle 		ieee754_setcx(IEEE754_INEXACT);
1411da177e4SLinus Torvalds 		if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
1429e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_UNDERFLOW);
1431da177e4SLinus Torvalds 		}
1441da177e4SLinus Torvalds 
1451da177e4SLinus Torvalds 		/* inexact must round of 3 bits
1461da177e4SLinus Torvalds 		 */
147de2fc342SRalf Baechle 		xm = ieee754dp_get_rounding(sn, xm);
1481da177e4SLinus Torvalds 		/* adjust exponent for rounding add overflowing
1491da177e4SLinus Torvalds 		 */
150ad8fb553SRalf Baechle 		if (xm >> (DP_FBITS + 3 + 1)) {
1511da177e4SLinus Torvalds 			/* add causes mantissa overflow */
1521da177e4SLinus Torvalds 			xm >>= 1;
1531da177e4SLinus Torvalds 			xe++;
1541da177e4SLinus Torvalds 		}
1551da177e4SLinus Torvalds 	}
1561da177e4SLinus Torvalds 	/* strip grs bits */
1571da177e4SLinus Torvalds 	xm >>= 3;
1581da177e4SLinus Torvalds 
15992a76f6dSAdam Buchbinder 	assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
1601da177e4SLinus Torvalds 	assert(xe >= DP_EMIN);
1611da177e4SLinus Torvalds 
1621da177e4SLinus Torvalds 	if (xe > DP_EMAX) {
1639e8bad1fSRalf Baechle 		ieee754_setcx(IEEE754_OVERFLOW);
1649e8bad1fSRalf Baechle 		ieee754_setcx(IEEE754_INEXACT);
1651da177e4SLinus Torvalds 		/* -O can be table indexed by (rm,sn) */
1661da177e4SLinus Torvalds 		switch (ieee754_csr.rm) {
16756a64733SRalf Baechle 		case FPU_CSR_RN:
1681da177e4SLinus Torvalds 			return ieee754dp_inf(sn);
16956a64733SRalf Baechle 		case FPU_CSR_RZ:
1701da177e4SLinus Torvalds 			return ieee754dp_max(sn);
17156a64733SRalf Baechle 		case FPU_CSR_RU:	/* toward +Infinity */
1721da177e4SLinus Torvalds 			if (sn == 0)
1731da177e4SLinus Torvalds 				return ieee754dp_inf(0);
1741da177e4SLinus Torvalds 			else
1751da177e4SLinus Torvalds 				return ieee754dp_max(1);
17656a64733SRalf Baechle 		case FPU_CSR_RD:	/* toward -Infinity */
1771da177e4SLinus Torvalds 			if (sn == 0)
1781da177e4SLinus Torvalds 				return ieee754dp_max(0);
1791da177e4SLinus Torvalds 			else
1801da177e4SLinus Torvalds 				return ieee754dp_inf(1);
1811da177e4SLinus Torvalds 		}
1821da177e4SLinus Torvalds 	}
1831da177e4SLinus Torvalds 	/* gen norm/denorm/zero */
1841da177e4SLinus Torvalds 
1851da177e4SLinus Torvalds 	if ((xm & DP_HIDDEN_BIT) == 0) {
1861da177e4SLinus Torvalds 		/* we underflow (tiny/zero) */
1871da177e4SLinus Torvalds 		assert(xe == DP_EMIN);
1881da177e4SLinus Torvalds 		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
1899e8bad1fSRalf Baechle 			ieee754_setcx(IEEE754_UNDERFLOW);
1901da177e4SLinus Torvalds 		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
1911da177e4SLinus Torvalds 	} else {
19292a76f6dSAdam Buchbinder 		assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
1931da177e4SLinus Torvalds 		assert(xm & DP_HIDDEN_BIT);
1941da177e4SLinus Torvalds 
1951da177e4SLinus Torvalds 		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
1961da177e4SLinus Torvalds 	}
1971da177e4SLinus Torvalds }
198