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

/* IEEE754 floating point arithmetic
 * single precision
 */
/*
 * 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 "ieee754sp.h"

ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y)
{
	COMPXSP;
	COMPYSP;

	EXPLODEXSP;
	EXPLODEYSP;

	CLEARCX;

	FLUSHXSP;
	FLUSHYSP;

	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 ieee754sp_nanxcpt(ieee754sp_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 ieee754sp_xcpt(ieee754sp_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 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;

	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);

	{
		int re = xe + ye;
		int rs = xs ^ ys;
		unsigned rm;

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

		/* multiply 32bits xm,ym to give high 32bits rm with stickness
		 */
		{
			unsigned short lxm = xm & 0xffff;
			unsigned short hxm = xm >> 16;
			unsigned short lym = ym & 0xffff;
			unsigned short hym = ym >> 16;
			unsigned lrm;
			unsigned hrm;

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

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

			{
				unsigned t = hxm * lym; /* 16 * 16 => 32 */
				{
					unsigned 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_MBITS + 1 + 3))) |
			    ((rm << (SP_MBITS + 1 + 3)) != 0);
			re++;
		} else {
			rm = (rm >> (32 - (SP_MBITS + 1 + 3 + 1))) |
			    ((rm << (SP_MBITS + 1 + 3 + 1)) != 0);
		}
		assert(rm & (SP_HIDDEN_BIT << 3));

		SPNORMRET2(rs, re, rm, "mul", x, y);
	}
}
Commit	Line	Data
1da177e4 LT	1	/* IEEE754 floating point arithmetic
	2	* single precision
	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 "ieee754sp.h"
	28
	29	ieee754sp ieee754sp_mul(ieee754sp x, ieee754sp y)
	30	{
	31	COMPXSP;
	32	COMPYSP;
	33
	34	EXPLODEXSP;
	35	EXPLODEYSP;
	36
	37	CLEARCX;
	38
	39	FLUSHXSP;
	40	FLUSHYSP;
	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 ieee754sp_nanxcpt(ieee754sp_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 ieee754sp_xcpt(ieee754sp_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 ieee754sp_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 ieee754sp_zero(xs ^ ys);
91
92
93	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_DNORM):
94	SPDNORMX;
95
96	case CLPAIR(IEEE754_CLASS_NORM, IEEE754_CLASS_DNORM):
97	SPDNORMY;
98	break;
99
100	case CLPAIR(IEEE754_CLASS_DNORM, IEEE754_CLASS_NORM):
101	SPDNORMX;
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 & SP_HIDDEN_BIT);
	109	assert(ym & SP_HIDDEN_BIT);
	110
	111	{
	112	int re = xe + ye;
	113	int rs = xs ^ ys;
	114	unsigned rm;
	115
	116	/* shunt to top of word */
	117	xm <<= 32 - (SP_MBITS + 1);
	118	ym <<= 32 - (SP_MBITS + 1);
	119
	120	/* multiply 32bits xm,ym to give high 32bits rm with stickness
	121	*/
	122	{
	123	unsigned short lxm = xm & 0xffff;
	124	unsigned short hxm = xm >> 16;
	125	unsigned short lym = ym & 0xffff;
	126	unsigned short hym = ym >> 16;
	127	unsigned lrm;
	128	unsigned hrm;
	129
	130	lrm = lxm * lym; /* 16 * 16 => 32 */
	131	hrm = hxm * hym; /* 16 * 16 => 32 */
	132
	133	{
70342287	134	unsigned t = lxm * hym; /* 16 * 16 => 32 */
1da177e4 LT	135	{
	136	unsigned at = lrm + (t << 16);
	137	hrm += at < lrm;
	138	lrm = at;
	139	}
	140	hrm = hrm + (t >> 16);
	141	}
	142
	143	{
70342287	144	unsigned t = hxm * lym; /* 16 * 16 => 32 */
1da177e4 LT	145	{
	146	unsigned at = lrm + (t << 16);
	147	hrm += at < lrm;
	148	lrm = at;
	149	}
	150	hrm = hrm + (t >> 16);
	151	}
	152	rm = hrm \| (lrm != 0);
	153	}
	154
	155	/*
	156	* sticky shift down to normal rounding precision
	157	*/
	158	if ((int) rm < 0) {
	159	rm = (rm >> (32 - (SP_MBITS + 1 + 3))) \|
	160	((rm << (SP_MBITS + 1 + 3)) != 0);
	161	re++;
	162	} else {
	163	rm = (rm >> (32 - (SP_MBITS + 1 + 3 + 1))) \|
	164	((rm << (SP_MBITS + 1 + 3 + 1)) != 0);
	165	}
	166	assert(rm & (SP_HIDDEN_BIT << 3));
	167
	168	SPNORMRET2(rs, re, rm, "mul", x, y);
	169	}
	170	}