Remove the implementation of gmp primops from the rts
[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 <math.h>
14
15 /*
16 * Encoding and decoding Doubles. Code based on the HBC code
17 * (lib/fltcode.c).
18 */
19
20 #if IEEE_FLOATING_POINT
21 #define MY_DMINEXP ((DBL_MIN_EXP) - (DBL_MANT_DIG) - 1)
22 /* DMINEXP is defined in values.h on Linux (for example) */
23 #define DHIGHBIT 0x00100000
24 #define DMSBIT 0x80000000
25
26 #define MY_FMINEXP ((FLT_MIN_EXP) - (FLT_MANT_DIG) - 1)
27 #define FHIGHBIT 0x00800000
28 #define FMSBIT 0x80000000
29 #endif
30
31 #if defined(WORDS_BIGENDIAN) || defined(FLOAT_WORDS_BIGENDIAN)
32 #define L 1
33 #define H 0
34 #else
35 #define L 0
36 #define H 1
37 #endif
38
39 #define __abs(a) (( (a) >= 0 ) ? (a) : (-(a)))
40
41 StgDouble
42 __2Int_encodeDouble (I_ j_high, I_ j_low, I_ e)
43 {
44 StgDouble r;
45
46 /* assuming 32 bit ints */
47 ASSERT(sizeof(int ) == 4 );
48
49 r = (StgDouble)((unsigned int)j_high);
50 r *= 4294967296.0; /* exp2f(32); */
51 r += (StgDouble)((unsigned int)j_low);
52
53 /* Now raise to the exponent */
54 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
55 r = ldexp(r, e);
56
57 /* sign is encoded in the size */
58 if (j_high < 0)
59 r = -r;
60
61 return r;
62 }
63
64 /* Special version for words */
65 StgDouble
66 __word_encodeDouble (W_ j, I_ e)
67 {
68 StgDouble r;
69
70 r = (StgDouble)j;
71
72 /* Now raise to the exponent */
73 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
74 r = ldexp(r, e);
75
76 return r;
77 }
78
79 /* Special version for small Integers */
80 StgDouble
81 __int_encodeDouble (I_ j, I_ e)
82 {
83 StgDouble r;
84
85 r = (StgDouble)__abs(j);
86
87 /* Now raise to the exponent */
88 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
89 r = ldexp(r, e);
90
91 /* sign is encoded in the size */
92 if (j < 0)
93 r = -r;
94
95 return r;
96 }
97
98 /* Special version for small Integers */
99 StgFloat
100 __int_encodeFloat (I_ j, I_ e)
101 {
102 StgFloat r;
103
104 r = (StgFloat)__abs(j);
105
106 /* Now raise to the exponent */
107 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
108 r = ldexp(r, e);
109
110 /* sign is encoded in the size */
111 if (j < 0)
112 r = -r;
113
114 return r;
115 }
116
117 /* Special version for small positive Integers */
118 StgFloat
119 __word_encodeFloat (W_ j, I_ e)
120 {
121 StgFloat r;
122
123 r = (StgFloat)j;
124
125 /* Now raise to the exponent */
126 if ( r != 0.0 ) /* Lennart suggests this avoids a bug in MIPS's ldexp */
127 r = ldexp(r, e);
128
129 return r;
130 }
131
132 /* This only supports IEEE floating point */
133
134 void
135 __decodeDouble_2Int (I_ *man_sign, W_ *man_high, W_ *man_low, I_ *exp, StgDouble dbl)
136 {
137 /* Do some bit fiddling on IEEE */
138 unsigned int low, high; /* assuming 32 bit ints */
139 int sign, iexp;
140 union { double d; unsigned int i[2]; } u; /* assuming 32 bit ints, 64 bit double */
141
142 ASSERT(sizeof(unsigned int ) == 4 );
143 ASSERT(sizeof(dbl ) == 8 );
144 ASSERT(sizeof(dbl ) == SIZEOF_DOUBLE);
145
146 u.d = dbl; /* grab chunks of the double */
147 low = u.i[L];
148 high = u.i[H];
149
150 if (low == 0 && (high & ~DMSBIT) == 0) {
151 *man_low = 0;
152 *man_high = 0;
153 *exp = 0L;
154 } else {
155 iexp = ((high >> 20) & 0x7ff) + MY_DMINEXP;
156 sign = high;
157
158 high &= DHIGHBIT-1;
159 if (iexp != MY_DMINEXP) /* don't add hidden bit to denorms */
160 high |= DHIGHBIT;
161 else {
162 iexp++;
163 /* A denorm, normalize the mantissa */
164 while (! (high & DHIGHBIT)) {
165 high <<= 1;
166 if (low & DMSBIT)
167 high++;
168 low <<= 1;
169 iexp--;
170 }
171 }
172 *exp = (I_) iexp;
173 *man_low = low;
174 *man_high = high;
175 *man_sign = (sign < 0) ? -1 : 1;
176 }
177 }
178
179 /* Convenient union types for checking the layout of IEEE 754 types -
180 based on defs in GNU libc <ieee754.h>
181 */
182
183 void
184 __decodeFloat_Int (I_ *man, I_ *exp, StgFloat flt)
185 {
186 /* Do some bit fiddling on IEEE */
187 int high, sign; /* assuming 32 bit ints */
188 union { float f; int i; } u; /* assuming 32 bit float and int */
189
190 ASSERT(sizeof(int ) == 4 );
191 ASSERT(sizeof(flt ) == 4 );
192 ASSERT(sizeof(flt ) == SIZEOF_FLOAT );
193
194 u.f = flt; /* grab the float */
195 high = u.i;
196
197 if ((high & ~FMSBIT) == 0) {
198 *man = 0;
199 *exp = 0;
200 } else {
201 *exp = ((high >> 23) & 0xff) + MY_FMINEXP;
202 sign = high;
203
204 high &= FHIGHBIT-1;
205 if (*exp != MY_FMINEXP) /* don't add hidden bit to denorms */
206 high |= FHIGHBIT;
207 else {
208 (*exp)++;
209 /* A denorm, normalize the mantissa */
210 while (! (high & FHIGHBIT)) {
211 high <<= 1;
212 (*exp)--;
213 }
214 }
215 *man = high;
216 if (sign < 0)
217 *man = - *man;
218 }
219 }
220
221 union stg_ieee754_flt
222 {
223 float f;
224 struct {
225
226 #if WORDS_BIGENDIAN
227 unsigned int negative:1;
228 unsigned int exponent:8;
229 unsigned int mantissa:23;
230 #else
231 unsigned int mantissa:23;
232 unsigned int exponent:8;
233 unsigned int negative:1;
234 #endif
235 } ieee;
236 struct {
237
238 #if WORDS_BIGENDIAN
239 unsigned int negative:1;
240 unsigned int exponent:8;
241 unsigned int quiet_nan:1;
242 unsigned int mantissa:22;
243 #else
244 unsigned int mantissa:22;
245 unsigned int quiet_nan:1;
246 unsigned int exponent:8;
247 unsigned int negative:1;
248 #endif
249 } ieee_nan;
250 };
251
252 /*
253
254 To recap, here's the representation of a double precision
255 IEEE floating point number:
256
257 sign 63 sign bit (0==positive, 1==negative)
258 exponent 62-52 exponent (biased by 1023)
259 fraction 51-0 fraction (bits to right of binary point)
260 */
261
262 union stg_ieee754_dbl
263 {
264 double d;
265 struct {
266
267 #if WORDS_BIGENDIAN
268 unsigned int negative:1;
269 unsigned int exponent:11;
270 unsigned int mantissa0:20;
271 unsigned int mantissa1:32;
272 #else
273 #if FLOAT_WORDS_BIGENDIAN
274 unsigned int mantissa0:20;
275 unsigned int exponent:11;
276 unsigned int negative:1;
277 unsigned int mantissa1:32;
278 #else
279 unsigned int mantissa1:32;
280 unsigned int mantissa0:20;
281 unsigned int exponent:11;
282 unsigned int negative:1;
283 #endif
284 #endif
285 } ieee;
286 /* This format makes it easier to see if a NaN is a signalling NaN. */
287 struct {
288
289 #if WORDS_BIGENDIAN
290 unsigned int negative:1;
291 unsigned int exponent:11;
292 unsigned int quiet_nan:1;
293 unsigned int mantissa0:19;
294 unsigned int mantissa1:32;
295 #else
296 #if FLOAT_WORDS_BIGENDIAN
297 unsigned int mantissa0:19;
298 unsigned int quiet_nan:1;
299 unsigned int exponent:11;
300 unsigned int negative:1;
301 unsigned int mantissa1:32;
302 #else
303 unsigned int mantissa1:32;
304 unsigned int mantissa0:19;
305 unsigned int quiet_nan:1;
306 unsigned int exponent:11;
307 unsigned int negative:1;
308 #endif
309 #endif
310 } ieee_nan;
311 };
312
313 /*
314 * Predicates for testing for extreme IEEE fp values. Used
315 * by the bytecode evaluator and the Prelude.
316 *
317 */
318
319 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
320 #ifdef IEEE_FLOATING_POINT
321
322 StgInt
323 isDoubleNaN(StgDouble d)
324 {
325 union stg_ieee754_dbl u;
326
327 u.d = d;
328
329 return (
330 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
331 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
332 /* and the mantissa non-zero? */
333 );
334 }
335
336 StgInt
337 isDoubleInfinite(StgDouble d)
338 {
339 union stg_ieee754_dbl u;
340
341 u.d = d;
342
343 /* Inf iff exponent is all ones, mantissa all zeros */
344 return (
345 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
346 u.ieee.mantissa0 == 0 &&
347 u.ieee.mantissa1 == 0
348 );
349 }
350
351 StgInt
352 isDoubleDenormalized(StgDouble d)
353 {
354 union stg_ieee754_dbl u;
355
356 u.d = d;
357
358 /* A (single/double/quad) precision floating point number
359 is denormalised iff:
360 - exponent is zero
361 - mantissa is non-zero.
362 - (don't care about setting of sign bit.)
363
364 */
365 return (
366 u.ieee.exponent == 0 &&
367 (u.ieee.mantissa0 != 0 ||
368 u.ieee.mantissa1 != 0)
369 );
370
371 }
372
373 StgInt
374 isDoubleNegativeZero(StgDouble d)
375 {
376 union stg_ieee754_dbl u;
377
378 u.d = d;
379 /* sign (bit 63) set (only) => negative zero */
380
381 return (
382 u.ieee.negative == 1 &&
383 u.ieee.exponent == 0 &&
384 u.ieee.mantissa0 == 0 &&
385 u.ieee.mantissa1 == 0);
386 }
387
388 /* Same tests, this time for StgFloats. */
389
390 /*
391 To recap, here's the representation of a single precision
392 IEEE floating point number:
393
394 sign 31 sign bit (0 == positive, 1 == negative)
395 exponent 30-23 exponent (biased by 127)
396 fraction 22-0 fraction (bits to right of binary point)
397 */
398
399
400 StgInt
401 isFloatNaN(StgFloat f)
402 {
403 union stg_ieee754_flt u;
404 u.f = f;
405
406 /* Floating point NaN iff exponent is all ones, mantissa is
407 non-zero (but see below.) */
408 return (
409 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
410 u.ieee.mantissa != 0);
411 }
412
413 StgInt
414 isFloatInfinite(StgFloat f)
415 {
416 union stg_ieee754_flt u;
417 u.f = f;
418
419 /* A float is Inf iff exponent is max (all ones),
420 and mantissa is min(all zeros.) */
421 return (
422 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
423 u.ieee.mantissa == 0);
424 }
425
426 StgInt
427 isFloatDenormalized(StgFloat f)
428 {
429 union stg_ieee754_flt u;
430 u.f = f;
431
432 /* A (single/double/quad) precision floating point number
433 is denormalised iff:
434 - exponent is zero
435 - mantissa is non-zero.
436 - (don't care about setting of sign bit.)
437
438 */
439 return (
440 u.ieee.exponent == 0 &&
441 u.ieee.mantissa != 0);
442 }
443
444 StgInt
445 isFloatNegativeZero(StgFloat f)
446 {
447 union stg_ieee754_flt u;
448 u.f = f;
449
450 /* sign (bit 31) set (only) => negative zero */
451 return (
452 u.ieee.negative &&
453 u.ieee.exponent == 0 &&
454 u.ieee.mantissa == 0);
455 }
456
457 #else /* ! IEEE_FLOATING_POINT */
458
459 /* Dummy definitions of predicates - they all return false */
460 StgInt isDoubleNaN(d) StgDouble d; { return 0; }
461 StgInt isDoubleInfinite(d) StgDouble d; { return 0; }
462 StgInt isDoubleDenormalized(d) StgDouble d; { return 0; }
463 StgInt isDoubleNegativeZero(d) StgDouble d; { return 0; }
464 StgInt isFloatNaN(f) StgFloat f; { return 0; }
465 StgInt isFloatInfinite(f) StgFloat f; { return 0; }
466 StgInt isFloatDenormalized(f) StgFloat f; { return 0; }
467 StgInt isFloatNegativeZero(f) StgFloat f; { return 0; }
468
469 #endif /* ! IEEE_FLOATING_POINT */