GNU Linux-libre 4.19.286-gnu1
[releases.git] / arch / mips / math-emu / ieee754sp.c
1 /* IEEE754 floating point arithmetic
2  * single precision
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 "ieee754sp.h"
25
26 int ieee754sp_class(union ieee754sp x)
27 {
28         COMPXSP;
29         EXPLODEXSP;
30         return xc;
31 }
32
33 static inline int ieee754sp_isnan(union ieee754sp x)
34 {
35         return ieee754_class_nan(ieee754sp_class(x));
36 }
37
38 static inline int ieee754sp_issnan(union ieee754sp x)
39 {
40         int qbit;
41
42         assert(ieee754sp_isnan(x));
43         qbit = (SPMANT(x) & SP_MBIT(SP_FBITS - 1)) == SP_MBIT(SP_FBITS - 1);
44         return ieee754_csr.nan2008 ^ qbit;
45 }
46
47
48 /*
49  * Raise the Invalid Operation IEEE 754 exception
50  * and convert the signaling NaN supplied to a quiet NaN.
51  */
52 union ieee754sp __cold ieee754sp_nanxcpt(union ieee754sp r)
53 {
54         assert(ieee754sp_issnan(r));
55
56         ieee754_setcx(IEEE754_INVALID_OPERATION);
57         if (ieee754_csr.nan2008) {
58                 SPMANT(r) |= SP_MBIT(SP_FBITS - 1);
59         } else {
60                 SPMANT(r) &= ~SP_MBIT(SP_FBITS - 1);
61                 if (!ieee754sp_isnan(r))
62                         SPMANT(r) |= SP_MBIT(SP_FBITS - 2);
63         }
64
65         return r;
66 }
67
68 static unsigned int ieee754sp_get_rounding(int sn, unsigned int xm)
69 {
70         /* inexact must round of 3 bits
71          */
72         if (xm & (SP_MBIT(3) - 1)) {
73                 switch (ieee754_csr.rm) {
74                 case FPU_CSR_RZ:
75                         break;
76                 case FPU_CSR_RN:
77                         xm += 0x3 + ((xm >> 3) & 1);
78                         /* xm += (xm&0x8)?0x4:0x3 */
79                         break;
80                 case FPU_CSR_RU:        /* toward +Infinity */
81                         if (!sn)        /* ?? */
82                                 xm += 0x8;
83                         break;
84                 case FPU_CSR_RD:        /* toward -Infinity */
85                         if (sn) /* ?? */
86                                 xm += 0x8;
87                         break;
88                 }
89         }
90         return xm;
91 }
92
93
94 /* generate a normal/denormal number with over,under handling
95  * sn is sign
96  * xe is an unbiased exponent
97  * xm is 3bit extended precision value.
98  */
99 union ieee754sp ieee754sp_format(int sn, int xe, unsigned int xm)
100 {
101         assert(xm);             /* we don't gen exact zeros (probably should) */
102
103         assert((xm >> (SP_FBITS + 1 + 3)) == 0);        /* no excess */
104         assert(xm & (SP_HIDDEN_BIT << 3));
105
106         if (xe < SP_EMIN) {
107                 /* strip lower bits */
108                 int es = SP_EMIN - xe;
109
110                 if (ieee754_csr.nod) {
111                         ieee754_setcx(IEEE754_UNDERFLOW);
112                         ieee754_setcx(IEEE754_INEXACT);
113
114                         switch(ieee754_csr.rm) {
115                         case FPU_CSR_RN:
116                         case FPU_CSR_RZ:
117                                 return ieee754sp_zero(sn);
118                         case FPU_CSR_RU:      /* toward +Infinity */
119                                 if (sn == 0)
120                                         return ieee754sp_min(0);
121                                 else
122                                         return ieee754sp_zero(1);
123                         case FPU_CSR_RD:      /* toward -Infinity */
124                                 if (sn == 0)
125                                         return ieee754sp_zero(0);
126                                 else
127                                         return ieee754sp_min(1);
128                         }
129                 }
130
131                 if (xe == SP_EMIN - 1 &&
132                     ieee754sp_get_rounding(sn, xm) >> (SP_FBITS + 1 + 3))
133                 {
134                         /* Not tiny after rounding */
135                         ieee754_setcx(IEEE754_INEXACT);
136                         xm = ieee754sp_get_rounding(sn, xm);
137                         xm >>= 1;
138                         /* Clear grs bits */
139                         xm &= ~(SP_MBIT(3) - 1);
140                         xe++;
141                 } else {
142                         /* sticky right shift es bits
143                          */
144                         xm = XSPSRS(xm, es);
145                         xe += es;
146                         assert((xm & (SP_HIDDEN_BIT << 3)) == 0);
147                         assert(xe == SP_EMIN);
148                 }
149         }
150         if (xm & (SP_MBIT(3) - 1)) {
151                 ieee754_setcx(IEEE754_INEXACT);
152                 if ((xm & (SP_HIDDEN_BIT << 3)) == 0) {
153                         ieee754_setcx(IEEE754_UNDERFLOW);
154                 }
155
156                 /* inexact must round of 3 bits
157                  */
158                 xm = ieee754sp_get_rounding(sn, xm);
159                 /* adjust exponent for rounding add overflowing
160                  */
161                 if (xm >> (SP_FBITS + 1 + 3)) {
162                         /* add causes mantissa overflow */
163                         xm >>= 1;
164                         xe++;
165                 }
166         }
167         /* strip grs bits */
168         xm >>= 3;
169
170         assert((xm >> (SP_FBITS + 1)) == 0);    /* no excess */
171         assert(xe >= SP_EMIN);
172
173         if (xe > SP_EMAX) {
174                 ieee754_setcx(IEEE754_OVERFLOW);
175                 ieee754_setcx(IEEE754_INEXACT);
176                 /* -O can be table indexed by (rm,sn) */
177                 switch (ieee754_csr.rm) {
178                 case FPU_CSR_RN:
179                         return ieee754sp_inf(sn);
180                 case FPU_CSR_RZ:
181                         return ieee754sp_max(sn);
182                 case FPU_CSR_RU:        /* toward +Infinity */
183                         if (sn == 0)
184                                 return ieee754sp_inf(0);
185                         else
186                                 return ieee754sp_max(1);
187                 case FPU_CSR_RD:        /* toward -Infinity */
188                         if (sn == 0)
189                                 return ieee754sp_max(0);
190                         else
191                                 return ieee754sp_inf(1);
192                 }
193         }
194         /* gen norm/denorm/zero */
195
196         if ((xm & SP_HIDDEN_BIT) == 0) {
197                 /* we underflow (tiny/zero) */
198                 assert(xe == SP_EMIN);
199                 if (ieee754_csr.mx & IEEE754_UNDERFLOW)
200                         ieee754_setcx(IEEE754_UNDERFLOW);
201                 return buildsp(sn, SP_EMIN - 1 + SP_EBIAS, xm);
202         } else {
203                 assert((xm >> (SP_FBITS + 1)) == 0);    /* no excess */
204                 assert(xm & SP_HIDDEN_BIT);
205
206                 return buildsp(sn, xe + SP_EBIAS, xm & ~SP_HIDDEN_BIT);
207         }
208 }