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 /* because variable sized structs are impossible, the internal
118 * value is represented as a pseudo-struct char array, addressed
121 * char sign; // 0 for positive, 1 for negative
122 * char exp[value_size];
123 * char mant[value_size];
127 #define _sign(a) (((char*)a)[SIGN_POS])
128 #define _exp(a) (&((char*)a)[EXPONENT_POS])
129 #define _mant(a) (&((char*)a)[MANTISSA_POS])
130 #define _desc(a) (*(descriptor_t *)&((char *)a)[DESCRIPTOR_POS])
132 #define _save_result(x) memcpy((x), sc_get_buffer(), value_size)
133 #define _shift_right(x, y, b) sc_shr((x), (y), value_size*4, 0, (b))
134 #define _shift_left(x, y, b) sc_shl((x), (y), value_size*4, 0, (b))
136 #define FC_DEFINE1(code) \
137 char *fc_##code(const void *a, void *result) { \
138 return _calc((const char*)a, NULL, FC_##code, (char*)result); \
141 #define FC_DEFINE2(code) \
142 char *fc_##code(const void *a, const void *b, void *result) { \
143 return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \
146 #define FUNC_PTR(code) fc_##code
149 # define DEBUGPRINTF(x) printf x
151 # define DEBUGPRINTF(x) ((void)0)
154 #if FLTCALC_TRACE_CALC
155 # define TRACEPRINTF(x) printf x
157 # define TRACEPRINTF(x) ((void)0)
160 static char *calc_buffer = NULL;
162 static fc_rounding_mode_t rounding_mode;
164 static int calc_buffer_size;
165 static int value_size;
167 static int EXPONENT_POS;
168 static int MANTISSA_POS;
169 static int DESCRIPTOR_POS;
171 static int max_precision;
176 static void _fail_char(const char *str, unsigned int len, int pos)
179 printf("ERROR: Unexpected character '%c'\n", *(str + pos));
181 printf("ERROR: Unexpected end of string\n");
182 while (len-- && *str) printf("%c", *str++); printf("\n");
183 while (pos--) printf(" "); printf("^\n");
184 /* the front end has to to check constant strings */
189 /** pack machine-like */
190 static char* _pack(const char *int_float, char *packed)
196 temp = alloca(value_size);
197 shift_val = alloca(value_size);
199 switch (_desc(int_float).clss) {
201 val_buffer = alloca(calc_buffer_size);
202 fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
203 int_float = val_buffer;
207 val_buffer = alloca(calc_buffer_size);
208 fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
209 _sign(val_buffer) = _sign(int_float);
210 int_float = val_buffer;
217 sc_val_from_ulong(_sign(int_float), temp);
219 sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
220 _shift_left(temp, sc_get_buffer(), packed);
222 /* extract exponent */
223 sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);
225 _shift_left(_exp(int_float), shift_val, temp);
227 sc_or(temp, packed, packed);
229 /* extract mantissa */
230 /* remove 2 rounding bits */
231 sc_val_from_ulong(2, shift_val);
232 _shift_right(_mant(int_float), shift_val, temp);
234 /* remove leading 1 (or 0 if denormalized) */
235 sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
236 sc_and(temp, shift_val, temp);
239 sc_or(temp, packed, packed);
244 char* _normalize(const char *in_val, char *out_val, int sticky)
247 char lsb, guard, round, round_dir = 0;
250 temp = alloca(value_size);
252 /* +2: save two rounding bits at the end */
253 hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
255 if (in_val != out_val)
257 _sign(out_val) = _sign(in_val);
258 memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
261 _desc(out_val).clss = NORMAL;
263 /* mantissa all zeros, so zero exponent (because of explicit one)*/
264 if (hsb == 2 + _desc(in_val).mantissa_size)
266 sc_val_from_ulong(0, _exp(out_val));
270 /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
274 sc_val_from_ulong(-hsb-1, temp);
276 _shift_right(_mant(in_val), temp, _mant(out_val));
278 /* remember if some bits were shifted away */
279 if (!sticky) sticky = sc_had_carry();
281 sc_add(_exp(in_val), temp, _exp(out_val));
286 sc_val_from_ulong(hsb+1, temp);
288 _shift_left(_mant(in_val), temp, _mant(out_val));
290 sc_sub(_exp(in_val), temp, _exp(out_val));
293 /* check for exponent underflow */
294 if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
295 DEBUGPRINTF(("Exponent underflow!\n"));
296 /* exponent underflow */
297 /* shift the mantissa right to have a zero exponent */
298 sc_val_from_ulong(1, temp);
299 sc_sub(temp, _exp(out_val), NULL);
301 _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
302 if (!sticky) sticky = sc_had_carry();
303 /* denormalized means exponent of zero */
304 sc_val_from_ulong(0, _exp(out_val));
306 _desc(out_val).clss = SUBNORMAL;
309 /* perform rounding by adding a value that clears the guard bit and the round bit
310 * and either causes a carry to round up or not */
311 /* get the last 3 bits of the value */
312 lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
313 guard = (lsb&0x2)>>1;
316 switch (rounding_mode)
319 /* round to nearest representable value, if in doubt choose the version
321 round_dir = guard && (sticky || round || lsb>>2);
324 /* if positive: round to one if the exact value is bigger, else to zero */
325 round_dir = (!_sign(out_val) && (guard || round || sticky));
328 /* if negative: round to one if the exact value is bigger, else to zero */
329 round_dir = (_sign(out_val) && (guard || round || sticky));
332 /* always round to 0 (chopping mode) */
336 DEBUGPRINTF(("Rounding (s%d, l%d, g%d, r%d, s%d) %s\n", _sign(out_val), lsb>>2, guard, round, sticky, (round_dir)?"up":"down"));
340 guard = (round^guard)<<1;
341 lsb = !(round || guard)<<2 | guard | round;
345 lsb = -((guard<<1) | round);
348 /* add the rounded value */
350 sc_val_from_long(lsb, temp);
351 sc_add(_mant(out_val), temp, _mant(out_val));
354 /* could have rounded down to zero */
355 if (sc_is_zero(_mant(out_val)) && (_desc(out_val).clss == SUBNORMAL))
356 _desc(out_val).clss = ZERO;
358 /* check for rounding overflow */
359 hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
360 if ((_desc(out_val).clss != SUBNORMAL) && (hsb < -1))
362 sc_val_from_ulong(1, temp);
363 _shift_right(_mant(out_val), temp, _mant(out_val));
365 sc_add(_exp(out_val), temp, _exp(out_val));
367 else if ((_desc(out_val).clss == SUBNORMAL) && (hsb == -1))
369 /* overflow caused the mantissa to be normal again,
370 * so adapt the exponent accordingly */
371 sc_val_from_ulong(1, temp);
372 sc_add(_exp(out_val), temp, _exp(out_val));
374 _desc(out_val).clss = NORMAL;
376 /* no further rounding is needed, because rounding overflow means
377 * the carry of the original rounding was propagated all the way
378 * up to the bit left of the radix point. This implies the bits
379 * to the right are all zeros (rounding is +1) */
381 /* check for exponent overflow */
382 sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
383 if (sc_comp(_exp(out_val), temp) != -1) {
384 DEBUGPRINTF(("Exponent overflow!\n"));
385 /* exponent overflow, reaction depends on rounding method:
387 * mode | sign of value | result
388 *--------------------------------------------------------------
389 * TO_NEAREST | + | +inf
391 *--------------------------------------------------------------
392 * TO_POSITIVE | + | +inf
393 * | - | smallest representable value
394 *--------------------------------------------------------------
395 * TO_NEAGTIVE | + | largest representable value
397 *--------------------------------------------------------------
398 * TO_ZERO | + | largest representable value
399 * | - | smallest representable value
400 *--------------------------------------------------------------*/
401 if (_sign(out_val) == 0)
403 /* value is positive */
404 switch (rounding_mode) {
407 _desc(out_val).clss = INF;
412 fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
415 /* value is negative */
416 switch (rounding_mode) {
419 _desc(out_val).clss = INF;
424 fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
433 * Operations involving NaN's must return NaN
435 #define handle_NAN(a, b, result) \
437 if (_desc(a).clss == NAN) { \
438 if (a != result) memcpy(result, a, calc_buffer_size); \
441 if (_desc(b).clss == NAN) { \
442 if (b != result) memcpy(result, b, calc_buffer_size); \
449 * calculate a + b, where a is the value with the bigger exponent
451 static char* _fadd(const char* a, const char* b, char* result)
459 handle_NAN(a, b, result);
461 /* make sure result has a descriptor */
462 if (result != a && result != b)
463 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
465 /* determine if this is an addition or subtraction */
466 sign = _sign(a) ^ _sign(b);
468 /* produce NaN on inf - inf */
469 if (sign && (_desc(a).clss == INF) && (_desc(b).clss == INF))
470 return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);
472 temp = alloca(value_size);
473 exp_diff = alloca(value_size);
475 /* get exponent difference */
476 sc_sub(_exp(a), _exp(b), exp_diff);
478 /* initially set sign to be the sign of a, special treatment of subtraction
479 * when exponents are equal is required though.
480 * Also special care about the sign is needed when the mantissas are equal
482 if (sign && sc_val_to_long(exp_diff) == 0) {
483 switch (sc_comp(_mant(a), _mant(b))) {
485 _sign(result) = _sign(a); /* abs(a) is bigger and a is negative */
488 _sign(result) = (rounding_mode == FC_TONEGATIVE);
491 _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
494 /* can't be reached */
499 _sign(result) = _sign(a);
501 /* sign has been taken care of, check for special cases */
502 if (_desc(a).clss == ZERO) {
503 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
506 if (_desc(b).clss == ZERO) {
507 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
511 if (_desc(a).clss == INF) {
512 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
515 if (_desc(b).clss == INF) {
516 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
520 /* shift the smaller value to the right to align the radix point */
521 /* subnormals have their radix point shifted to the right,
522 * take care of this first */
523 if ((_desc(b).clss == SUBNORMAL) && (_desc(a).clss != SUBNORMAL))
525 sc_val_from_ulong(1, temp);
526 sc_sub(exp_diff, temp, exp_diff);
529 _shift_right(_mant(b), exp_diff, temp);
530 sticky = sc_had_carry();
534 /* if subtracting a little more than the represented value or adding a little
535 * more than the represented value to a negative value this, in addition to the
536 * still set sticky bit, takes account of the 'little more' */
537 char *temp1 = alloca(calc_buffer_size);
538 sc_val_from_ulong(1, temp1);
539 sc_add(temp, temp1, temp);
543 if (sc_comp(_mant(a), temp) == -1)
544 sc_sub(temp, _mant(a), _mant(result));
546 sc_sub(_mant(a), temp, _mant(result));
548 sc_add(_mant(a), temp, _mant(result));
551 /* _normalize expects a 'normal' radix point, adding two subnormals
552 * results in a subnormal radix point -> shifting before normalizing */
553 if ((_desc(a).clss == SUBNORMAL) && (_desc(b).clss == SUBNORMAL))
555 sc_val_from_ulong(1, NULL);
556 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
559 /* resulting exponent is the bigger one */
560 memmove(_exp(result), _exp(a), value_size);
562 return _normalize(result, result, sticky);
568 static char* _fmul(const char* a, const char* b, char* result)
572 handle_NAN(a, b, result);
574 temp = alloca(value_size);
576 if (result != a && result != b)
577 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
579 _sign(result) = _sign(a) ^ _sign(b);
581 /* produce NaN on 0 * inf */
582 if (_desc(a).clss == ZERO) {
583 if (_desc(b).clss == INF)
584 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
586 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
589 if (_desc(b).clss == ZERO) {
590 if (_desc(a).clss == INF)
591 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
593 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
597 if (_desc(a).clss == INF) {
598 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
601 if (_desc(b).clss == INF) {
602 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
606 /* exp = exp(a) + exp(b) - excess */
607 sc_add(_exp(a), _exp(b), _exp(result));
609 sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
610 sc_sub(_exp(result), temp, _exp(result));
612 /* mixed normal, subnormal values introduce an error of 1, correct it */
613 if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
615 sc_val_from_ulong(1, temp);
616 sc_add(_exp(result), temp, _exp(result));
619 sc_mul(_mant(a), _mant(b), _mant(result));
621 /* realign result: after a multiplication the digits right of the radix
622 * point are the sum of the factors' digits after the radix point. As all
623 * values are normalized they both have the same amount of these digits,
624 * which has to be restored by proper shifting
625 * +2 because of the two rounding bits */
626 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
628 _shift_right(_mant(result), temp, _mant(result));
630 return _normalize(result, result, sc_had_carry());
636 static char* _fdiv(const char* a, const char* b, char* result)
638 char *temp, *dividend;
640 handle_NAN(a, b, result);
642 temp = alloca(value_size);
643 dividend = alloca(value_size);
645 if (result != a && result != b)
646 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
648 _sign(result) = _sign(a) ^ _sign(b);
650 /* produce nan on 0/0 and inf/inf */
651 if (_desc(a).clss == ZERO) {
652 if (_desc(b).clss == ZERO)
654 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
657 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
661 if (_desc(b).clss == INF) {
662 if (_desc(a).clss == INF)
664 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
667 sc_val_from_ulong(0, NULL);
668 _save_result(_exp(result));
669 _save_result(_mant(result));
670 _desc(result).clss = ZERO;
675 if (_desc(a).clss == INF) {
677 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
680 if (_desc(b).clss == ZERO) {
681 /* division by zero */
683 fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
685 fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
689 /* exp = exp(a) - exp(b) + excess - 1*/
690 sc_sub(_exp(a), _exp(b), _exp(result));
691 sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
692 sc_add(_exp(result), temp, _exp(result));
694 /* mixed normal, subnormal values introduce an error of 1, correct it */
695 if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
697 sc_val_from_ulong(1, temp);
698 sc_add(_exp(result), temp, _exp(result));
701 /* mant(res) = mant(a) / 1/2mant(b) */
702 /* to gain more bits of precision in the result the dividend could be
703 * shifted left, as this operation does not loose bits. This would not
704 * fit into the integer precision, but due to the rounding bits (which
705 * are always zero because the values are all normalized) the divisor
706 * can be shifted right instead to achieve the same result */
707 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
709 _shift_left(_mant(a), temp, dividend);
712 char *divisor = alloca(calc_buffer_size);
713 sc_val_from_ulong(1, divisor);
714 _shift_right(_mant(b), divisor, divisor);
715 sc_div(dividend, divisor, _mant(result));
718 return _normalize(result, result, sc_had_carry());
722 static void _power_of_ten(int exp, descriptor_t *desc, char *result)
730 /* set new descriptor (else result is supposed to already have one) */
732 memcpy(&_desc(result), desc, sizeof(descriptor_t));
734 build = alloca(value_size);
735 temp = alloca(value_size);
737 sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));
741 /* temp is value of ten now */
742 sc_val_from_ulong(10, NULL);
745 for (exp--; exp > 0; exp--) {
747 sc_mul(build, temp, NULL);
751 /* temp is amount of left shift needed to put the value left of the radix point */
752 sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
754 _shift_left(build, temp, _mant(result));
756 _normalize(result, result, 0);
762 * Truncate the fractional part away.
764 * This does not clip to any integer rang.
766 static char* _trunc(const char *a, char *result)
769 * When exponent == 0 all bits left of the radix point
770 * are the integral part of the value. For 15bit exp_size
771 * this would require a left shift of max. 16383 bits which
773 * But it is enough to ensure that no bit right of the radix
774 * point remains set. This restricts the interesting
775 * exponents to the interval [0, mant_size-1].
776 * Outside this interval the truncated value is either 0 or
777 * it does not have fractional parts.
780 int exp_bias, exp_val;
783 temp = alloca(value_size);
786 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
788 exp_bias = (1<<_desc(a).exponent_size)/2-1;
789 exp_val = sc_val_to_long(_exp(a)) - exp_bias;
792 sc_val_from_ulong(0, NULL);
793 _save_result(_exp(result));
794 _save_result(_mant(result));
795 _desc(result).clss = ZERO;
800 if (exp_val > _desc(a).mantissa_size) {
802 memcpy(result, a, calc_buffer_size);
807 /* set up a proper mask to delete all bits right of the
808 * radix point if the mantissa had been shifted until exp == 0 */
809 sc_max_from_bits(1 + exp_val, 0, temp);
810 sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
811 _shift_left(temp, sc_get_buffer(), temp);
813 /* and the mask and return the result */
814 sc_and(_mant(a), temp, _mant(result));
816 if (a != result) memcpy(_exp(result), _exp(a), value_size);
822 * This does value sanity checking(or should do it), sets up any prerequisites,
823 * calls the proper internal functions, clears up and returns
826 char* _calc(const char *a, const char *b, int opcode, char *result)
829 #ifdef FLTCALC_TRACE_CALC
832 buffer = alloca(100);
835 if (result == NULL) result = calc_buffer;
837 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
841 /* make the value with the bigger exponent the first one */
842 TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
843 if (sc_comp(_exp(a), _exp(b)) == -1)
849 TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
850 temp = alloca(calc_buffer_size);
851 memcpy(temp, b, calc_buffer_size);
852 _sign(temp) = !_sign(b);
853 if (sc_comp(_exp(a), _exp(temp)) == -1)
854 _fadd(temp, a, result);
856 _fadd(a, temp, result);
859 TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
863 TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
867 TRACEPRINTF(("negated "));
868 if (a != result) memcpy(result, a, calc_buffer_size);
869 _sign(result) = !_sign(a);
872 TRACEPRINTF(("truncated to integer "));
876 TRACEPRINTF(("rounded to integer "));
881 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
886 * functions defined in fltcalc.h
888 const void *fc_get_buffer(void)
893 int fc_get_buffer_length(void)
895 return calc_buffer_size;
898 char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
911 int exp_int, hsb, state;
916 char *mant_str, *exp_val, *power_val;
919 if (result == NULL) result = calc_buffer;
921 exp_val = alloca(value_size);
922 power_val = alloca(calc_buffer_size);
923 mant_str = alloca((len)?(len):(strlen(str)));
925 _desc(result).exponent_size = exp_size;
926 _desc(result).mantissa_size = mant_size;
927 _desc(result).clss = NORMAL;
934 while (len == 0 || str-old_str < len)
951 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
958 state = RIGHT_OF_DOT;
969 _fail_char(old_str, len, str - old_str);
975 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
976 mant_str[pos++] = *(str++);
980 state = RIGHT_OF_DOT;
991 mant_str[pos] = '\0';
995 _fail_char(old_str, len, str - old_str);
1001 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1002 mant_str[pos++] = *(str++);
1013 mant_str[pos] = '\0';
1017 _fail_char(old_str, len, str - old_str);
1027 if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
1031 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1032 mant_str[pos] = '\0';
1039 _fail_char(old_str, len, str - old_str);
1045 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1050 case '\0': goto done;
1053 _fail_char(old_str, len, str - old_str);
1056 } /* switch(state) */
1059 sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1061 /* shift to put value left of radix point */
1062 sc_val_from_ulong(mant_size + 2, exp_val);
1064 _shift_left(_mant(result), exp_val, _mant(result));
1066 sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
1068 _normalize(result, result, 0);
1070 if (state == EXPONENT) {
1071 exp_int -= atoi(str-pos);
1074 _power_of_ten(exp_int, &_desc(result), power_val);
1076 _fdiv(result, power_val, result);
1080 /* XXX excuse of an implementation to make things work */
1082 char *tmp = alloca(calc_buffer_size);
1085 #ifdef HAVE_LONG_DOUBLE
1086 val = strtold(str, NULL);
1087 DEBUGPRINTF(("val_from_str(%s)\n", str));
1088 fc_val_from_float(val, 15, 64, tmp);
1090 val = strtod(str, NULL);
1091 DEBUGPRINTF(("val_from_str(%s)\n", str));
1092 fc_val_from_float(val, 11, 52, tmp);
1093 #endif /* HAVE_LONG_DOUBLE */
1094 return fc_cast(tmp, exp_size, mant_size, result);
1098 char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
1101 int bias_res, bias_val, mant_val;
1103 UINT32 sign, exponent, mantissa0, mantissa1;
1106 bias_res = ((1<<exp_size)/2-1);
1108 #ifdef HAVE_LONG_DOUBLE
1111 sign = (srcval.val.high & 0x00008000) != 0;
1112 exponent = (srcval.val.high & 0x00007FFF) ;
1113 mantissa0 = srcval.val.mid;
1114 mantissa1 = srcval.val.low;
1115 #else /* no long double */
1118 sign = (srcval.val.high & 0x80000000) != 0;
1119 exponent = (srcval.val.high & 0x7FF00000) >> 20;
1120 mantissa0 = srcval.val.high & 0x000FFFFF;
1121 mantissa1 = srcval.val.low;
1124 #ifdef HAVE_LONG_DOUBLE
1125 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)); */
1126 DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1128 TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1129 DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1132 if (result == NULL) result = calc_buffer;
1133 temp = alloca(value_size);
1135 _desc(result).exponent_size = exp_size;
1136 _desc(result).mantissa_size = mant_size;
1139 _sign(result) = sign;
1141 /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1142 * encoding is needed. the function can return immediately in these cases */
1144 _desc(result).clss = NAN;
1145 TRACEPRINTF(("val_from_float resulted in NAN\n"));
1148 else if (isinf(l)) {
1149 _desc(result).clss = INF;
1150 TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1154 /* build exponent, because input and output exponent and mantissa sizes may differ
1155 * this looks more complicated than it is: unbiased input exponent + output bias,
1156 * minus the mantissa difference which is added again later when the output float
1157 * becomes normalized */
1158 #ifdef HAVE_EXPLICIT_ONE
1159 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1161 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1164 /* build mantissa representation */
1165 #ifndef HAVE_EXPLICIT_ONE
1168 /* insert the hidden bit */
1169 sc_val_from_ulong(1, temp);
1170 sc_val_from_ulong(mant_val + 2, NULL);
1171 _shift_left(temp, sc_get_buffer(), NULL);
1176 sc_val_from_ulong(0, NULL);
1179 _save_result(_mant(result));
1181 /* bits from the upper word */
1182 sc_val_from_ulong(mantissa0, temp);
1183 sc_val_from_ulong(34, NULL);
1184 _shift_left(temp, sc_get_buffer(), temp);
1185 sc_or(_mant(result), temp, _mant(result));
1187 /* bits from the lower word */
1188 sc_val_from_ulong(mantissa1, temp);
1189 sc_val_from_ulong(2, NULL);
1190 _shift_left(temp, sc_get_buffer(), temp);
1191 sc_or(_mant(result), temp, _mant(result));
1193 /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1194 * origin one to the left */
1197 sc_val_from_ulong(1, NULL);
1198 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1201 _normalize(result, result, 0);
1203 TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, calc_buffer_size, FC_PACKED)));
1208 LLDBL fc_val_to_float(const void *val)
1222 #ifdef HAVE_LONG_DOUBLE
1223 char result_exponent = 15;
1224 char result_mantissa = 64;
1226 char result_exponent = 11;
1227 char result_mantissa = 52;
1230 temp = alloca(calc_buffer_size);
1231 #ifdef HAVE_EXPLICIT_ONE
1232 value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1234 value = fc_cast(val, result_exponent, result_mantissa, temp);
1237 sign = _sign(value);
1239 /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1240 * lead to wrong results */
1241 exponent = sc_val_to_long(_exp(value)) ;
1243 sc_val_from_ulong(2, NULL);
1244 _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1249 for (byte_offset = 0; byte_offset < 4; byte_offset++)
1250 mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1252 for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1253 mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1255 #ifdef HAVE_LONG_DOUBLE
1256 buildval.val.high = sign << 15;
1257 buildval.val.high |= exponent;
1258 buildval.val.mid = mantissa0;
1259 buildval.val.low = mantissa1;
1260 #else /* no long double */
1261 mantissa0 &= 0x000FFFFF; /* get rid of garbage */
1262 buildval.val.high = sign << 31;
1263 buildval.val.high |= exponent << 20;
1264 buildval.val.high |= mantissa0;
1265 buildval.val.low = mantissa1;
1268 TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1272 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1274 const char *value = (const char*) val;
1276 int exp_offset, val_bias, res_bias;
1278 if (result == NULL) result = calc_buffer;
1279 temp = alloca(value_size);
1281 if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1283 if (value != result) memcpy(result, value, calc_buffer_size);
1287 /* set the descriptor of the new value */
1288 _desc(result).exponent_size = exp_size;
1289 _desc(result).mantissa_size = mant_size;
1290 _desc(result).clss = _desc(value).clss;
1292 _sign(result) = _sign(value);
1294 /* when the mantissa sizes differ normalizing has to shift to align it.
1295 * this would change the exponent, which is unwanted. So calculate this
1296 * offset and add it */
1297 val_bias = (1<<_desc(value).exponent_size)/2-1;
1298 res_bias = (1<<exp_size)/2-1;
1300 exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1301 sc_val_from_long(exp_offset, temp);
1302 sc_add(_exp(value), temp, _exp(result));
1304 /* _normalize expects normalized radix point */
1305 if (_desc(val).clss == SUBNORMAL) {
1306 sc_val_from_ulong(1, NULL);
1307 _shift_left(_mant(val), sc_get_buffer(), _mant(result));
1308 } else if (value != result) {
1309 memcpy(_mant(result), _mant(value), value_size);
1311 memmove(_mant(result), _mant(value), value_size);
1314 _normalize(result, result, 0);
1315 TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, value_size, FC_PACKED)));
1319 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1321 if (result == NULL) result = calc_buffer;
1323 _desc(result).exponent_size = exponent_size;
1324 _desc(result).mantissa_size = mantissa_size;
1325 _desc(result).clss = NORMAL;
1329 sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1331 sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1332 sc_val_from_ulong(2, NULL);
1333 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1338 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1340 if (result == NULL) result = calc_buffer;
1342 fc_get_max(exponent_size, mantissa_size, result);
1348 char* fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1350 if (result == NULL) result = calc_buffer;
1352 _desc(result).exponent_size = exponent_size;
1353 _desc(result).mantissa_size = mantissa_size;
1354 _desc(result).clss = NAN;
1358 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1360 /* signaling NaN has non-zero mantissa with msb not set */
1361 sc_val_from_ulong(1, _mant(result));
1366 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1368 if (result == NULL) result = calc_buffer;
1370 _desc(result).exponent_size = exponent_size;
1371 _desc(result).mantissa_size = mantissa_size;
1372 _desc(result).clss = NAN;
1376 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1378 /* quiet NaN has the msb of the mantissa set, so shift one there */
1379 sc_val_from_ulong(1, _mant(result));
1380 /* mantissa_size >+< 1 because of two extra rounding bits */
1381 sc_val_from_ulong(mantissa_size + 1, NULL);
1382 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1387 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1389 if (result == NULL) result = calc_buffer;
1391 _desc(result).exponent_size = exponent_size;
1392 _desc(result).mantissa_size = mantissa_size;
1393 _desc(result).clss = NORMAL;
1397 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1399 sc_val_from_ulong(0, _mant(result));
1404 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1406 if (result == NULL) result = calc_buffer;
1408 fc_get_plusinf(exponent_size, mantissa_size, result);
1414 int fc_comp(const void *a, const void *b)
1416 const char *val_a = (const char*)a;
1417 const char *val_b = (const char*)b;
1421 * shortcut: if both values are identical, they are either
1422 * Unordered if NaN or equal
1425 return _desc(val_a).clss == NAN ? 2 : 0;
1427 /* unordered if one is a NaN */
1428 if (_desc(val_a).clss == NAN || _desc(val_b).clss == NAN)
1431 /* zero is equal independent of sign */
1432 if ((_desc(val_a).clss == ZERO) && (_desc(val_b).clss == ZERO))
1435 /* different signs make compare easy */
1436 if (_sign(val_a) != _sign(val_b))
1437 return (_sign(val_a)==0)?(1):(-1);
1439 mul = _sign(a) ? -1 : 1;
1441 /* both infinity means equality */
1442 if ((_desc(val_a).clss == INF) && (_desc(val_b).clss == INF))
1445 /* infinity is bigger than the rest */
1446 if (_desc(val_a).clss == INF)
1448 if (_desc(val_b).clss == INF)
1451 /* check first exponent, that mantissa if equal */
1452 switch (sc_comp(_exp(val_a), _exp(val_b))) {
1458 return sc_comp(_mant(val_a), _mant(val_b)) * mul;
1464 int fc_is_zero(const void *a)
1466 return _desc(a).clss == ZERO;
1469 int fc_is_negative(const void *a)
1474 int fc_is_inf(const void *a)
1476 return _desc(a).clss == INF;
1479 int fc_is_nan(const void *a)
1481 return _desc(a).clss == NAN;
1484 int fc_is_subnormal(const void *a)
1486 return _desc(a).clss == SUBNORMAL;
1489 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1494 val = (const char*)a;
1496 mul_1 = alloca(calc_buffer_size);
1500 switch (_desc(val).clss) {
1502 if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1503 else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1506 snprintf(buf, buflen, "NAN");
1509 snprintf(buf, buflen, "0.0");
1512 /* XXX to be implemented */
1513 #ifdef HAVE_LONG_DOUBLE
1514 /* XXX 30 is arbitrary */
1515 snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1517 snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1523 switch (_desc(val).clss) {
1525 if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1526 else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1529 snprintf(buf, buflen, "NAN");
1532 snprintf(buf, buflen, "0.0");
1535 #ifdef HAVE_LONG_DOUBLE
1536 snprintf(buf, buflen, "%LA", fc_val_to_float(val));
1538 snprintf(buf, buflen, "%A", fc_val_to_float(val));
1545 snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), value_size*4, SC_HEX, 0));
1546 buf[buflen - 1] = '\0';
1552 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1554 /* this is used to cache the packed version of the value */
1555 static char *pack = NULL;
1557 if (pack == NULL) pack = xmalloc(value_size);
1560 _pack((const char*)value, pack);
1562 return sc_sub_bits(pack, num_bits, byte_ofs);
1565 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1567 if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1568 rounding_mode = mode;
1570 return rounding_mode;
1573 fc_rounding_mode_t fc_get_rounding_mode(void)
1575 return rounding_mode;
1578 void init_fltcalc(int precision)
1580 if (calc_buffer == NULL) {
1581 /* does nothing if already init */
1582 if (precision == 0) precision = FC_DEFAULT_PRECISION;
1584 init_strcalc(precision + 4);
1586 /* needs additionally two bits to round, a bit as explicit 1., and one for
1587 * addition overflow */
1588 max_precision = sc_get_precision() - 4;
1589 if (max_precision < precision)
1590 printf("WARING: not enough precision available, using %d\n", max_precision);
1592 rounding_mode = FC_TONEAREST;
1593 value_size = sc_get_buffer_length();
1595 EXPONENT_POS = SIGN_POS + sizeof(char);
1596 MANTISSA_POS = EXPONENT_POS + value_size;
1597 DESCRIPTOR_POS = MANTISSA_POS + value_size;
1598 calc_buffer_size = DESCRIPTOR_POS + sizeof(descriptor_t);
1600 calc_buffer = xmalloc(calc_buffer_size);
1601 memset(calc_buffer, 0, calc_buffer_size);
1602 DEBUGPRINTF(("init fltcalc:\n\tVALUE_SIZE = %d\n\tSIGN_POS = %d\n\tEXPONENT_POS = %d\n\tMANTISSA_POS = %d\n\tDESCRIPTOR_POS = %d\n\tCALC_BUFFER_SIZE = %d\n\tcalc_buffer = %p\n\n", value_size, SIGN_POS, EXPONENT_POS, MANTISSA_POS, DESCRIPTOR_POS, calc_buffer_size, calc_buffer));
1603 #ifdef HAVE_LONG_DOUBLE
1604 DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1606 DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1608 #ifdef WORDS_BIGENDIAN
1609 DEBUGPRINTF(("\tWord order is big endian\n\n"));
1611 DEBUGPRINTF(("\tWord order is little endian\n\n"));
1616 void finish_fltcalc (void) {
1617 free(calc_buffer); calc_buffer = NULL;
1620 /* definition of interface functions */