3 * File name: ir/tv/fltcalc.c
9 * Copyright: (c) 2003 Universität Karlsruhe
10 * Licence: This file protected by GPL - GNU GENERAL PUBLIC LICENSE.
22 #include <math.h> /* need isnan() and isinf() (will be changed)*/
23 /* undef some reused constants defined by math.h */
28 #ifdef HAVE_INTTYPES_H
43 typedef uint32_t UINT32;
45 #ifdef HAVE_LONG_DOUBLE
46 #ifdef WORDS_BIGENDIAN
53 volatile long double d;
62 volatile long double d;
66 #ifdef WORDS_BIGENDIAN
103 #define CLEAR_BUFFER(buffer) memset(buffer, 0, CALC_BUFFER_SIZE)
105 /* because variable sized structs are impossible, the internal
106 * value is represented as a pseudo-struct char array, addressed
109 * char sign; // 0 for positive, 1 for negative
110 * char exp[VALUE_SIZE];
111 * char mant[VALUE_SIZE];
115 #define _sign(a) (((char*)a)[SIGN_POS])
116 #define _exp(a) (&((char*)a)[EXPONENT_POS])
117 #define _mant(a) (&((char*)a)[MANTISSA_POS])
118 #define _desc(a) (*(descriptor_t *)&((char*)a)[DESCRIPTOR_POS])
120 #define _save_result(x) memcpy((x), sc_get_buffer(), VALUE_SIZE)
121 #define _shift_right(x, y, b) sc_shr((x), (y), VALUE_SIZE*4, 0, (b))
122 #define _shift_left(x, y, b) sc_shl((x), (y), VALUE_SIZE*4, 0, (b))
124 #define FC_DEFINE1(code) char* fc_##code(const void *a, void *result) \
126 return _calc((const char*)a, NULL, FC_##code, (char*)result); \
129 #define FC_DEFINE2(code) char* fc_##code(const void *a, const void *b, void *result) \
131 return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \
134 #define FUNC_PTR(code) fc_##code
137 # define DEBUGPRINTF(x) printf x
139 # define DEBUGPRINTF(x) ((void)0)
142 #if FLTCALC_TRACE_CALC
143 # define TRACEPRINTF(x) printf x
145 # define TRACEPRINTF(x) ((void)0)
148 static char *calc_buffer = NULL;
150 static fc_rounding_mode_t ROUNDING_MODE;
152 static int CALC_BUFFER_SIZE;
153 static int VALUE_SIZE;
155 static int EXPONENT_POS;
156 static int MANTISSA_POS;
157 static int DESCRIPTOR_POS;
159 static int max_precision;
164 static void _fail_char(const char *str, unsigned int len, int pos)
167 printf("ERROR: Unexpected character '%c'\n", *(str + pos));
169 printf("ERROR: Unexpected end of string\n");
170 while (len-- && *str) printf("%c", *str++); printf("\n");
171 while (pos--) printf(" "); printf("^\n");
172 /* the front end has to to check constant strings */
177 /* pack machine-like */
178 static char* _pack(const char *int_float, char *packed)
184 temp = alloca(VALUE_SIZE);
185 shift_val = alloca(VALUE_SIZE);
187 switch (_desc(int_float).class) {
189 val_buffer = alloca(CALC_BUFFER_SIZE);
190 fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
191 int_float = val_buffer;
195 val_buffer = alloca(CALC_BUFFER_SIZE);
196 fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
197 _sign(val_buffer) = _sign(int_float);
198 int_float = val_buffer;
205 sc_val_from_ulong(_sign(int_float), temp);
207 sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
208 _shift_left(temp, sc_get_buffer(), packed);
210 /* extract exponent */
211 sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);
213 _shift_left(_exp(int_float), shift_val, temp);
215 sc_or(temp, packed, packed);
217 /* extract mantissa */
218 /* remove 2 rounding bits */
219 sc_val_from_ulong(2, shift_val);
220 _shift_right(_mant(int_float), shift_val, temp);
222 /* remove leading 1 (or 0 if denormalized) */
223 sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
224 sc_and(temp, shift_val, temp);
227 sc_or(temp, packed, packed);
232 char* _normalize(const char *in_val, char *out_val, int sticky)
235 char lsb, guard, round, round_dir = 0;
238 temp = alloca(VALUE_SIZE);
240 /* +2: save two rounding bits at the end */
241 hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
243 if (in_val != out_val)
245 _sign(out_val) = _sign(in_val);
246 memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
249 _desc(out_val).class = NORMAL;
251 /* mantissa all zeroes, so zero exponent (because of explicit one)*/
252 if (hsb == 2 + _desc(in_val).mantissa_size)
254 sc_val_from_ulong(0, _exp(out_val));
258 /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
262 sc_val_from_ulong(-hsb-1, temp);
264 _shift_right(_mant(in_val), temp, _mant(out_val));
266 /* remember if some bits were shifted away */
267 if (!sticky) sticky = sc_had_carry();
269 sc_add(_exp(in_val), temp, _exp(out_val));
274 sc_val_from_ulong(hsb+1, temp);
276 _shift_left(_mant(in_val), temp, _mant(out_val));
278 sc_sub(_exp(in_val), temp, _exp(out_val));
281 /* check for exponent underflow */
282 if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
283 DEBUGPRINTF(("Exponent underflow!\n"));
284 /* exponent underflow */
285 /* shift the mantissa right to have a zero exponent */
286 sc_val_from_ulong(1, temp);
287 sc_sub(temp, _exp(out_val), NULL);
289 _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
290 if (!sticky) sticky = sc_had_carry();
291 /* denormalized means exponent of zero */
292 sc_val_from_ulong(0, _exp(out_val));
294 _desc(out_val).class = SUBNORMAL;
297 /* perform rounding by adding a value that clears the guard bit and the round bit
298 * and either causes a carry to round up or not */
299 /* get the last 3 bits of the value */
300 lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
301 guard = (lsb&0x2)>>1;
304 switch (ROUNDING_MODE)
307 /* round to nearest representable value, if in doubt choose the version
309 round_dir = guard && (sticky || round || lsb>>2);
312 /* if positive: round to one if the exact value is bigger, else to zero */
313 round_dir = (!_sign(out_val) && (guard || round || sticky));
316 /* if negative: round to one if the exact value is bigger, else to zero */
317 round_dir = (_sign(out_val) && (guard || round || sticky));
320 /* always round to 0 (chopping mode) */
324 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"));
328 guard = (round^guard)<<1;
329 lsb = !(round || guard)<<2 | guard | round;
333 lsb = -((guard<<1) | round);
336 /* add the rounded value */
338 sc_val_from_long(lsb, temp);
339 sc_add(_mant(out_val), temp, _mant(out_val));
342 /* could have rounded down to zero */
343 if (sc_is_zero(_mant(out_val)) && (_desc(out_val).class == SUBNORMAL))
344 _desc(out_val).class = ZERO;
346 /* check for rounding overflow */
347 hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
348 if ((_desc(out_val).class != SUBNORMAL) && (hsb < -1))
350 sc_val_from_ulong(1, temp);
351 _shift_right(_mant(out_val), temp, _mant(out_val));
353 sc_add(_exp(out_val), temp, _exp(out_val));
355 else if ((_desc(out_val).class == SUBNORMAL) && (hsb == -1))
357 /* overflow caused the matissa to be normal again,
358 * so adapt the exponent accordingly */
359 sc_val_from_ulong(1, temp);
360 sc_add(_exp(out_val), temp, _exp(out_val));
362 _desc(out_val).class = NORMAL;
364 /* no further rounding is needed, because rounding overflow means
365 * the carry of the original rounding was propagated all the way
366 * up to the bit left of the radix point. This implies the bits
367 * to the right are all zeros (rounding is +1) */
369 /* check for exponent overflow */
370 sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
371 if (sc_comp(_exp(out_val), temp) != -1) {
372 DEBUGPRINTF(("Exponent overflow!\n"));
373 /* exponent overflow, reaction depends on rounding method:
375 * mode | sign of value | result
376 *--------------------------------------------------------------
377 * TO_NEAREST | + | +inf
379 *--------------------------------------------------------------
380 * TO_POSITIVE | + | +inf
381 * | - | smallest representable value
382 *--------------------------------------------------------------
383 * TO_NEAGTIVE | + | largest representable value
385 *--------------------------------------------------------------
386 * TO_ZERO | + | largest representable value
387 * | - | smallest representable value
388 *--------------------------------------------------------------*/
389 if (_sign(out_val) == 0)
391 /* value is positive */
392 switch (ROUNDING_MODE) {
395 _desc(out_val).class = INF;
400 fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
403 /* value is negative */
404 switch (ROUNDING_MODE) {
407 _desc(out_val).class = INF;
412 fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
421 * calculate a + b, where a is the value with the bigger exponent
423 static char* _add(const char* a, const char* b, char* result)
431 if (_desc(a).class == NAN) {
432 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
435 if (_desc(b).class == NAN) {
436 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
440 /* make sure result has a descriptor */
441 if (result != a && result != b)
442 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
444 /* determine if this is an addition or subtraction */
445 sign = _sign(a) ^ _sign(b);
447 /* produce nan on inf - inf */
448 if (sign && (_desc(a).class == INF) && (_desc(b).class == INF))
449 return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);
451 temp = alloca(VALUE_SIZE);
452 exp_diff = alloca(VALUE_SIZE);
454 /* get exponent difference */
455 sc_sub(_exp(a), _exp(b), exp_diff);
457 /* initially set sign to be the sign of a, special treatment of subtraction
458 * when exponents are equal is required though.
459 * Also special care about the sign is needed when the mantissas are equal
461 if (sign && sc_val_to_long(exp_diff) == 0) {
462 switch (sc_comp(_mant(a), _mant(b))) {
464 if (_sign(a)) _sign(result) = 1; /* abs(a) is bigger and a is negative */
465 else _sign(result) = 0;
468 if (ROUNDING_MODE == FC_TONEGATIVE)
474 if (_sign(b)) _sign(result) = 1; /* abs(b) is bigger and b is negative */
475 else _sign(result) = 0;
478 /* can't be reached */
482 _sign(result) = _sign(a);
485 /* sign has been taken care of, check for special cases */
486 if (_desc(a).class == ZERO) {
487 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
490 if (_desc(b).class == ZERO) {
491 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
495 if (_desc(a).class == INF) {
496 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
499 if (_desc(b).class == INF) {
500 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
504 /* shift the smaller value to the right to align the radix point */
505 /* subnormals have their radix point shifted to the right,
506 * take care of this first */
507 if ((_desc(b).class == SUBNORMAL) && (_desc(a).class != SUBNORMAL))
509 sc_val_from_ulong(1, temp);
510 sc_sub(exp_diff, temp, exp_diff);
513 _shift_right(_mant(b), exp_diff, temp);
514 sticky = sc_had_carry();
518 /* if subtracting a little more than the represented value or adding a little
519 * more than the represented value to a negative value this, in addition to the
520 * still set sticky bit, takes account of the 'little more' */
521 char *temp1 = alloca(CALC_BUFFER_SIZE);
522 sc_val_from_ulong(1, temp1);
523 sc_add(temp, temp1, temp);
527 if (sc_comp(_mant(a), temp) == -1)
528 sc_sub(temp, _mant(a), _mant(result));
530 sc_sub(_mant(a), temp, _mant(result));
532 sc_add(_mant(a), temp, _mant(result));
535 /* _normalize expects a 'normal' radix point, adding two subnormals
536 * results in a subnormal radix point -> shifting before normalizing */
537 if ((_desc(a).class == SUBNORMAL) && (_desc(b).class == SUBNORMAL))
539 sc_val_from_ulong(1, NULL);
540 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
543 /* resulting exponent is the bigger one */
544 memmove(_exp(result), _exp(a), VALUE_SIZE);
546 return _normalize(result, result, sticky);
549 static char* _mul(const char* a, const char* b, char* result)
553 if (_desc(a).class == NAN) {
554 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
557 if (_desc(b).class == NAN) {
558 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
562 temp = alloca(VALUE_SIZE);
564 if (result != a && result != b)
565 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
567 _sign(result) = _sign(a) ^ _sign(b);
569 /* produce nan on 0 * inf */
570 if (_desc(a).class == ZERO) {
571 if (_desc(b).class == INF)
572 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
574 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
577 if (_desc(b).class == ZERO) {
578 if (_desc(a).class == INF)
579 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
581 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1);
585 if (_desc(a).class == INF) {
586 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
589 if (_desc(b).class == INF) {
590 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1);
594 /* exp = exp(a) + exp(b) - excess */
595 sc_add(_exp(a), _exp(b), _exp(result));
597 sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
598 sc_sub(_exp(result), temp, _exp(result));
600 /* mixed normal, subnormal values introduce an error of 1, correct it */
601 if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL))
603 sc_val_from_ulong(1, temp);
604 sc_add(_exp(result), temp, _exp(result));
607 sc_mul(_mant(a), _mant(b), _mant(result));
609 /* realign result: after a multiplication the digits right of the radix
610 * point are the sum of the factors' digits after the radix point. As all
611 * values are normalized they both have the same amount of these digits,
612 * which has to be restored by proper shifting
613 * +2 because of the two rounding bits */
614 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
616 _shift_right(_mant(result), temp, _mant(result));
618 return _normalize(result, result, sc_had_carry());
621 static char* _div(const char* a, const char* b, char* result)
623 char *temp, *dividend;
625 if (_desc(a).class == NAN) {
626 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
629 if (_desc(b).class == NAN) {
630 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
634 temp = alloca(VALUE_SIZE);
635 dividend = alloca(VALUE_SIZE);
637 if (result != a && result != b)
638 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
640 _sign(result) = _sign(a) ^ _sign(b);
642 /* produce nan on 0/0 and inf/inf */
643 if (_desc(a).class == ZERO) {
644 if (_desc(b).class == ZERO)
646 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
649 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
653 if (_desc(b).class == INF) {
654 if (_desc(a).class == INF)
656 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
659 sc_val_from_ulong(0, NULL);
660 _save_result(_exp(result));
661 _save_result(_mant(result));
662 _desc(result).class = ZERO;
667 if (_desc(a).class == INF) {
669 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
672 if (_desc(b).class == ZERO) {
673 /* division by zero */
675 fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
677 fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
681 /* exp = exp(a) - exp(b) + excess - 1*/
682 sc_sub(_exp(a), _exp(b), _exp(result));
683 sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
684 sc_add(_exp(result), temp, _exp(result));
686 /* mixed normal, subnormal values introduce an error of 1, correct it */
687 if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL))
689 sc_val_from_ulong(1, temp);
690 sc_add(_exp(result), temp, _exp(result));
693 /* mant(res) = mant(a) / 1/2mant(b) */
694 /* to gain more bits of precision in the result the dividend could be
695 * shifted left, as this operation does not loose bits. This would not
696 * fit into the integer precision, but due to the rounding bits (which
697 * are always zero because the values are all normalized) the divisor
698 * can be shifted right instead to achieve the same result */
699 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
701 _shift_left(_mant(a), temp, dividend);
704 char *divisor = alloca(CALC_BUFFER_SIZE);
705 sc_val_from_ulong(1, divisor);
706 _shift_right(_mant(b), divisor, divisor);
707 sc_div(dividend, divisor, _mant(result));
710 return _normalize(result, result, sc_had_carry());
713 void _power_of_ten(int exp, descriptor_t *desc, char *result)
721 /* set new descriptor (else result is supposed to already have one) */
723 memcpy(&_desc(result), desc, sizeof(descriptor_t));
725 build = alloca(VALUE_SIZE);
726 temp = alloca(VALUE_SIZE);
728 sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));
732 /* temp is value of ten now */
733 sc_val_from_ulong(10, NULL);
736 for (exp--; exp > 0; exp--) {
738 sc_mul(build, temp, NULL);
742 /* temp is amount of leftshift needed to put the value left of the radix point */
743 sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
745 _shift_left(build, temp, _mant(result));
747 _normalize(result, result, 0);
751 static char* _trunc(const char *a, char *result)
753 /* when exponent == 0 all bits left of the radix point
754 * are the integral part of the value. For 15bit exp_size
755 * this would require a leftshift of max. 16383 bits which
757 * But it is enough to ensure that no bit right of the radix
758 * point remains set. This restricts the interesting
759 * exponents to the interval [0, mant_size-1].
760 * Outside this interval the truncated value is either 0 or
761 * it is are already truncated */
763 int exp_bias, exp_val;
766 temp = alloca(VALUE_SIZE);
769 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
771 exp_bias = (1<<_desc(a).exponent_size)/2-1;
772 exp_val = sc_val_to_long(_exp(a)) - exp_bias;
775 sc_val_from_ulong(0, NULL);
776 _save_result(_exp(result));
777 _save_result(_mant(result));
778 _desc(result).class = ZERO;
783 if (exp_val > _desc(a).mantissa_size) {
785 memcpy(result, a, CALC_BUFFER_SIZE);
790 /* set up a proper mask to delete all bits right of the
791 * radix point if the mantissa had been shifted until exp == 0 */
792 sc_max_from_bits(1 + exp_val, 0, temp);
793 sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
794 _shift_left(temp, sc_get_buffer(), temp);
796 /* and the mask and return the result */
797 sc_and(_mant(a), temp, _mant(result));
799 if (a != result) memcpy(_exp(result), _exp(a), VALUE_SIZE);
805 * This does value sanity checking(or should do it), sets up any prerequisites,
806 * calls the proper internal functions, clears up and returns
808 char* _calc(const char *a, const char *b, int opcode, char *result)
811 #ifdef FLTCALC_TRACE_CALC
814 buffer = alloca(100);
817 if (result == NULL) result = calc_buffer;
819 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_PACKED)));
823 /* make the value with the bigger exponent the first one */
824 TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_PACKED)));
825 if (sc_comp(_exp(a), _exp(b)) == -1)
831 TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_PACKED)));
832 temp = alloca(CALC_BUFFER_SIZE);
833 memcpy(temp, b, CALC_BUFFER_SIZE);
834 _sign(temp) = !_sign(b);
835 if (sc_comp(_exp(a), _exp(temp)) == -1)
836 _add(temp, a, result);
838 _add(a, temp, result);
841 TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_PACKED)));
845 TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_PACKED)));
849 TRACEPRINTF(("negated "));
850 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
851 _sign(result) = !_sign(a);
860 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_PACKED)));
865 * functions defined in fltcalc.h
867 const void *fc_get_buffer(void)
872 const int fc_get_buffer_length(void)
874 return CALC_BUFFER_SIZE;
877 char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
890 int exp_int, hsb, state;
895 char *mant_str, *exp_val, *power_val;
897 if (result == NULL) result = calc_buffer;
899 exp_val = alloca(VALUE_SIZE);
900 power_val = alloca(CALC_BUFFER_SIZE);
901 mant_str = alloca((len)?(len):(strlen(str)));
903 _desc(result).exponent_size = exp_size;
904 _desc(result).mantissa_size = mant_size;
905 _desc(result).class = NORMAL;
912 while (len == 0 || str-old_str < len)
929 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
936 state = RIGHT_OF_DOT;
947 _fail_char(old_str, len, str - old_str);
953 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
954 mant_str[pos++] = *(str++);
958 state = RIGHT_OF_DOT;
969 mant_str[pos] = '\0';
973 _fail_char(old_str, len, str - old_str);
979 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
980 mant_str[pos++] = *(str++);
991 mant_str[pos] = '\0';
995 _fail_char(old_str, len, str - old_str);
1005 if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
1009 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1010 mant_str[pos] = '\0';
1017 _fail_char(old_str, len, str - old_str);
1023 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
1028 case '\0': goto done;
1031 _fail_char(old_str, len, str - old_str);
1034 } /* switch(state) */
1037 sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1039 /* shift to put value left of radix point */
1040 sc_val_from_ulong(mant_size + 2, exp_val);
1042 _shift_left(_mant(result), exp_val, _mant(result));
1044 sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
1046 _normalize(result, result, 0);
1048 if (state == EXPONENT) {
1049 exp_int -= atoi(str-pos);
1052 _power_of_ten(exp_int, &_desc(result), power_val);
1054 _div(result, power_val, result);
1059 /* XXX excuse of an implementation to make things work */
1061 #ifdef HAVE_LONG_DOUBLE
1062 val = strtold(str, NULL);
1064 val = strtod(str, NULL);
1067 DEBUGPRINTF(("val_from_str(%s)\n", str));
1068 return fc_val_from_float(val, exp_size, mant_size, result);
1072 char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
1075 int bias_res, bias_val, mant_val;
1077 UINT32 sign, exponent, mantissa0, mantissa1;
1080 bias_res = ((1<<exp_size)/2-1);
1082 #ifdef HAVE_LONG_DOUBLE
1085 sign = (srcval.val.high & 0x00008000) != 0;
1086 exponent = (srcval.val.high & 0x00007FFF) ;
1087 mantissa0 = srcval.val.mid;
1088 mantissa1 = srcval.val.low;
1089 #else /* no long double */
1092 sign = (srcval.val.high & 0x80000000) != 0;
1093 exponent = (srcval.val.high & 0x7FF00000) >> 20;
1094 mantissa0 = srcval.val.high & 0x000FFFFF;
1095 mantissa1 = srcval.val.low;
1098 #ifdef HAVE_LONG_DOUBLE
1099 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)); */
1100 DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1102 TRACEPRINTF(("val_from_float(%.8X%.8X)\n", srcval.val.high, srcval.val.low));
1103 DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1106 if (result == NULL) result = calc_buffer;
1107 temp = alloca(VALUE_SIZE);
1109 _desc(result).exponent_size = exp_size;
1110 _desc(result).mantissa_size = mant_size;
1113 _sign(result) = sign;
1115 /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1116 * encoding is needed. the function can return immediately in these cases */
1118 _desc(result).class = NAN;
1119 TRACEPRINTF(("val_from_float resulted in NAN\n"));
1122 else if (isinf(l)) {
1123 _desc(result).class = INF;
1124 TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1128 /* build exponent, because input and output exponent and mantissa sizes may differ
1129 * this looks more complicated than it is: unbiased input exponent + output bias,
1130 * minus the mantissa difference which is added again later when the output float
1131 * becomes normalized */
1132 #ifdef HAVE_EXPLICIT_ONE
1133 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size-1), _exp(result));
1135 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1138 /* build mantissa representation */
1139 #ifndef HAVE_EXPLICIT_ONE
1142 /* insert the hidden bit */
1143 sc_val_from_ulong(1, temp);
1144 sc_val_from_ulong(mant_val + 2, NULL);
1145 _shift_left(temp, sc_get_buffer(), NULL);
1150 sc_val_from_ulong(0, NULL);
1153 _save_result(_mant(result));
1155 /* bits from the upper word */
1156 sc_val_from_ulong(mantissa0, temp);
1157 sc_val_from_ulong(34, NULL);
1158 _shift_left(temp, sc_get_buffer(), temp);
1159 sc_or(_mant(result), temp, _mant(result));
1161 /* bits from the lower word */
1162 sc_val_from_ulong(mantissa1, temp);
1163 sc_val_from_ulong(2, NULL);
1164 _shift_left(temp, sc_get_buffer(), temp);
1165 sc_or(_mant(result), temp, _mant(result));
1167 /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1168 * origin one to the left */
1171 sc_val_from_ulong(1, NULL);
1172 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1175 _normalize(result, result, 0);
1177 TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, CALC_BUFFER_SIZE, FC_PACKED)));
1182 LLDBL fc_val_to_float(const void *val)
1196 #ifdef HAVE_LONG_DOUBLE
1197 char result_exponent = 15;
1198 char result_mantissa = 64;
1200 char result_exponent = 11;
1201 char result_mantissa = 52;
1204 temp = alloca(CALC_BUFFER_SIZE);
1205 #ifdef HAVE_EXPLICIT_ONE
1206 value = fc_cast(val, result_exponent, result_mantissa-1, temp);
1208 value = fc_cast(val, result_exponent, result_mantissa, temp);
1211 sign = _sign(value);
1213 /* @@@ long double exponent is 15bit, so the use of sc_val_to_long should not
1214 * lead to wrong results */
1215 exponent = sc_val_to_long(_exp(value)) ;
1217 sc_val_from_ulong(2, NULL);
1218 _shift_right(_mant(value), sc_get_buffer(), _mant(value));
1223 for (byte_offset = 0; byte_offset < 4; byte_offset++)
1224 mantissa1 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << (byte_offset<<3);
1226 for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1227 mantissa0 |= sc_sub_bits(_mant(value), result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1229 #ifdef HAVE_LONG_DOUBLE
1230 buildval.val.high = sign << 15;
1231 buildval.val.high |= exponent;
1232 buildval.val.mid = mantissa0;
1233 buildval.val.low = mantissa1;
1234 #else /* no long double */
1235 mantissa0 &= 0x000FFFFF; /* get rid of garbage */
1236 buildval.val.high = sign << 31;
1237 buildval.val.high |= exponent << 20;
1238 buildval.val.high |= mantissa0;
1239 buildval.val.low = mantissa1;
1242 TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1246 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1248 const char *value = (const char*) val;
1250 int exp_offset, val_bias, res_bias;
1252 if (result == NULL) result = calc_buffer;
1253 temp = alloca(VALUE_SIZE);
1255 if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1257 if (value != result) memcpy(result, value, CALC_BUFFER_SIZE);
1261 /* set the descriptor of the new value */
1262 _desc(result).exponent_size = exp_size;
1263 _desc(result).mantissa_size = mant_size;
1264 _desc(result).class = _desc(value).class;
1266 _sign(result) = _sign(value);
1268 /* when the mantissa sizes differ normalizing has to shift to align it.
1269 * this would change the exponent, which is unwanted. So calculate this
1270 * offset and add it */
1271 val_bias = (1<<_desc(value).exponent_size)/2-1;
1272 res_bias = (1<<exp_size)/2-1;
1274 exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1275 sc_val_from_long(exp_offset, temp);
1276 sc_add(_exp(value), temp, _exp(result));
1278 /* _normalize expects normalized radix point */
1279 if (_desc(val).class == SUBNORMAL) {
1280 sc_val_from_ulong(1, NULL);
1281 _shift_left(_mant(val), sc_get_buffer(), _mant(result));
1282 } else if (value != result) {
1283 memcpy(_mant(result), _mant(value), VALUE_SIZE);
1285 memmove(_mant(result), _mant(value), VALUE_SIZE);
1288 _normalize(result, result, 0);
1289 TRACEPRINTF(("Cast results in %s\n", fc_print(result, temp, VALUE_SIZE, FC_PACKED)));
1293 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1295 if (result == NULL) result = calc_buffer;
1297 _desc(result).exponent_size = exponent_size;
1298 _desc(result).mantissa_size = mantissa_size;
1299 _desc(result).class = NORMAL;
1303 sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1305 sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1306 sc_val_from_ulong(2, NULL);
1307 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1312 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1314 if (result == NULL) result = calc_buffer;
1316 fc_get_max(exponent_size, mantissa_size, result);
1322 char* fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1324 if (result == NULL) result = calc_buffer;
1326 _desc(result).exponent_size = exponent_size;
1327 _desc(result).mantissa_size = mantissa_size;
1328 _desc(result).class = NAN;
1332 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1334 /* signalling nan has non-zero mantissa with msb not set */
1335 sc_val_from_ulong(1, _mant(result));
1340 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1342 if (result == NULL) result = calc_buffer;
1344 _desc(result).exponent_size = exponent_size;
1345 _desc(result).mantissa_size = mantissa_size;
1346 _desc(result).class = NAN;
1350 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1352 /* quiet nan has the msb of the mantissa set, so shift one there */
1353 sc_val_from_ulong(1, _mant(result));
1354 /* mantissa_size >+< 1 because of two extra rounding bits */
1355 sc_val_from_ulong(mantissa_size + 1, NULL);
1356 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1361 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1363 if (result == NULL) result = calc_buffer;
1365 _desc(result).exponent_size = exponent_size;
1366 _desc(result).mantissa_size = mantissa_size;
1367 _desc(result).class = NORMAL;
1371 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1373 sc_val_from_ulong(0, _mant(result));
1378 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1380 if (result == NULL) result = calc_buffer;
1382 fc_get_plusinf(exponent_size, mantissa_size, result);
1388 int fc_comp(const void *a, const void *b)
1390 const char *val_a = (const char*)a;
1391 const char *val_b = (const char*)b;
1394 if (_desc(val_a).class == NAN || _desc(val_b).class == NAN) return 2;
1395 /* zero is equal independent of sign */
1396 if ((_desc(val_a).class == ZERO) && (_desc(val_b).class == ZERO)) return 0;
1397 /* different signs make compare easy */
1398 if (_sign(val_a) != _sign(val_b)) return (_sign(val_a)==0)?(1):(-1);
1399 /* both infinity means equality */
1400 if ((_desc(val_a).class == INF) && (_desc(val_b).class == INF)) return 0;
1401 /* infinity is bigger than the rest */
1402 if (_desc(val_a).class == INF) return _sign(val_a)?(-1):(1);
1403 if (_desc(val_b).class == INF) return _sign(val_b)?(1):(-1);
1405 switch (sc_comp(_exp(val_a), _exp(val_b))) {
1411 return sc_comp(_mant(val_a), _mant(val_b));
1417 int fc_is_zero(const void *a)
1419 return _desc((const char*)a).class == ZERO;
1422 int fc_is_negative(const void *a)
1424 return _sign((const char*)a);
1427 int fc_is_inf(const void *a)
1429 return _desc(a).class == INF;
1432 int fc_is_nan(const void *a)
1434 return _desc(a).class == NAN;
1437 int fc_is_subnormal(const void *a)
1439 return _desc(a).class == SUBNORMAL;
1442 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1447 val = (const char*)a;
1449 mul_1 = alloca(CALC_BUFFER_SIZE);
1453 switch (_desc(val).class) {
1455 if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1456 else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1459 snprintf(buf, buflen, "NAN");
1462 snprintf(buf, buflen, "0.0");
1465 /* XXX to be implemented */
1466 #ifdef HAVE_LONG_DOUBLE
1467 /* XXX 30 is arbitrary */
1468 snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1470 snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1476 switch (_desc(val).class) {
1478 if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1479 else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1482 snprintf(buf, buflen, "NAN");
1485 snprintf(buf, buflen, "0.0");
1488 #ifdef HAVE_LONG_DOUBLE
1489 snprintf(buf, buflen, "%LA", fc_val_to_float(val));
1491 snprintf(buf, buflen, "%A", fc_val_to_float(val));
1498 snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), VALUE_SIZE*4, SC_HEX));
1504 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1506 /* this is used to cache the packed version of the value */
1507 static char *pack = NULL;
1509 if (pack == NULL) pack = malloc(VALUE_SIZE);
1512 _pack((const char*)value, pack);
1514 return sc_sub_bits(pack, num_bits, byte_ofs);
1517 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1519 if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1520 ROUNDING_MODE = mode;
1522 return ROUNDING_MODE;
1525 fc_rounding_mode_t fc_get_rounding_mode(void)
1527 return ROUNDING_MODE;
1530 void init_fltcalc(int precision)
1532 if (calc_buffer == NULL) {
1533 /* does nothing if already init */
1534 if (precision == 0) precision = FC_DEFAULT_PRECISION;
1536 init_strcalc(precision + 4);
1538 /* needs additionally two bits to round, a bit as explicit 1., and one for
1539 * addition overflow */
1540 max_precision = sc_get_precision() - 4;
1541 if (max_precision < precision)
1542 printf("WARING: not enough precision available, using %d\n", max_precision);
1544 ROUNDING_MODE = FC_TONEAREST;
1545 VALUE_SIZE = sc_get_buffer_length();
1547 EXPONENT_POS = SIGN_POS + sizeof(char);
1548 MANTISSA_POS = EXPONENT_POS + VALUE_SIZE;
1549 DESCRIPTOR_POS = MANTISSA_POS + VALUE_SIZE;
1550 CALC_BUFFER_SIZE = DESCRIPTOR_POS + sizeof(descriptor_t);
1552 calc_buffer = malloc(CALC_BUFFER_SIZE);
1553 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));
1554 #ifdef HAVE_LONG_DOUBLE
1555 DEBUGPRINTF(("\tUsing long double (1-15-64) interface\n"));
1557 DEBUGPRINTF(("\tUsing double (1-11-52) interface\n"));
1559 #ifdef WORDS_BIGENDIAN
1560 DEBUGPRINTF(("\tWord order is big endian\n\n"));
1562 DEBUGPRINTF(("\tWord order is little endian\n\n"));
1567 void finish_fltcalc (void) {
1568 free(calc_buffer); calc_buffer = NULL;
1571 /* definition of interface functions */