a8f4803772f40489099ec18c1e26a8d97fa8d8a1
[packages/base.git] / cbits / primFloat.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 "HsFFI.h"
11 #include "Rts.h" // XXX wrong (for IEEE_FLOATING_POINT and WORDS_BIGENDIAN)
12
13 #define IEEE_FLOATING_POINT 1
14
15 union stg_ieee754_flt
16 {
17 float f;
18 struct {
19
20 #if WORDS_BIGENDIAN
21 unsigned int negative:1;
22 unsigned int exponent:8;
23 unsigned int mantissa:23;
24 #else
25 unsigned int mantissa:23;
26 unsigned int exponent:8;
27 unsigned int negative:1;
28 #endif
29 } ieee;
30 struct {
31
32 #if WORDS_BIGENDIAN
33 unsigned int negative:1;
34 unsigned int exponent:8;
35 unsigned int quiet_nan:1;
36 unsigned int mantissa:22;
37 #else
38 unsigned int mantissa:22;
39 unsigned int quiet_nan:1;
40 unsigned int exponent:8;
41 unsigned int negative:1;
42 #endif
43 } ieee_nan;
44 };
45
46 /*
47
48 To recap, here's the representation of a double precision
49 IEEE floating point number:
50
51 sign 63 sign bit (0==positive, 1==negative)
52 exponent 62-52 exponent (biased by 1023)
53 fraction 51-0 fraction (bits to right of binary point)
54 */
55
56 union stg_ieee754_dbl
57 {
58 double d;
59 struct {
60
61 #if WORDS_BIGENDIAN
62 unsigned int negative:1;
63 unsigned int exponent:11;
64 unsigned int mantissa0:20;
65 unsigned int mantissa1:32;
66 #else
67 #if FLOAT_WORDS_BIGENDIAN
68 unsigned int mantissa0:20;
69 unsigned int exponent:11;
70 unsigned int negative:1;
71 unsigned int mantissa1:32;
72 #else
73 unsigned int mantissa1:32;
74 unsigned int mantissa0:20;
75 unsigned int exponent:11;
76 unsigned int negative:1;
77 #endif
78 #endif
79 } ieee;
80 /* This format makes it easier to see if a NaN is a signalling NaN. */
81 struct {
82
83 #if WORDS_BIGENDIAN
84 unsigned int negative:1;
85 unsigned int exponent:11;
86 unsigned int quiet_nan:1;
87 unsigned int mantissa0:19;
88 unsigned int mantissa1:32;
89 #else
90 #if FLOAT_WORDS_BIGENDIAN
91 unsigned int mantissa0:19;
92 unsigned int quiet_nan:1;
93 unsigned int exponent:11;
94 unsigned int negative:1;
95 unsigned int mantissa1:32;
96 #else
97 unsigned int mantissa1:32;
98 unsigned int mantissa0:19;
99 unsigned int quiet_nan:1;
100 unsigned int exponent:11;
101 unsigned int negative:1;
102 #endif
103 #endif
104 } ieee_nan;
105 };
106
107 /*
108 * Predicates for testing for extreme IEEE fp values.
109 */
110
111 /* In case you don't suppport IEEE, you'll just get dummy defs.. */
112 #ifdef IEEE_FLOATING_POINT
113
114 HsInt
115 isDoubleNaN(HsDouble d)
116 {
117 union stg_ieee754_dbl u;
118
119 u.d = d;
120
121 return (
122 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
123 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
124 /* and the mantissa non-zero? */
125 );
126 }
127
128 HsInt
129 isDoubleInfinite(HsDouble d)
130 {
131 union stg_ieee754_dbl u;
132
133 u.d = d;
134
135 /* Inf iff exponent is all ones, mantissa all zeros */
136 return (
137 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
138 u.ieee.mantissa0 == 0 &&
139 u.ieee.mantissa1 == 0
140 );
141 }
142
143 HsInt
144 isDoubleDenormalized(HsDouble d)
145 {
146 union stg_ieee754_dbl u;
147
148 u.d = d;
149
150 /* A (single/double/quad) precision floating point number
151 is denormalised iff:
152 - exponent is zero
153 - mantissa is non-zero.
154 - (don't care about setting of sign bit.)
155
156 */
157 return (
158 u.ieee.exponent == 0 &&
159 (u.ieee.mantissa0 != 0 ||
160 u.ieee.mantissa1 != 0)
161 );
162
163 }
164
165 HsInt
166 isDoubleNegativeZero(HsDouble d)
167 {
168 union stg_ieee754_dbl u;
169
170 u.d = d;
171 /* sign (bit 63) set (only) => negative zero */
172
173 return (
174 u.ieee.negative == 1 &&
175 u.ieee.exponent == 0 &&
176 u.ieee.mantissa0 == 0 &&
177 u.ieee.mantissa1 == 0);
178 }
179
180 /* Same tests, this time for HsFloats. */
181
182 /*
183 To recap, here's the representation of a single precision
184 IEEE floating point number:
185
186 sign 31 sign bit (0 == positive, 1 == negative)
187 exponent 30-23 exponent (biased by 127)
188 fraction 22-0 fraction (bits to right of binary point)
189 */
190
191
192 HsInt
193 isFloatNaN(HsFloat f)
194 {
195 union stg_ieee754_flt u;
196 u.f = f;
197
198 /* Floating point NaN iff exponent is all ones, mantissa is
199 non-zero (but see below.) */
200 return (
201 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
202 u.ieee.mantissa != 0);
203 }
204
205 HsInt
206 isFloatInfinite(HsFloat f)
207 {
208 union stg_ieee754_flt u;
209 u.f = f;
210
211 /* A float is Inf iff exponent is max (all ones),
212 and mantissa is min(all zeros.) */
213 return (
214 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
215 u.ieee.mantissa == 0);
216 }
217
218 HsInt
219 isFloatDenormalized(HsFloat f)
220 {
221 union stg_ieee754_flt u;
222 u.f = f;
223
224 /* A (single/double/quad) precision floating point number
225 is denormalised iff:
226 - exponent is zero
227 - mantissa is non-zero.
228 - (don't care about setting of sign bit.)
229
230 */
231 return (
232 u.ieee.exponent == 0 &&
233 u.ieee.mantissa != 0);
234 }
235
236 HsInt
237 isFloatNegativeZero(HsFloat f)
238 {
239 union stg_ieee754_flt u;
240 u.f = f;
241
242 /* sign (bit 31) set (only) => negative zero */
243 return (
244 u.ieee.negative &&
245 u.ieee.exponent == 0 &&
246 u.ieee.mantissa == 0);
247 }
248
249 /*
250 There are glibc versions around with buggy rintf or rint, hence we
251 provide our own. We always round ties to even, so we can be simpler.
252 */
253
254 #define FLT_HIDDEN 0x800000
255 #define FLT_POWER2 0x1000000
256
257 HsFloat
258 rintFloat(HsFloat f)
259 {
260 union stg_ieee754_flt u;
261 u.f = f;
262 /* if real exponent > 22, it's already integral, infinite or nan */
263 if (u.ieee.exponent > 149) /* 22 + 127 */
264 {
265 return u.f;
266 }
267 if (u.ieee.exponent < 126) /* (-1) + 127, abs(f) < 0.5 */
268 {
269 /* only used for rounding to Integral a, so don't care about -0.0 */
270 return 0.0;
271 }
272 /* 0.5 <= abs(f) < 2^23 */
273 unsigned int half, mask, mant, frac;
274 half = 1 << (149 - u.ieee.exponent); /* bit for 0.5 */
275 mask = 2*half - 1; /* fraction bits */
276 mant = u.ieee.mantissa | FLT_HIDDEN; /* add hidden bit */
277 frac = mant & mask; /* get fraction */
278 mant ^= frac; /* truncate mantissa */
279 if ((frac < half) || ((frac == half) && ((mant & (2*half)) == 0)))
280 {
281 /* this means we have to truncate */
282 if (mant == 0)
283 {
284 /* f == ±0.5, return 0.0 */
285 return 0.0;
286 }
287 else
288 {
289 /* remove hidden bit and set mantissa */
290 u.ieee.mantissa = mant ^ FLT_HIDDEN;
291 return u.f;
292 }
293 }
294 else
295 {
296 /* round away from zero, increment mantissa */
297 mant += 2*half;
298 if (mant == FLT_POWER2)
299 {
300 /* next power of 2, increase exponent an set mantissa to 0 */
301 u.ieee.mantissa = 0;
302 u.ieee.exponent += 1;
303 return u.f;
304 }
305 else
306 {
307 /* remove hidden bit and set mantissa */
308 u.ieee.mantissa = mant ^ FLT_HIDDEN;
309 return u.f;
310 }
311 }
312 }
313
314 #define DBL_HIDDEN 0x100000
315 #define DBL_POWER2 0x200000
316 #define LTOP_BIT 0x80000000
317
318 HsDouble
319 rintDouble(HsDouble d)
320 {
321 union stg_ieee754_dbl u;
322 u.d = d;
323 /* if real exponent > 51, it's already integral, infinite or nan */
324 if (u.ieee.exponent > 1074) /* 51 + 1023 */
325 {
326 return u.d;
327 }
328 if (u.ieee.exponent < 1022) /* (-1) + 1023, abs(d) < 0.5 */
329 {
330 /* only used for rounding to Integral a, so don't care about -0.0 */
331 return 0.0;
332 }
333 unsigned int half, mask, mant, frac;
334 if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
335 {
336 /* the fractional part meets the higher part of the mantissa */
337 half = 1 << (1042 - u.ieee.exponent); /* bit for 0.5 */
338 mask = 2*half - 1; /* fraction bits */
339 mant = u.ieee.mantissa0 | DBL_HIDDEN; /* add hidden bit */
340 frac = mant & mask; /* get fraction */
341 mant ^= frac; /* truncate mantissa */
342 if ((frac < half) ||
343 ((frac == half) && (u.ieee.mantissa1 == 0) /* a tie */
344 && ((mant & (2*half)) == 0)))
345 {
346 /* truncate */
347 if (mant == 0)
348 {
349 /* d = ±0.5, return 0.0 */
350 return 0.0;
351 }
352 /* remove hidden bit and set mantissa */
353 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
354 u.ieee.mantissa1 = 0;
355 return u.d;
356 }
357 else /* round away from zero */
358 {
359 /* zero low mantissa bits */
360 u.ieee.mantissa1 = 0;
361 /* increment integer part of mantissa */
362 mant += 2*half;
363 if (mant == DBL_POWER2)
364 {
365 /* power of 2, increment exponent and zero mantissa */
366 u.ieee.mantissa0 = 0;
367 u.ieee.exponent += 1;
368 return u.d;
369 }
370 /* remove hidden bit */
371 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
372 return u.d;
373 }
374 }
375 else
376 {
377 /* 20 <= real exponent < 52, fractional part entirely in mantissa1 */
378 half = 1 << (1074 - u.ieee.exponent); /* bit for 0.5 */
379 mask = 2*half - 1; /* fraction bits */
380 mant = u.ieee.mantissa1; /* no hidden bit here */
381 frac = mant & mask; /* get fraction */
382 mant ^= frac; /* truncate mantissa */
383 if ((frac < half) ||
384 ((frac == half) && /* tie */
385 (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1) /* yuck */
386 : (mant & (2*half)))
387 == 0)))
388 {
389 /* truncate */
390 u.ieee.mantissa1 = mant;
391 return u.d;
392 }
393 else
394 {
395 /* round away from zero */
396 /* increment mantissa */
397 mant += 2*half;
398 u.ieee.mantissa1 = mant;
399 if (mant == 0)
400 {
401 /* low part of mantissa overflowed */
402 /* increment high part of mantissa */
403 mant = u.ieee.mantissa0 + 1;
404 if (mant == DBL_HIDDEN)
405 {
406 /* hit power of 2 */
407 /* zero mantissa */
408 u.ieee.mantissa0 = 0;
409 /* and increment exponent */
410 u.ieee.exponent += 1;
411 return u.d;
412 }
413 else
414 {
415 u.ieee.mantissa0 = mant;
416 return u.d;
417 }
418 }
419 else
420 {
421 return u.d;
422 }
423 }
424 }
425 }
426
427 #else /* ! IEEE_FLOATING_POINT */
428
429 /* Dummy definitions of predicates - they all return false */
430 HsInt isDoubleNaN(d) HsDouble d; { return 0; }
431 HsInt isDoubleInfinite(d) HsDouble d; { return 0; }
432 HsInt isDoubleDenormalized(d) HsDouble d; { return 0; }
433 HsInt isDoubleNegativeZero(d) HsDouble d; { return 0; }
434 HsInt isFloatNaN(f) HsFloat f; { return 0; }
435 HsInt isFloatInfinite(f) HsFloat f; { return 0; }
436 HsInt isFloatDenormalized(f) HsFloat f; { return 0; }
437 HsInt isFloatNegativeZero(f) HsFloat f; { return 0; }
438
439
440 /* For exotic floating point formats, we can't do much */
441 /* We suppose the format has not too many bits */
442 /* I hope nobody tries to build GHC where this is wrong */
443
444 #define FLT_UPP 536870912.0
445
446 HsFloat
447 rintFloat(HsFloat f)
448 {
449 if ((f > FLT_UPP) || (f < (-FLT_UPP)))
450 {
451 return f;
452 }
453 else
454 {
455 int i = (int)f;
456 float g = i;
457 float d = f - g;
458 if (d > 0.5)
459 {
460 return g + 1.0;
461 }
462 if (d == 0.5)
463 {
464 return (i & 1) ? (g + 1.0) : g;
465 }
466 if (d == -0.5)
467 {
468 return (i & 1) ? (g - 1.0) : g;
469 }
470 if (d < -0.5)
471 {
472 return g - 1.0;
473 }
474 return g;
475 }
476 }
477
478 #define DBL_UPP 2305843009213693952.0
479
480 HsDouble
481 rintDouble(HsDouble d)
482 {
483 if ((d > DBL_UPP) || (d < (-DBL_UPP)))
484 {
485 return d;
486 }
487 else
488 {
489 HsInt64 i = (HsInt64)d;
490 double e = i;
491 double r = d - e;
492 if (r > 0.5)
493 {
494 return e + 1.0;
495 }
496 if (r == 0.5)
497 {
498 return (i & 1) ? (e + 1.0) : e;
499 }
500 if (r == -0.5)
501 {
502 return (i & 1) ? (e - 1.0) : e;
503 }
504 if (r < -0.5)
505 {
506 return e - 1.0;
507 }
508 return e;
509 }
510 }
511
512 #endif /* ! IEEE_FLOATING_POINT */