[mirror_ubuntu-bionic-kernel.git] / arch / mips / math-emu / dp_mul.c

/* IEEE754 floating point arithmetic
 * double precision: common utilities
 */
/*
 * MIPS floating point support
 * Copyright (C) 1994-2000 Algorithmics Ltd.
 *
 * ########################################################################
 *
 *  This program is free software; you can distribute it and/or modify it
 *  under the terms of the GNU General Public License (Version 2) as
 *  published by the Free Software Foundation.
 *
 *  This program is distributed in the hope it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 *  for more details.
 *
 *  You should have received a copy of the GNU General Public License along
 *  with this program; if not, write to the Free Software Foundation, Inc.,
 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
 *
 * ########################################################################
 */


#include "ieee754dp.h"

ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y)
{
	COMPXDP;
	COMPYDP;

	EXPLODEXDP;
	EXPLODEYDP;

	CLEARCX;

	FLUSHXDP;
	FLUSHYDP;

	switch (CLPAIR(xc, yc)) {
	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	case CLPAIR(IEEE754_CLASS_SNAN, 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):
	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y);

	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):
		SETCX(IEEE754_INVALID_OPERATION);
		return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y);

	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 ieee754dp_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 ieee754dp_zero(xs ^ ys);


	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
		DPDNORMX;

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

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

	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
		break;
	}
	/* rm = xm * ym, re = xe+ye basically */
	assert(xm & DP_HIDDEN_BIT);
	assert(ym & DP_HIDDEN_BIT);
	{
		int re = xe + ye;
		int rs = xs ^ ys;
		u64 rm;

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

		/* multiply 32bits xm,ym to give high 32bits rm with stickness
		 */

		/* 32 * 32 => 64 */
#define DPXMULT(x, y)	((u64)(x) * (u64)y)

		{
			unsigned lxm = xm;
			unsigned hxm = xm >> 32;
			unsigned lym = ym;
			unsigned hym = ym >> 32;
			u64 lrm;
			u64 hrm;

			lrm = DPXMULT(lxm, lym);
			hrm = DPXMULT(hxm, hym);

			{
				u64 t = DPXMULT(lxm, hym);
				{
					u64 at =
					    lrm + (t << 32);
					hrm += at < lrm;
					lrm = at;
				}
				hrm = hrm + (t >> 32);
			}

			{
				u64 t = DPXMULT(hxm, lym);
				{
					u64 at =
					    lrm + (t << 32);
					hrm += at < lrm;
					lrm = at;
				}
				hrm = hrm + (t >> 32);
			}
			rm = hrm | (lrm != 0);
		}

		/*
		 * sticky shift down to normal rounding precision
		 */
		if ((s64) rm < 0) {
			rm =
			    (rm >> (64 - (DP_MBITS + 1 + 3))) |
			    ((rm << (DP_MBITS + 1 + 3)) != 0);
			re++;
		} else {
			rm =
			    (rm >> (64 - (DP_MBITS + 1 + 3 + 1))) |
			    ((rm << (DP_MBITS + 1 + 3 + 1)) != 0);
		}
		assert(rm & (DP_HIDDEN_BIT << 3));
		DPNORMRET2(rs, re, rm, "mul", x, y);
	}
}
Commit	Line	Data
1da177e4 LT	1	/* IEEE754 floating point arithmetic
	2	* double precision: common utilities
	3	*/
	4	/*
	5	* MIPS floating point support
	6	* Copyright (C) 1994-2000 Algorithmics Ltd.
1da177e4 LT	7	*
	8	* ########################################################################
	9	*
	10	* This program is free software; you can distribute it and/or modify it
	11	* under the terms of the GNU General Public License (Version 2) as
	12	* published by the Free Software Foundation.
	13	*
	14	* This program is distributed in the hope it will be useful, but WITHOUT
	15	* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
	16	* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
	17	* for more details.
	18	*
	19	* You should have received a copy of the GNU General Public License along
	20	* with this program; if not, write to the Free Software Foundation, Inc.,
	21	* 59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
	22	*
	23	* ########################################################################
	24	*/
	25
	26
	27	#include "ieee754dp.h"
	28
	29	ieee754dp ieee754dp_mul(ieee754dp x, ieee754dp y)
	30	{
	31	COMPXDP;
	32	COMPYDP;
	33
	34	EXPLODEXDP;
	35	EXPLODEYDP;
	36
	37	CLEARCX;
	38
	39	FLUSHXDP;
	40	FLUSHYDP;
	41
	42	switch (CLPAIR(xc, yc)) {
	43	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_QNAN):
	44	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_SNAN):
	45	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_SNAN):
	46	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_SNAN):
	47	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_SNAN):
	48	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_SNAN):
	49	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_SNAN):
	50	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_ZERO):
	51	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_NORM):
	52	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_DNORM):
	53	case CLPAIR(IEEE754_CLASS_SNAN, IEEE754_CLASS_INF):
	54	SETCX(IEEE754_INVALID_OPERATION);
	55	return ieee754dp_nanxcpt(ieee754dp_indef(), "mul", x, y);
	56
	57	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_QNAN):
	58	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_QNAN):
	59	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_QNAN):
	60	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_QNAN):
	61	return y;
	62
	63	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_QNAN):
	64	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_ZERO):
	65	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_NORM):
	66	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_DNORM):
	67	case CLPAIR(IEEE754_CLASS_QNAN, IEEE754_CLASS_INF):
	68	return x;
	69
	70
71	/* Infinity handling */
72
73	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_ZERO):
74	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_INF):
75	SETCX(IEEE754_INVALID_OPERATION);
76	return ieee754dp_xcpt(ieee754dp_indef(), "mul", x, y);
77
78	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_INF):
79	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_INF):
80	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_NORM):
81	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_DNORM):
82	case CLPAIR(IEEE754_CLASS_INF, IEEE754_CLASS_INF):
83	return ieee754dp_inf(xs ^ ys);
84
85	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_ZERO):
86	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_NORM):
87	case CLPAIR(IEEE754_CLASS_ZERO, IEEE754_CLASS_DNORM):
88	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_ZERO):
89	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_ZERO):
90	return ieee754dp_zero(xs ^ ys);
91
92
93	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
94	DPDNORMX;
95
96	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
97	DPDNORMY;
98	break;
99
100	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
101	DPDNORMX;
102	break;
103
104	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_NORM):
105	break;
106	}
25985edc	107	/* rm = xm * ym, re = xe+ye basically */
1da177e4 LT	108	assert(xm & DP_HIDDEN_BIT);
	109	assert(ym & DP_HIDDEN_BIT);
	110	{
	111	int re = xe + ye;
	112	int rs = xs ^ ys;
	113	u64 rm;
	114
	115	/* shunt to top of word */
	116	xm <<= 64 - (DP_MBITS + 1);
	117	ym <<= 64 - (DP_MBITS + 1);
	118
	119	/* multiply 32bits xm,ym to give high 32bits rm with stickness
	120	*/
	121
	122	/* 32 * 32 => 64 */
21a151d8	123	#define DPXMULT(x, y) ((u64)(x) * (u64)y)
1da177e4 LT	124
	125	{
	126	unsigned lxm = xm;
	127	unsigned hxm = xm >> 32;
	128	unsigned lym = ym;
	129	unsigned hym = ym >> 32;
	130	u64 lrm;
	131	u64 hrm;
	132
	133	lrm = DPXMULT(lxm, lym);
	134	hrm = DPXMULT(hxm, hym);
	135
	136	{
	137	u64 t = DPXMULT(lxm, hym);
	138	{
	139	u64 at =
	140	lrm + (t << 32);
	141	hrm += at < lrm;
	142	lrm = at;
	143	}
	144	hrm = hrm + (t >> 32);
	145	}
	146
	147	{
	148	u64 t = DPXMULT(hxm, lym);
	149	{
	150	u64 at =
	151	lrm + (t << 32);
	152	hrm += at < lrm;
	153	lrm = at;
	154	}
	155	hrm = hrm + (t >> 32);
	156	}
	157	rm = hrm \| (lrm != 0);
	158	}
	159
	160	/*
	161	* sticky shift down to normal rounding precision
	162	*/
	163	if ((s64) rm < 0) {
	164	rm =
	165	(rm >> (64 - (DP_MBITS + 1 + 3))) \|
	166	((rm << (DP_MBITS + 1 + 3)) != 0);
	167	re++;
	168	} else {
	169	rm =
	170	(rm >> (64 - (DP_MBITS + 1 + 3 + 1))) \|
	171	((rm << (DP_MBITS + 1 + 3 + 1)) != 0);
	172	}
	173	assert(rm & (DP_HIDDEN_BIT << 3));
	174	DPNORMRET2(rs, re, rm, "mul", x, y);
	175	}
	176	}