trunk/src/float_emul.c

/*
 *  Copyright (C) 2004-2007  Anders Gavare.  All rights reserved.
 *
 *  Redistribution and use in source and binary forms, with or without
 *  modification, are permitted provided that the following conditions are met:
 *
 *  1. Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *  2. Redistributions in binary form must reproduce the above copyright  
 *     notice, this list of conditions and the following disclaimer in the 
 *     documentation and/or other materials provided with the distribution.
 *  3. The name of the author may not be used to endorse or promote products
 *     derived from this software without specific prior written permission.
 *
 *  THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 *  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 *  ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE   
 *  FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 *  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 *  OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 *  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 *  OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 *  SUCH DAMAGE.
 *
 *  $Id: float_emul.c,v 1.9 2006/12/30 13:30:52 debug Exp $
 *
 *  Floating point emulation routines.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "float_emul.h"
#include "misc.h"


/*  #define IEEE_DEBUG  */


/*
 *  ieee_interpret_float_value():
 *
 *  Interprets a float value from binary IEEE format into an ieee_float_value
 *  struct.
 */
void ieee_interpret_float_value(uint64_t x, struct ieee_float_value *fvp,
        int fmt)
{
        int n_frac = 0, n_exp = 0;
        int i, nan, sign = 0, exponent;
        double fraction;

        memset(fvp, 0, sizeof(struct ieee_float_value));

        /*  n_frac and n_exp:  */
        switch (fmt) {
        case IEEE_FMT_S:        n_frac = 23; n_exp = 8; break;
        case IEEE_FMT_W:        n_frac = 31; n_exp = 0; break;
        case IEEE_FMT_D:        n_frac = 52; n_exp = 11; break;
        case IEEE_FMT_L:        n_frac = 63; n_exp = 0; break;
        default:fatal("ieee_interpret_float_value(): "
                    "unimplemented format %i\n", fmt);
        }

        /*  exponent:  */
        exponent = 0;
        switch (fmt) {
        case IEEE_FMT_W:
                x &= 0xffffffffULL;
        case IEEE_FMT_L:
                break;
        case IEEE_FMT_S:
                x &= 0xffffffffULL;
        case IEEE_FMT_D:
                exponent = (x >> n_frac) & ((1 << n_exp) - 1);
                exponent -= (1 << (n_exp-1)) - 1;
                break;
        default:fatal("ieee_interpret_float_value(): unimplemented "
                    "format %i\n", fmt);
        }

        /*  nan:  */
        nan = 0;
        switch (fmt) {
        case IEEE_FMT_S:
                if (x == 0x7fffffffULL || x == 0x7fbfffffULL)
                        nan = 1;
                break;
        case IEEE_FMT_D:
                if (x == 0x7fffffffffffffffULL ||
                    x == 0x7ff7ffffffffffffULL)
                        nan = 1;
                break;
        }

        if (nan) {
                fvp->f = 1.0;
                goto no_reasonable_result;
        }

        /*  fraction:  */
        fraction = 0.0;
        switch (fmt) {
        case IEEE_FMT_W:
                {
                        int32_t r_int = x;
                        fraction = r_int;
                }
                break;
        case IEEE_FMT_L:
                {
                        int64_t r_int = x;
                        fraction = r_int;
                }
                break;
        case IEEE_FMT_S:
        case IEEE_FMT_D:
                /*  sign:  */
                sign = (x >> 31) & 1;
                if (fmt == IEEE_FMT_D)
                        sign = (x >> 63) & 1;

                fraction = 0.0;
                for (i=0; i<n_frac; i++) {
                        int bit = (x >> i) & 1;
                        fraction /= 2.0;
                        if (bit)
                                fraction += 1.0;
                }
                /*  Add implicit bit 0:  */
                fraction = (fraction / 2.0) + 1.0;
                break;
        default:fatal("ieee_interpret_float_value(): "
                    "unimplemented format %i\n", fmt);
        }

        /*  form the value:  */
        fvp->f = fraction;

#ifdef IEEE_DEBUG
        fatal("{ ieee: x=%016"PRIx64" sign=%i exponent=%i frac=%f ",
            (uint64_t) x, sign, exponent, fraction);
#endif

        /*  TODO: this is awful for exponents of large magnitude.  */
        if (exponent > 0) {
                /*
                 *  NOTE / TODO:
                 *
                 *  This is an ulgy workaround on Alpha, where it seems that
                 *  multiplying by 2, 1024 times causes a floating point
                 *  exception. (Triggered by running for example NetBSD/pmax
                 *  2.0 emulated on an Alpha host.)
                 */
                if (exponent == 1024)
                        exponent = 1023;

                while (exponent-- > 0)
                        fvp->f *= 2.0;
        } else if (exponent < 0) {
                while (exponent++ < 0)
                        fvp->f /= 2.0;
        }

        if (sign)
                fvp->f = -fvp->f;

no_reasonable_result:
        fvp->nan = nan;

#ifdef IEEE_DEBUG
        fatal("f=%f }\n", fvp->f);
#endif
}


/*
 *  ieee_store_float_value():
 *
 *  Generates a 64-bit IEEE-formated value in a specific format.
 */
uint64_t ieee_store_float_value(double nf, int fmt, int nan)
{
        int n_frac = 0, n_exp = 0, signofs=0;
        int i, exponent;
        uint64_t r = 0, r2;
        int64_t r3;

        /*  n_frac and n_exp:  */
        switch (fmt) {
        case IEEE_FMT_S:        n_frac = 23; n_exp = 8; signofs = 31; break;
        case IEEE_FMT_W:        n_frac = 31; n_exp = 0; signofs = 31; break;
        case IEEE_FMT_D:        n_frac = 52; n_exp = 11; signofs = 63; break;
        case IEEE_FMT_L:        n_frac = 63; n_exp = 0; signofs = 63; break;
        default:fatal("ieee_store_float_value(): unimplemented format"
                    " %i\n", fmt);
        }

        if ((fmt == IEEE_FMT_S || fmt == IEEE_FMT_D) && nan)
                goto store_nan;

        /*  fraction:  */
        switch (fmt) {
        case IEEE_FMT_W:
        case IEEE_FMT_L:
                /*
                 *  This causes an implicit conversion of double to integer.
                 *  If nf < 0.0, then r2 will begin with a sequence of binary
                 *  1's, which is ok.
                 */
                r3 = nf;
                r2 = r3;
                r |= r2;

                if (fmt == IEEE_FMT_W)
                        r &= 0xffffffffULL;
                break;
        case IEEE_FMT_S:
        case IEEE_FMT_D:
#ifdef IEEE_DEBUG
                fatal("{ ieee store f=%f ", nf);
#endif
                /*  sign bit:  */
                if (nf < 0.0) {
                        r |= ((uint64_t)1 << signofs);
                        nf = -nf;
                }

                /*
                 *  How to convert back from double to exponent + fraction:
                 *  The fraction should be 1.xxx, that is
                 *  1.0 <= fraction < 2.0
                 *
                 *  This method is very slow but should work:
                 *  (TODO: Fix the performance problem!)
                 */
                exponent = 0;
                while (nf < 1.0 && exponent > -1023) {
                        nf *= 2.0;
                        exponent --;
                }
                while (nf >= 2.0 && exponent < 1023) {
                        nf /= 2.0;
                        exponent ++;
                }

                /*  Here:   1.0 <= nf < 2.0  */
#ifdef IEEE_DEBUG
                fatal(" nf=%f", nf);
#endif
                nf -= 1.0;      /*  remove implicit first bit  */
                for (i=n_frac-1; i>=0; i--) {
                        nf *= 2.0;
                        if (nf >= 1.0) {
                                r |= ((uint64_t)1 << i);
                                nf -= 1.0;
                        }
                }

                /*  Insert the exponent into the resulting word:  */
                /*  (First bias, then make sure it's within range)  */
                exponent += (((uint64_t)1 << (n_exp-1)) - 1);
                if (exponent < 0)
                        exponent = 0;
                if (exponent >= ((int64_t)1 << n_exp))
                        exponent = ((int64_t)1 << n_exp) - 1;
                r |= (uint64_t)exponent << n_frac;

                /*  Special case for 0.0:  */
                if (exponent == 0)
                        r = 0;

#ifdef IEEE_DEBUG
                fatal(" exp=%i, r = %016"PRIx64" }\n", exponent, (uint64_t) r);
#endif
                break;
        default:/*  TODO  */
                fatal("ieee_store_float_value(): unimplemented format %i\n",
                    fmt);
        }

store_nan:
        if (nan) {
                if (fmt == IEEE_FMT_S)
                        r = 0x7fffffffULL;
                else if (fmt == IEEE_FMT_D)
                        r = 0x7fffffffffffffffULL;
                else
                        r = 0x7fffffffULL;
        }

        if (fmt == IEEE_FMT_S || fmt == IEEE_FMT_W)
                r &= 0xffffffff;

        return r;
}

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