2 * Copyright (C) 1995-2007 University of Karlsruhe. All right reserved.
4 * This file is part of libFirm.
6 * This file may be distributed and/or modified under the terms of the
7 * GNU General Public License version 2 as published by the Free Software
8 * Foundation and appearing in the file LICENSE.GPL included in the
9 * packaging of this file.
11 * Licensees holding valid libFirm Professional Edition licenses may use
12 * this file in accordance with the libFirm Commercial License.
13 * Agreement provided with the Software.
15 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
16 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * @brief tarval floating point calculations
24 * @author Mathias Heil
35 #include <math.h> /* need isnan() and isinf() (will be changed)*/
36 /* undef some reused constants defined by math.h */
41 #ifdef HAVE_INTTYPES_H
42 # include <inttypes.h>
55 typedef uint32_t UINT32;
57 #ifdef HAVE_LONG_DOUBLE
58 #ifdef WORDS_BIGENDIAN
65 volatile long double d;
74 volatile long double d;
78 #ifdef WORDS_BIGENDIAN
98 * possible float states
101 NORMAL, /**< normal representation, implicit 1 */
103 SUBNORMAL, /**< denormals, implicit 0 */
105 NAN, /**< Not A Number */
108 /** A descriptor for an IEEE float value. */
110 unsigned char exponent_size; /**< size of exponent in bits */
111 unsigned char mantissa_size; /**< size of mantissa in bits */
112 value_class_t clss; /**< state of this float */
115 #define CLEAR_BUFFER(buffer) memset(buffer, 0, calc_buffer_size)
117 /* our floating point value */
121 char value[1]; /* exp[value_size] + mant[value_size] */
124 #define _exp(a) &((a)->value[0])
125 #define _mant(a) &((a)->value[value_size])
127 #define _save_result(x) memcpy((x), sc_get_buffer(), value_size)
128 #define _shift_right(x, y, b) sc_shr((x), (y), value_size*4, 0, (b))
129 #define _shift_left(x, y, b) sc_shl((x), (y), value_size*4, 0, (b))
132 #define FUNC_PTR(code) fc_##code
135 # define DEBUGPRINTF(x) printf x
137 # define DEBUGPRINTF(x) ((void)0)
140 #if FLTCALC_TRACE_CALC
141 # define TRACEPRINTF(x) printf x
143 # define TRACEPRINTF(x) ((void)0)
146 static fp_value *calc_buffer = NULL;
148 static fc_rounding_mode_t rounding_mode;
150 static int calc_buffer_size;
151 static int value_size;
152 static int max_precision;
157 static void _fail_char(const char *str, unsigned int len, int pos) {
159 printf("ERROR: Unexpected character '%c'\n", *(str + pos));
161 printf("ERROR: Unexpected end of string\n");
162 while (len-- && *str) printf("%c", *str++); printf("\n");
163 while (pos--) printf(" "); printf("^\n");
164 /* the front end has to to check constant strings */
169 /** pack machine-like */
170 static void *_pack(const fp_value *int_float, void *packed) {
173 fp_value *val_buffer;
175 temp = alloca(value_size);
176 shift_val = alloca(value_size);
178 switch (int_float->desc.clss) {
180 val_buffer = alloca(calc_buffer_size);
181 fc_get_qnan(int_float->desc.exponent_size, int_float->desc.mantissa_size, val_buffer);
182 int_float = val_buffer;
186 val_buffer = alloca(calc_buffer_size);
187 fc_get_plusinf(int_float->desc.exponent_size, int_float->desc.mantissa_size, val_buffer);
188 val_buffer->sign = int_float->sign;
189 int_float = val_buffer;
196 sc_val_from_ulong(int_float->sign, temp);
198 sc_val_from_ulong(int_float->desc.exponent_size + int_float->desc.mantissa_size, NULL);
199 _shift_left(temp, sc_get_buffer(), packed);
201 /* extract exponent */
202 sc_val_from_ulong(int_float->desc.mantissa_size, shift_val);
204 _shift_left(_exp(int_float), shift_val, temp);
206 sc_or(temp, packed, packed);
208 /* extract mantissa */
209 /* remove 2 rounding bits */
210 sc_val_from_ulong(2, shift_val);
211 _shift_right(_mant(int_float), shift_val, temp);
213 /* remove leading 1 (or 0 if denormalized) */
214 sc_max_from_bits(int_float->desc.mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
215 sc_and(temp, shift_val, temp);
218 sc_or(temp, packed, packed);
223 static void _normalize(const fp_value *in_val, fp_value *out_val, int sticky) {
225 char lsb, guard, round, round_dir = 0;
226 char *temp = alloca(value_size);
228 /* +2: save two rounding bits at the end */
229 hsb = 2 + in_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
231 if (in_val != out_val) {
232 out_val->sign = in_val->sign;
233 memcpy(&out_val->desc, &in_val->desc, sizeof(out_val->desc));
236 out_val->desc.clss = NORMAL;
238 /* mantissa all zeros, so zero exponent (because of explicit one)*/
239 if (hsb == 2 + in_val->desc.mantissa_size) {
240 sc_val_from_ulong(0, _exp(out_val));
244 /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
247 sc_val_from_ulong(-hsb-1, temp);
249 _shift_right(_mant(in_val), temp, _mant(out_val));
251 /* remember if some bits were shifted away */
252 if (!sticky) sticky = sc_had_carry();
254 sc_add(_exp(in_val), temp, _exp(out_val));
255 } else if (hsb > -1) {
257 sc_val_from_ulong(hsb+1, temp);
259 _shift_left(_mant(in_val), temp, _mant(out_val));
261 sc_sub(_exp(in_val), temp, _exp(out_val));
264 /* check for exponent underflow */
265 if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
266 DEBUGPRINTF(("Exponent underflow!\n"));
267 /* exponent underflow */
268 /* shift the mantissa right to have a zero exponent */
269 sc_val_from_ulong(1, temp);
270 sc_sub(temp, _exp(out_val), NULL);
272 _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
273 if (!sticky) sticky = sc_had_carry();
274 /* denormalized means exponent of zero */
275 sc_val_from_ulong(0, _exp(out_val));
277 out_val->desc.clss = SUBNORMAL;
280 /* perform rounding by adding a value that clears the guard bit and the round bit
281 * and either causes a carry to round up or not */
282 /* get the last 3 bits of the value */
283 lsb = sc_sub_bits(_mant(out_val), out_val->desc.mantissa_size + 2, 0) & 0x7;
284 guard = (lsb&0x2)>>1;
287 switch (rounding_mode) {
289 /* round to nearest representable value, if in doubt choose the version
291 round_dir = guard && (sticky || round || lsb>>2);
294 /* if positive: round to one if the exact value is bigger, else to zero */
295 round_dir = (!out_val->sign && (guard || round || sticky));
298 /* if negative: round to one if the exact value is bigger, else to zero */
299 round_dir = (out_val->sign && (guard || round || sticky));
302 /* always round to 0 (chopping mode) */
306 DEBUGPRINTF(("Rounding (s%d, l%d, g%d, r%d, s%d) %s\n", out_val->sign, lsb>>2, guard, round, sticky, (round_dir)?"up":"down"));
308 if (round_dir == 1) {
309 guard = (round^guard)<<1;
310 lsb = !(round || guard)<<2 | guard | round;
312 lsb = -((guard<<1) | round);
315 /* add the rounded value */
317 sc_val_from_long(lsb, temp);
318 sc_add(_mant(out_val), temp, _mant(out_val));
321 /* could have rounded down to zero */
322 if (sc_is_zero(_mant(out_val)) && (out_val->desc.clss == SUBNORMAL))
323 out_val->desc.clss = ZERO;
325 /* check for rounding overflow */
326 hsb = 2 + out_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
327 if ((out_val->desc.clss != SUBNORMAL) && (hsb < -1)) {
328 sc_val_from_ulong(1, temp);
329 _shift_right(_mant(out_val), temp, _mant(out_val));
331 sc_add(_exp(out_val), temp, _exp(out_val));
332 } else if ((out_val->desc.clss == SUBNORMAL) && (hsb == -1)) {
333 /* overflow caused the mantissa to be normal again,
334 * so adapt the exponent accordingly */
335 sc_val_from_ulong(1, temp);
336 sc_add(_exp(out_val), temp, _exp(out_val));
338 out_val->desc.clss = NORMAL;
340 /* no further rounding is needed, because rounding overflow means
341 * the carry of the original rounding was propagated all the way
342 * up to the bit left of the radix point. This implies the bits
343 * to the right are all zeros (rounding is +1) */
345 /* check for exponent overflow */
346 sc_val_from_ulong((1 << out_val->desc.exponent_size) - 1, temp);
347 if (sc_comp(_exp(out_val), temp) != -1) {
348 DEBUGPRINTF(("Exponent overflow!\n"));
349 /* exponent overflow, reaction depends on rounding method:
351 * mode | sign of value | result
352 *--------------------------------------------------------------
353 * TO_NEAREST | + | +inf
355 *--------------------------------------------------------------
356 * TO_POSITIVE | + | +inf
357 * | - | smallest representable value
358 *--------------------------------------------------------------
359 * TO_NEAGTIVE | + | largest representable value
361 *--------------------------------------------------------------
362 * TO_ZERO | + | largest representable value
363 * | - | smallest representable value
364 *--------------------------------------------------------------*/
365 if (out_val->sign == 0) {
366 /* value is positive */
367 switch (rounding_mode) {
370 out_val->desc.clss = INF;
375 fc_get_max(out_val->desc.exponent_size, out_val->desc.mantissa_size, out_val);
378 /* value is negative */
379 switch (rounding_mode) {
382 out_val->desc.clss = INF;
387 fc_get_min(out_val->desc.exponent_size, out_val->desc.mantissa_size, out_val);
394 * Operations involving NaN's must return NaN
396 #define handle_NAN(a, b, result) \
398 if (a->desc.clss == NAN) { \
399 if (a != result) memcpy(result, a, calc_buffer_size); \
402 if (b->desc.clss == NAN) { \
403 if (b != result) memcpy(result, b, calc_buffer_size); \
410 * calculate a + b, where a is the value with the bigger exponent
412 static void _fadd(const fp_value *a, const fp_value *b, fp_value *result) {
419 handle_NAN(a, b, result);
421 /* make sure result has a descriptor */
422 if (result != a && result != b)
423 memcpy(&result->desc, &a->desc, sizeof(descriptor_t));
425 /* determine if this is an addition or subtraction */
426 sign = a->sign ^ b->sign;
428 /* produce NaN on inf - inf */
429 if (sign && (a->desc.clss == INF) && (b->desc.clss == INF)) {
430 fc_get_qnan(a->desc.exponent_size, b->desc.mantissa_size, result);
434 temp = alloca(value_size);
435 exp_diff = alloca(value_size);
437 /* get exponent difference */
438 sc_sub(_exp(a), _exp(b), exp_diff);
440 /* initially set sign to be the sign of a, special treatment of subtraction
441 * when exponents are equal is required though.
442 * Also special care about the sign is needed when the mantissas are equal
444 if (sign && sc_val_to_long(exp_diff) == 0) {
445 switch (sc_comp(_mant(a), _mant(b))) {
447 res_sign = a->sign; /* abs(a) is bigger and a is negative */
450 res_sign = (rounding_mode == FC_TONEGATIVE);
453 res_sign = b->sign; /* abs(b) is bigger and b is negative */
456 /* can't be reached */
463 result->sign = res_sign;
465 /* sign has been taken care of, check for special cases */
466 if (a->desc.clss == ZERO || b->desc.clss == INF) {
468 memcpy(result, b, calc_buffer_size);
469 result->sign = res_sign;
472 if (b->desc.clss == ZERO || a->desc.clss == INF) {
474 memcpy(result, a, calc_buffer_size);
475 result->sign = res_sign;
479 /* shift the smaller value to the right to align the radix point */
480 /* subnormals have their radix point shifted to the right,
481 * take care of this first */
482 if ((b->desc.clss == SUBNORMAL) && (a->desc.clss != SUBNORMAL)) {
483 sc_val_from_ulong(1, temp);
484 sc_sub(exp_diff, temp, exp_diff);
487 _shift_right(_mant(b), exp_diff, temp);
488 sticky = sc_had_carry();
490 if (sticky && sign) {
491 /* if subtracting a little more than the represented value or adding a little
492 * more than the represented value to a negative value this, in addition to the
493 * still set sticky bit, takes account of the 'little more' */
494 char *temp1 = alloca(calc_buffer_size);
495 sc_val_from_ulong(1, temp1);
496 sc_add(temp, temp1, temp);
500 if (sc_comp(_mant(a), temp) == -1)
501 sc_sub(temp, _mant(a), _mant(result));
503 sc_sub(_mant(a), temp, _mant(result));
505 sc_add(_mant(a), temp, _mant(result));
508 /* _normalize expects a 'normal' radix point, adding two subnormals
509 * results in a subnormal radix point -> shifting before normalizing */
510 if ((a->desc.clss == SUBNORMAL) && (b->desc.clss == SUBNORMAL)) {
511 sc_val_from_ulong(1, NULL);
512 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
515 /* resulting exponent is the bigger one */
516 memmove(_exp(result), _exp(a), value_size);
518 _normalize(result, result, sticky);
524 static void _fmul(const fp_value *a, const fp_value *b, fp_value *result) {
528 handle_NAN(a, b, result);
530 temp = alloca(value_size);
532 if (result != a && result != b)
533 memcpy(&result->desc, &a->desc, sizeof(descriptor_t));
535 result->sign = res_sign = a->sign ^ b->sign;
537 /* produce NaN on 0 * inf */
538 if (a->desc.clss == ZERO) {
539 if (b->desc.clss == INF)
540 fc_get_qnan(a->desc.exponent_size, a->desc.mantissa_size, result);
543 memcpy(result, a, calc_buffer_size);
544 result->sign = res_sign;
548 if (b->desc.clss == ZERO) {
549 if (a->desc.clss == INF)
550 fc_get_qnan(a->desc.exponent_size, a->desc.mantissa_size, result);
553 memcpy(result, b, calc_buffer_size);
554 result->sign = res_sign;
559 if (a->desc.clss == INF) {
561 memcpy(result, a, calc_buffer_size);
562 result->sign = res_sign;
565 if (b->desc.clss == INF) {
567 memcpy(result, b, calc_buffer_size);
568 result->sign = res_sign;
572 /* exp = exp(a) + exp(b) - excess */
573 sc_add(_exp(a), _exp(b), _exp(result));
575 sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 1, temp);
576 sc_sub(_exp(result), temp, _exp(result));
578 /* mixed normal, subnormal values introduce an error of 1, correct it */
579 if ((a->desc.clss == SUBNORMAL) ^ (b->desc.clss == SUBNORMAL)) {
580 sc_val_from_ulong(1, temp);
581 sc_add(_exp(result), temp, _exp(result));
584 sc_mul(_mant(a), _mant(b), _mant(result));
586 /* realign result: after a multiplication the digits right of the radix
587 * point are the sum of the factors' digits after the radix point. As all
588 * values are normalized they both have the same amount of these digits,
589 * which has to be restored by proper shifting
590 * +2 because of the two rounding bits */
591 sc_val_from_ulong(2 + result->desc.mantissa_size, temp);
593 _shift_right(_mant(result), temp, _mant(result));
595 _normalize(result, result, sc_had_carry());
601 static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result) {
602 char *temp, *dividend;
605 handle_NAN(a, b, result);
607 temp = alloca(value_size);
608 dividend = alloca(value_size);
610 if (result != a && result != b)
611 memcpy(&result->desc, &a->desc, sizeof(descriptor_t));
613 result->sign = res_sign = a->sign ^ b->sign;
615 /* produce nan on 0/0 and inf/inf */
616 if (a->desc.clss == ZERO) {
617 if (b->desc.clss == ZERO)
619 fc_get_qnan(a->desc.exponent_size, a->desc.mantissa_size, result);
623 memcpy(result, a, calc_buffer_size);
624 result->sign = res_sign;
629 if (b->desc.clss == INF) {
630 if (a->desc.clss == INF)
632 fc_get_qnan(a->desc.exponent_size, a->desc.mantissa_size, result);
635 sc_val_from_ulong(0, NULL);
636 _save_result(_exp(result));
637 _save_result(_mant(result));
638 result->desc.clss = ZERO;
643 if (a->desc.clss == INF) {
646 memcpy(result, a, calc_buffer_size);
647 result->sign = res_sign;
650 if (b->desc.clss == ZERO) {
651 /* division by zero */
653 fc_get_minusinf(a->desc.exponent_size, a->desc.mantissa_size, result);
655 fc_get_plusinf(a->desc.exponent_size, a->desc.mantissa_size, result);
659 /* exp = exp(a) - exp(b) + excess - 1*/
660 sc_sub(_exp(a), _exp(b), _exp(result));
661 sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 2, temp);
662 sc_add(_exp(result), temp, _exp(result));
664 /* mixed normal, subnormal values introduce an error of 1, correct it */
665 if ((a->desc.clss == SUBNORMAL) ^ (b->desc.clss == SUBNORMAL)) {
666 sc_val_from_ulong(1, temp);
667 sc_add(_exp(result), temp, _exp(result));
670 /* mant(res) = mant(a) / 1/2mant(b) */
671 /* to gain more bits of precision in the result the dividend could be
672 * shifted left, as this operation does not loose bits. This would not
673 * fit into the integer precision, but due to the rounding bits (which
674 * are always zero because the values are all normalized) the divisor
675 * can be shifted right instead to achieve the same result */
676 sc_val_from_ulong(2 + result->desc.mantissa_size, temp);
678 _shift_left(_mant(a), temp, dividend);
681 char *divisor = alloca(calc_buffer_size);
682 sc_val_from_ulong(1, divisor);
683 _shift_right(_mant(b), divisor, divisor);
684 sc_div(dividend, divisor, _mant(result));
687 _normalize(result, result, sc_had_carry());
691 static void _power_of_ten(int exp, descriptor_t *desc, char *result) {
698 /* set new descriptor (else result is supposed to already have one) */
700 memcpy(&result->desc, desc, sizeof(descriptor_t));
702 build = alloca(value_size);
703 temp = alloca(value_size);
705 sc_val_from_ulong((1 << result->desc.exponent_size)/2-1, _exp(result));
708 /* temp is value of ten now */
709 sc_val_from_ulong(10, NULL);
712 for (exp--; exp > 0; exp--) {
714 sc_mul(build, temp, NULL);
718 /* temp is amount of left shift needed to put the value left of the radix point */
719 sc_val_from_ulong(result->desc.mantissa_size + 2, temp);
721 _shift_left(build, temp, _mant(result));
723 _normalize(result, result, 0);
729 * Truncate the fractional part away.
731 * This does not clip to any integer rang.
733 static void _trunc(const fp_value *a, fp_value *result) {
735 * When exponent == 0 all bits left of the radix point
736 * are the integral part of the value. For 15bit exp_size
737 * this would require a left shift of max. 16383 bits which
739 * But it is enough to ensure that no bit right of the radix
740 * point remains set. This restricts the interesting
741 * exponents to the interval [0, mant_size-1].
742 * Outside this interval the truncated value is either 0 or
743 * it does not have fractional parts.
746 int exp_bias, exp_val;
749 temp = alloca(value_size);
752 memcpy(&result->desc, &a->desc, sizeof(descriptor_t));
754 exp_bias = (1 << (a->desc.exponent_size - 1)) - 1;
755 exp_val = sc_val_to_long(_exp(a)) - exp_bias;
758 sc_val_from_ulong(0, NULL);
759 _save_result(_exp(result));
760 _save_result(_mant(result));
761 result->desc.clss = ZERO;
766 if (exp_val > a->desc.mantissa_size) {
768 memcpy(result, a, calc_buffer_size);
773 /* set up a proper mask to delete all bits right of the
774 * radix point if the mantissa had been shifted until exp == 0 */
775 sc_max_from_bits(1 + exp_val, 0, temp);
776 sc_val_from_long(a->desc.mantissa_size - exp_val + 2, NULL);
777 _shift_left(temp, sc_get_buffer(), temp);
779 /* and the mask and return the result */
780 sc_and(_mant(a), temp, _mant(result));
782 if (a != result) memcpy(_exp(result), _exp(a), value_size);
788 * functions defined in fltcalc.h
790 const void *fc_get_buffer(void) {
794 int fc_get_buffer_length(void) {
795 return calc_buffer_size;
798 void *fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, void *result) {
810 int exp_int, hsb, state;
815 char *mant_str, *exp_val, *power_val;
818 if (result == NULL) result = calc_buffer;
820 exp_val = alloca(value_size);
821 power_val = alloca(calc_buffer_size);
822 mant_str = alloca((len)?(len):(strlen(str)));
824 result->desc.exponent_size = exp_size;
825 result->desc.mantissa_size = mant_size;
826 result->desc.clss = NORMAL;
833 while (len == 0 || str-old_str < len) {
849 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
856 state = RIGHT_OF_DOT;
867 _fail_char(old_str, len, str - old_str);
873 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
874 mant_str[pos++] = *(str++);
878 state = RIGHT_OF_DOT;
889 mant_str[pos] = '\0';
893 _fail_char(old_str, len, str - old_str);
899 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
900 mant_str[pos++] = *(str++);
911 mant_str[pos] = '\0';
915 _fail_char(old_str, len, str - old_str);
925 if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
929 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
930 mant_str[pos] = '\0';
937 _fail_char(old_str, len, str - old_str);
943 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
948 case '\0': goto done;
951 _fail_char(old_str, len, str - old_str);
954 } /* switch(state) */
957 sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
959 /* shift to put value left of radix point */
960 sc_val_from_ulong(mant_size + 2, exp_val);
962 _shift_left(_mant(result), exp_val, _mant(result));
964 sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
966 _normalize(result, result, 0);
968 if (state == EXPONENT) {
969 exp_int -= atoi(str-pos);
972 _power_of_ten(exp_int, &result->desc, power_val);
974 _fdiv(result, power_val, result);
978 /* XXX excuse of an implementation to make things work */
980 fp_value *tmp = alloca(calc_buffer_size);
983 #ifdef HAVE_LONG_DOUBLE
984 val = strtold(str, NULL);
985 DEBUGPRINTF(("val_from_str(%s)\n", str));
986 fc_val_from_ieee754(val, 15, 64, tmp);
988 val = strtod(str, NULL);
989 DEBUGPRINTF(("val_from_str(%s)\n", str));
990 fc_val_from_ieee754(val, 11, 52, tmp);
991 #endif /* HAVE_LONG_DOUBLE */
992 return fc_cast(tmp, exp_size, mant_size, result);
996 fp_value *fc_val_from_ieee754(LLDBL l, char exp_size, char mant_size, fp_value *result) {
998 int bias_res, bias_val, mant_val;
1000 UINT32 sign, exponent, mantissa0, mantissa1;
1003 bias_res = ((1<<exp_size)/2-1);
1005 #ifdef HAVE_LONG_DOUBLE
1008 sign = (srcval.val.high & 0x00008000) != 0;
1009 exponent = (srcval.val.high & 0x00007FFF) ;
1010 mantissa0 = srcval.val.mid;
1011 mantissa1 = srcval.val.low;
1012 #else /* no long double */
1015 sign = (srcval.val.high & 0x80000000) != 0;
1016 exponent = (srcval.val.high & 0x7FF00000) >> 20;
1017 mantissa0 = srcval.val.high & 0x000FFFFF;
1018 mantissa1 = srcval.val.low;
1021 #ifdef HAVE_LONG_DOUBLE
1022 TRACEPRINTF(("val_from_float(%.8X%.8X%.8X)\n", ((int*)&l)[2], ((int*)&l)[1], ((int*)&l)[0]));/* srcval.val.high, srcval.val.mid, srcval.val.low)); */
1023 DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1025 TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1026 DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1029 if (result == NULL) result = calc_buffer;
1030 temp = alloca(value_size);
1032 result->desc.exponent_size = exp_size;
1033 result->desc.mantissa_size = mant_size;
1036 result->sign = sign;
1038 /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1039 * encoding is needed. the function can return immediately in these cases */
1041 result->desc.clss = NAN;
1042 TRACEPRINTF(("val_from_float resulted in NAN\n"));
1045 else if (isinf(l)) {
1046 result->desc.clss = INF;
1047 TRACEPRINTF(("val_from_float resulted in %sINF\n", (result->sign == 1) ? "-" : ""));
1051 /* build exponent, because input and output exponent and mantissa sizes may differ
1052 * this looks more complicated than it is: unbiased input exponent + output bias,
1053 * minus the mantissa difference which is added again later when the output float
1054 * becomes normalized */
1055 #ifdef HAVE_EXPLICIT_ONE
1056 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1058 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1061 /* build mantissa representation */
1062 #ifndef HAVE_EXPLICIT_ONE
1063 if (exponent != 0) {
1064 /* insert the hidden bit */
1065 sc_val_from_ulong(1, temp);
1066 sc_val_from_ulong(mant_val + 2, NULL);
1067 _shift_left(temp, sc_get_buffer(), NULL);
1072 sc_val_from_ulong(0, NULL);
1075 _save_result(_mant(result));
1077 /* bits from the upper word */
1078 sc_val_from_ulong(mantissa0, temp);
1079 sc_val_from_ulong(34, NULL);
1080 _shift_left(temp, sc_get_buffer(), temp);
1081 sc_or(_mant(result), temp, _mant(result));
1083 /* bits from the lower word */
1084 sc_val_from_ulong(mantissa1, temp);
1085 sc_val_from_ulong(2, NULL);
1086 _shift_left(temp, sc_get_buffer(), temp);
1087 sc_or(_mant(result), temp, _mant(result));
1089 /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1090 * origin one to the left */
1091 if (exponent == 0) {
1092 sc_val_from_ulong(1, NULL);
1093 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1096 _normalize(result, result, 0);
1098 TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1103 LLDBL fc_val_to_ieee754(const fp_value *val) {
1105 fp_value *temp = NULL;
1116 #ifdef HAVE_LONG_DOUBLE
1117 char result_exponent = 15;
1118 char result_mantissa = 64;
1120 char result_exponent = 11;
1121 char result_mantissa = 52;
1124 temp = alloca(calc_buffer_size);
1125 #ifdef HAVE_EXPLICIT_ONE
1126 value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1128 value = fc_cast(val, result_exponent, result_mantissa, temp);
1133 /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1134 * lead to wrong results */
1135 exponent = sc_val_to_long(_exp(value)) ;
1137 sc_val_from_ulong(2, NULL);
1138 _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1143 for (byte_offset = 0; byte_offset < 4; byte_offset++)
1144 mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1146 for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1147 mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1149 #ifdef HAVE_LONG_DOUBLE
1150 buildval.val.high = sign << 15;
1151 buildval.val.high |= exponent;
1152 buildval.val.mid = mantissa0;
1153 buildval.val.low = mantissa1;
1154 #else /* no long double */
1155 mantissa0 &= 0x000FFFFF; /* get rid of garbage */
1156 buildval.val.high = sign << 31;
1157 buildval.val.high |= exponent << 20;
1158 buildval.val.high |= mantissa0;
1159 buildval.val.low = mantissa1;
1162 TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1166 fp_value *fc_cast(const fp_value *value, char exp_size, char mant_size, fp_value *result) {
1168 int exp_offset, val_bias, res_bias;
1170 if (result == NULL) result = calc_buffer;
1171 temp = alloca(value_size);
1173 if (value->desc.exponent_size == exp_size && value->desc.mantissa_size == mant_size) {
1174 if (value != result)
1175 memcpy(result, value, calc_buffer_size);
1179 /* set the descriptor of the new value */
1180 result->desc.exponent_size = exp_size;
1181 result->desc.mantissa_size = mant_size;
1182 result->desc.clss = value->desc.clss;
1184 result->sign = value->sign;
1186 /* when the mantissa sizes differ normalizing has to shift to align it.
1187 * this would change the exponent, which is unwanted. So calculate this
1188 * offset and add it */
1189 val_bias = (1 << (value->desc.exponent_size - 1)) - 1;
1190 res_bias = (1 << (exp_size - 1)) - 1;
1192 exp_offset = (res_bias - val_bias) - (value->desc.mantissa_size - mant_size);
1193 sc_val_from_long(exp_offset, temp);
1194 sc_add(_exp(value), temp, _exp(result));
1196 /* _normalize expects normalized radix point */
1197 if (value->desc.clss == SUBNORMAL) {
1198 sc_val_from_ulong(1, NULL);
1199 _shift_left(_mant(value), sc_get_buffer(), _mant(result));
1200 } else if (value != result) {
1201 memcpy(_mant(result), _mant(value), value_size);
1203 memmove(_mant(result), _mant(value), value_size);
1206 _normalize(result, result, 0);
1207 TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, value_size, FC_PACKED)));
1211 fp_value *fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1212 if (result == NULL) result = calc_buffer;
1214 result->desc.exponent_size = exponent_size;
1215 result->desc.mantissa_size = mantissa_size;
1216 result->desc.clss = NORMAL;
1220 sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1222 sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1223 sc_val_from_ulong(2, NULL);
1224 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1229 fp_value *fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1230 if (result == NULL) result = calc_buffer;
1232 fc_get_max(exponent_size, mantissa_size, result);
1238 fp_value *fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1239 if (result == NULL) result = calc_buffer;
1241 result->desc.exponent_size = exponent_size;
1242 result->desc.mantissa_size = mantissa_size;
1243 result->desc.clss = NAN;
1247 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1249 /* signaling NaN has non-zero mantissa with msb not set */
1250 sc_val_from_ulong(1, _mant(result));
1255 fp_value *fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1256 if (result == NULL) result = calc_buffer;
1258 result->desc.exponent_size = exponent_size;
1259 result->desc.mantissa_size = mantissa_size;
1260 result->desc.clss = NAN;
1264 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1266 /* quiet NaN has the msb of the mantissa set, so shift one there */
1267 sc_val_from_ulong(1, _mant(result));
1268 /* mantissa_size >+< 1 because of two extra rounding bits */
1269 sc_val_from_ulong(mantissa_size + 1, NULL);
1270 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1275 fp_value *fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1276 if (result == NULL) result = calc_buffer;
1278 result->desc.exponent_size = exponent_size;
1279 result->desc.mantissa_size = mantissa_size;
1280 result->desc.clss = NORMAL;
1284 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1286 sc_val_from_ulong(0, _mant(result));
1291 fp_value *fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, fp_value *result) {
1292 if (result == NULL) result = calc_buffer;
1294 fc_get_plusinf(exponent_size, mantissa_size, result);
1300 int fc_comp(const fp_value *val_a, const fp_value *val_b) {
1304 * shortcut: if both values are identical, they are either
1305 * Unordered if NaN or equal
1308 return val_a->desc.clss == NAN ? 2 : 0;
1310 /* unordered if one is a NaN */
1311 if (val_a->desc.clss == NAN || val_b->desc.clss == NAN)
1314 /* zero is equal independent of sign */
1315 if ((val_a->desc.clss == ZERO) && (val_b->desc.clss == ZERO))
1318 /* different signs make compare easy */
1319 if (val_a->sign != val_b->sign)
1320 return (val_a->sign == 0) ? (1) : (-1);
1322 mul = val_a->sign ? -1 : 1;
1324 /* both infinity means equality */
1325 if ((val_a->desc.clss == INF) && (val_b->desc.clss == INF))
1328 /* infinity is bigger than the rest */
1329 if (val_a->desc.clss == INF)
1331 if (val_b->desc.clss == INF)
1334 /* check first exponent, that mantissa if equal */
1335 switch (sc_comp(_exp(val_a), _exp(val_b))) {
1341 return sc_comp(_mant(val_a), _mant(val_b)) * mul;
1347 int fc_is_zero(const fp_value *a) {
1348 return a->desc.clss == ZERO;
1351 int fc_is_negative(const fp_value *a) {
1355 int fc_is_inf(const fp_value *a) {
1356 return a->desc.clss == INF;
1359 int fc_is_nan(const fp_value *a) {
1360 return a->desc.clss == NAN;
1363 int fc_is_subnormal(const fp_value *a) {
1364 return a->desc.clss == SUBNORMAL;
1367 char *fc_print(const fp_value *val, char *buf, int buflen, unsigned base) {
1370 mul_1 = alloca(calc_buffer_size);
1374 switch (val->desc.clss) {
1376 if (buflen >= 8 + val->sign) sprintf(buf, "%sINFINITY", val->sign ? "-":"");
1377 else snprintf(buf, buflen, "%sINF", val->sign ? "-":NULL);
1380 snprintf(buf, buflen, "NAN");
1383 snprintf(buf, buflen, "0.0");
1386 /* XXX to be implemented */
1387 #ifdef HAVE_LONG_DOUBLE
1388 /* XXX 30 is arbitrary */
1389 snprintf(buf, buflen, "%.30LE", fc_val_to_ieee754(val));
1391 snprintf(buf, buflen, "%.18E", fc_val_to_ieee754(val));
1397 switch (val->desc.clss) {
1399 if (buflen >= 8+val->sign) sprintf(buf, "%sINFINITY", val->sign?"-":"");
1400 else snprintf(buf, buflen, "%sINF", val->sign?"-":NULL);
1403 snprintf(buf, buflen, "NAN");
1406 snprintf(buf, buflen, "0.0");
1409 #ifdef HAVE_LONG_DOUBLE
1410 snprintf(buf, buflen, "%LA", fc_val_to_ieee754(val));
1412 snprintf(buf, buflen, "%A", fc_val_to_ieee754(val));
1419 snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), value_size*4, SC_HEX, 0));
1420 buf[buflen - 1] = '\0';
1426 unsigned char fc_sub_bits(const fp_value *value, unsigned num_bits, unsigned byte_ofs) {
1427 /* this is used to cache the packed version of the value */
1428 static char *pack = NULL;
1430 if (pack == NULL) pack = xmalloc(value_size);
1435 return sc_sub_bits(pack, num_bits, byte_ofs);
1438 int fc_zero_mantissa(const fp_value *value) {
1439 return sc_is_zero(_mant(value));
1442 int fc_get_exponent(const fp_value *value) {
1443 int exp_bias = (1 << (value->desc.exponent_size - 1)) - 1;
1444 return sc_val_to_long(_exp(value)) - exp_bias;
1448 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode) {
1449 if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1450 rounding_mode = mode;
1452 return rounding_mode;
1455 fc_rounding_mode_t fc_get_rounding_mode(void) {
1456 return rounding_mode;
1459 void init_fltcalc(int precision) {
1460 if (calc_buffer == NULL) {
1461 /* does nothing if already init */
1462 if (precision == 0) precision = FC_DEFAULT_PRECISION;
1464 init_strcalc(precision + 4);
1466 /* needs additionally two bits to round, a bit as explicit 1., and one for
1467 * addition overflow */
1468 max_precision = sc_get_precision() - 4;
1469 if (max_precision < precision)
1470 printf("WARING: not enough precision available, using %d\n", max_precision);
1472 rounding_mode = FC_TONEAREST;
1473 value_size = sc_get_buffer_length();
1474 calc_buffer_size = sizeof(fp_value) + 2*value_size - 1;
1476 calc_buffer = xmalloc(calc_buffer_size);
1477 memset(calc_buffer, 0, calc_buffer_size);
1478 DEBUGPRINTF(("init fltcalc:\n\tVALUE_SIZE = %d\ntCALC_BUFFER_SIZE = %d\n\tcalc_buffer = %p\n\n", value_size, calc_buffer_size, calc_buffer));
1479 #ifdef HAVE_LONG_DOUBLE
1480 DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1482 DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1484 #ifdef WORDS_BIGENDIAN
1485 DEBUGPRINTF(("\tWord order is big endian\n\n"));
1487 DEBUGPRINTF(("\tWord order is little endian\n\n"));
1492 void finish_fltcalc (void) {
1493 free(calc_buffer); calc_buffer = NULL;
1496 /* definition of interface functions */
1497 fp_value *fc_add(const fp_value *a, const fp_value *b, fp_value *result) {
1498 if (result == NULL) result = calc_buffer;
1500 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1501 TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
1503 /* make the value with the bigger exponent the first one */
1504 if (sc_comp(_exp(a), _exp(b)) == -1)
1505 _fadd(b, a, result);
1507 _fadd(a, b, result);
1509 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1513 fp_value *fc_sub(const fp_value *a, const fp_value *b, fp_value *result) {
1516 if (result == NULL) result = calc_buffer;
1518 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1519 TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
1521 temp = alloca(calc_buffer_size);
1522 memcpy(temp, b, calc_buffer_size);
1523 temp->sign = !b->sign;
1524 if (sc_comp(_exp(a), _exp(temp)) == -1)
1525 _fadd(temp, a, result);
1527 _fadd(a, temp, result);
1529 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1533 fp_value *fc_mul(const fp_value *a, const fp_value *b, fp_value *result) {
1534 if (result == NULL) result = calc_buffer;
1536 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1537 TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
1539 _fmul(a, b, result);
1541 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1545 fp_value *fc_div(const fp_value *a, const fp_value *b, fp_value *result) {
1546 if (result == NULL) result = calc_buffer;
1548 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1549 TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
1551 _fdiv(a, b, result);
1553 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1557 fp_value *fc_neg(const fp_value *a, fp_value *result) {
1558 if (result == NULL) result = calc_buffer;
1560 TRACEPRINTF(("- %s ", fc_print(a, buffer, 100, FC_PACKED)));
1563 memcpy(result, a, calc_buffer_size);
1564 result->sign = !a->sign;
1566 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1570 fp_value *fc_int(const fp_value *a, fp_value *result) {
1571 if (result == NULL) result = calc_buffer;
1573 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1574 TRACEPRINTF(("truncated to integer "));
1578 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
1582 fp_value *fc_rnd(const fp_value *a, fp_value *result) {
1583 if (result == NULL) result = calc_buffer;
1585 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
1586 TRACEPRINTF(("rounded to integer "));
1588 assert(!"fc_rnd() not yet implemented");
1590 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));