[linux-2.6-block.git] / arch / mips / math-emu / ieee754dp.c

// SPDX-License-Identifier: GPL-2.0-only
/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.
 */

#include <linux/compiler.h>

#include "ieee754dp.h"

int ieee754dp_class(union ieee754dp x)
{
	COMPXDP;
	EXPLODEXDP;
	return xc;
}

static inline int ieee754dp_isnan(union ieee754dp x)
{
	return ieee754_class_nan(ieee754dp_class(x));
}

static inline int ieee754dp_issnan(union ieee754dp x)
{
	int qbit;

	assert(ieee754dp_isnan(x));
	qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1);
	return ieee754_csr.nan2008 ^ qbit;
}


/*
 * Raise the Invalid Operation IEEE 754 exception
 * and convert the signaling NaN supplied to a quiet NaN.
 */
union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
{
	assert(ieee754dp_issnan(r));

	ieee754_setcx(IEEE754_INVALID_OPERATION);
	if (ieee754_csr.nan2008) {
		DPMANT(r) |= DP_MBIT(DP_FBITS - 1);
	} else {
		DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1);
		if (!ieee754dp_isnan(r))
			DPMANT(r) |= DP_MBIT(DP_FBITS - 2);
	}

	return r;
}

static u64 ieee754dp_get_rounding(int sn, u64 xm)
{
	/* inexact must round of 3 bits
	 */
	if (xm & (DP_MBIT(3) - 1)) {
		switch (ieee754_csr.rm) {
		case FPU_CSR_RZ:
			break;
		case FPU_CSR_RN:
			xm += 0x3 + ((xm >> 3) & 1);
			/* xm += (xm&0x8)?0x4:0x3 */
			break;
		case FPU_CSR_RU:	/* toward +Infinity */
			if (!sn)	/* ?? */
				xm += 0x8;
			break;
		case FPU_CSR_RD:	/* toward -Infinity */
			if (sn) /* ?? */
				xm += 0x8;
			break;
		}
	}
	return xm;
}


/* generate a normal/denormal number with over,under handling
 * sn is sign
 * xe is an unbiased exponent
 * xm is 3bit extended precision value.
 */
union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
{
	assert(xm);		/* we don't gen exact zeros (probably should) */

	assert((xm >> (DP_FBITS + 1 + 3)) == 0);	/* no excess */
	assert(xm & (DP_HIDDEN_BIT << 3));

	if (xe < DP_EMIN) {
		/* strip lower bits */
		int es = DP_EMIN - xe;

		if (ieee754_csr.nod) {
			ieee754_setcx(IEEE754_UNDERFLOW);
			ieee754_setcx(IEEE754_INEXACT);

			switch(ieee754_csr.rm) {
			case FPU_CSR_RN:
			case FPU_CSR_RZ:
				return ieee754dp_zero(sn);
			case FPU_CSR_RU:    /* toward +Infinity */
				if (sn == 0)
					return ieee754dp_min(0);
				else
					return ieee754dp_zero(1);
			case FPU_CSR_RD:    /* toward -Infinity */
				if (sn == 0)
					return ieee754dp_zero(0);
				else
					return ieee754dp_min(1);
			}
		}

		if (xe == DP_EMIN - 1 &&
		    ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
		{
			/* Not tiny after rounding */
			ieee754_setcx(IEEE754_INEXACT);
			xm = ieee754dp_get_rounding(sn, xm);
			xm >>= 1;
			/* Clear grs bits */
			xm &= ~(DP_MBIT(3) - 1);
			xe++;
		}
		else {
			/* sticky right shift es bits
			 */
			xm = XDPSRS(xm, es);
			xe += es;
			assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
			assert(xe == DP_EMIN);
		}
	}
	if (xm & (DP_MBIT(3) - 1)) {
		ieee754_setcx(IEEE754_INEXACT);
		if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
			ieee754_setcx(IEEE754_UNDERFLOW);
		}

		/* inexact must round of 3 bits
		 */
		xm = ieee754dp_get_rounding(sn, xm);
		/* adjust exponent for rounding add overflowing
		 */
		if (xm >> (DP_FBITS + 3 + 1)) {
			/* add causes mantissa overflow */
			xm >>= 1;
			xe++;
		}
	}
	/* strip grs bits */
	xm >>= 3;

	assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
	assert(xe >= DP_EMIN);

	if (xe > DP_EMAX) {
		ieee754_setcx(IEEE754_OVERFLOW);
		ieee754_setcx(IEEE754_INEXACT);
		/* -O can be table indexed by (rm,sn) */
		switch (ieee754_csr.rm) {
		case FPU_CSR_RN:
			return ieee754dp_inf(sn);
		case FPU_CSR_RZ:
			return ieee754dp_max(sn);
		case FPU_CSR_RU:	/* toward +Infinity */
			if (sn == 0)
				return ieee754dp_inf(0);
			else
				return ieee754dp_max(1);
		case FPU_CSR_RD:	/* toward -Infinity */
			if (sn == 0)
				return ieee754dp_max(0);
			else
				return ieee754dp_inf(1);
		}
	}
	/* gen norm/denorm/zero */

	if ((xm & DP_HIDDEN_BIT) == 0) {
		/* we underflow (tiny/zero) */
		assert(xe == DP_EMIN);
		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
			ieee754_setcx(IEEE754_UNDERFLOW);
		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
	} else {
		assert((xm >> (DP_FBITS + 1)) == 0);	/* no excess */
		assert(xm & DP_HIDDEN_BIT);

		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
	}
}
Commit	Line	Data
9d5a6349	1	// SPDX-License-Identifier: GPL-2.0-only
1da177e4 LT	2	/* IEEE754 floating point arithmetic
	3	* double precision: common utilities
	4	*/
	5	/*
	6	* MIPS floating point support
	7	* Copyright (C) 1994-2000 Algorithmics Ltd.
1da177e4 LT	8	*/
1da177e4 LT	9
cae55066	10	#include <linux/compiler.h>
1da177e4 LT	11
	12	#include "ieee754dp.h"
	13
2209bcb1	14	int ieee754dp_class(union ieee754dp x)
1da177e4 LT	15	{
	16	COMPXDP;
	17	EXPLODEXDP;
	18	return xc;
	19	}
	20
e06b530b	21	static inline int ieee754dp_isnan(union ieee754dp x)
1da177e4	22	{
c9a10845	23	return ieee754_class_nan(ieee754dp_class(x));
1da177e4 LT	24	}
1da177e4 LT	25
f71baa11	26	static inline int ieee754dp_issnan(union ieee754dp x)
1da177e4	27	{
90d53a91 MR	28	int qbit;
90d53a91 MR	29
1da177e4	30	assert(ieee754dp_isnan(x));
90d53a91 MR	31	qbit = (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1);
90d53a91 MR	32	return ieee754_csr.nan2008 ^ qbit;
1da177e4 LT	33	}
	34
	35
d5afa7e9 MR	36	/*
	37	* Raise the Invalid Operation IEEE 754 exception
	38	* and convert the signaling NaN supplied to a quiet NaN.
	39	*/
90efba36	40	union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
1da177e4	41	{
d5afa7e9	42	assert(ieee754dp_issnan(r));
1da177e4	43
d5afa7e9	44	ieee754_setcx(IEEE754_INVALID_OPERATION);
acd9e20c	45	if (ieee754_csr.nan2008) {
90d53a91	46	DPMANT(r) \|= DP_MBIT(DP_FBITS - 1);
acd9e20c MR	47	} else {
	48	DPMANT(r) &= ~DP_MBIT(DP_FBITS - 1);
	49	if (!ieee754dp_isnan(r))
	50	DPMANT(r) \|= DP_MBIT(DP_FBITS - 2);
	51	}
90d53a91 MR	52
90d53a91 MR	53	return r;
1da177e4 LT	54	}
1da177e4 LT	55
de2fc342	56	static u64 ieee754dp_get_rounding(int sn, u64 xm)
1da177e4 LT	57	{
	58	/* inexact must round of 3 bits
	59	*/
	60	if (xm & (DP_MBIT(3) - 1)) {
	61	switch (ieee754_csr.rm) {
56a64733	62	case FPU_CSR_RZ:
1da177e4	63	break;
56a64733	64	case FPU_CSR_RN:
1da177e4 LT	65	xm += 0x3 + ((xm >> 3) & 1);
	66	/* xm += (xm&0x8)?0x4:0x3 */
	67	break;
56a64733	68	case FPU_CSR_RU: /* toward +Infinity */
1da177e4 LT	69	if (!sn) /* ?? */
	70	xm += 0x8;
	71	break;
56a64733	72	case FPU_CSR_RD: /* toward -Infinity */
70342287	73	if (sn) /* ?? */
1da177e4 LT	74	xm += 0x8;
	75	break;
	76	}
	77	}
	78	return xm;
	79	}
	80
	81
	82	/* generate a normal/denormal number with over,under handling
	83	* sn is sign
	84	* xe is an unbiased exponent
	85	* xm is 3bit extended precision value.
	86	*/
2209bcb1	87	union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
1da177e4 LT	88	{
	89	assert(xm); /* we don't gen exact zeros (probably should) */
	90
92a76f6d	91	assert((xm >> (DP_FBITS + 1 + 3)) == 0); /* no excess */
1da177e4 LT	92	assert(xm & (DP_HIDDEN_BIT << 3));
	93
	94	if (xe < DP_EMIN) {
	95	/* strip lower bits */
	96	int es = DP_EMIN - xe;
	97
	98	if (ieee754_csr.nod) {
9e8bad1f RB	99	ieee754_setcx(IEEE754_UNDERFLOW);
9e8bad1f RB	100	ieee754_setcx(IEEE754_INEXACT);
1da177e4 LT	101
1da177e4 LT	102	switch(ieee754_csr.rm) {
56a64733 RB	103	case FPU_CSR_RN:
56a64733 RB	104	case FPU_CSR_RZ:
1da177e4	105	return ieee754dp_zero(sn);
56a64733	106	case FPU_CSR_RU: /* toward +Infinity */
47fa0c02	107	if (sn == 0)
1da177e4 LT	108	return ieee754dp_min(0);
	109	else
	110	return ieee754dp_zero(1);
56a64733	111	case FPU_CSR_RD: /* toward -Infinity */
47fa0c02	112	if (sn == 0)
1da177e4 LT	113	return ieee754dp_zero(0);
	114	else
	115	return ieee754dp_min(1);
	116	}
	117	}
	118
de2fc342 RB	119	if (xe == DP_EMIN - 1 &&
de2fc342 RB	120	ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
1da177e4 LT	121	{
1da177e4 LT	122	/* Not tiny after rounding */
9e8bad1f	123	ieee754_setcx(IEEE754_INEXACT);
de2fc342	124	xm = ieee754dp_get_rounding(sn, xm);
1da177e4 LT	125	xm >>= 1;
	126	/* Clear grs bits */
	127	xm &= ~(DP_MBIT(3) - 1);
	128	xe++;
	129	}
	130	else {
	131	/* sticky right shift es bits
	132	*/
	133	xm = XDPSRS(xm, es);
	134	xe += es;
	135	assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
	136	assert(xe == DP_EMIN);
	137	}
	138	}
	139	if (xm & (DP_MBIT(3) - 1)) {
9e8bad1f	140	ieee754_setcx(IEEE754_INEXACT);
1da177e4	141	if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
9e8bad1f	142	ieee754_setcx(IEEE754_UNDERFLOW);
1da177e4 LT	143	}
	144
	145	/* inexact must round of 3 bits
	146	*/
de2fc342	147	xm = ieee754dp_get_rounding(sn, xm);
1da177e4 LT	148	/* adjust exponent for rounding add overflowing
1da177e4 LT	149	*/
ad8fb553	150	if (xm >> (DP_FBITS + 3 + 1)) {
1da177e4 LT	151	/* add causes mantissa overflow */
	152	xm >>= 1;
	153	xe++;
	154	}
	155	}
	156	/* strip grs bits */
	157	xm >>= 3;
	158
92a76f6d	159	assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */
1da177e4 LT	160	assert(xe >= DP_EMIN);
	161
	162	if (xe > DP_EMAX) {
9e8bad1f RB	163	ieee754_setcx(IEEE754_OVERFLOW);
9e8bad1f RB	164	ieee754_setcx(IEEE754_INEXACT);
1da177e4 LT	165	/* -O can be table indexed by (rm,sn) */
1da177e4 LT	166	switch (ieee754_csr.rm) {
56a64733	167	case FPU_CSR_RN:
1da177e4	168	return ieee754dp_inf(sn);
56a64733	169	case FPU_CSR_RZ:
1da177e4	170	return ieee754dp_max(sn);
56a64733	171	case FPU_CSR_RU: /* toward +Infinity */
1da177e4 LT	172	if (sn == 0)
	173	return ieee754dp_inf(0);
	174	else
	175	return ieee754dp_max(1);
56a64733	176	case FPU_CSR_RD: /* toward -Infinity */
1da177e4 LT	177	if (sn == 0)
	178	return ieee754dp_max(0);
	179	else
	180	return ieee754dp_inf(1);
	181	}
	182	}
	183	/* gen norm/denorm/zero */
	184
	185	if ((xm & DP_HIDDEN_BIT) == 0) {
	186	/* we underflow (tiny/zero) */
	187	assert(xe == DP_EMIN);
	188	if (ieee754_csr.mx & IEEE754_UNDERFLOW)
9e8bad1f	189	ieee754_setcx(IEEE754_UNDERFLOW);
1da177e4 LT	190	return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
1da177e4 LT	191	} else {
92a76f6d	192	assert((xm >> (DP_FBITS + 1)) == 0); /* no excess */
1da177e4 LT	193	assert(xm & DP_HIDDEN_BIT);
	194
	195	return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
	196	}
	197	}