Fix the behaviour of scaleFloat; part of #3898
[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 isDoubleFinite(HsDouble d)
116 {
117 union stg_ieee754_dbl u;
118
119 u.d = d;
120
121 return u.ieee.exponent != 2047;
122 }
123
124 HsInt
125 isDoubleNaN(HsDouble d)
126 {
127 union stg_ieee754_dbl u;
128
129 u.d = d;
130
131 return (
132 u.ieee.exponent == 2047 /* 2^11 - 1 */ && /* Is the exponent all ones? */
133 (u.ieee.mantissa0 != 0 || u.ieee.mantissa1 != 0)
134 /* and the mantissa non-zero? */
135 );
136 }
137
138 HsInt
139 isDoubleInfinite(HsDouble d)
140 {
141 union stg_ieee754_dbl u;
142
143 u.d = d;
144
145 /* Inf iff exponent is all ones, mantissa all zeros */
146 return (
147 u.ieee.exponent == 2047 /* 2^11 - 1 */ &&
148 u.ieee.mantissa0 == 0 &&
149 u.ieee.mantissa1 == 0
150 );
151 }
152
153 HsInt
154 isDoubleDenormalized(HsDouble d)
155 {
156 union stg_ieee754_dbl u;
157
158 u.d = d;
159
160 /* A (single/double/quad) precision floating point number
161 is denormalised iff:
162 - exponent is zero
163 - mantissa is non-zero.
164 - (don't care about setting of sign bit.)
165
166 */
167 return (
168 u.ieee.exponent == 0 &&
169 (u.ieee.mantissa0 != 0 ||
170 u.ieee.mantissa1 != 0)
171 );
172
173 }
174
175 HsInt
176 isDoubleNegativeZero(HsDouble d)
177 {
178 union stg_ieee754_dbl u;
179
180 u.d = d;
181 /* sign (bit 63) set (only) => negative zero */
182
183 return (
184 u.ieee.negative == 1 &&
185 u.ieee.exponent == 0 &&
186 u.ieee.mantissa0 == 0 &&
187 u.ieee.mantissa1 == 0);
188 }
189
190 /* Same tests, this time for HsFloats. */
191
192 /*
193 To recap, here's the representation of a single precision
194 IEEE floating point number:
195
196 sign 31 sign bit (0 == positive, 1 == negative)
197 exponent 30-23 exponent (biased by 127)
198 fraction 22-0 fraction (bits to right of binary point)
199 */
200
201
202 HsInt
203 isFloatFinite(HsFloat f)
204 {
205 union stg_ieee754_flt u;
206 u.f = f;
207 return u.ieee.exponent != 255;
208 }
209
210 HsInt
211 isFloatNaN(HsFloat f)
212 {
213 union stg_ieee754_flt u;
214 u.f = f;
215
216 /* Floating point NaN iff exponent is all ones, mantissa is
217 non-zero (but see below.) */
218 return (
219 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
220 u.ieee.mantissa != 0);
221 }
222
223 HsInt
224 isFloatInfinite(HsFloat f)
225 {
226 union stg_ieee754_flt u;
227 u.f = f;
228
229 /* A float is Inf iff exponent is max (all ones),
230 and mantissa is min(all zeros.) */
231 return (
232 u.ieee.exponent == 255 /* 2^8 - 1 */ &&
233 u.ieee.mantissa == 0);
234 }
235
236 HsInt
237 isFloatDenormalized(HsFloat f)
238 {
239 union stg_ieee754_flt u;
240 u.f = f;
241
242 /* A (single/double/quad) precision floating point number
243 is denormalised iff:
244 - exponent is zero
245 - mantissa is non-zero.
246 - (don't care about setting of sign bit.)
247
248 */
249 return (
250 u.ieee.exponent == 0 &&
251 u.ieee.mantissa != 0);
252 }
253
254 HsInt
255 isFloatNegativeZero(HsFloat f)
256 {
257 union stg_ieee754_flt u;
258 u.f = f;
259
260 /* sign (bit 31) set (only) => negative zero */
261 return (
262 u.ieee.negative &&
263 u.ieee.exponent == 0 &&
264 u.ieee.mantissa == 0);
265 }
266
267 /*
268 There are glibc versions around with buggy rintf or rint, hence we
269 provide our own. We always round ties to even, so we can be simpler.
270 */
271
272 #define FLT_HIDDEN 0x800000
273 #define FLT_POWER2 0x1000000
274
275 HsFloat
276 rintFloat(HsFloat f)
277 {
278 union stg_ieee754_flt u;
279 u.f = f;
280 /* if real exponent > 22, it's already integral, infinite or nan */
281 if (u.ieee.exponent > 149) /* 22 + 127 */
282 {
283 return u.f;
284 }
285 if (u.ieee.exponent < 126) /* (-1) + 127, abs(f) < 0.5 */
286 {
287 /* only used for rounding to Integral a, so don't care about -0.0 */
288 return 0.0;
289 }
290 /* 0.5 <= abs(f) < 2^23 */
291 unsigned int half, mask, mant, frac;
292 half = 1 << (149 - u.ieee.exponent); /* bit for 0.5 */
293 mask = 2*half - 1; /* fraction bits */
294 mant = u.ieee.mantissa | FLT_HIDDEN; /* add hidden bit */
295 frac = mant & mask; /* get fraction */
296 mant ^= frac; /* truncate mantissa */
297 if ((frac < half) || ((frac == half) && ((mant & (2*half)) == 0)))
298 {
299 /* this means we have to truncate */
300 if (mant == 0)
301 {
302 /* f == ±0.5, return 0.0 */
303 return 0.0;
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 else
313 {
314 /* round away from zero, increment mantissa */
315 mant += 2*half;
316 if (mant == FLT_POWER2)
317 {
318 /* next power of 2, increase exponent an set mantissa to 0 */
319 u.ieee.mantissa = 0;
320 u.ieee.exponent += 1;
321 return u.f;
322 }
323 else
324 {
325 /* remove hidden bit and set mantissa */
326 u.ieee.mantissa = mant ^ FLT_HIDDEN;
327 return u.f;
328 }
329 }
330 }
331
332 #define DBL_HIDDEN 0x100000
333 #define DBL_POWER2 0x200000
334 #define LTOP_BIT 0x80000000
335
336 HsDouble
337 rintDouble(HsDouble d)
338 {
339 union stg_ieee754_dbl u;
340 u.d = d;
341 /* if real exponent > 51, it's already integral, infinite or nan */
342 if (u.ieee.exponent > 1074) /* 51 + 1023 */
343 {
344 return u.d;
345 }
346 if (u.ieee.exponent < 1022) /* (-1) + 1023, abs(d) < 0.5 */
347 {
348 /* only used for rounding to Integral a, so don't care about -0.0 */
349 return 0.0;
350 }
351 unsigned int half, mask, mant, frac;
352 if (u.ieee.exponent < 1043) /* 20 + 1023, real exponent < 20 */
353 {
354 /* the fractional part meets the higher part of the mantissa */
355 half = 1 << (1042 - u.ieee.exponent); /* bit for 0.5 */
356 mask = 2*half - 1; /* fraction bits */
357 mant = u.ieee.mantissa0 | DBL_HIDDEN; /* add hidden bit */
358 frac = mant & mask; /* get fraction */
359 mant ^= frac; /* truncate mantissa */
360 if ((frac < half) ||
361 ((frac == half) && (u.ieee.mantissa1 == 0) /* a tie */
362 && ((mant & (2*half)) == 0)))
363 {
364 /* truncate */
365 if (mant == 0)
366 {
367 /* d = ±0.5, return 0.0 */
368 return 0.0;
369 }
370 /* remove hidden bit and set mantissa */
371 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
372 u.ieee.mantissa1 = 0;
373 return u.d;
374 }
375 else /* round away from zero */
376 {
377 /* zero low mantissa bits */
378 u.ieee.mantissa1 = 0;
379 /* increment integer part of mantissa */
380 mant += 2*half;
381 if (mant == DBL_POWER2)
382 {
383 /* power of 2, increment exponent and zero mantissa */
384 u.ieee.mantissa0 = 0;
385 u.ieee.exponent += 1;
386 return u.d;
387 }
388 /* remove hidden bit */
389 u.ieee.mantissa0 = mant ^ DBL_HIDDEN;
390 return u.d;
391 }
392 }
393 else
394 {
395 /* 20 <= real exponent < 52, fractional part entirely in mantissa1 */
396 half = 1 << (1074 - u.ieee.exponent); /* bit for 0.5 */
397 mask = 2*half - 1; /* fraction bits */
398 mant = u.ieee.mantissa1; /* no hidden bit here */
399 frac = mant & mask; /* get fraction */
400 mant ^= frac; /* truncate mantissa */
401 if ((frac < half) ||
402 ((frac == half) && /* tie */
403 (((half == LTOP_BIT) ? (u.ieee.mantissa0 & 1) /* yuck */
404 : (mant & (2*half)))
405 == 0)))
406 {
407 /* truncate */
408 u.ieee.mantissa1 = mant;
409 return u.d;
410 }
411 else
412 {
413 /* round away from zero */
414 /* increment mantissa */
415 mant += 2*half;
416 u.ieee.mantissa1 = mant;
417 if (mant == 0)
418 {
419 /* low part of mantissa overflowed */
420 /* increment high part of mantissa */
421 mant = u.ieee.mantissa0 + 1;
422 if (mant == DBL_HIDDEN)
423 {
424 /* hit power of 2 */
425 /* zero mantissa */
426 u.ieee.mantissa0 = 0;
427 /* and increment exponent */
428 u.ieee.exponent += 1;
429 return u.d;
430 }
431 else
432 {
433 u.ieee.mantissa0 = mant;
434 return u.d;
435 }
436 }
437 else
438 {
439 return u.d;
440 }
441 }
442 }
443 }
444
445 #else /* ! IEEE_FLOATING_POINT */
446
447 /* Dummy definitions of predicates - they all return "normal" values */
448 HsInt isDoubleFinite(d) HsDouble d; { return 1;}
449 HsInt isDoubleNaN(d) HsDouble d; { return 0; }
450 HsInt isDoubleInfinite(d) HsDouble d; { return 0; }
451 HsInt isDoubleDenormalized(d) HsDouble d; { return 0; }
452 HsInt isDoubleNegativeZero(d) HsDouble d; { return 0; }
453 HsInt isFloatFinite(f) HsFloat f; { return 1; }
454 HsInt isFloatNaN(f) HsFloat f; { return 0; }
455 HsInt isFloatInfinite(f) HsFloat f; { return 0; }
456 HsInt isFloatDenormalized(f) HsFloat f; { return 0; }
457 HsInt isFloatNegativeZero(f) HsFloat f; { return 0; }
458
459
460 /* For exotic floating point formats, we can't do much */
461 /* We suppose the format has not too many bits */
462 /* I hope nobody tries to build GHC where this is wrong */
463
464 #define FLT_UPP 536870912.0
465
466 HsFloat
467 rintFloat(HsFloat f)
468 {
469 if ((f > FLT_UPP) || (f < (-FLT_UPP)))
470 {
471 return f;
472 }
473 else
474 {
475 int i = (int)f;
476 float g = i;
477 float d = f - g;
478 if (d > 0.5)
479 {
480 return g + 1.0;
481 }
482 if (d == 0.5)
483 {
484 return (i & 1) ? (g + 1.0) : g;
485 }
486 if (d == -0.5)
487 {
488 return (i & 1) ? (g - 1.0) : g;
489 }
490 if (d < -0.5)
491 {
492 return g - 1.0;
493 }
494 return g;
495 }
496 }
497
498 #define DBL_UPP 2305843009213693952.0
499
500 HsDouble
501 rintDouble(HsDouble d)
502 {
503 if ((d > DBL_UPP) || (d < (-DBL_UPP)))
504 {
505 return d;
506 }
507 else
508 {
509 HsInt64 i = (HsInt64)d;
510 double e = i;
511 double r = d - e;
512 if (r > 0.5)
513 {
514 return e + 1.0;
515 }
516 if (r == 0.5)
517 {
518 return (i & 1) ? (e + 1.0) : e;
519 }
520 if (r == -0.5)
521 {
522 return (i & 1) ? (e - 1.0) : e;
523 }
524 if (r < -0.5)
525 {
526 return e - 1.0;
527 }
528 return e;
529 }
530 }
531
532 #endif /* ! IEEE_FLOATING_POINT */