- int hsb;
- char lsb, guard, round, round_dir = 0;
- char *temp;
-
- temp = alloca(VALUE_SIZE);
-
- /* +2: save two rounding bits at the end */
- hsb = 2 + _desc(in_val).mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
-
- if (in_val != out_val)
- {
- _sign(out_val) = _sign(in_val);
- memcpy(&_desc(out_val), &_desc(in_val), sizeof(descriptor_t));
- }
-
- _desc(out_val).class = NORMAL;
-
- /* mantissa all zeroes, so zero exponent (because of explicit one)*/
- if (hsb == 2 + _desc(in_val).mantissa_size)
- {
- sc_val_from_ulong(0, _exp(out_val));
- hsb = -1;
- }
-
- /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
- if (hsb < -1)
- {
- /* shift right */
- sc_val_from_ulong(-hsb-1, temp);
-
- _shift_right(_mant(in_val), temp, _mant(out_val));
-
- /* remember if some bits were shifted away */
- if (!sticky) sticky = sc_had_carry();
-
- sc_add(_exp(in_val), temp, _exp(out_val));
- }
- else if (hsb > -1)
- {
- /* shift left */
- sc_val_from_ulong(hsb+1, temp);
-
- _shift_left(_mant(in_val), temp, _mant(out_val));
-
- sc_sub(_exp(in_val), temp, _exp(out_val));
- }
-
- /* check for exponent underflow */
- if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
- DEBUGPRINTF(("Exponent underflow!\n"));
- /* exponent underflow */
- /* shift the mantissa right to have a zero exponent */
- sc_val_from_ulong(1, temp);
- sc_sub(temp, _exp(out_val), NULL);
-
- _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
- if (!sticky) sticky = sc_had_carry();
- /* denormalized means exponent of zero */
- sc_val_from_ulong(0, _exp(out_val));
-
- _desc(out_val).class = SUBNORMAL;
- }
-
- /* perform rounding by adding a value that clears the guard bit and the round bit
- * and either causes a carry to round up or not */
- /* get the last 3 bits of the value */
- lsb = sc_sub_bits(_mant(out_val), _desc(out_val).mantissa_size + 2, 0) & 0x7;
- guard = (lsb&0x2)>>1;
- round = lsb&0x1;
-
- switch (ROUNDING_MODE)
- {
- case FC_TONEAREST:
- /* round to nearest representable value, if in doubt choose the version
- * with lsb == 0 */
- round_dir = guard && (sticky || round || lsb>>2);
- break;
- case FC_TOPOSITIVE:
- /* if positive: round to one if the exact value is bigger, else to zero */
- round_dir = (!_sign(out_val) && (guard || round || sticky));
- break;
- case FC_TONEGATIVE:
- /* if negative: round to one if the exact value is bigger, else to zero */
- round_dir = (_sign(out_val) && (guard || round || sticky));
- break;
- case FC_TOZERO:
- /* always round to 0 (chopping mode) */
- round_dir = 0;
- break;
- }
- 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"));
-
- if (round_dir == 1)
- {
- guard = (round^guard)<<1;
- lsb = !(round || guard)<<2 | guard | round;
- }
- else
- {
- lsb = -((guard<<1) | round);
- }
-
- /* add the rounded value */
- if (lsb != 0) {
- sc_val_from_long(lsb, temp);
- sc_add(_mant(out_val), temp, _mant(out_val));
- }
-
- /* could have rounded down to zero */
- if (sc_is_zero(_mant(out_val)) && (_desc(out_val).class == SUBNORMAL))
- _desc(out_val).class = ZERO;
-
- /* check for rounding overflow */
- hsb = 2 + _desc(out_val).mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
- if ((_desc(out_val).class != SUBNORMAL) && (hsb < -1))
- {
- sc_val_from_ulong(1, temp);
- _shift_right(_mant(out_val), temp, _mant(out_val));
-
- sc_add(_exp(out_val), temp, _exp(out_val));
- }
- else if ((_desc(out_val).class == SUBNORMAL) && (hsb == -1))
- {
- /* overflow caused the matissa to be normal again,
- * so adapt the exponent accordingly */
- sc_val_from_ulong(1, temp);
- sc_add(_exp(out_val), temp, _exp(out_val));
-
- _desc(out_val).class = NORMAL;
- }
- /* no further rounding is needed, because rounding overflow means
- * the carry of the original rounding was propagated all the way
- * up to the bit left of the radix point. This implies the bits
- * to the right are all zeros (rounding is +1) */
-
- /* check for exponent overflow */
- sc_val_from_ulong((1 << _desc(out_val).exponent_size) - 1, temp);
- if (sc_comp(_exp(out_val), temp) != -1) {
- DEBUGPRINTF(("Exponent overflow!\n"));
- /* exponent overflow, reaction depends on rounding method:
- *
- * mode | sign of value | result
- *--------------------------------------------------------------
- * TO_NEAREST | + | +inf
- * | - | -inf
- *--------------------------------------------------------------
- * TO_POSITIVE | + | +inf
- * | - | smallest representable value
- *--------------------------------------------------------------
- * TO_NEAGTIVE | + | largest representable value
- * | - | -inf
- *--------------------------------------------------------------
- * TO_ZERO | + | largest representable value
- * | - | smallest representable value
- *--------------------------------------------------------------*/
- if (_sign(out_val) == 0)
- {
- /* value is positive */
- switch (ROUNDING_MODE) {
- case FC_TONEAREST:
- case FC_TOPOSITIVE:
- _desc(out_val).class = INF;
- break;
-
- case FC_TONEGATIVE:
- case FC_TOZERO:
- fc_get_max(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
- }
- } else {
- /* value is negative */
- switch (ROUNDING_MODE) {
- case FC_TONEAREST:
- case FC_TONEGATIVE:
- _desc(out_val).class = INF;
- break;
-
- case FC_TOPOSITIVE:
- case FC_TOZERO:
- fc_get_min(_desc(out_val).exponent_size, _desc(out_val).mantissa_size, out_val);
- }
- }
- }
-
- return out_val;
+ int exact = 1;
+ int hsb;
+ char lsb, guard, round, round_dir = 0;
+ char *temp = (char*) alloca(value_size);
+
+ /* save rounding bits at the end */
+ hsb = ROUNDING_BITS + in_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(in_val)) - 1;
+
+ if (in_val != out_val) {
+ out_val->sign = in_val->sign;
+ out_val->desc = in_val->desc;
+ }
+
+ out_val->desc.clss = NORMAL;
+
+ /* mantissa all zeros, so zero exponent (because of explicit one) */
+ if (hsb == ROUNDING_BITS + in_val->desc.mantissa_size) {
+ sc_val_from_ulong(0, _exp(out_val));
+ hsb = -1;
+ }
+
+ /* shift the first 1 into the left of the radix point (i.e. hsb == -1) */
+ if (hsb < -1) {
+ /* shift right */
+ sc_val_from_ulong(-hsb-1, temp);
+
+ _shift_right(_mant(in_val), temp, _mant(out_val));
+
+ /* remember if some bits were shifted away */
+ if (sc_had_carry()) {
+ exact = 0;
+ sticky = 1;
+ }
+ sc_add(_exp(in_val), temp, _exp(out_val));
+ } else if (hsb > -1) {
+ /* shift left */
+ sc_val_from_ulong(hsb+1, temp);
+
+ _shift_left(_mant(in_val), temp, _mant(out_val));
+
+ sc_sub(_exp(in_val), temp, _exp(out_val));
+ }
+
+ /* check for exponent underflow */
+ if (sc_is_negative(_exp(out_val)) || sc_is_zero(_exp(out_val))) {
+ DEBUGPRINTF(("Exponent underflow!\n"));
+ /* exponent underflow */
+ /* shift the mantissa right to have a zero exponent */
+ sc_val_from_ulong(1, temp);
+ sc_sub(temp, _exp(out_val), NULL);
+
+ _shift_right(_mant(out_val), sc_get_buffer(), _mant(out_val));
+ if (sc_had_carry()) {
+ exact = 0;
+ sticky = 1;
+ }
+ /* denormalized means exponent of zero */
+ sc_val_from_ulong(0, _exp(out_val));
+
+ out_val->desc.clss = SUBNORMAL;
+ }
+
+ /* perform rounding by adding a value that clears the guard bit and the round bit
+ * and either causes a carry to round up or not */
+ /* get the last 3 bits of the value */
+ lsb = sc_sub_bits(_mant(out_val), out_val->desc.mantissa_size + ROUNDING_BITS, 0) & 0x7;
+ guard = (lsb&0x2)>>1;
+ round = lsb&0x1;
+
+ switch (rounding_mode) {
+ case FC_TONEAREST:
+ /* round to nearest representable value, if in doubt choose the version
+ * with lsb == 0 */
+ round_dir = guard && (sticky || round || lsb>>2);
+ break;
+ case FC_TOPOSITIVE:
+ /* if positive: round to one if the exact value is bigger, else to zero */
+ round_dir = (!out_val->sign && (guard || round || sticky));
+ break;
+ case FC_TONEGATIVE:
+ /* if negative: round to one if the exact value is bigger, else to zero */
+ round_dir = (out_val->sign && (guard || round || sticky));
+ break;
+ case FC_TOZERO:
+ /* always round to 0 (chopping mode) */
+ round_dir = 0;
+ break;
+ }
+ DEBUGPRINTF(("Rounding (s%d, l%d, g%d, r%d, s%d) %s\n", out_val->sign, lsb>>2, guard, round, sticky, (round_dir)?"up":"down"));
+
+ if (round_dir == 1) {
+ guard = (round^guard)<<1;
+ lsb = !(round || guard)<<2 | guard | round;
+ } else {
+ lsb = -((guard<<1) | round);
+ }
+
+ /* add the rounded value */
+ if (lsb != 0) {
+ sc_val_from_long(lsb, temp);
+ sc_add(_mant(out_val), temp, _mant(out_val));
+ exact = 0;
+ }
+
+ /* could have rounded down to zero */
+ if (sc_is_zero(_mant(out_val)) && (out_val->desc.clss == SUBNORMAL))
+ out_val->desc.clss = ZERO;
+
+ /* check for rounding overflow */
+ hsb = ROUNDING_BITS + out_val->desc.mantissa_size - sc_get_highest_set_bit(_mant(out_val)) - 1;
+ if ((out_val->desc.clss != SUBNORMAL) && (hsb < -1)) {
+ sc_val_from_ulong(1, temp);
+ _shift_right(_mant(out_val), temp, _mant(out_val));
+ if (exact && sc_had_carry())
+ exact = 0;
+ sc_add(_exp(out_val), temp, _exp(out_val));
+ } else if ((out_val->desc.clss == SUBNORMAL) && (hsb == -1)) {
+ /* overflow caused the mantissa to be normal again,
+ * so adapt the exponent accordingly */
+ sc_val_from_ulong(1, temp);
+ sc_add(_exp(out_val), temp, _exp(out_val));
+
+ out_val->desc.clss = NORMAL;
+ }
+ /* no further rounding is needed, because rounding overflow means
+ * the carry of the original rounding was propagated all the way
+ * up to the bit left of the radix point. This implies the bits
+ * to the right are all zeros (rounding is +1) */
+
+ /* check for exponent overflow */
+ sc_val_from_ulong((1 << out_val->desc.exponent_size) - 1, temp);
+ if (sc_comp(_exp(out_val), temp) != -1) {
+ DEBUGPRINTF(("Exponent overflow!\n"));
+ /* exponent overflow, reaction depends on rounding method:
+ *
+ * mode | sign of value | result
+ *--------------------------------------------------------------
+ * TO_NEAREST | + | +inf
+ * | - | -inf
+ *--------------------------------------------------------------
+ * TO_POSITIVE | + | +inf
+ * | - | smallest representable value
+ *--------------------------------------------------------------
+ * TO_NEAGTIVE | + | largest representable value
+ * | - | -inf
+ *--------------------------------------------------------------
+ * TO_ZERO | + | largest representable value
+ * | - | smallest representable value
+ *--------------------------------------------------------------*/
+ if (out_val->sign == 0) {
+ /* value is positive */
+ switch (rounding_mode) {
+ case FC_TONEAREST:
+ case FC_TOPOSITIVE:
+ out_val->desc.clss = INF;
+ break;
+
+ case FC_TONEGATIVE:
+ case FC_TOZERO:
+ fc_get_max(&out_val->desc, out_val);
+ }
+ } else {
+ /* value is negative */
+ switch (rounding_mode) {
+ case FC_TONEAREST:
+ case FC_TONEGATIVE:
+ out_val->desc.clss = INF;
+ break;
+
+ case FC_TOPOSITIVE:
+ case FC_TOZERO:
+ fc_get_min(&out_val->desc, out_val);
+ }
+ }
+ }
+ return exact;