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 *
8 *  This program is free software; you can distribute it and/or modify it
9 *  under the terms of the GNU General Public License (Version 2) as
10 *  published by the Free Software Foundation.
11 *
12 *  This program is distributed in the hope it will be useful, but WITHOUT
13 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 *  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15 *  for more details.
16 *
17 *  You should have received a copy of the GNU General Public License along
18 *  with this program; if not, write to the Free Software Foundation, Inc.,
19 *  51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA.
20 */
21
22#include <linux/compiler.h>
23
24#include "ieee754dp.h"
25
26int ieee754dp_class(union ieee754dp x)
27{
28	COMPXDP;
29	EXPLODEXDP;
30	return xc;
31}
32
33static inline int ieee754dp_isnan(union ieee754dp x)
34{
35	return ieee754_class_nan(ieee754dp_class(x));
36}
37
38static inline int ieee754dp_issnan(union ieee754dp x)
39{
40	assert(ieee754dp_isnan(x));
41	return (DPMANT(x) & DP_MBIT(DP_FBITS - 1)) == DP_MBIT(DP_FBITS - 1);
42}
43
44
45/*
46 * Raise the Invalid Operation IEEE 754 exception
47 * and convert the signaling NaN supplied to a quiet NaN.
48 */
49union ieee754dp __cold ieee754dp_nanxcpt(union ieee754dp r)
50{
51	assert(ieee754dp_issnan(r));
52
53	ieee754_setcx(IEEE754_INVALID_OPERATION);
54	return ieee754dp_indef();
55}
56
57static u64 ieee754dp_get_rounding(int sn, u64 xm)
58{
59	/* inexact must round of 3 bits
60	 */
61	if (xm & (DP_MBIT(3) - 1)) {
62		switch (ieee754_csr.rm) {
63		case FPU_CSR_RZ:
64			break;
65		case FPU_CSR_RN:
66			xm += 0x3 + ((xm >> 3) & 1);
67			/* xm += (xm&0x8)?0x4:0x3 */
68			break;
69		case FPU_CSR_RU:	/* toward +Infinity */
70			if (!sn)	/* ?? */
71				xm += 0x8;
72			break;
73		case FPU_CSR_RD:	/* toward -Infinity */
74			if (sn) /* ?? */
75				xm += 0x8;
76			break;
77		}
78	}
79	return xm;
80}
81
82
83/* generate a normal/denormal number with over,under handling
84 * sn is sign
85 * xe is an unbiased exponent
86 * xm is 3bit extended precision value.
87 */
88union ieee754dp ieee754dp_format(int sn, int xe, u64 xm)
89{
90	assert(xm);		/* we don't gen exact zeros (probably should) */
91
92	assert((xm >> (DP_FBITS + 1 + 3)) == 0);	/* no execess */
93	assert(xm & (DP_HIDDEN_BIT << 3));
94
95	if (xe < DP_EMIN) {
96		/* strip lower bits */
97		int es = DP_EMIN - xe;
98
99		if (ieee754_csr.nod) {
100			ieee754_setcx(IEEE754_UNDERFLOW);
101			ieee754_setcx(IEEE754_INEXACT);
102
103			switch(ieee754_csr.rm) {
104			case FPU_CSR_RN:
105			case FPU_CSR_RZ:
106				return ieee754dp_zero(sn);
107			case FPU_CSR_RU:    /* toward +Infinity */
108				if (sn == 0)
109					return ieee754dp_min(0);
110				else
111					return ieee754dp_zero(1);
112			case FPU_CSR_RD:    /* toward -Infinity */
113				if (sn == 0)
114					return ieee754dp_zero(0);
115				else
116					return ieee754dp_min(1);
117			}
118		}
119
120		if (xe == DP_EMIN - 1 &&
121		    ieee754dp_get_rounding(sn, xm) >> (DP_FBITS + 1 + 3))
122		{
123			/* Not tiny after rounding */
124			ieee754_setcx(IEEE754_INEXACT);
125			xm = ieee754dp_get_rounding(sn, xm);
126			xm >>= 1;
127			/* Clear grs bits */
128			xm &= ~(DP_MBIT(3) - 1);
129			xe++;
130		}
131		else {
132			/* sticky right shift es bits
133			 */
134			xm = XDPSRS(xm, es);
135			xe += es;
136			assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
137			assert(xe == DP_EMIN);
138		}
139	}
140	if (xm & (DP_MBIT(3) - 1)) {
141		ieee754_setcx(IEEE754_INEXACT);
142		if ((xm & (DP_HIDDEN_BIT << 3)) == 0) {
143			ieee754_setcx(IEEE754_UNDERFLOW);
144		}
145
146		/* inexact must round of 3 bits
147		 */
148		xm = ieee754dp_get_rounding(sn, xm);
149		/* adjust exponent for rounding add overflowing
150		 */
151		if (xm >> (DP_FBITS + 3 + 1)) {
152			/* add causes mantissa overflow */
153			xm >>= 1;
154			xe++;
155		}
156	}
157	/* strip grs bits */
158	xm >>= 3;
159
160	assert((xm >> (DP_FBITS + 1)) == 0);	/* no execess */
161	assert(xe >= DP_EMIN);
162
163	if (xe > DP_EMAX) {
164		ieee754_setcx(IEEE754_OVERFLOW);
165		ieee754_setcx(IEEE754_INEXACT);
166		/* -O can be table indexed by (rm,sn) */
167		switch (ieee754_csr.rm) {
168		case FPU_CSR_RN:
169			return ieee754dp_inf(sn);
170		case FPU_CSR_RZ:
171			return ieee754dp_max(sn);
172		case FPU_CSR_RU:	/* toward +Infinity */
173			if (sn == 0)
174				return ieee754dp_inf(0);
175			else
176				return ieee754dp_max(1);
177		case FPU_CSR_RD:	/* toward -Infinity */
178			if (sn == 0)
179				return ieee754dp_max(0);
180			else
181				return ieee754dp_inf(1);
182		}
183	}
184	/* gen norm/denorm/zero */
185
186	if ((xm & DP_HIDDEN_BIT) == 0) {
187		/* we underflow (tiny/zero) */
188		assert(xe == DP_EMIN);
189		if (ieee754_csr.mx & IEEE754_UNDERFLOW)
190			ieee754_setcx(IEEE754_UNDERFLOW);
191		return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
192	} else {
193		assert((xm >> (DP_FBITS + 1)) == 0);	/* no execess */
194		assert(xm & DP_HIDDEN_BIT);
195
196		return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
197	}
198}
199