__fmodx

pub fn __fmodx(a: f80, b: f80) callconv(.c) f80

File

Code

pub fn __fmodx(a: f80, b: f80) callconv(.c) f80 {
    const T = f80;
    const Z = @Int(.unsigned, @bitSizeOf(T));

    const significandBits = math.floatMantissaBits(T);
    const fractionalBits = math.floatFractionalBits(T);
    const exponentBits = math.floatExponentBits(T);

    const signBit = (@as(Z, 1) << (significandBits + exponentBits));
    const maxExponent = ((1 << exponentBits) - 1);

    var aRep: Z = @bitCast(a);
    var bRep: Z = @bitCast(b);

    const signA = aRep & signBit;
    var expA: i32 = @intCast((@as(Z, @bitCast(a)) >> significandBits) & maxExponent);
    var expB: i32 = @intCast((@as(Z, @bitCast(b)) >> significandBits) & maxExponent);

    // There are 3 cases where the answer is undefined, check for:
    //   - fmodx(val, 0)
    //   - fmodx(val, NaN)
    //   - fmodx(inf, val)
    // The sign on checked values does not matter.
    // Doing (a * b) / (a * b) produces undefined results
    // because the three cases always produce undefined calculations:
    //   - 0 / 0
    //   - val * NaN
    //   - inf / inf
    if (b == 0 or math.isNan(b) or expA == maxExponent) {
        return (a * b) / (a * b);
    }

    // Remove the sign from both
    aRep &= ~signBit;
    bRep &= ~signBit;
    if (aRep <= bRep) {
        if (aRep == bRep) {
            return 0 * a;
        }
        return a;
    }

    if (expA == 0) expA = normalize(f80, &aRep);
    if (expB == 0) expB = normalize(f80, &bRep);

    var highA: u64 = 0;
    const highB: u64 = 0;
    var lowA: u64 = @truncate(aRep);
    const lowB: u64 = @truncate(bRep);

    while (expA > expB) : (expA -= 1) {
        var high = highA -% highB;
        const low = lowA -% lowB;
        if (lowA < lowB) {
            high -%= 1;
        }
        if (high >> 63 == 0) {
            if ((high | low) == 0) {
                return 0 * a;
            }
            highA = 2 *% high + (low >> 63);
            lowA = 2 *% low;
        } else {
            highA = 2 *% highA + (lowA >> 63);
            lowA = 2 *% lowA;
        }
    }

    var high = highA -% highB;
    const low = lowA -% lowB;
    if (lowA < lowB) {
        high -%= 1;
    }
    if (high >> 63 == 0) {
        if ((high | low) == 0) {
            return 0 * a;
        }
        highA = high;
        lowA = low;
    }

    while ((lowA >> fractionalBits) == 0) {
        lowA = 2 *% lowA;
        expA = expA - 1;
    }

    // Combine the exponent with the sign and significand, normalize if happened to be denormalized
    if (expA < -fractionalBits) {
        return @bitCast(signA);
    } else if (expA <= 0) {
        return @bitCast((lowA >> @intCast(1 - expA)) | signA);
    } else {
        return @bitCast(lowA | (@as(Z, @as(u16, @intCast(expA))) << significandBits) | signA);
    }
}