0.4.4/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	dpavlin	20	/*
2	dpavlin	34	* Copyright (C) 2004-2007 Anders Gavare. All rights reserved.
3	dpavlin	20	*
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	dpavlin	34	* $Id: float_emul.c,v 1.9 2006/12/30 13:30:52 debug Exp $
28	dpavlin	20	*
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	dpavlin	24	fatal("{ ieee: x=%016"PRIx64" sign=%i exponent=%i frac=%f ",
146			(uint64_t) x, sign, exponent, fraction);
147	dpavlin	20	#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	dpavlin	32	* The fraction should be 1.xxx, that is
236	dpavlin	20	* 1.0 <= fraction < 2.0
237			*
238			* This method is very slow but should work:
239	dpavlin	32	* (TODO: Fix the performance problem!)
240	dpavlin	20	*/
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	dpavlin	24	fatal(" exp=%i, r = %016"PRIx64" }\n", exponent, (uint64_t) r);
279	dpavlin	20	#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