6
$\begingroup$

IEEE Std 754-2008 is the modern definition of Floating-Point Arithmetic. It requires that division (among other operations) performs

as if it first produced an intermediate result correct to infinite precision (..), and then rounded that intermediate result, if necessary, to fit in the destination’s format.

I ask for a simple algorithm performing that operation, tailored to 8-bit CPUs with 8x8->16 bit multiplier; speed is a second(ary) criteria.

For simplicity I restrict to binary32 type

binary32
Image credit: wikipedia

and positive "normal" input with exponents such that no overflow or underflow occurs (so that we can ignore any sign or exponent consideration beyond some limited final exponent adjustment according to the result of the division), and rounding to nearest even:

the floating-point number nearest to the infinitely precise result shall be delivered; if the two nearest floating-point numbers bracketing an unrepresentable infinitely precise result are equally near, the one with an even least significant digit shall be delivered.

I ask because I strongly suspect that the AVR libc implementation of float division (typically used on the Arduino Uno) deviates from the standard beyond being locked to roundTiesToEven, and I wonder how hard it would be to fix that. The current code seems to be divsf3, which invokes divsf3x to perform the division to apparently 40-bit mantissa, then rounding with fp_round. Isn't that very sketch doomed, BTW?

$\endgroup$
  • 1
    $\begingroup$ The ultimate tiebreaker can be remainder > numerator/2 after the integer division done on the mantissas $\endgroup$ – ratchet freak Aug 31 '17 at 11:32
4
$\begingroup$

When using properly-rounded higher-precision floating-point arithmetic to compute operations in lower-precision arithmetic, there is the potential issue of double rounding: The result is first rounded to the higher precision, then rounded again during conversion to the lower precision. This can lead to incorrect results, but that is not always the case.

Samuel Figueroa, "When is double rounding innocuous?", ACM SIGNUM Newsletter, Vol. 30, No. 3, July 1995, pp. 21-26

shows that in the case of floating-point division, the higher-precision format needs to have twice the precision of the lower-precision format to avoid issues due to double rounding. This conditions is satisfied were we to use IEEE-754 binary64 to perform the division, then round the result to IEEE-754 binary32. On the other hand, if the higher precision uses 40 mantissa bits, its precision is less than twice that of binary32, so some incorrect results due to double rounding would be expected.

However, from the description in the question it doesn't sound like this is the case. Rather, it seems that there is an internal computation that computes many more quotient bits than are strictly necessary, and this un-rounded result is rounded once to get the binary32 result. After normalization of dividend and divisor a correctly rounded binary32 division requires a total of 25 quotient bits (24 bits for the significand including the integer bit that is implicit in the stored format, plus a round bit), plus a sticky bit that is the OR of all other potential quotient bits.

Note that we don't have to compute the potentially infinite number of quotient bits in order to determine the value of the sticky bit. It suffices to know whether the remainder after computation of 25 quotient bits is zero or not. A zero remainder means the sticky bit is zero, a non-zero remainder means the sticky bit is one. The rounding in nearest-or-even mode requires three pieces of information: the least significant bit of the significand, the round bit, and the sticky bit. The significand needs to be increased when either (1) the round bit is set and the sticky bit is set or (2; tie case) the round bit is set, the sticky bit is clear, but the LSB of the significand is set (meaning it is odd).

Conceptually, for division in a semi-logarithmic floating-point format one needs to subtract the exponents and divide the significand. The conceptually simplest way to perform the latter operation is to perform binary long-hand division. In each step, one compares the current remainder with the divisor. If it is greater than or equal to the divisor, the next quotient bits is one and one can subtract the divisor from the remainder. Then double the remainder for the next step.

Using a small 8x8->16 bit multiply would conceptually be much more involved. One would either build a narrow division operation from this building block and perform a high-radix long-hand division (e.g. base 256), or build a much wider multiplication operation with is then used to compute the quotient with a Newton-Raphson iteration for the reciprocal, using fixed-point computation. In the latter case, the rounding is more involved as an accurate remainder does not fall naturally out of the computation. In addition, a table of starting approximations, say 8 bits, for the reciprocal is typically needed. The following research report discusses relevant trade-offs in the context of the 8-bit Kestrel processor:

Eric Rice and Richard Hughey, "Multiprecision division on small-word parallel processors: Expanded Version". University of California Santa Cruz, April 1998 (online)

I incorporated the above algorithmic sketch into the (tested) ISO-C code below. This code is written for clarity and does not incorporate various obvious optimizations, however I have tried to minimize the size of variables under the assumption that the asker's platform has support for 8-bit operations and limited support for 16-bit operations, similar to a Z80.

#include <stdint.h>

#define SIGN_MASK     (0x80000000)  /* sign bit of IEEE-754 binary32 format */
#define EXPO_MASK     (0x7f800000)  /* exponent bits of IEEE-754 binary32 */
#define MANT_MASK     (0x007fffff)  /* significand bits of IEEE-754 binary32 */
#define MANT_INTBIT   (0x00800000)  /* where integer bit would be if stored */
#define QNAN_BIT      (0x00400000)  /* distinguishes QNaN from SNaN */
#define INFTY         (0x7f800000)  /* infinity */
#define INDEFINITE    (0xffc00000)  /* special QNaN */
#define MAX_NORM_EXPO (254)         /* maximum biased exponent for normals */
#define MIN_NORM_EXPO (1)           /* minimum biased exponent for normals */
#define EXPO_BIAS     (127)         /* exponent bias IEEE-754 binary32 */
#define MANT_BITS     (24)          /* effective # of mantissa bits binary32 */
#define QUOT_BITS     (MANT_BITS+1) /* one extra bit for rounding */
#define EXPO_SHFT     (MANT_BITS-1) /* integer bit is implicit */
#define ROUND_BIT     (1)           /* mask for round bit in 25-bit quotient */
#define MANT_LSB      (2)           /* mask for significand lsb in quotient */

uint32_t fp32_div_kernel (uint32_t a, uint32_t b)
{
    uint32_t r, x, y, sign;
    uint16_t expo_x, expo_y, expo_res;
    uint8_t i, odd, rnd, sticky;

    /* extract biased exponents and sign bits */
    expo_x = (a & EXPO_MASK) >> EXPO_SHFT;
    expo_y = (b & EXPO_MASK) >> EXPO_SHFT;
    sign = (a ^ b) & SIGN_MASK;

    if ((expo_x >= MIN_NORM_EXPO) && (expo_x <= MAX_NORM_EXPO) &&
        (expo_y >= MIN_NORM_EXPO) && (expo_y <= MIN_NORM_EXPO)) { /* fast path */
divide:
        /* extract significands */
        x = (a & MANT_MASK) | MANT_INTBIT;
        y = (b & MANT_MASK) | MANT_INTBIT;
        /* compute exponent of result */
        expo_res = expo_x - expo_y + EXPO_BIAS;
        /* dividend may not be smaller than divisor: normalize */
        if (x < y) {
            x = x << 1;
            expo_res--;
        }
        /* generate quotient one bit at at time */
        r = 0;
        for (i = 0; i < QUOT_BITS; i++) {
            r = r << 1;
            if (x >= y) {
                x = x - y;
                r = r | 1;
            }
            x = x << 1;
        }
        /* OR remainder bits into sticky bit */
        sticky = (x != 0);
        if ((expo_res >= MIN_NORM_EXPO) && 
            (expo_res <= MAX_NORM_EXPO)) { /* normal, may overflow to infinity*/
            /* extract round and lsb bits */
            rnd = (r & ROUND_BIT);
            odd = (r & MANT_LSB) != 0;
            /* remove round bit from quotient and round to-nearest-even */
            r = (r >> 1) + (rnd & (sticky | odd));
            /* combine exponent and significand */
            r = ((uint32_t)expo_res << EXPO_SHFT) + (r - MANT_INTBIT);
        } else if ((int16_t)expo_res > MAX_NORM_EXPO) { // overflow: infinity
            r = INFTY;
        } else { /* underflow: result is zero, subnormal, or smallest normal */
            uint8_t shift = (uint8_t)(1 - expo_res);
            /* clamp shift count */
            if (shift > QUOT_BITS) shift = QUOT_BITS;
            /* OR shifted-off bits of significand into sticky bit */
            sticky = sticky | ((r & ~(~0 << shift)) != 0);
            /* denormalize significand */
            r = r >> shift;
            /* extract round and lsb bits */
            rnd = (r & ROUND_BIT);
            odd = (r & MANT_LSB) != 0;
            /* remove round bit from quotient and round to-nearest-even */
            r = (r >> 1) + (rnd & (sticky | odd));
        }
        /* combine sign bit with combo of exponent and significand */
        r = r | sign;
    } else { /* slow path */
        /* take absolute value of arguments */
        x = a & ~SIGN_MASK;
        y = b & ~SIGN_MASK;
        /* if dividend is a NaN, convert that NaN into a QNaN and return it */
        if (x > INFTY) return a | QNAN_BIT;
        /* if divisor is a NaN, convert that NaN into a QNaN and return it */
        if (y > INFTY) return b | QNAN_BIT;
        /* dividend and divisor are both zero or infinity: invalid operation */
        if (((x == 0) && (y == 0)) || 
            ((x == INFTY) && (y == INFTY))) return INDEFINITE;
        /* 0/y or x/INF -> 0 */
        if ((x == 0) || (y == INFTY)) return sign;
        /* x/0 or INF/y -> INF */
        if ((y == 0) || (x == INFTY)) return sign | INFTY;
        /* if dividend is a subnormal, normalize it */
        if (expo_x == 0) {
            expo_x++;
            do {
                a = a + a;
                expo_x--;
            } while (!(a & MANT_INTBIT));
        }
        /* if divisor is a subnormal, normalize it */
        if (expo_y == 0) {
            expo_y++;
            do {
                b = b + b;
                expo_y--;
            } while (!(b & MANT_INTBIT));
        }
        /* now that dividend and divisor are normalized, do the division */
        goto divide;
    }
    return r;
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.