+
+/**
+ * calculate a + b, where a is the value with the bigger exponent
+ */
+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->clss == FC_INF) && (b->clss == FC_INF)) {
+ fc_exact = 0;
+ fc_get_qnan(&a->desc, result);
+ return;
+ }
+
+ temp = (char*) alloca(value_size);
+ exp_diff = (char*) 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->clss == FC_ZERO || b->clss == FC_INF) {
+ if (b != result)
+ memcpy(result, b, calc_buffer_size);
+ fc_exact = b->clss == FC_NORMAL;
+ result->sign = res_sign;
+ return;
+ }
+ if (b->clss == FC_ZERO || a->clss == FC_INF) {
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ fc_exact = a->clss == FC_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->clss == FC_SUBNORMAL) && (a->clss != FC_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 = (char*) 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->clss == FC_SUBNORMAL) && (b->clss == FC_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);
+}
+
+/**
+ * calculate a * b
+ */
+static void _fmul(const fp_value *a, const fp_value *b, fp_value *result)
+{
+ int sticky;
+ char *temp;
+ char res_sign;
+
+ fc_exact = 1;
+
+ handle_NAN(a, b, result);
+
+ temp = (char*) alloca(value_size);
+
+ if (result != a && result != b)
+ result->desc = a->desc;
+
+ result->sign = res_sign = a->sign ^ b->sign;
+
+ /* produce NaN on 0 * inf */
+ if (a->clss == FC_ZERO) {
+ if (b->clss == FC_INF) {
+ fc_get_qnan(&a->desc, result);
+ fc_exact = 0;
+ } else {
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ result->sign = res_sign;
+ }
+ return;
+ }
+ if (b->clss == FC_ZERO) {
+ if (a->clss == FC_INF) {
+ fc_get_qnan(&a->desc, result);
+ fc_exact = 0;
+ } else {
+ if (b != result)
+ memcpy(result, b, calc_buffer_size);
+ result->sign = res_sign;
+ }
+ return;
+ }
+
+ if (a->clss == FC_INF) {
+ fc_exact = 0;
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ result->sign = res_sign;
+ return;
+ }
+ if (b->clss == FC_INF) {
+ fc_exact = 0;
+ if (b != result)
+ memcpy(result, b, calc_buffer_size);
+ result->sign = res_sign;
+ return;
+ }
+
+ /* exp = exp(a) + exp(b) - excess */
+ sc_add(_exp(a), _exp(b), _exp(result));
+
+ sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 1, temp);
+ sc_sub(_exp(result), temp, _exp(result));
+
+ /* mixed normal, subnormal values introduce an error of 1, correct it */
+ if ((a->clss == FC_SUBNORMAL) ^ (b->clss == FC_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
+ * because of the rounding bits */
+ sc_val_from_ulong(ROUNDING_BITS + result->desc.mantissa_size, temp);
+
+ _shift_right(_mant(result), temp, _mant(result));
+ sticky = sc_had_carry();
+ fc_exact &= !sticky;
+
+ fc_exact &= normalize(result, result, sticky);
+}
+
+/**
+ * calculate a / b
+ */
+static void _fdiv(const fp_value *a, const fp_value *b, fp_value *result)
+{
+ int sticky;
+ char *temp, *dividend;
+ char res_sign;
+
+ fc_exact = 1;
+
+ handle_NAN(a, b, result);
+
+ temp = (char*) alloca(value_size);
+ dividend = (char*) alloca(value_size);
+
+ if (result != a && result != b)
+ result->desc = a->desc;
+
+ result->sign = res_sign = a->sign ^ b->sign;
+
+ /* produce FC_NAN on 0/0 and inf/inf */
+ if (a->clss == FC_ZERO) {
+ if (b->clss == FC_ZERO) {
+ /* 0/0 -> NaN */
+ fc_get_qnan(&a->desc, result);
+ fc_exact = 0;
+ } else {
+ /* 0/x -> a */
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ result->sign = res_sign;
+ }
+ return;
+ }
+
+ if (b->clss == FC_INF) {
+ fc_exact = 0;
+ if (a->clss == FC_INF) {
+ /* inf/inf -> NaN */
+ fc_get_qnan(&a->desc, result);
+ } else {
+ /* x/inf -> 0 */
+ sc_val_from_ulong(0, NULL);
+ _save_result(_exp(result));
+ _save_result(_mant(result));
+ result->clss = FC_ZERO;
+ }
+ return;
+ }
+
+ if (a->clss == FC_INF) {
+ fc_exact = 0;
+ /* inf/x -> inf */
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+ result->sign = res_sign;
+ return;
+ }
+ if (b->clss == FC_ZERO) {
+ fc_exact = 0;
+ /* division by zero */
+ if (result->sign)
+ fc_get_minusinf(&a->desc, result);
+ else
+ fc_get_plusinf(&a->desc, result);
+ return;
+ }
+
+ /* exp = exp(a) - exp(b) + excess - 1*/
+ sc_sub(_exp(a), _exp(b), _exp(result));
+ sc_val_from_ulong((1 << (a->desc.exponent_size - 1)) - 2, temp);
+ sc_add(_exp(result), temp, _exp(result));
+
+ /* mixed normal, subnormal values introduce an error of 1, correct it */
+ if ((a->clss == FC_SUBNORMAL) ^ (b->clss == FC_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(ROUNDING_BITS + result->desc.mantissa_size, temp);
+
+ _shift_left(_mant(a), temp, dividend);
+
+ {
+ char *divisor = (char*) alloca(calc_buffer_size);
+ sc_val_from_ulong(1, divisor);
+ _shift_right(_mant(b), divisor, divisor);
+ sc_div(dividend, divisor, _mant(result));
+ sticky = sc_had_carry();
+ fc_exact &= !sticky;
+ }
+
+ fc_exact &= normalize(result, result, sticky);
+}
+
+#if 0
+static void _power_of_ten(int exp, float_descriptor_t *desc, char *result)
+{
+ char *build;
+ char *temp;
+
+ /* positive sign */
+ result->sign = 0;
+
+ /* set new descriptor (else result is supposed to already have one) */
+ if (desc != NULL)
+ result->desc = *desc;
+
+ build = alloca(value_size);
+ temp = alloca(value_size);
+
+ sc_val_from_ulong((1 << (result->desc.exponent_size - 1)) - 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(result->desc.mantissa_size + ROUNDING_BITS, temp);
+
+ _shift_left(build, temp, _mant(result));
+
+ _normalize(result, result, 0);
+ }
+}
+#endif
+
+/**
+ * Truncate the fractional part away.
+ *
+ * This does not clip to any integer range.
+ */
+static void _trunc(const fp_value *a, fp_value *result)
+{
+ /*
+ * When exponent == 0 all bits left of the radix point
+ * are the integral part of the value. For 15bit exp_size
+ * this would require a left shift of max. 16383 bits which
+ * is too much.
+ * But it is enough to ensure that no bit right of the radix
+ * point remains set. This restricts the interesting
+ * exponents to the interval [0, mant_size-1].
+ * Outside this interval the truncated value is either 0 or
+ * it does not have fractional parts.
+ */
+
+ int exp_bias, exp_val;
+ char *temp;
+
+ /* fixme: can be exact */
+ fc_exact = 0;
+
+ temp = (char*) alloca(value_size);
+
+ if (a != result) {
+ result->desc = a->desc;
+ result->clss = a->clss;
+ }
+
+ exp_bias = (1 << (a->desc.exponent_size - 1)) - 1;
+ exp_val = sc_val_to_long(_exp(a)) - exp_bias;
+
+ if (exp_val < 0) {
+ sc_val_from_ulong(0, NULL);
+ _save_result(_exp(result));
+ _save_result(_mant(result));
+ result->clss = FC_ZERO;
+
+ return;
+ }
+
+ if (exp_val > (long)a->desc.mantissa_size) {
+ if (a != result)
+ memcpy(result, a, calc_buffer_size);
+
+ return;
+ }
+
+ /* set up a proper mask to delete all bits right of the
+ * radix point if the mantissa had been shifted until exp == 0 */
+ sc_max_from_bits(1 + exp_val, 0, temp);
+ sc_val_from_long(a->desc.mantissa_size - exp_val + 2, NULL);
+ _shift_left(temp, sc_get_buffer(), temp);
+
+ /* and the mask and return the result */
+ sc_and(_mant(a), temp, _mant(result));
+
+ if (a != result) {
+ memcpy(_exp(result), _exp(a), value_size);
+ result->sign = a->sign;
+ }
+}