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

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

#include "ieee754sp.h"

union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y)
{
	int re;
	int rs;
	unsigned int rm;
	unsigned short lxm;
	unsigned short hxm;
	unsigned short lym;
	unsigned short hym;
	unsigned int lrm;
	unsigned int hrm;
	unsigned int t;
	unsigned int at;

	COMPXSP;
	COMPYSP;

	EXPLODEXSP;
	EXPLODEYSP;

	ieee754_clearcx();

	FLUSHXSP;
	FLUSHYSP;

	switch (CLPAIR(xc, yc)) {
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
		return ieee754sp_nanxcpt(y);

	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
		return ieee754sp_nanxcpt(x);

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
		return y;

	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
		return x;


	/*
	 * Infinity handling
	 */
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
		ieee754_setcx(IEEE754_INVALID_OPERATION);
		return ieee754sp_indef();

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
		return ieee754sp_inf(xs ^ ys);

	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
		return ieee754sp_zero(xs ^ ys);


	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
		SPDNORMX;
		fallthrough;
	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
		SPDNORMY;
		break;

	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
		SPDNORMX;
		break;

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
		break;
	}
	/* rm = xm * ym, re = xe+ye basically */
	assert(xm & SP_HIDDEN_BIT);
	assert(ym & SP_HIDDEN_BIT);

	re = xe + ye;
	rs = xs ^ ys;

	/* shunt to top of word */
	xm <<= 32 - (SP_FBITS + 1);
	ym <<= 32 - (SP_FBITS + 1);

	/*
	 * Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
	 */
	lxm = xm & 0xffff;
	hxm = xm >> 16;
	lym = ym & 0xffff;
	hym = ym >> 16;

	lrm = lxm * lym;	/* 16 * 16 => 32 */
	hrm = hxm * hym;	/* 16 * 16 => 32 */

	t = lxm * hym; /* 16 * 16 => 32 */
	at = lrm + (t << 16);
	hrm += at < lrm;
	lrm = at;
	hrm = hrm + (t >> 16);

	t = hxm * lym; /* 16 * 16 => 32 */
	at = lrm + (t << 16);
	hrm += at < lrm;
	lrm = at;
	hrm = hrm + (t >> 16);

	rm = hrm | (lrm != 0);

	/*
	 * Sticky shift down to normal rounding precision.
	 */
	if ((int) rm < 0) {
		rm = (rm >> (32 - (SP_FBITS + 1 + 3))) |
		    ((rm << (SP_FBITS + 1 + 3)) != 0);
		re++;
	} else {
		rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) |
		     ((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
	}
	assert(rm & (SP_HIDDEN_BIT << 3));

	return ieee754sp_format(rs, re, rm);
}
Commit	Line	Data
9d5a6349	1	// SPDX-License-Identifier: GPL-2.0-only
1da177e4 LT	2	/* IEEE754 floating point arithmetic
	3	* single precision
	4	*/
	5	/*
	6	* MIPS floating point support
	7	* Copyright (C) 1994-2000 Algorithmics Ltd.
1da177e4 LT	8	*/
1da177e4 LT	9
1da177e4 LT	10	#include "ieee754sp.h"
1da177e4 LT	11
2209bcb1	12	union ieee754sp ieee754sp_mul(union ieee754sp x, union ieee754sp y)
1da177e4	13	{
3f7cac41 RB	14	int re;
3f7cac41 RB	15	int rs;
a58f85b5	16	unsigned int rm;
3f7cac41 RB	17	unsigned short lxm;
	18	unsigned short hxm;
	19	unsigned short lym;
	20	unsigned short hym;
a58f85b5 AM	21	unsigned int lrm;
	22	unsigned int hrm;
	23	unsigned int t;
	24	unsigned int at;
3f7cac41	25
1da177e4 LT	26	COMPXSP;
	27	COMPYSP;
	28
	29	EXPLODEXSP;
	30	EXPLODEYSP;
	31
9e8bad1f	32	ieee754_clearcx();
1da177e4 LT	33
	34	FLUSHXSP;
	35	FLUSHYSP;
	36
	37	switch (CLPAIR(xc, yc)) {
1da177e4	38	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
1da177e4 LT	39	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	40	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	41	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	42	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
d5afa7e9 MR	43	return ieee754sp_nanxcpt(y);
	44
	45	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	46	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
1da177e4 LT	47	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
	48	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
	49	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
	50	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
d5afa7e9	51	return ieee754sp_nanxcpt(x);
1da177e4 LT	52
	53	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	54	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	55	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	56	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
	57	return y;
	58
	59	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	60	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	61	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	62	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	63	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
	64	return x;
	65
	66
3f7cac41 RB	67	/*
	68	* Infinity handling
	69	*/
1da177e4 LT	70	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
1da177e4 LT	71	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
9e8bad1f	72	ieee754_setcx(IEEE754_INVALID_OPERATION);
90efba36	73	return ieee754sp_indef();
1da177e4 LT	74
	75	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
	76	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
	77	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
	78	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
	79	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
	80	return ieee754sp_inf(xs ^ ys);
	81
	82	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
	83	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
	84	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
	85	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
	86	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
	87	return ieee754sp_zero(xs ^ ys);
	88
	89
	90	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
	91	SPDNORMX;
c9b02990	92	fallthrough;
1da177e4 LT	93	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
	94	SPDNORMY;
	95	break;
	96
	97	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
	98	SPDNORMX;
	99	break;
	100
	101	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
	102	break;
	103	}
25985edc	104	/* rm = xm * ym, re = xe+ye basically */
1da177e4 LT	105	assert(xm & SP_HIDDEN_BIT);
	106	assert(ym & SP_HIDDEN_BIT);
	107
3f7cac41 RB	108	re = xe + ye;
	109	rs = xs ^ ys;
	110
	111	/* shunt to top of word */
	112	xm <<= 32 - (SP_FBITS + 1);
	113	ym <<= 32 - (SP_FBITS + 1);
	114
	115	/*
	116	* Multiply 32 bits xm, ym to give high 32 bits rm with stickness.
	117	*/
	118	lxm = xm & 0xffff;
	119	hxm = xm >> 16;
	120	lym = ym & 0xffff;
	121	hym = ym >> 16;
	122
	123	lrm = lxm * lym; /* 16 * 16 => 32 */
	124	hrm = hxm * hym; /* 16 * 16 => 32 */
	125
	126	t = lxm * hym; /* 16 * 16 => 32 */
	127	at = lrm + (t << 16);
	128	hrm += at < lrm;
	129	lrm = at;
	130	hrm = hrm + (t >> 16);
	131
	132	t = hxm * lym; /* 16 * 16 => 32 */
	133	at = lrm + (t << 16);
	134	hrm += at < lrm;
	135	lrm = at;
	136	hrm = hrm + (t >> 16);
	137
	138	rm = hrm \| (lrm != 0);
	139
	140	/*
	141	* Sticky shift down to normal rounding precision.
	142	*/
	143	if ((int) rm < 0) {
	144	rm = (rm >> (32 - (SP_FBITS + 1 + 3))) \|
	145	((rm << (SP_FBITS + 1 + 3)) != 0);
	146	re++;
	147	} else {
	148	rm = (rm >> (32 - (SP_FBITS + 1 + 3 + 1))) \|
	149	((rm << (SP_FBITS + 1 + 3 + 1)) != 0);
1da177e4	150	}
3f7cac41 RB	151	assert(rm & (SP_HIDDEN_BIT << 3));
	152
	153	return ieee754sp_format(rs, re, rm);
1da177e4	154	}