+ /* 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 */
+ _sign(result) = _sign(a); /* abs(a) is bigger and a is negative */
+ break;
+ case 0: /* a == b */
+ _sign(result) = (rounding_mode == FC_TONEGATIVE);
+ break;
+ case -1: /* a < b */
+ _sign(result) = _sign(b); /* abs(b) is bigger and b is negative */
+ 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).clss == ZERO) {
+ if (b != result) memcpy(result+SIGN_POS+1, b+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
+ return result;
+ }
+ if (_desc(b).clss == ZERO) {
+ if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
+ return result;
+ }
+
+ if (_desc(a).clss == INF) {
+ if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-SIGN_POS-1);
+ return result;
+ }
+ if (_desc(b).clss == 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).clss == SUBNORMAL) && (_desc(a).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();
+
+ 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).clss == SUBNORMAL) && (_desc(b).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);
+
+ return _normalize(result, result, sticky);
+}
+
+/**
+ * calculate a * b
+ */
+static char* _fmul(const char* a, const char* b, char* result)
+{
+ char *temp;
+
+ handle_NAN(a, b, 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).clss == ZERO) {
+ if (_desc(b).clss == 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).clss == ZERO) {
+ if (_desc(a).clss == 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).clss == INF) {
+ if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
+ return result;
+ }
+ if (_desc(b).clss == 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).clss == SUBNORMAL) ^ (_desc(b).clss == 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());
+}
+
+/**
+ * calculate a / b
+ */
+static char* _fdiv(const char* a, const char* b, char* result)
+{
+ char *temp, *dividend;
+
+ handle_NAN(a, b, result);
+
+ temp = alloca(value_size);
+ dividend = 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/0 and inf/inf */
+ if (_desc(a).clss == ZERO) {
+ if (_desc(b).clss == ZERO)
+ /* 0/0 -> nan */
+ fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
+ else
+ /* 0/x -> a */
+ if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
+ return result;
+ }
+
+ if (_desc(b).clss == INF) {
+ if (_desc(a).clss == INF)
+ /* inf/inf -> nan */
+ fc_get_qnan(_desc(a).exponent_size, _desc(a).mantissa_size, result);
+ else {
+ /* x/inf -> 0 */
+ sc_val_from_ulong(0, NULL);
+ _save_result(_exp(result));
+ _save_result(_mant(result));
+ _desc(result).clss = ZERO;
+ }
+ return result;
+ }
+
+ if (_desc(a).clss == INF) {
+ /* inf/x -> inf */
+ if (a != result) memcpy(result+SIGN_POS+1, a+SIGN_POS+1, calc_buffer_size-1);
+ return result;
+ }
+ if (_desc(b).clss == ZERO) {
+ /* division by zero */
+ if (_sign(result))
+ fc_get_minusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
+ else
+ fc_get_plusinf(_desc(a).exponent_size, _desc(a).mantissa_size, result);
+ return result;
+ }
+
+ /* exp = exp(a) - exp(b) + excess - 1*/
+ sc_sub(_exp(a), _exp(b), _exp(result));
+ sc_val_from_ulong((1 << _desc(a).exponent_size)/2-2, temp);
+ sc_add(_exp(result), temp, _exp(result));
+
+ /* mixed normal, subnormal values introduce an error of 1, correct it */
+ if ((_desc(a).clss == SUBNORMAL) ^ (_desc(b).clss == SUBNORMAL))
+ {
+ sc_val_from_ulong(1, temp);
+ sc_add(_exp(result), temp, _exp(result));
+ }
+
+ /* mant(res) = mant(a) / 1/2mant(b) */
+ /* to gain more bits of precision in the result the dividend could be
+ * shifted left, as this operation does not loose bits. This would not
+ * fit into the integer precision, but due to the rounding bits (which
+ * are always zero because the values are all normalized) the divisor
+ * can be shifted right instead to achieve the same result */
+ sc_val_from_ulong(2 + _desc(result).mantissa_size, temp);
+
+ _shift_left(_mant(a), temp, dividend);
+
+ {
+ char *divisor = alloca(calc_buffer_size);
+ sc_val_from_ulong(1, divisor);
+ _shift_right(_mant(b), divisor, divisor);
+ sc_div(dividend, divisor, _mant(result));
+ }
+
+ return _normalize(result, result, sc_had_carry());
+}
+
+#if 0
+static void _power_of_ten(int exp, descriptor_t *desc, char *result)
+{
+ char *build;
+ char *temp;
+
+ /* positive sign */
+ _sign(result) = 0;
+
+ /* set new descriptor (else result is supposed to already have one) */
+ if (desc != NULL)
+ memcpy(&_desc(result), desc, sizeof(descriptor_t));
+
+ build = alloca(value_size);
+ temp = alloca(value_size);
+
+ sc_val_from_ulong((1 << _desc(result).exponent_size)/2-1, _exp(result));
+
+ if (exp > 0)
+ {
+ /* temp is value of ten now */
+ sc_val_from_ulong(10, NULL);
+ _save_result(temp);
+
+ for (exp--; exp > 0; exp--) {
+ _save_result(build);
+ sc_mul(build, temp, NULL);
+ }
+ _save_result(build);
+
+ /* temp is amount of left shift needed to put the value left of the radix point */
+ sc_val_from_ulong(_desc(result).mantissa_size + 2, temp);
+
+ _shift_left(build, temp, _mant(result));
+
+ _normalize(result, result, 0);
+ }