-static char* _add(const char* a, const char* b, char* result)
-{
- char *temp;
- char *exp_diff;
-
- char sign;
- char sticky;
-
- if (_desc(a).class == NAN) {
- if (a != result) memcpy(result, a, calc_buffer_size);
- return result;
- }
- if (_desc(b).class == NAN) {
- if (b != result) memcpy(result, b, calc_buffer_size);
- return result;
- }
-
- /* make sure result has a descriptor */
- if (result != a && result != b)
- memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
-
- /* determine if this is an addition or subtraction */
- sign = _sign(a) ^ _sign(b);
-
- /* produce nan on inf - inf */
- if (sign && (_desc(a).class == INF) && (_desc(b).class == INF))
- return fc_get_qnan(_desc(a).exponent_size, _desc(b).mantissa_size, result);
-
- temp = alloca(value_size);
- exp_diff = alloca(value_size);
-
- /* get exponent difference */
- sc_sub(_exp(a), _exp(b), exp_diff);
-
- /* initially set sign to be the sign of a, special treatment of subtraction
- * when exponents are equal is required though.
- * Also special care about the sign is needed when the mantissas are equal
- * (+/- 0 ?) */
- if (sign && sc_val_to_long(exp_diff) == 0) {
- switch (sc_comp(_mant(a), _mant(b))) {
- case 1: /* a > b */
- if (_sign(a)) _sign(result) = 1; /* abs(a) is bigger and a is negative */
- else _sign(result) = 0;
- break;
- case 0: /* a == b */
- if (rounding_mode == FC_TONEGATIVE)
- _sign(result) = 1;
- else
- _sign(result) = 0;
- break;
- case -1: /* a < b */
- if (_sign(b)) _sign(result) = 1; /* abs(b) is bigger and b is negative */
- else _sign(result) = 0;
- break;
- default:
- /* can't be reached */
- break;
- }
- } else {
- _sign(result) = _sign(a);
- }
-
- /* sign has been taken care of, check for special cases */
- if (_desc(a).class == ZERO) {
- if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
- return result;
- }
- if (_desc(b).class == ZERO) {
- if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
- return result;
- }
-
- if (_desc(a).class == INF) {
- if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
- return result;
- }
- if (_desc(b).class == INF) {
- if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
- return result;
- }
-
- /* shift the smaller value to the right to align the radix point */
- /* subnormals have their radix point shifted to the right,
- * take care of this first */
- if ((_desc(b).class == SUBNORMAL) && (_desc(a).class != SUBNORMAL))
- {
- sc_val_from_ulong(1, temp);
- sc_sub(exp_diff, temp, exp_diff);
- }
-
- _shift_right(_mant(b), exp_diff, temp);
- sticky = sc_had_carry();
-
- if (sticky && sign)
- {
- /* if subtracting a little more than the represented value or adding a little
- * more than the represented value to a negative value this, in addition to the
- * still set sticky bit, takes account of the 'little more' */
- char *temp1 = alloca(calc_buffer_size);
- sc_val_from_ulong(1, temp1);
- sc_add(temp, temp1, temp);
- }
-
- if (sign) {
- if (sc_comp(_mant(a), temp) == -1)
- sc_sub(temp, _mant(a), _mant(result));
- else
- sc_sub(_mant(a), temp, _mant(result));
- } else {
- sc_add(_mant(a), temp, _mant(result));
- }
-
- /* _normalize expects a 'normal' radix point, adding two subnormals
- * results in a subnormal radix point -> shifting before normalizing */
- if ((_desc(a).class == SUBNORMAL) && (_desc(b).class == SUBNORMAL))
- {
- sc_val_from_ulong(1, NULL);
- _shift_left(_mant(result), sc_get_buffer(), _mant(result));
- }
-
- /* resulting exponent is the bigger one */
- memmove(_exp(result), _exp(a), value_size);
-
- return _normalize(result, result, sticky);
-}
-
-static char* _mul(const char* a, const char* b, char* result)
-{
- char *temp;
-
- if (_desc(a).class == NAN) {
- if (a != result) memcpy(result, a, calc_buffer_size);
- return result;
- }
- if (_desc(b).class == NAN) {
- if (b != result) memcpy(result, b, calc_buffer_size);
- return result;
- }
-
- temp = alloca(value_size);
-
- if (result != a && result != b)
- memcpy(&_desc(result), &_desc(a), sizeof(descriptor_t));
-
- _sign(result) = _sign(a) ^ _sign(b);
-
- /* produce nan on 0 * inf */
- if (_desc(a).class == ZERO) {
- if (_desc(b).class == INF)
- fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
- else
- if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
- return result;
- }
- if (_desc(b).class == ZERO) {
- if (_desc(a).class == INF)
- fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
- else
- if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
- return result;
- }
-
- if (_desc(a).class == INF) {
- if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
- return result;
- }
- if (_desc(b).class == INF) {
- if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-1);
- return result;
- }
-
- /* exp = exp(a) + exp(b) - excess */
- sc_add(_exp(a), _exp(b), _exp(result));
-
- sc_val_from_ulong((1<<_desc(a).exponent_size)/2-1, temp);
- sc_sub(_exp(result), temp, _exp(result));
-
- /* mixed normal, subnormal values introduce an error of 1, correct it */
- if ((_desc(a).class == SUBNORMAL) ^ (_desc(b).class == SUBNORMAL))
- {
- sc_val_from_ulong(1, temp);
- sc_add(_exp(result), temp, _exp(result));
- }
-
- sc_mul(_mant(a), _mant(b), _mant(result));
-
- /* realign result: after a multiplication the digits right of the radix
- * point are the sum of the factors' digits after the radix point. As all
- * values are normalized they both have the same amount of these digits,
- * which has to be restored by proper shifting
- * +2 because of the two rounding bits */
- sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
-
- _shift_right(_mant(result), temp, _mant(result));
-
- return _normalize(result, result, sc_had_carry());
+static void _fadd(const fp_value *a, const fp_value *b, fp_value *result) {
+ char *temp;
+ char *exp_diff;
+
+ char sign, res_sign;
+ char sticky;
+
+ fc_exact = 1;
+
+ handle_NAN(a, b, result);
+
+ /* make sure result has a descriptor */
+ if (result != a && result != b)
+ result->desc = a->desc;
+
+ /* determine if this is an addition or subtraction */
+ sign = a->sign ^ b->sign;
+
+ /* produce NaN on inf - inf */
+ if (sign && (a->desc.clss == INF) && (b->desc.clss == INF)) {
+ fc_exact = 0;
+ fc_get_qnan(&a->desc, result);
+ return;
+ }
+
+ temp = alloca(value_size);
+ exp_diff = alloca(value_size);
+
+ /* get exponent difference */
+ sc_sub(_exp(a), _exp(b), exp_diff);
+
+ /* initially set sign to be the sign of a, special treatment of subtraction
+ * when exponents are equal is required though.
+ * Also special care about the sign is needed when the mantissas are equal
+ * (+/- 0 ?) */
+ if (sign && sc_val_to_long(exp_diff) == 0) {
+ switch (sc_comp(_mant(a), _mant(b))) {
+ case 1: /* a > b */
+ res_sign = a->sign; /* abs(a) is bigger and a is negative */
+ break;
+ case 0: /* a == b */
+ res_sign = (rounding_mode == FC_TONEGATIVE);
+ break;
+ case -1: /* a < b */
+ res_sign = b->sign; /* abs(b) is bigger and b is negative */
+ break;
+ default:
+ /* can't be reached */
+ res_sign = 0;
+ break;
+ }
+ }
+ else
+ res_sign = a->sign;
+ result->sign = res_sign;
+
+ /* sign has been taken care of, check for special cases */
+ if (a->desc.clss == ZERO || b->desc.clss == INF) {
+ if (b != result)
+ memcpy(result, b, calc_buffer_size);
+ fc_exact = b->desc.clss == NORMAL;
+ result->sign = res_sign;
+ return;
+ }
+ if (b->desc.clss == ZERO || a->desc.clss == INF) {
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ fc_exact = a->desc.clss == NORMAL;
+ result->sign = res_sign;
+ return;
+ }
+
+ /* shift the smaller value to the right to align the radix point */
+ /* subnormals have their radix point shifted to the right,
+ * take care of this first */
+ if ((b->desc.clss == SUBNORMAL) && (a->desc.clss != SUBNORMAL)) {
+ sc_val_from_ulong(1, temp);
+ sc_sub(exp_diff, temp, exp_diff);
+ }
+
+ _shift_right(_mant(b), exp_diff, temp);
+ sticky = sc_had_carry();
+ fc_exact &= !sticky;
+
+ if (sticky && sign) {
+ /* if subtracting a little more than the represented value or adding a little
+ * more than the represented value to a negative value this, in addition to the
+ * still set sticky bit, takes account of the 'little more' */
+ char *temp1 = alloca(calc_buffer_size);
+ sc_val_from_ulong(1, temp1);
+ sc_add(temp, temp1, temp);
+ }
+
+ if (sign) {
+ if (sc_comp(_mant(a), temp) == -1)
+ sc_sub(temp, _mant(a), _mant(result));
+ else
+ sc_sub(_mant(a), temp, _mant(result));
+ } else {
+ sc_add(_mant(a), temp, _mant(result));
+ }
+
+ /* _normalize expects a 'normal' radix point, adding two subnormals
+ * results in a subnormal radix point -> shifting before normalizing */
+ if ((a->desc.clss == SUBNORMAL) && (b->desc.clss == SUBNORMAL)) {
+ sc_val_from_ulong(1, NULL);
+ _shift_left(_mant(result), sc_get_buffer(), _mant(result));
+ }
+
+ /* resulting exponent is the bigger one */
+ memmove(_exp(result), _exp(a), value_size);
+
+ fc_exact &= normalize(result, result, sticky);