1 /* IEEE754 floating point arithmetic
2 * double precision: common utilities
5 * MIPS floating point support
6 * Copyright (C) 1994-2000 Algorithmics Ltd.
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.
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
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.
22 #include <linux/compiler.h>
24 #include "ieee754dp.h"
26 int ieee754dp_class(union ieee754dp x
)
33 static inline int ieee754dp_isnan(union ieee754dp x
)
35 return ieee754_class_nan(ieee754dp_class(x
));
38 static inline int ieee754dp_issnan(union ieee754dp x
)
42 assert(ieee754dp_isnan(x
));
43 qbit
= (DPMANT(x
) & DP_MBIT(DP_FBITS
- 1)) == DP_MBIT(DP_FBITS
- 1);
44 return ieee754_csr
.nan2008
^ qbit
;
49 * Raise the Invalid Operation IEEE 754 exception
50 * and convert the signaling NaN supplied to a quiet NaN.
52 union ieee754dp __cold
ieee754dp_nanxcpt(union ieee754dp r
)
54 assert(ieee754dp_issnan(r
));
56 ieee754_setcx(IEEE754_INVALID_OPERATION
);
57 if (ieee754_csr
.nan2008
)
58 DPMANT(r
) |= DP_MBIT(DP_FBITS
- 1);
60 r
= ieee754dp_indef();
65 static u64
ieee754dp_get_rounding(int sn
, u64 xm
)
67 /* inexact must round of 3 bits
69 if (xm
& (DP_MBIT(3) - 1)) {
70 switch (ieee754_csr
.rm
) {
74 xm
+= 0x3 + ((xm
>> 3) & 1);
75 /* xm += (xm&0x8)?0x4:0x3 */
77 case FPU_CSR_RU
: /* toward +Infinity */
81 case FPU_CSR_RD
: /* toward -Infinity */
91 /* generate a normal/denormal number with over,under handling
93 * xe is an unbiased exponent
94 * xm is 3bit extended precision value.
96 union ieee754dp
ieee754dp_format(int sn
, int xe
, u64 xm
)
98 assert(xm
); /* we don't gen exact zeros (probably should) */
100 assert((xm
>> (DP_FBITS
+ 1 + 3)) == 0); /* no excess */
101 assert(xm
& (DP_HIDDEN_BIT
<< 3));
104 /* strip lower bits */
105 int es
= DP_EMIN
- xe
;
107 if (ieee754_csr
.nod
) {
108 ieee754_setcx(IEEE754_UNDERFLOW
);
109 ieee754_setcx(IEEE754_INEXACT
);
111 switch(ieee754_csr
.rm
) {
114 return ieee754dp_zero(sn
);
115 case FPU_CSR_RU
: /* toward +Infinity */
117 return ieee754dp_min(0);
119 return ieee754dp_zero(1);
120 case FPU_CSR_RD
: /* toward -Infinity */
122 return ieee754dp_zero(0);
124 return ieee754dp_min(1);
128 if (xe
== DP_EMIN
- 1 &&
129 ieee754dp_get_rounding(sn
, xm
) >> (DP_FBITS
+ 1 + 3))
131 /* Not tiny after rounding */
132 ieee754_setcx(IEEE754_INEXACT
);
133 xm
= ieee754dp_get_rounding(sn
, xm
);
136 xm
&= ~(DP_MBIT(3) - 1);
140 /* sticky right shift es bits
144 assert((xm
& (DP_HIDDEN_BIT
<< 3)) == 0);
145 assert(xe
== DP_EMIN
);
148 if (xm
& (DP_MBIT(3) - 1)) {
149 ieee754_setcx(IEEE754_INEXACT
);
150 if ((xm
& (DP_HIDDEN_BIT
<< 3)) == 0) {
151 ieee754_setcx(IEEE754_UNDERFLOW
);
154 /* inexact must round of 3 bits
156 xm
= ieee754dp_get_rounding(sn
, xm
);
157 /* adjust exponent for rounding add overflowing
159 if (xm
>> (DP_FBITS
+ 3 + 1)) {
160 /* add causes mantissa overflow */
168 assert((xm
>> (DP_FBITS
+ 1)) == 0); /* no excess */
169 assert(xe
>= DP_EMIN
);
172 ieee754_setcx(IEEE754_OVERFLOW
);
173 ieee754_setcx(IEEE754_INEXACT
);
174 /* -O can be table indexed by (rm,sn) */
175 switch (ieee754_csr
.rm
) {
177 return ieee754dp_inf(sn
);
179 return ieee754dp_max(sn
);
180 case FPU_CSR_RU
: /* toward +Infinity */
182 return ieee754dp_inf(0);
184 return ieee754dp_max(1);
185 case FPU_CSR_RD
: /* toward -Infinity */
187 return ieee754dp_max(0);
189 return ieee754dp_inf(1);
192 /* gen norm/denorm/zero */
194 if ((xm
& DP_HIDDEN_BIT
) == 0) {
195 /* we underflow (tiny/zero) */
196 assert(xe
== DP_EMIN
);
197 if (ieee754_csr
.mx
& IEEE754_UNDERFLOW
)
198 ieee754_setcx(IEEE754_UNDERFLOW
);
199 return builddp(sn
, DP_EMIN
- 1 + DP_EBIAS
, xm
);
201 assert((xm
>> (DP_FBITS
+ 1)) == 0); /* no excess */
202 assert(xm
& DP_HIDDEN_BIT
);
204 return builddp(sn
, xe
+ DP_EBIAS
, xm
& ~DP_HIDDEN_BIT
);