Add 'hadrian/' from commit '45f3bff7016a2a0cd9a5455a882ced984655e90b'
[ghc.git] / rts / StgPrimFloat.c
1 /* -----------------------------------------------------------------------------
2 *
3 * (c) Lennart Augustsson
4 * (c) The GHC Team, 1998-2000
5 *
6 * Miscellaneous support for floating-point primitives
7 *
8 * ---------------------------------------------------------------------------*/
9
10 #include "PosixSource.h"
11 #include "Rts.h"
12
13 #include "StgPrimFloat.h"
14
15 #include <math.h>
16 #include <float.h>
17
18 #define IEEE_FLOATING_POINT 1
19
20 #if FLT_RADIX != 2
21 # error FLT_RADIX != 2 not supported
22 #endif
23
24 /*
25 * Encoding and decoding Doubles. Code based on the HBC code
26 * (lib/fltcode.c).
27 */
28
29 #if IEEE_FLOATING_POINT
30 #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
31 /* DMINEXP is defined in values.h on Linux (for example) */
32 #define DHIGHBIT 0x00100000
33 #define DMSBIT 0x80000000
34
35 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
36 #define FHIGHBIT 0x00800000
37 #define FMSBIT 0x80000000
38 #endif
39
40 #if defined(WORDS_BIGENDIAN) || defined(FLOAT_WORDS_BIGENDIAN)
41 #define L 1
42 #define H 0
43 #else
44 #define L 0
45 #define H 1
46 #endif
47
48 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
49
50 /* Special version for words */
51 StgDouble
52 __word_encodeDouble (W_ j, I_ e)
53 {
54 StgDouble r;
55
56 r = (StgDouble)j;
57
58 /* Now raise to the exponent */
59 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
60 r = ldexp(r, e);
61
62 return r;
63 }
64
65 /* Special version for small Integers */
66 StgDouble
67 __int_encodeDouble (I_ j, I_ e)
68 {
69 StgDouble r;
70
71 r = (StgDouble)__abs(j);
72
73 /* Now raise to the exponent */
74 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
75 r = ldexp(r, e);
76
77 /* sign is encoded in the size */
78 if (j < 0)
79 r = -r;
80
81 return r;
82 }
83
84 /* Special version for small Integers */
85 StgFloat
86 __int_encodeFloat (I_ j, I_ e)
87 {
88 StgFloat r;
89
90 r = (StgFloat)__abs(j);
91
92 /* Now raise to the exponent */
93 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
94 r = ldexp(r, e);
95
96 /* sign is encoded in the size */
97 if (j < 0)
98 r = -r;
99
100 return r;
101 }
102
103 /* Special version for small positive Integers */
104 StgFloat
105 __word_encodeFloat (W_ j, I_ e)
106 {
107 StgFloat r;
108
109 r = (StgFloat)j;
110
111 /* Now raise to the exponent */
112 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
113 r = ldexp(r, e);
114
115 return r;
116 }
117
118 /* This only supports IEEE floating point */
119
120 void
121 __decodeDouble_2Int (I_ *man_sign, W_ *man_high, W_ *man_low, I_ *exp, StgDouble dbl)
122 {
123 /* Do some bit fiddling on IEEE */
124 unsigned int low, high; /* assuming 32 bit ints */
125 int sign, iexp;
126 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
127
128 ASSERT(sizeof(unsigned int ) == 4 );
129 ASSERT(sizeof(dbl ) == 8 );
130 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
131
132 u.d = dbl; /* grab chunks of the double */
133 low = u.i[L];
134 high = u.i[H];
135
136 if (low == 0 && (high & ~DMSBIT) == 0) {
137 *man_low = 0;
138 *man_high = 0;
139 *exp = 0L;
140 } else {
141 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
142 sign = high;
143
144 high &= DHIGHBIT-1;
145 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
146 high |= DHIGHBIT;
147 else {
148 iexp++;
149 /* A denorm, normalize the mantissa */
150 while (! (high & DHIGHBIT)) {
151 high <<= 1;
152 if (low & DMSBIT)
153 high++;
154 low <<= 1;
155 iexp--;
156 }
157 }
158 *exp = (I_) iexp;
159 *man_low = low;
160 *man_high = high;
161 *man_sign = (sign < 0) ? -1 : 1;
162 }
163 }
164
165 /* This is expected to replace uses of __decodeDouble_2Int() in the long run */
166 StgInt
167 __decodeDouble_Int64 (StgInt64 *const mantissa, const StgDouble dbl)
168 {
169 #if 0
170 // We can't use this yet as-is, see ticket #9810
171 if (dbl) {
172 int exp = 0;
173 *mantissa = (StgInt64)scalbn(frexp(dbl, &exp), DBL_MANT_DIG);
174 return exp-DBL_MANT_DIG;
175 } else {
176 *mantissa = 0;
177 return 0;
178 }
179 #else
180 I_ man_sign = 0;
181 W_ man_high = 0, man_low = 0;
182 I_ exp = 0;
183
184 __decodeDouble_2Int (&man_sign, &man_high, &man_low, &exp, dbl);
185 ASSIGN_Int64((W_*)mantissa, ((((StgInt64)man_high << 32)
186 | (StgInt64)man_low)
187 * (StgInt64)man_sign));
188 return exp;
189 #endif
190 }
191
192 /* Convenient union types for checking the layout of IEEE 754 types -
193 based on defs in GNU libc <ieee754.h>
194 */
195
196 void
197 __decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt)
198 {
199 /* Do some bit fiddling on IEEE */
200 int high, sign; /* assuming 32 bit ints */
201 union { float f; int i; } u; /* assuming 32 bit float and int */
202
203 ASSERT(sizeof(int ) == 4 );
204 ASSERT(sizeof(flt ) == 4 );
205 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
206
207 u.f = flt; /* grab the float */
208 high = u.i;
209
210 if ((high & ~FMSBIT) == 0) {
211 *man = 0;
212 *exp = 0;
213 } else {
214 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
215 sign = high;
216
217 high &= FHIGHBIT-1;
218 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
219 high |= FHIGHBIT;
220 else {
221 (*exp)++;
222 /* A denorm, normalize the mantissa */
223 while (! (high & FHIGHBIT)) {
224 high <<= 1;
225 (*exp)--;
226 }
227 }
228 *man = high;
229 if (sign < 0)
230 *man = - *man;
231 }
232 }
233