fadd.c

// Function that adds 64-bit floats implemented using bit magic

double fadd(double xf, double yf) {
    const uint64_t exp_mask = ((UINT64_C(1)<<11)-1);
    const int64_t  exp_min = -1023;
    const int64_t  exp_max = 1022;
    const uint64_t mant_mask = ((UINT64_C(1)<<52)-1);
    // Extract fields from floats
    union{uint64_t i; double f;} v;
    v.f = xf;
    uint64_t x = v.i;
    int64_t  xexp = (int64_t)((x >> 52) & exp_mask) - 1023;
    uint64_t xmant = (x & mant_mask);
    if(xexp != exp_min) {
        xmant |= (uint64_t)1 << 52;
    }
    v.f = yf;
    uint64_t y = v.i;
    int64_t  yexp = (int64_t)((y >> 52) & exp_mask) - 1023;
    uint64_t ymant = (y & mant_mask);
    if(yexp != exp_min) {
        ymant |= (uint64_t)1 << 52;
    }
    // Handle special cases: infinities, NaN and zero
    if(xexp == exp_min && xmant == 0) return yf;
    if(yexp == exp_min && ymant == 0) return xf;
    if(xexp == exp_max) return xf;
    if(yexp == exp_max) return yf;
    // Now we treat floats like 2^exp * man, for this we need to subtract
    // 52 bits that we used to make mantissa an integer
    xexp -= 52;
    yexp -= 52;
    uint64_t rmant;
    int64_t rexp;
    // If exponents aren't equal we normalize one of those bois
    if(xexp != yexp) {
        // Make sure xexp < yexp (swap values around)
        if(xexp > yexp) {
            uint64_t temp_mant = xmant;
            int64_t temp_exp = xexp;
            xmant = ymant;
            xexp = yexp;
            ymant = temp_mant;
            yexp = temp_exp;
        }
        // Figure out whether we need to normalize
        int diff = yexp - xexp;
        if(diff == 63) {
            // Special case: if difference is 64 we round the lower  bit of
            // x's mantissa based on the two highest bit of y's mantissa
            // I'm not going to handle rounding modes because I'm too lazy
            if((ymant>>51)>=2) {
                xmant += 1;
                if(xmant == 0) {
                    xexp += 1;
                }
            }
            rexp = xexp;
            rmant = xmant;
            goto finish;
        }
        else if(diff > 63) {
            rexp = xexp;
            rmant = xmant;
            goto finish;
        }
        else {
            // Normalize x until xexp == yexp
            xmant >>= diff;
            xexp += diff;
        }
    }
    // Now add exponents and mantissas
    rexp = xexp;
    rmant = xmant + ymant;
finish:
    // If rmant is zero we're done, otherwise we normalize mantissa
    if(rmant != 0) {
        // Find the highest set bit
        int highest_set_bit = 0;
        int bit;
        for(bit = 0; bit != 64; ++bit) {
            if((rmant >> (bit+1)) == 0) {
                break;
            }
        }
        highest_set_bit = bit;
        // Restore the mantissa
        rmant <<= 52 - highest_set_bit;
        rexp  +=  highest_set_bit; //idk if ok
    }
    uint64_t rbits;
    if(rexp < exp_min) {
        rbits = 0;
    }
    else if(rexp > exp_max) {
        rbits = UINT64_C(0x7ff0000000000000);
    }
    else {
        rbits = (rmant & mant_mask) | ((uint64_t)(rexp - exp_min)<<52);
    }
    v.i = rbits;
    return v.f;
}