File: /usr/src/linux/arch/mips64/math-emu/ieee754dp.c
1 /* IEEE754 floating point arithmetic
2 * double precision: common utilities
3 */
4 /*
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd. All rights reserved.
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
30 int ieee754dp_class(ieee754dp x)
31 {
32 COMPXDP;
33 EXPLODEXDP;
34 return xc;
35 }
36
37 int ieee754dp_isnan(ieee754dp x)
38 {
39 return ieee754dp_class(x) >= IEEE754_CLASS_SNAN;
40 }
41
42 int ieee754dp_issnan(ieee754dp x)
43 {
44 assert(ieee754dp_isnan(x));
45 if (ieee754_csr.noq)
46 return 1;
47 return !(DPMANT(x) & DP_MBIT(DP_MBITS - 1));
48 }
49
50
51 ieee754dp ieee754dp_xcpt(ieee754dp r, const char *op, ...)
52 {
53 struct ieee754xctx ax;
54 if (!TSTX())
55 return r;
56
57 ax.op = op;
58 ax.rt = IEEE754_RT_DP;
59 ax.rv.dp = r;
60 va_start(ax.ap, op);
61 ieee754_xcpt(&ax);
62 return ax.rv.dp;
63 }
64
65 ieee754dp ieee754dp_nanxcpt(ieee754dp r, const char *op, ...)
66 {
67 struct ieee754xctx ax;
68
69 assert(ieee754dp_isnan(r));
70
71 if (!ieee754dp_issnan(r)) /* QNAN does not cause invalid op !! */
72 return r;
73
74 if (!SETCX(IEEE754_INVALID_OPERATION)) {
75 /* not enabled convert to a quiet NaN */
76 if (ieee754_csr.noq)
77 return r;
78 DPMANT(r) |= DP_MBIT(DP_MBITS - 1);
79 return r;
80 }
81
82 ax.op = op;
83 ax.rt = 0;
84 ax.rv.dp = r;
85 va_start(ax.ap, op);
86 ieee754_xcpt(&ax);
87 return ax.rv.dp;
88 }
89
90 ieee754dp ieee754dp_bestnan(ieee754dp x, ieee754dp y)
91 {
92 assert(ieee754dp_isnan(x));
93 assert(ieee754dp_isnan(y));
94
95 if (DPMANT(x) > DPMANT(y))
96 return x;
97 else
98 return y;
99 }
100
101
102 /* generate a normal/denormal number with over,under handeling
103 * sn is sign
104 * xe is an unbiased exponent
105 * xm is 3bit extended precision value.
106 */
107 ieee754dp ieee754dp_format(int sn, int xe, unsigned long long xm)
108 {
109 assert(xm); /* we dont gen exact zeros (probably should) */
110
111 assert((xm >> (DP_MBITS + 1 + 3)) == 0); /* no execess */
112 assert(xm & (DP_HIDDEN_BIT << 3));
113
114 if (xe < DP_EMIN) {
115 /* strip lower bits */
116 int es = DP_EMIN - xe;
117
118 if (ieee754_csr.nod) {
119 SETCX(IEEE754_UNDERFLOW);
120 return ieee754dp_zero(sn);
121 }
122
123 /* sticky right shift es bits
124 */
125 xm = XDPSRS(xm, es);
126 xe += es;
127
128 assert((xm & (DP_HIDDEN_BIT << 3)) == 0);
129 assert(xe == DP_EMIN);
130 }
131 if (xm & (DP_MBIT(3) - 1)) {
132 SETCX(IEEE754_INEXACT);
133 /* inexact must round of 3 bits
134 */
135 switch (ieee754_csr.rm) {
136 case IEEE754_RZ:
137 break;
138 case IEEE754_RN:
139 xm += 0x3 + ((xm >> 3) & 1);
140 /* xm += (xm&0x8)?0x4:0x3 */
141 break;
142 case IEEE754_RU: /* toward +Infinity */
143 if (!sn) /* ?? */
144 xm += 0x8;
145 break;
146 case IEEE754_RD: /* toward -Infinity */
147 if (sn) /* ?? */
148 xm += 0x8;
149 break;
150 }
151 /* adjust exponent for rounding add overflowing
152 */
153 if (xm >> (DP_MBITS + 3 + 1)) { /* add causes mantissa overflow */
154 xm >>= 1;
155 xe++;
156 }
157 }
158 /* strip grs bits */
159 xm >>= 3;
160
161 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */
162 assert(xe >= DP_EMIN);
163
164 if (xe > DP_EMAX) {
165 SETCX(IEEE754_OVERFLOW);
166 /* -O can be table indexed by (rm,sn) */
167 switch (ieee754_csr.rm) {
168 case IEEE754_RN:
169 return ieee754dp_inf(sn);
170 case IEEE754_RZ:
171 return ieee754dp_max(sn);
172 case IEEE754_RU: /* toward +Infinity */
173 if (sn == 0)
174 return ieee754dp_inf(0);
175 else
176 return ieee754dp_max(1);
177 case IEEE754_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 SETCX(IEEE754_UNDERFLOW);
190 return builddp(sn, DP_EMIN - 1 + DP_EBIAS, xm);
191 } else {
192 assert((xm >> (DP_MBITS + 1)) == 0); /* no execess */
193 assert(xm & DP_HIDDEN_BIT);
194
195 return builddp(sn, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT);
196 }
197 }
198