1/* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4/*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
7 * http://www.algor.co.uk
8 *
9 * ########################################################################
10 *
11 *  This program is free software; you can distribute it and/or modify it
12 *  under the terms of the GNU General Public License (Version 2) as
13 *  published by the Free Software Foundation.
14 *
15 *  This program is distributed in the hope it will be useful, but WITHOUT
16 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
18 *  for more details.
19 *
20 *  You should have received a copy of the GNU General Public License along
21 *  with this program; if not, write to the Free Software Foundation, Inc.,
22 *  59 Temple Place - Suite 330, Boston MA 02111-1307, USA.
23 *
24 * ########################################################################
25 */
26
27
28#include "ieee754dp.h"
29
30int ieee754dp_class(ieee754dp x)
31{
32	COMPXDP;
33	EXPLODEXDP;
34	return xc;
35}
36
37int ieee754dp_isnan(ieee754dp x)
38{
39	return ieee754dp_class(x) >= IEEE754_CLASS_SNAN;
40}
41
42int ieee754dp_issnan(ieee754dp x)
43{
44	assert(ieee754dp_isnan(x));
45	return ((DPMANT(x) & DP_MBIT(DP_MBITS-1)) == DP_MBIT(DP_MBITS-1));
46}
47
48
49ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...)
50{
51	struct ieee754xctx ax;
52	if (!TSTX())
53		return r;
54
55	ax.op = op;
56	ax.rt = IEEE754_RT_DP;
57	ax.rv.dp = r;
58	va_start(ax.ap, op);
59	ieee754_xcpt(&ax);
60	return ax.rv.dp;
61}
62
63ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...)
64{
65	struct ieee754xctx ax;
66
67	assert(ieee754dp_isnan(r));
68
69	if (!ieee754dp_issnan(r))	/* QNAN does not cause invalid op !! */
70		return r;
71
72	if (!SETANDTESTCX(IEEE754_INVALID_OPERATION)) {
73		/* not enabled convert to a quiet NaN */
74		DPMANT(r) &= (~DP_MBIT(DP_MBITS-1));
75		if (ieee754dp_isnan(r))
76			return r;
77		else
78			return ieee754dp_indef();
79	}
80
81	ax.op = op;
82	ax.rt = 0;
83	ax.rv.dp = r;
84	va_start(ax.ap, op);
85	ieee754_xcpt(&ax);
86	return ax.rv.dp;
87}
88
89ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y)
90{
91	assert(ieee754dp_isnan(x));
92	assert(ieee754dp_isnan(y));
93
94	if (DPMANT(x) > DPMANT(y))
95		return x;
96	else
97		return y;
98}
99
100
101static u64 get_rounding(int sn, u64 xm)
102{
103	/* inexact must round of 3 bits
104	 */
105	if (xm & (DP_MBIT(3) - 1)) {
106		switch (ieee754_csr.rm) {
107		case IEEE754_RZ:
108			break;
109		case IEEE754_RN:
110			xm += 0x3 + ((xm >> 3) & 1);
111			/* xm += (xm&0x8)?0x4:0x3 */
112			break;
113		case IEEE754_RU:	/* toward +Infinity */
114			if (!sn)	/* ?? */
115				xm += 0x8;
116			break;
117		case IEEE754_RD:	/* toward -Infinity */
118			if (sn)	/* ?? */
119				xm += 0x8;
120			break;
121		}
122	}
123	return xm;
124}
125
126
127/* generate a normal/denormal number with over,under handling
128 * sn is sign
129 * xe is an unbiased exponent
130 * xm is 3bit extended precision value.
131 */
132ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
133{
134	assert(xm);		/* we don't gen exact zeros (probably should) */
135
136	assert((xm >> (DP_MBITS + 1 + 3)) == 0);	/* no execess */
137	assert(xm & (DP_HIDDEN_BIT << 3));
138
139	if (xe < DP_EMIN) {
140		/* strip lower bits */
141		int es = DP_EMIN - xe;
142
143		if (ieee754_csr.nod) {
144			SETCX(IEEE754_UNDERFLOW);
145			SETCX(IEEE754_INEXACT);
146
147			switch(ieee754_csr.rm) {
148			case IEEE754_RN:
149				return ieee754dp_zero(sn);
150			case IEEE754_RZ:
151				return ieee754dp_zero(sn);
152			case IEEE754_RU:    /* toward +Infinity */
153				if(sn == 0)
154					return ieee754dp_min(0);
155				else
156					return ieee754dp_zero(1);
157			case IEEE754_RD:    /* toward -Infinity */
158				if(sn == 0)
159					return ieee754dp_zero(0);
160				else
161					return ieee754dp_min(1);
162			}
163		}
164
165		if (xe == DP_EMIN - 1
166				&& get_rounding(sn, xm) >> (DP_MBITS + 1 + 3))
167		{
168			/* Not tiny after rounding */
169			SETCX(IEEE754_INEXACT);
170			xm = get_rounding(sn, xm);
171			xm >>= 1;
172			/* Clear grs bits */
173			xm &= ~(DP_MBIT(3) - 1);
174			xe++;
175		}
176		else {
177			/* sticky right shift es bits
178			 */
179			xm = XDPSRS(xm, es);
180			xe += es;
181			assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
182			assert(xe == DP_EMIN);
183		}
184	}
185	if (xm & (DP_MBIT(3) - 1)) {
186		SETCX(IEEE754_INEXACT);
187		if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
188			SETCX(IEEE754_UNDERFLOW);
189		}
190
191		/* inexact must round of 3 bits
192		 */
193		xm = get_rounding(sn, xm);
194		/* adjust exponent for rounding add overflowing
195		 */
196		if (xm >> (DP_MBITS + 3 + 1)) {
197			/* add causes mantissa overflow */
198			xm >>= 1;
199			xe++;
200		}
201	}
202	/* strip grs bits */
203	xm >>= 3;
204
205	assert((xm >> (DP_MBITS + 1)) == 0);	/* no execess */
206	assert(xe >= DP_EMIN);
207
208	if (xe > DP_EMAX) {
209		SETCX(IEEE754_OVERFLOW);
210		SETCX(IEEE754_INEXACT);
211		/* -O can be table indexed by (rm,sn) */
212		switch (ieee754_csr.rm) {
213		case IEEE754_RN:
214			return ieee754dp_inf(sn);
215		case IEEE754_RZ:
216			return ieee754dp_max(sn);
217		case IEEE754_RU:	/* toward +Infinity */
218			if (sn == 0)
219				return ieee754dp_inf(0);
220			else
221				return ieee754dp_max(1);
222		case IEEE754_RD:	/* toward -Infinity */
223			if (sn == 0)
224				return ieee754dp_max(0);
225			else
226				return ieee754dp_inf(1);
227		}
228	}
229	/* gen norm/denorm/zero */
230
231	if ((xm & DP_HIDDEN_BIT) == 0) {
232		/* we underflow (tiny/zero) */
233		assert(xe == DP_EMIN);
234		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
235			SETCX(IEEE754_UNDERFLOW);
236		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
237	} else {
238		assert((xm >> (DP_MBITS + 1)) == 0);	/* no execess */
239		assert(xm & DP_HIDDEN_BIT);
240
241		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
242	}
243}
244