1/* IEEE754 floating point arithmetic
2 * single precision
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
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
29ieee754sp 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	}
107	/* rm = xm * ym, re = xe+ye basically */
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			{
134				unsigned t = lxm * hym;	/* 16 * 16 => 32 */
135				{
136					unsigned at = lrm + (t << 16);
137					hrm += at < lrm;
138					lrm = at;
139				}
140				hrm = hrm + (t >> 16);
141			}
142
143			{
144				unsigned t = hxm * lym;	/* 16 * 16 => 32 */
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}
171