float_emul.cc Source File

Back to the index.

float_emul.cc
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2004-2018 Anders Gavare. All rights reserved.
3  *
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions are met:
6  *
7  * 1. Redistributions of source code must retain the above copyright
8  * notice, this list of conditions and the following disclaimer.
9  * 2. Redistributions in binary form must reproduce the above copyright
10  * notice, this list of conditions and the following disclaimer in the
11  * documentation and/or other materials provided with the distribution.
12  * 3. The name of the author may not be used to endorse or promote products
13  * derived from this software without specific prior written permission.
14  *
15  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
16  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
17  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
18  * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
19  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
21  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
22  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25  * SUCH DAMAGE.
26  *
27  *
28  * Floating point emulation routines.
29  */
30 
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string.h>
34 #include <math.h>
35 
36 #include "float_emul.h"
37 #include "misc.h"
38 
39 
40 /* #define IEEE_DEBUG */
41 
42 
43 /*
44  * ieee_interpret_float_value():
45  *
46  * Interprets a float value from binary IEEE format into an ieee_float_value
47  * struct.
48  */
49 void ieee_interpret_float_value(uint64_t x, struct ieee_float_value *fvp,
50  int fmt)
51 {
52  memset(fvp, 0, sizeof(struct ieee_float_value));
53 
54 #if 0
55  // HACK: Use the host's float/double representation:
56  switch (fmt) {
57  case IEEE_FMT_S:
58  {
59  uint32_t x2 = x;
60  void* p = (void*) &x2;
61  float *pf = (float*) p;
62  fvp->f = *pf;
63  }
64  break;
65 
66  case IEEE_FMT_D:
67  {
68  void* p = (void*) &x;
69  double *pf = (double*) p;
70  fvp->f = *pf;
71  }
72  break;
73 
74  case IEEE_FMT_W:
75  case IEEE_FMT_L:
76  {
77  fvp->f = x;
78  }
79  break;
80 
81  default:fatal("ieee_interpret_float_value(): "
82  "unimplemented format %i\n", fmt);
83  }
84 
85  fvp->nan = isnan(fvp->f);
86 #else
87  int n_frac = 0, n_exp = 0;
88  int i, nan, sign = 0, exponent;
89  double fraction;
90 
91  /* n_frac and n_exp: */
92  switch (fmt) {
93  case IEEE_FMT_S: n_frac = 23; n_exp = 8; break;
94  case IEEE_FMT_W: n_frac = 31; n_exp = 0; break;
95  case IEEE_FMT_D: n_frac = 52; n_exp = 11; break;
96  case IEEE_FMT_L: n_frac = 63; n_exp = 0; break;
97  default:fatal("ieee_interpret_float_value(): "
98  "unimplemented format %i\n", fmt);
99  }
100 
101  /* Get the Exponent: */
102  exponent = 0;
103  switch (fmt) {
104  case IEEE_FMT_W:
105  x &= 0xffffffffULL;
106  case IEEE_FMT_L:
107  break;
108  case IEEE_FMT_S:
109  x &= 0xffffffffULL;
110  case IEEE_FMT_D:
111  exponent = (x >> n_frac) & ((1 << n_exp) - 1);
112  exponent -= (1 << (n_exp-1)) - 1;
113  break;
114  default:fatal("ieee_interpret_float_value(): unimplemented "
115  "format %i\n", fmt);
116  }
117 
118  /* Is this a Not-A-Number? */
119  nan = 0;
120  switch (fmt) {
121  case IEEE_FMT_S:
122  sign = (x >> 31) & 1;
123  if ((x & ~0x80000000ULL) == 0x7f800000ULL) {
124  fvp->f = 1.0 / 0.0;
125  goto zero_or_no_reasonable_result;
126  }
127  if ((x & 0x7f800000ULL) == 0x7f800000ULL)
128  nan = 1;
129  break;
130  case IEEE_FMT_D:
131  sign = (x >> 63) & 1;
132  if ((x & ~0x8000000000000000ULL) == 0x7ff0000000000000ULL) {
133  fvp->f = 1.0 / 0.0;
134  goto zero_or_no_reasonable_result;
135  }
136  if ((x & 0x7ff0000000000000ULL) == 0x7ff0000000000000ULL)
137  nan = 1;
138  break;
139  }
140 
141  if (nan) {
142  fvp->f = NAN;
143  goto no_reasonable_result;
144  }
145 
146  /* Calculate the fraction: */
147  fraction = 0.0;
148 
149  switch (fmt) {
150 
151  case IEEE_FMT_W:
152  {
153  int32_t r_int = x;
154  fraction = r_int;
155  }
156  break;
157 
158  case IEEE_FMT_L:
159  {
160  int64_t r_int = x;
161  fraction = r_int;
162  }
163  break;
164 
165  case IEEE_FMT_S:
166  case IEEE_FMT_D:
167  if (x == 0 ||
168  (fmt == IEEE_FMT_D && x == 0x8000000000000000ULL) ||
169  (fmt == IEEE_FMT_S && x == 0x80000000ULL)) {
170  fvp->f = 0.0;
171  goto zero_or_no_reasonable_result;
172  }
173 
174  fraction = 0.0;
175  for (i=0; i<n_frac; i++) {
176  int bit = (x >> i) & 1;
177  fraction /= 2.0;
178  if (bit)
179  fraction += 1.0;
180  }
181 
182  /* Add implicit bit 0: */
183  fraction = (fraction / 2.0) + 1.0;
184  break;
185 
186  default:fatal("ieee_interpret_float_value(): "
187  "unimplemented format %i\n", fmt);
188  }
189 
190  /* form the value: */
191  fvp->f = fraction;
192 
193 #ifdef IEEE_DEBUG
194  fatal("{ ieee: x=%016"PRIx64" => sign=%i exponent=%i frac=%f ",
195  (uint64_t) x, sign, exponent, fraction);
196 #endif
197 
198  /* TODO: this is awful for exponents of large magnitude. */
199  if (exponent > 0) {
200  /*
201  * NOTE / TODO:
202  *
203  * This is an ulgy workaround on Alpha, where it seems that
204  * multiplying by 2, 1024 times causes a floating point
205  * exception. (Triggered by running for example NetBSD/pmax
206  * 2.0 emulated on an Alpha host.)
207  */
208  if (exponent == 1024)
209  exponent = 1023;
210 
211  while (exponent-- > 0)
212  fvp->f *= 2.0;
213  } else if (exponent < 0) {
214  while (exponent++ < 0)
215  fvp->f /= 2.0;
216  }
217 
218 zero_or_no_reasonable_result:
219  if (sign)
220  fvp->f = -fvp->f;
221 
222 no_reasonable_result:
223  fvp->nan = nan;
224 
225 #ifdef IEEE_DEBUG
226  fatal("nan=%i (f=%f) }\n", nan, fvp->f);
227 #endif
228 
229 #endif
230 }
231 
232 
233 /*
234  * ieee_store_float_value():
235  *
236  * Generates a 64-bit IEEE-formated value in a specific format.
237  */
238 uint64_t ieee_store_float_value(double nf, int fmt)
239 {
240  int n_frac = 0, n_exp = 0, signofs = 0, i, exponent;
241  uint64_t r = 0, r2;
242  int64_t r3;
243 
244  /* n_frac and n_exp: */
245  switch (fmt) {
246  case IEEE_FMT_S: n_frac = 23; n_exp = 8; signofs = 31; break;
247  case IEEE_FMT_W: n_frac = 31; n_exp = 0; signofs = 31; break;
248  case IEEE_FMT_D: n_frac = 52; n_exp = 11; signofs = 63; break;
249  case IEEE_FMT_L: n_frac = 63; n_exp = 0; signofs = 63; break;
250  default:fatal("ieee_store_float_value(): unimplemented format"
251  " %i\n", fmt);
252  }
253 
254  switch (fmt) {
255  case IEEE_FMT_W:
256  case IEEE_FMT_L:
257  /*
258  * This causes an implicit conversion of double to integer.
259  * If nf < 0.0, then r2 will begin with a sequence of binary
260  * 1's, which is ok.
261  */
262  r3 = (int64_t) nf;
263  r2 = r3;
264  r |= r2;
265  break;
266  case IEEE_FMT_S:
267  case IEEE_FMT_D:
268  /* sign bit: */
269  if (signbit(nf)) {
270  r |= ((uint64_t)1 << signofs);
271  }
272 
273  // printf("fpclassify(nf) = %i\n", fpclassify(nf));
274  switch (fpclassify(nf)) {
275  case FP_INFINITE:
276  if (fmt == IEEE_FMT_D)
277  r |= 0x7ff0000000000000ULL;
278  else
279  r |= 0x7f800000ULL;
280  break;
281  case FP_NAN:
282  if (fmt == IEEE_FMT_D)
283  r |= 0x7fffffffffffffffULL;
284  else
285  r |= 0x7fffffffULL;
286  break;
287  case FP_NORMAL:
288  if (signbit(nf))
289  nf = -nf;
290 
291  /*
292  * How to convert back from double to exponent + fraction:
293  * The fraction should be 1.xxx, that is
294  * 1.0 <= fraction < 2.0
295  *
296  * This method is very slow but should work:
297  * (TODO: Fix the performance problem!)
298  */
299  exponent = 0;
300  while (nf < 1.0 && exponent > -1023) {
301  nf *= 2.0;
302  exponent --;
303  }
304  while (nf >= 2.0 && exponent < 1023) {
305  nf /= 2.0;
306  exponent ++;
307  }
308 
309  /* Here: 1.0 <= nf < 2.0 */
310  nf -= 1.0; /* remove implicit first bit */
311  for (i=n_frac-1; i>=0; i--) {
312  nf *= 2.0;
313  if (nf >= 1.0) {
314  r |= ((uint64_t)1 << i);
315  nf -= 1.0;
316  }
317  }
318 
319  /* Insert the exponent into the resulting word: */
320  /* (First bias, then make sure it's within range) */
321  exponent += (((uint64_t)1 << (n_exp-1)) - 1);
322  if (exponent < 0)
323  exponent = 0;
324  if (exponent >= ((int64_t)1 << n_exp))
325  exponent = ((int64_t)1 << n_exp) - 1;
326  r |= (uint64_t)exponent << n_frac;
327 
328  /* Special case for 0.0: */
329  if (exponent == 0)
330  r = 0;
331  break;
332  case FP_SUBNORMAL:
333  // TODO
334  break;
335  case FP_ZERO:
336  // r already has zeros in the lowest bits. Done.
337  break;
338  }
339  break;
340  default:/* TODO */
341  fatal("ieee_store_float_value(): unimplemented format %i\n",
342  fmt);
343  }
344 
345  if (fmt == IEEE_FMT_S || fmt == IEEE_FMT_W)
346  r = (uint32_t) r;
347 
348  return r;
349 }
350 
IEEE_FMT_D
#define IEEE_FMT_D
Definition: float_emul.h:44
ieee_interpret_float_value
void ieee_interpret_float_value(uint64_t x, struct ieee_float_value *fvp, int fmt)
Definition: float_emul.cc:49
fatal
void fatal(const char *fmt,...)
Definition: main.cc:152
misc.h
ieee_float_value::nan
int nan
Definition: float_emul.h:40
ieee_float_value
Definition: float_emul.h:38
float_emul.h
ieee_float_value::f
double f
Definition: float_emul.h:39
IEEE_FMT_S
#define IEEE_FMT_S
Definition: float_emul.h:43
IEEE_FMT_W
#define IEEE_FMT_W
Definition: float_emul.h:45
IEEE_FMT_L
#define IEEE_FMT_L
Definition: float_emul.h:46
ieee_store_float_value
uint64_t ieee_store_float_value(double nf, int fmt)
Definition: float_emul.cc:238

Generated on Tue Mar 24 2020 14:04:48 for GXemul by doxygen 1.8.17