float_emul.cc Source File

Back to the index.

float_emul.cc
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2004-2009 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  int n_frac = 0, n_exp = 0;
53  int i, nan, sign = 0, exponent;
54  double fraction;
55 
56  memset(fvp, 0, sizeof(struct ieee_float_value));
57 
58  /* n_frac and n_exp: */
59  switch (fmt) {
60  case IEEE_FMT_S: n_frac = 23; n_exp = 8; break;
61  case IEEE_FMT_W: n_frac = 31; n_exp = 0; break;
62  case IEEE_FMT_D: n_frac = 52; n_exp = 11; break;
63  case IEEE_FMT_L: n_frac = 63; n_exp = 0; break;
64  default:fatal("ieee_interpret_float_value(): "
65  "unimplemented format %i\n", fmt);
66  }
67 
68  /* exponent: */
69  exponent = 0;
70  switch (fmt) {
71  case IEEE_FMT_W:
72  x &= 0xffffffffULL;
73  case IEEE_FMT_L:
74  break;
75  case IEEE_FMT_S:
76  x &= 0xffffffffULL;
77  case IEEE_FMT_D:
78  exponent = (x >> n_frac) & ((1 << n_exp) - 1);
79  exponent -= (1 << (n_exp-1)) - 1;
80  break;
81  default:fatal("ieee_interpret_float_value(): unimplemented "
82  "format %i\n", fmt);
83  }
84 
85  /* nan: */
86  nan = 0;
87  switch (fmt) {
88  case IEEE_FMT_S:
89  if (x == 0x7fffffffULL || x == 0x7fbfffffULL)
90  nan = 1;
91  break;
92  case IEEE_FMT_D:
93  if (x == 0x7fffffffffffffffULL ||
94  x == 0x7ff7ffffffffffffULL)
95  nan = 1;
96  break;
97  }
98 
99  if (nan) {
100  fvp->f = 1.0;
101  goto no_reasonable_result;
102  }
103 
104  /* fraction: */
105  fraction = 0.0;
106  switch (fmt) {
107  case IEEE_FMT_W:
108  {
109  int32_t r_int = x;
110  fraction = r_int;
111  }
112  break;
113  case IEEE_FMT_L:
114  {
115  int64_t r_int = x;
116  fraction = r_int;
117  }
118  break;
119  case IEEE_FMT_S:
120  case IEEE_FMT_D:
121  /* sign: */
122  sign = (x >> 31) & 1;
123  if (fmt == IEEE_FMT_D)
124  sign = (x >> 63) & 1;
125 
126  fraction = 0.0;
127  for (i=0; i<n_frac; i++) {
128  int bit = (x >> i) & 1;
129  fraction /= 2.0;
130  if (bit)
131  fraction += 1.0;
132  }
133  /* Add implicit bit 0: */
134  fraction = (fraction / 2.0) + 1.0;
135  break;
136  default:fatal("ieee_interpret_float_value(): "
137  "unimplemented format %i\n", fmt);
138  }
139 
140  /* form the value: */
141  fvp->f = fraction;
142 
143 #ifdef IEEE_DEBUG
144  fatal("{ ieee: x=%016"PRIx64" sign=%i exponent=%i frac=%f ",
145  (uint64_t) x, sign, exponent, fraction);
146 #endif
147 
148  /* TODO: this is awful for exponents of large magnitude. */
149  if (exponent > 0) {
150  /*
151  * NOTE / TODO:
152  *
153  * This is an ulgy workaround on Alpha, where it seems that
154  * multiplying by 2, 1024 times causes a floating point
155  * exception. (Triggered by running for example NetBSD/pmax
156  * 2.0 emulated on an Alpha host.)
157  */
158  if (exponent == 1024)
159  exponent = 1023;
160 
161  while (exponent-- > 0)
162  fvp->f *= 2.0;
163  } else if (exponent < 0) {
164  while (exponent++ < 0)
165  fvp->f /= 2.0;
166  }
167 
168  if (sign)
169  fvp->f = -fvp->f;
170 
171 no_reasonable_result:
172  fvp->nan = nan;
173 
174 #ifdef IEEE_DEBUG
175  fatal("f=%f }\n", fvp->f);
176 #endif
177 }
178 
179 
180 /*
181  * ieee_store_float_value():
182  *
183  * Generates a 64-bit IEEE-formated value in a specific format.
184  */
185 uint64_t ieee_store_float_value(double nf, int fmt, int nan)
186 {
187  int n_frac = 0, n_exp = 0, signofs=0;
188  int i, exponent;
189  uint64_t r = 0, r2;
190  int64_t r3;
191 
192  /* n_frac and n_exp: */
193  switch (fmt) {
194  case IEEE_FMT_S: n_frac = 23; n_exp = 8; signofs = 31; break;
195  case IEEE_FMT_W: n_frac = 31; n_exp = 0; signofs = 31; break;
196  case IEEE_FMT_D: n_frac = 52; n_exp = 11; signofs = 63; break;
197  case IEEE_FMT_L: n_frac = 63; n_exp = 0; signofs = 63; break;
198  default:fatal("ieee_store_float_value(): unimplemented format"
199  " %i\n", fmt);
200  }
201 
202  if ((fmt == IEEE_FMT_S || fmt == IEEE_FMT_D) && nan)
203  goto store_nan;
204 
205  /* fraction: */
206  switch (fmt) {
207  case IEEE_FMT_W:
208  case IEEE_FMT_L:
209  /*
210  * This causes an implicit conversion of double to integer.
211  * If nf < 0.0, then r2 will begin with a sequence of binary
212  * 1's, which is ok.
213  */
214  r3 = (int64_t) nf;
215  r2 = r3;
216  r |= r2;
217 
218  if (fmt == IEEE_FMT_W)
219  r &= 0xffffffffULL;
220  break;
221  case IEEE_FMT_S:
222  case IEEE_FMT_D:
223 #ifdef IEEE_DEBUG
224  fatal("{ ieee store f=%f ", nf);
225 #endif
226  /* sign bit: */
227  if (nf < 0.0) {
228  r |= ((uint64_t)1 << signofs);
229  nf = -nf;
230  }
231 
232  /*
233  * How to convert back from double to exponent + fraction:
234  * The fraction should be 1.xxx, that is
235  * 1.0 <= fraction < 2.0
236  *
237  * This method is very slow but should work:
238  * (TODO: Fix the performance problem!)
239  */
240  exponent = 0;
241  while (nf < 1.0 && exponent > -1023) {
242  nf *= 2.0;
243  exponent --;
244  }
245  while (nf >= 2.0 && exponent < 1023) {
246  nf /= 2.0;
247  exponent ++;
248  }
249 
250  /* Here: 1.0 <= nf < 2.0 */
251 #ifdef IEEE_DEBUG
252  fatal(" nf=%f", nf);
253 #endif
254  nf -= 1.0; /* remove implicit first bit */
255  for (i=n_frac-1; i>=0; i--) {
256  nf *= 2.0;
257  if (nf >= 1.0) {
258  r |= ((uint64_t)1 << i);
259  nf -= 1.0;
260  }
261  }
262 
263  /* Insert the exponent into the resulting word: */
264  /* (First bias, then make sure it's within range) */
265  exponent += (((uint64_t)1 << (n_exp-1)) - 1);
266  if (exponent < 0)
267  exponent = 0;
268  if (exponent >= ((int64_t)1 << n_exp))
269  exponent = ((int64_t)1 << n_exp) - 1;
270  r |= (uint64_t)exponent << n_frac;
271 
272  /* Special case for 0.0: */
273  if (exponent == 0)
274  r = 0;
275 
276 #ifdef IEEE_DEBUG
277  fatal(" exp=%i, r = %016"PRIx64" }\n", exponent, (uint64_t) r);
278 #endif
279  break;
280  default:/* TODO */
281  fatal("ieee_store_float_value(): unimplemented format %i\n",
282  fmt);
283  }
284 
285 store_nan:
286  if (nan) {
287  if (fmt == IEEE_FMT_S)
288  r = 0x7fffffffULL;
289  else if (fmt == IEEE_FMT_D)
290  r = 0x7fffffffffffffffULL;
291  else
292  r = 0x7fffffffULL;
293  }
294 
295  if (fmt == IEEE_FMT_S || fmt == IEEE_FMT_W)
296  r &= 0xffffffff;
297 
298  return r;
299 }
300 
#define IEEE_FMT_S
Definition: float_emul.h:43
void fatal(const char *fmt,...)
Definition: main.cc:152
#define IEEE_FMT_W
Definition: float_emul.h:45
#define IEEE_FMT_L
Definition: float_emul.h:46
#define IEEE_FMT_D
Definition: float_emul.h:44
uint64_t ieee_store_float_value(double nf, int fmt, int nan)
Definition: float_emul.cc:185
void ieee_interpret_float_value(uint64_t x, struct ieee_float_value *fvp, int fmt)
Definition: float_emul.cc:49

Generated on Sun Sep 30 2018 16:05:18 for GXemul by doxygen 1.8.13