2 * Authors: Matthias Heil
8 #include <math.h> /* need isnan() and isinf() (will be changed)*/
9 /* undef some reused constants defined by math.h */
21 typedef uint32_t UINT32;
22 #ifdef WORD_LITTLE_ENDIAN
23 #undef WORD_LITTLE_ENDIAN
25 #ifdef WORD_BIG_ENDIAN
26 #undef WORD_BIG_ENDIAN
28 #define WORD_LITTLE_ENDIAN
30 #ifdef HAVE_LONG_DOUBLE
31 #ifdef WORD_LITTLE_ENDIAN
38 volatile long double d;
47 volatile long double d;
51 #ifdef WORD_LITTLE_ENDIAN
88 #define CLEAR_BUFFER(buffer) memset(buffer, 0, CALC_BUFFER_SIZE)
90 /* because variable sized structs are impossible, the internal
91 * value is represented as a pseudo-struct char array, addressed
94 * char sign; // 0 for positive, 1 for negative
95 * char exp[VALUE_SIZE];
96 * char mant[VALUE_SIZE];
100 #define _sign(a) (((char*)a)[SIGN_POS])
101 #define _exp(a) (&((char*)a)[EXPONENT_POS])
102 #define _mant(a) (&((char*)a)[MANTISSA_POS])
103 #define _desc(a) (*(descriptor_t *)&((char*)a)[DESCRIPTOR_POS])
105 #define _save_result(x) memcpy((x), sc_get_buffer(), VALUE_SIZE)
106 #define _shift_right(x, y, b) sc_shr((x), (y), VALUE_SIZE*4, 0, (b))
107 #define _shift_left(x, y, b) sc_shl((x), (y), VALUE_SIZE*4, 0, (b))
109 #define FC_DEFINE1(code) char* fc_##code(const void *a, void *result) \
111 return _calc((const char*)a, NULL, FC_##code, (char*)result); \
114 #define FC_DEFINE2(code) char* fc_##code(const void *a, const void *b, void *result) \
116 return _calc((const char*)a, (const char*)b, FC_##code, (char*)result); \
119 #define FUNC_PTR(code) fc_##code
122 # define DEBUGPRINTF(x) printf x
124 # define DEBUGPRINTF(x) ((void)0)
127 #if FLTCALC_TRACE_CALC
128 # define TRACEPRINTF(x) printf x
130 # define TRACEPRINTF(x) ((void)0)
133 static char *calc_buffer = NULL;
135 static fc_rounding_mode_t ROUNDING_MODE;
137 static int CALC_BUFFER_SIZE;
138 static int VALUE_SIZE;
140 static int EXPONENT_POS;
141 static int MANTISSA_POS;
142 static int DESCRIPTOR_POS;
144 static int max_precision;
149 static void _fail_char(const char *str, unsigned int len, int pos)
152 printf("ERROR: Unexpected character '%c'\n", *(str + pos));
154 printf("ERROR: Unexpected end of string\n");
155 while (len-- && *str) printf("%c", *str++); printf("\n");
156 while (pos--) printf(" "); printf("^\n");
157 /* the front end has to to check constant strings */
163 static char* _pack(const char *int_float, char *packed)
169 temp = alloca(VALUE_SIZE);
170 shift_val = alloca(VALUE_SIZE);
172 switch (_desc(int_float).class) {
174 val_buffer = alloca(CALC_BUFFER_SIZE);
175 fc_get_qnan(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
176 int_float = val_buffer;
180 val_buffer = alloca(CALC_BUFFER_SIZE);
181 fc_get_plusinf(_desc(int_float).exponent_size, _desc(int_float).mantissa_size, val_buffer);
182 _sign(val_buffer) = _sign(int_float);
183 int_float = val_buffer;
190 sc_val_from_ulong(_sign(int_float), temp);
192 sc_val_from_ulong(_desc(int_float).exponent_size + _desc(int_float).mantissa_size, NULL);
193 _shift_left(temp, sc_get_buffer(), packed);
195 /* extract exponent */
196 sc_val_from_ulong(_desc(int_float).mantissa_size, shift_val);
198 _shift_left(_exp(int_float), shift_val, temp);
200 sc_or(temp, packed, packed);
202 /* extract mantissa */
203 /* remove 2 rounding bits */
204 sc_val_from_ulong(2, shift_val);
205 _shift_right(_mant(int_float), shift_val, temp);
207 /* remove leading 1 (or 0 if denormalized) */
208 sc_max_from_bits(_desc(int_float).mantissa_size, 0, shift_val); /* all mantissa bits are 1's */
209 sc_and(temp, shift_val, temp);
212 sc_or(temp, packed, packed);
217 char* _normalize(const char *in_val, char *out_val, int sticky)
220 char lsb, guard, round, round_dir = 0;
223 temp = alloca(VALUE_SIZE);
225 /* +2: save two rounding bits at the end */
226 hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
228 if (in_val != out_val)
230 _sign(out_val) = _sign(in_val);
231 memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
234 _desc(out_val).class = NORMAL;
236 /* mantissa all zeroes, so zero exponent (because of explicit one)*/
237 if (hsb == 2 + _desc(in_val).mantissa_size)
239 sc_val_from_ulong(0, _exp(out_val));
243 /* shift the first 1 ito 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));
259 sc_val_from_ulong(hsb+1, temp);
261 _shift_left(_mant(in_val), temp, _mant(out_val));
263 sc_sub(_exp(in_val), temp, _exp(out_val));
266 /* check for exponent underflow */
267 if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
268 DEBUGPRINTF(("Exponent underflow!\n"));
269 /* exponent underflow */
270 /* shift the mantissa right to have a zero exponent */
271 sc_val_from_ulong(1, temp);
272 sc_sub(temp, _exp(out_val), NULL);
274 _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
275 if (!sticky) sticky = sc_had_carry();
276 /* denormalized means exponent of zero */
277 sc_val_from_ulong(0, _exp(out_val));
279 _desc(out_val).class = SUBNORMAL;
282 /* perform rounding by adding a value that clears the guard bit and the round bit
283 * and either causes a carry to round up or not */
284 /* get the last 3 bits of the value */
285 lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
286 guard = (lsb&0x2)>>1;
289 switch (ROUNDING_MODE)
292 /* round to nearest representable value, if in doubt choose the version
294 round_dir = guard && (sticky || round || lsb>>2);
297 /* if positive: round to one if the exact value is bigger, else to zero */
298 round_dir = (!_sign(out_val) && (guard || round || sticky));
301 /* if negative: round to one if the exact value is bigger, else to zero */
302 round_dir = (_sign(out_val) && (guard || round || sticky));
305 /* always round to 0 (chopping mode) */
309 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"));
313 guard = (round^guard)<<1;
314 lsb = !(round || guard)<<2 | guard | round;
318 lsb = -((guard<<1) | round);
321 /* add the rounded value */
323 sc_val_from_long(lsb, temp);
324 sc_add(_mant(out_val), temp, _mant(out_val));
327 /* could have rounded down to zero */
328 if (sc_is_zero(_mant(out_val)) && (_desc(out_val).class == SUBNORMAL))
329 _desc(out_val).class = ZERO;
331 /* check for rounding overflow */
332 hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
333 if ((_desc(out_val).class != SUBNORMAL) && (hsb < -1))
335 sc_val_from_ulong(1, temp);
336 _shift_right(_mant(out_val), temp, _mant(out_val));
338 sc_add(_exp(out_val), temp, _exp(out_val));
340 else if ((_desc(out_val).class == SUBNORMAL) && (hsb == -1))
342 /* overflow caused the matissa to be normal again,
343 * so adapt the exponent accordingly */
344 sc_val_from_ulong(1, temp);
345 sc_add(_exp(out_val), temp, _exp(out_val));
347 _desc(out_val).class = NORMAL;
349 /* no further rounding is needed, because rounding overflow means
350 * the carry of the original rounding was propagated all the way
351 * up to the bit left of the radix point. This implies the bits
352 * to the right are all zeros (rounding is +1) */
354 /* check for exponent overflow */
355 sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
356 if (sc_comp(_exp(out_val), temp) != -1) {
357 DEBUGPRINTF(("Exponent overflow!\n"));
358 /* exponent overflow, reaction depends on rounding method:
360 * mode | sign of value | result
361 *--------------------------------------------------------------
362 * TO_NEAREST | + | +inf
364 *--------------------------------------------------------------
365 * TO_POSITIVE | + | +inf
366 * | - | smallest representable value
367 *--------------------------------------------------------------
368 * TO_NEAGTIVE | + | largest representable value
370 *--------------------------------------------------------------
371 * TO_ZERO | + | largest representable value
372 * | - | smallest representable value
373 *--------------------------------------------------------------*/
374 if (_sign(out_val) == 0)
376 /* value is positive */
377 switch (ROUNDING_MODE) {
380 _desc(out_val).class = INF;
385 fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
388 /* value is negative */
389 switch (ROUNDING_MODE) {
392 _desc(out_val).class = INF;
397 fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
406 * calculate a + b, where a is the value with the bigger exponent
408 static char* _add(const char* a, const char* b, char* result)
416 if (_desc(a).class == NAN) {
417 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
420 if (_desc(b).class == NAN) {
421 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
425 /* make sure result has a descriptor */
426 if (result != a && result != b)
427 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
429 /* determine if this is an addition or subtraction */
430 sign = _sign(a) ^ _sign(b);
431 DEBUGPRINTF(("sign a: %d, sign b: %d -> %s\n", _sign(a), _sign(b), sign?"sub":"add"));
433 /* produce nan on inf - inf */
434 if (sign && (_desc(a).class == INF) && (_desc(b).class == INF))
435 return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);
437 temp = alloca(VALUE_SIZE);
438 exp_diff = alloca(VALUE_SIZE);
440 /* get exponent difference */
441 sc_sub(_exp(a), _exp(b), exp_diff);
443 /* initially set sign to be the sign of a, special treatment of subtraction
444 * when exponents are equal is required though.
445 * Also special care about the sign is needed when the mantissas are equal
447 if (sign && sc_val_to_long(exp_diff) == 0) {
448 switch (sc_comp(_mant(a), _mant(b))) {
450 if (_sign(a)) _sign(result) = 1; /* abs(a) is bigger and a is negative */
451 else _sign(result) = 0;
454 if (ROUNDING_MODE == FC_TONEGATIVE)
460 if (_sign(b)) _sign(result) = 1; /* abs(b) is bigger and b is negative */
461 else _sign(result) = 0;
464 /* can't be reached */
468 _sign(result) = _sign(a);
471 /* sign has been taken care of, check for special cases */
472 if (_desc(a).class == ZERO) {
473 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
476 if (_desc(b).class == ZERO) {
477 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
481 if (_desc(a).class == INF) {
482 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
485 if (_desc(b).class == INF) {
486 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-SIGN_POS-1);
490 /* shift the smaller value to the right to align the radix point */
491 /* subnormals have their radix point shifted to the right,
492 * take care of this first */
493 if ((_desc(b).class == SUBNORMAL) && (_desc(a).class != SUBNORMAL))
495 sc_val_from_ulong(1, temp);
496 sc_sub(exp_diff, temp, exp_diff);
499 _shift_right(_mant(b), exp_diff, temp);
500 sticky = sc_had_carry();
504 /* if subtracting a little more than the represented value or adding a little
505 * more than the represented value to a negative value this, in addition to the
506 * still set sticky bit, takes account of the 'little more' */
507 char *temp1 = alloca(CALC_BUFFER_SIZE);
508 sc_val_from_ulong(1, temp1);
509 sc_add(temp, temp1, temp);
513 if (sc_comp(_mant(a), temp) == -1)
514 sc_sub(temp, _mant(a), _mant(result));
516 sc_sub(_mant(a), temp, _mant(result));
518 sc_add(_mant(a), temp, _mant(result));
521 /* _normalize expects a 'normal' radix point, adding two subnormals
522 * results in a subnormal radix point -> shifting before normalizing */
523 if ((_desc(a).class == SUBNORMAL) && (_desc(b).class == SUBNORMAL))
525 sc_val_from_ulong(1, NULL);
526 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
529 /* resulting exponent is the bigger one */
530 memmove(_exp(result), _exp(a), VALUE_SIZE);
532 return _normalize(result, result, sticky);
535 static char* _mul(const char* a, const char* b, char* result)
539 if (_desc(a).class == NAN) {
540 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
543 if (_desc(b).class == NAN) {
544 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
548 temp = alloca(VALUE_SIZE);
550 if (result != a && result != b)
551 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
553 _sign(result) = _sign(a) ^ _sign(b);
555 /* produce nan on 0 * inf */
556 if (_desc(a).class == ZERO) {
557 if (_desc(b).class == INF)
558 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
560 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
563 if (_desc(b).class == ZERO) {
564 if (_desc(a).class == INF)
565 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
567 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1);
571 if (_desc(a).class == INF) {
572 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
575 if (_desc(b).class == INF) {
576 if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, CALC_BUFFER_SIZE-1);
580 /* exp = exp(a) + exp(b) - excess */
581 sc_add(_exp(a), _exp(b), _exp(result));
583 sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
584 sc_sub(_exp(result), temp, _exp(result));
586 /* mixed normal, subnormal values introduce an error of 1, correct it */
587 if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL))
589 sc_val_from_ulong(1, temp);
590 sc_add(_exp(result), temp, _exp(result));
593 sc_mul(_mant(a), _mant(b), _mant(result));
595 /* realign result: after a multiplication the digits right of the radix
596 * point are the sum of the factors' digits after the radix point. As all
597 * values are normalized they both have the same amount of these digits,
598 * which has to be restored by proper shifting
599 * +2 because of the two rounding bits */
600 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
602 _shift_right(_mant(result), temp, _mant(result));
604 return _normalize(result, result, sc_had_carry());
607 static char* _div(const char* a, const char* b, char* result)
609 char *temp, *dividend;
611 if (_desc(a).class == NAN) {
612 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
615 if (_desc(b).class == NAN) {
616 if (b != result) memcpy(result, b, CALC_BUFFER_SIZE);
620 temp = alloca(VALUE_SIZE);
621 dividend = alloca(VALUE_SIZE);
623 if (result != a && result != b)
624 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
626 _sign(result) = _sign(a) ^ _sign(b);
628 /* produce nan on 0/0 and inf/inf */
629 if (_desc(a).class == ZERO) {
630 if (_desc(b).class == ZERO)
632 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
635 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
639 if (_desc(b).class == INF) {
640 if (_desc(a).class == INF)
642 fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
645 sc_val_from_ulong(0, NULL);
646 _save_result(_exp(result));
647 _save_result(_mant(result));
648 _desc(result).class = ZERO;
653 if (_desc(a).class == INF) {
655 if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, CALC_BUFFER_SIZE-1);
658 if (_desc(b).class == ZERO) {
659 /* division by zero */
661 fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
663 fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
667 /* exp = exp(a) - exp(b) + excess - 1*/
668 sc_sub(_exp(a), _exp(b), _exp(result));
669 sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
670 sc_add(_exp(result), temp, _exp(result));
672 /* mixed normal, subnormal values introduce an error of 1, correct it */
673 if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL))
675 sc_val_from_ulong(1, temp);
676 sc_add(_exp(result), temp, _exp(result));
679 /* mant(res) = mant(a) / 1/2mant(b) */
680 /* to gain more bits of precision in the result the dividend could be
681 * shifted left, as this operation does not loose bits. This would not
682 * fit into the integer precision, but due to the rounding bits (which
683 * are always zero because the values are all normalized) the divisor
684 * can be shifted right instead to achieve the same result */
685 sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
687 _shift_left(_mant(a), temp, dividend);
690 char *divisor = alloca(CALC_BUFFER_SIZE);
691 sc_val_from_ulong(1, divisor);
692 _shift_right(_mant(b), divisor, divisor);
693 sc_div(dividend, divisor, _mant(result));
696 return _normalize(result, result, sc_had_carry());
699 void _power_of_ten(int exp, descriptor_t *desc, char *result)
707 /* set new descriptor (else result is supposed to already have one) */
709 memcpy(&_desc(result), desc, sizeof(descriptor_t));
711 build = alloca(VALUE_SIZE);
712 temp = alloca(VALUE_SIZE);
714 sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));
718 /* temp is value of ten now */
719 sc_val_from_ulong(10, NULL);
722 for (exp--; exp > 0; exp--) {
724 sc_mul(build, temp, NULL);
728 /* temp is amount of leftshift needed to put the value left of the radix point */
729 sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
731 _shift_left(build, temp, _mant(result));
733 _normalize(result, result, 0);
737 static char* _trunc(const char *a, char *result)
739 /* when exponent == 0 all bits left of the radix point
740 * are the integral part of the value. For 15bit exp_size
741 * this would require a leftshift of max. 16383 bits which
743 * But it is enough to ensure that no bit right of the radix
744 * point remains set. This restricts the interesting
745 * exponents to the interval [0, mant_size-1].
746 * Outside this interval the truncated value is either 0 or
747 * it is are already truncated */
749 int exp_bias, exp_val;
752 temp = alloca(VALUE_SIZE);
755 memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
757 exp_bias = (1<<_desc(a).exponent_size)/2-1;
758 exp_val = sc_val_to_long(_exp(a)) - exp_bias;
761 sc_val_from_ulong(0, NULL);
762 _save_result(_exp(result));
763 _save_result(_mant(result));
764 _desc(result).class = ZERO;
769 if (exp_val > _desc(a).mantissa_size) {
771 memcpy(result, a, CALC_BUFFER_SIZE);
776 /* set up a proper mask to delete all bits right of the
777 * radix point if the mantissa had been shifted until exp == 0 */
778 sc_max_from_bits(1 + exp_val, 0, temp);
779 sc_val_from_long(_desc(a).mantissa_size - exp_val + 2, NULL);
780 _shift_left(temp, sc_get_buffer(), temp);
782 /* and the mask and return the result */
783 sc_and(_mant(a), temp, _mant(result));
785 if (a != result) memcpy(_exp(result), _exp(a), VALUE_SIZE);
791 * This does value sanity checking(or should do it), sets up any prerequisites,
792 * calls the proper internal functions, clears up and returns
794 char* _calc(const char *a, const char *b, int opcode, char *result)
797 #ifdef FLTCALC_TRACE_CALC
800 buffer = alloca(100);
803 if (result == NULL) result = calc_buffer;
805 TRACEPRINTF(("%s ", fc_print(a, buffer, 100, FC_HEX)));
809 /* make the value with the bigger exponent the first one */
810 TRACEPRINTF(("+ %s ", fc_print(b, buffer, 100, FC_HEX)));
811 if (sc_comp(_exp(a), _exp(b)) == -1)
817 TRACEPRINTF(("- %s ", fc_print(b, buffer, 100, FC_HEX)));
818 temp = alloca(CALC_BUFFER_SIZE);
819 memcpy(temp, b, CALC_BUFFER_SIZE);
820 _sign(temp) = !_sign(b);
821 if (sc_comp(_exp(a), _exp(temp)) == -1)
822 _add(temp, a, result);
824 _add(a, temp, result);
827 TRACEPRINTF(("* %s ", fc_print(b, buffer, 100, FC_HEX)));
831 TRACEPRINTF(("/ %s ", fc_print(b, buffer, 100, FC_HEX)));
835 TRACEPRINTF(("negated "));
836 if (a != result) memcpy(result, a, CALC_BUFFER_SIZE);
837 _sign(result) = !_sign(a);
846 TRACEPRINTF(("= %s\n", fc_print(result, buffer, 100, FC_HEX)));
851 * functions defined in fltcalc.h
853 const void *fc_get_buffer(void)
858 const int fc_get_buffer_length(void)
860 return CALC_BUFFER_SIZE;
863 char* fc_val_from_str(const char *str, unsigned int len, char exp_size, char mant_size, char *result)
876 int exp_int, hsb, state;
881 char *mant_str, *exp_val, *power_val;
883 if (result == NULL) result = calc_buffer;
885 exp_val = alloca(VALUE_SIZE);
886 power_val = alloca(CALC_BUFFER_SIZE);
887 mant_str = alloca((len)?(len):(strlen(str)));
889 _desc(result).exponent_size = exp_size;
890 _desc(result).mantissa_size = mant_size;
891 _desc(result).class = NORMAL;
898 while (len == 0 || str-old_str < len)
915 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
922 state = RIGHT_OF_DOT;
933 _fail_char(old_str, len, str - old_str);
939 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
940 mant_str[pos++] = *(str++);
944 state = RIGHT_OF_DOT;
955 mant_str[pos] = '\0';
959 _fail_char(old_str, len, str - old_str);
965 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
966 mant_str[pos++] = *(str++);
977 mant_str[pos] = '\0';
981 _fail_char(old_str, len, str - old_str);
991 if (*(str-1) != 'e' && *(str-1) != 'E') _fail_char(old_str, len, str - old_str);
995 case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9':
996 mant_str[pos] = '\0';
1003 _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':
1014 case '\0': goto done;
1017 _fail_char(old_str, len, str - old_str);
1023 sc_val_from_str(mant_str, strlen(mant_str), _mant(result));
1025 /* shift to put value left of radix point */
1026 sc_val_from_ulong(mant_size + 2, exp_val);
1028 _shift_left(_mant(result), exp_val, _mant(result));
1030 sc_val_from_ulong((1 << exp_size)/2-1, _exp(result));
1032 _normalize(result, result, 0);
1034 if (state == EXPONENT) {
1035 exp_int -= atoi(str-pos);
1038 _power_of_ten(exp_int, &_desc(result), power_val);
1040 _div(result, power_val, result);
1045 /* XXX excuse of an implementation to make things work */
1047 #ifdef HAVE_LONG_DOUBLE
1048 val = strtold(str, NULL);
1050 val = strtod(str, NULL);
1053 DEBUGPRINTF(("val_from_str(%s)\n", str));
1054 return fc_val_from_float(val, exp_size, mant_size, result);
1058 char* fc_val_from_float(LLDBL l, char exp_size, char mant_size, char* result)
1061 int bias_res = ((1<<exp_size)/2-1);
1068 #ifdef HAVE_LONG_DOUBLE
1071 UINT32 sign = (srcval.val.high & 0x80000000) != 0;
1072 UINT32 exponent = (srcval.val.high & 0x7FFF0000) >> 16;
1073 UINT32 mantissa1 = srcval.val.mid;
1074 UINT32 mantissa0 = srcval.val.low;
1075 #else /* no long double */
1078 UINT32 sign = (srcval.val.high & 0x80000000) != 0;
1079 UINT32 exponent = (srcval.val.high & 0x7FF00000) >> 20;
1080 UINT32 mantissa0 = srcval.val.high & 0x000FFFFF;
1081 UINT32 mantissa1 = srcval.val.low;
1084 #ifdef HAVE_LONG_DOUBLE
1085 TRACEPRINTF(("val_from_float(%.8X%.8X%.8X)\n", ((int*)&l)[2], ((int*)&l)[1], ((int*)&l)[0]));
1086 DEBUGPRINTF(("(%d-%.4X-%.8X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1088 TRACEPRINTF(("val_from_float(%.8X%.8X)\n", ((int*)&l)[1], ((int*)&l)[0]));
1089 DEBUGPRINTF(("(%d-%.3X-%.5X%.8X)\n", sign, exponent, mantissa0, mantissa1));
1092 if (result == NULL) result = calc_buffer;
1093 temp = alloca(VALUE_SIZE);
1095 _desc(result).exponent_size = exp_size;
1096 _desc(result).mantissa_size = mant_size;
1099 _sign(result) = sign;
1101 /* sign and flag suffice to identify nan or inf, no exponent/mantissa
1102 * encoding is needed. the function can return immediately in these cases */
1104 _desc(result).class = NAN;
1105 TRACEPRINTF(("val_from_float resulted in NAN\n"));
1108 else if (isinf(l)) {
1109 _desc(result).class = INF;
1110 TRACEPRINTF(("val_from_float resulted in %sINF\n", (_sign(result)==1)?"-":""));
1114 /* extract exponent */
1115 sc_val_from_long((exponent-bias_val+bias_res)-(mant_val-mant_size), _exp(result));
1117 /* extract mantissa */
1120 /* insert the hidden bit */
1121 sc_val_from_ulong(1, temp);
1122 sc_val_from_ulong(mant_val + 2, NULL);
1123 _shift_left(temp, sc_get_buffer(), NULL);
1127 sc_val_from_ulong(0, NULL);
1130 _save_result(_mant(result));
1132 /* bits from the upper word */
1133 sc_val_from_ulong(mantissa0, temp);
1134 sc_val_from_ulong(34, NULL);
1135 _shift_left(temp, sc_get_buffer(), temp);
1136 sc_or(_mant(result), temp, _mant(result));
1138 /* bits from the lower word */
1139 sc_val_from_ulong(mantissa1, temp);
1140 sc_val_from_ulong(2, NULL);
1141 _shift_left(temp, sc_get_buffer(), temp);
1142 sc_or(_mant(result), temp, _mant(result));
1144 /* _normalize expects the radix point to be normal, so shift mantissa of subnormal
1145 * origin one to the left */
1148 sc_val_from_ulong(1, NULL);
1149 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1152 _normalize(result, result, 0);
1154 TRACEPRINTF(("val_from_float results in %s\n", fc_print(result, temp, CALC_BUFFER_SIZE, FC_HEX)));
1159 LLDBL fc_val_to_float(const void *val)
1174 #ifdef HAVE_LONG_DOUBLE
1175 char result_mantissa = 64;
1176 char result_exponent = 15;
1178 char result_mantissa = 52;
1179 char result_exponent = 11;
1182 temp = alloca(CALC_BUFFER_SIZE);
1183 pack = alloca(VALUE_SIZE);
1184 value = fc_cast(val, result_exponent, result_mantissa, temp);
1186 sign = _sign(value);
1187 /* long double exponent is 15bit, so the use of sc_val_to_long should not
1188 * lead to wrong results */
1189 exponent = sc_val_to_long(_exp(value)) ;
1196 for (byte_offset = 0; byte_offset < 4; byte_offset++)
1197 mantissa1 |= sc_sub_bits(pack, result_mantissa, byte_offset) << (byte_offset<<3);
1199 for (; (byte_offset<<3) < result_mantissa; byte_offset++)
1200 mantissa0 |= sc_sub_bits(pack, result_mantissa, byte_offset) << ((byte_offset-4)<<3);
1202 #ifndef HAVE_LONG_DOUBLE
1203 mantissa0 &= 0x000FFFFF;
1206 #ifdef HAVE_LONG_DOUBLE
1207 buildval.val.high = sign << 31;
1208 buildval.val.high |= exponent << 16;
1209 buildval.val.mid = mantissa1;
1210 buildval.val.low = mantissa0;
1211 #else /* no long double */
1212 buildval.val.high = sign << 31;
1213 buildval.val.high |= exponent << 20;
1214 buildval.val.high |= mantissa0;
1215 buildval.val.low = mantissa1;
1218 TRACEPRINTF(("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1));
1219 //printf("val_to_float: %d-%x-%x%x\n", sign, exponent, mantissa0, mantissa1);
1223 char* fc_cast(const void *val, char exp_size, char mant_size, char *result)
1225 const char *value = (const char*) val;
1227 int exp_offset, val_bias, res_bias;
1229 if (result == NULL) result = calc_buffer;
1230 temp = alloca(CALC_BUFFER_SIZE);
1232 if (_desc(value).exponent_size == exp_size && _desc(value).mantissa_size == mant_size)
1234 if (value != result) memcpy(result, value, CALC_BUFFER_SIZE);
1238 /* set the descriptor of the new value */
1239 _desc(result).exponent_size = exp_size;
1240 _desc(result).mantissa_size = mant_size;
1241 _desc(result).class = _desc(value).class;
1243 _sign(result) = _sign(value);
1245 /* when the mantissa sizes differ normalizing has to shift to align it.
1246 * this would change the exponent, which is unwanted. So calculate this
1247 * offset and add it */
1248 val_bias = (1<<_desc(value).exponent_size)/2-1;
1249 res_bias = (1<<exp_size)/2-1;
1250 exp_offset = (res_bias - val_bias) - (_desc(value).mantissa_size - mant_size);
1252 sc_val_from_long(exp_offset, temp);
1253 sc_add(_exp(value), temp, _exp(result));
1255 if (value != result) memcpy(_mant(result), _mant(value), VALUE_SIZE);
1256 else memmove(_mant(result), _mant(value), VALUE_SIZE);
1258 return _normalize(result, result, 0);
1261 char* fc_get_max(unsigned int exponent_size, unsigned int mantissa_size, char* result)
1263 if (result == NULL) result = calc_buffer;
1265 _desc(result).exponent_size = exponent_size;
1266 _desc(result).mantissa_size = mantissa_size;
1267 _desc(result).class = NORMAL;
1271 sc_val_from_ulong((1<<exponent_size) - 2, _exp(result));
1273 sc_max_from_bits(mantissa_size + 1, 0, _mant(result));
1274 sc_val_from_ulong(2, NULL);
1275 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1280 char* fc_get_min(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1282 if (result == NULL) result = calc_buffer;
1284 fc_get_max(exponent_size, mantissa_size, result);
1290 char* fc_get_snan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1292 if (result == NULL) result = calc_buffer;
1294 _desc(result).exponent_size = exponent_size;
1295 _desc(result).mantissa_size = mantissa_size;
1296 _desc(result).class = NAN;
1300 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1302 /* signalling nan has non-zero mantissa with msb not set */
1303 sc_val_from_ulong(1, _mant(result));
1308 char* fc_get_qnan(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1310 if (result == NULL) result = calc_buffer;
1312 _desc(result).exponent_size = exponent_size;
1313 _desc(result).mantissa_size = mantissa_size;
1314 _desc(result).class = NAN;
1318 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1320 /* quiet nan has the msb of the mantissa set, so shift one there */
1321 sc_val_from_ulong(1, _mant(result));
1322 /* mantissa_size >+< 1 because of two extra rounding bits */
1323 sc_val_from_ulong(mantissa_size + 1, NULL);
1324 _shift_left(_mant(result), sc_get_buffer(), _mant(result));
1329 char* fc_get_plusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1331 if (result == NULL) result = calc_buffer;
1333 _desc(result).exponent_size = exponent_size;
1334 _desc(result).mantissa_size = mantissa_size;
1335 _desc(result).class = NORMAL;
1339 sc_val_from_ulong((1<<exponent_size)-1, _exp(result));
1341 sc_val_from_ulong(0, _mant(result));
1346 char* fc_get_minusinf(unsigned int exponent_size, unsigned int mantissa_size, char *result)
1348 if (result == NULL) result = calc_buffer;
1350 fc_get_plusinf(exponent_size, mantissa_size, result);
1356 int fc_comp(const void *a, const void *b)
1358 const char *val_a = (const char*)a;
1359 const char *val_b = (const char*)b;
1362 if (_desc(val_a).class == NAN || _desc(val_b).class == NAN) return 2;
1363 /* zero is equal independent of sign */
1364 if ((_desc(val_a).class == ZERO) && (_desc(val_b).class == ZERO)) return 0;
1365 /* different signs make compare easy */
1366 if (_sign(val_a) != _sign(val_b)) return (_sign(val_a)==0)?(1):(-1);
1367 /* both infinity means equality */
1368 if ((_desc(val_a).class == INF) && (_desc(val_b).class == INF)) return 0;
1369 /* infinity is bigger than the rest */
1370 if (_desc(val_a).class == INF) return _sign(val_a)?(-1):(1);
1371 if (_desc(val_b).class == INF) return _sign(val_b)?(1):(-1);
1373 switch (sc_comp(_exp(val_a), _exp(val_b))) {
1379 return sc_comp(_mant(val_a), _mant(val_b));
1385 int fc_is_zero(const void *a)
1387 return _desc((const char*)a).class == ZERO;
1390 int fc_is_negative(const void *a)
1392 return _sign((const char*)a);
1395 int fc_is_inf(const void *a)
1397 return _desc(a).class == INF;
1400 int fc_is_nan(const void *a)
1402 return _desc(a).class == NAN;
1405 int fc_is_subnormal(const void *a)
1407 return _desc(a).class == SUBNORMAL;
1410 char *fc_print(const void *a, char *buf, int buflen, unsigned base)
1415 val = (const char*)a;
1417 mul_1 = alloca(CALC_BUFFER_SIZE);
1421 switch (_desc(val).class) {
1423 if (buflen >= 8+_sign(val)) sprintf(buf, "%sINFINITY", _sign(val)?"-":"");
1424 else snprintf(buf, buflen, "%sINF", _sign(val)?"-":NULL);
1427 snprintf(buf, buflen, "NAN");
1430 snprintf(buf, buflen, "0.0");
1433 /* XXX to be implemented */
1434 #ifdef HAVE_LONG_DOUBLE
1435 /* XXX 30 is arbitrary */
1436 snprintf(buf, buflen, "%.30LE", fc_val_to_float(val));
1438 snprintf(buf, buflen, "%.18E", fc_val_to_float(val));
1443 snprintf(buf, buflen, "%s", sc_print(_pack(val, mul_1), 0, SC_HEX));
1449 unsigned char fc_sub_bits(const void *value, unsigned num_bits, unsigned byte_ofs)
1451 /* this is used to cache the packed version of the value */
1452 static char *pack = NULL;
1454 if (pack == NULL) pack = malloc(VALUE_SIZE);
1457 _pack((const char*)value, pack);
1459 return sc_sub_bits(pack, num_bits, byte_ofs);
1462 fc_rounding_mode_t fc_set_rounding_mode(fc_rounding_mode_t mode)
1464 if (mode == FC_TONEAREST || mode == FC_TOPOSITIVE || mode == FC_TONEGATIVE || mode == FC_TOZERO)
1465 ROUNDING_MODE = mode;
1467 return ROUNDING_MODE;
1470 fc_rounding_mode_t fc_get_rounding_mode(void)
1472 return ROUNDING_MODE;
1475 void init_fltcalc(int precision)
1477 if (calc_buffer == NULL) {
1478 /* does nothing if already init */
1479 if (precision == 0) precision = FC_DEFAULT_PRECISION;
1481 init_strcalc(precision + 4);
1483 /* needs additionally two bits to round, a bit as explicit 1., and one for
1484 * addition overflow */
1485 max_precision = sc_get_precision() - 4;
1486 if (max_precision < precision)
1487 printf("WARING: not enough precision available, using %d\n", max_precision);
1489 ROUNDING_MODE = FC_TONEAREST;
1490 VALUE_SIZE = sc_get_buffer_length();
1492 EXPONENT_POS = SIGN_POS + sizeof(char);
1493 MANTISSA_POS = EXPONENT_POS + VALUE_SIZE;
1494 DESCRIPTOR_POS = MANTISSA_POS + VALUE_SIZE;
1495 CALC_BUFFER_SIZE = DESCRIPTOR_POS + sizeof(descriptor_t);
1497 calc_buffer = malloc(CALC_BUFFER_SIZE);
1498 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));
1502 /* definition of interface functions */