A bitwise-magic fractional part function

frac.rs

/// A "fractional part" function on f32s that does bitwise magic
/// on the representation of its argument.
pub fn frac(f: f32) -> f32 {
    f32::from_bits(_frac(f.to_bits()))
}

/// A "fractional part" function that operates on a signed integer
/// representation of an IEEE single-precision floating point number.
fn _frac(i: u32) -> u32 {
    // A floating point number consists of a sign bit,
    // an 8-bit exponent, and a 23-bit "mantissa".
    //
    //   s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm
    //
    let s = i & 0x8000_0000;
    let mut e = (i >> 23) & 0xff;
    let mut m = i & 0x7f_ffff;  // I keep thinking "Mantissa" would
                                // make for a cute given name.

    // The top bit of m is the suggested "is-quiet" bit for NaNs.
    let is_quiet = 1 << 22;

    // OK, first, let's deal with the pathological cases.
    //
    //   * When e == 0xff and m == 0, the value is ±inf.
    //     We'll return a quiet NaN in this case.
    //     If i represents ±inf, then i | is_quiet works.
    //
    //   * When e == 0xff and m != 0, the value is NaN.
    //     We'll propagate the NaN value, quieted.
    //     So we return i | is_quiet in this case, too.
    //
    if e == 0xff { return i | is_quiet; }

    // Otherwise, the value is:
    //
    //     (−1)^s * 2^(e−127) * 0b1.mmmmmmmm...   if e == 0
    //     (−1)^s * 2^(1−127) * 0b0.mmmmmmmm...   otherwise
    //
    // We don't have any work to do when e < 127: those numbers are
    // so small in absolute value they only have a fractional part.
    if e < 127 { return i }

    // Since e >= 127, we know we're dealing with the first case.
    // The mantissa has a leading "1." that we need to keep in mind.
    //
    // We wish to turn e.g.
    //
    //   e.g. 101010.001000111  (e == 132, m == 01010001000111000000000)
    //        ^~~~~~~~~~~~~~~~         ^        ~~~~~~~~~~~~~~
    //   into      0.001000111  (e == 124, m == 00011100000000000000000)
    //                 ^~~~~~~         ^        ~~~~~~
    //
    // The exponent controls "which 1 the ^ points at", and m is all
    // of the bits after it (~~~). Note what happened to the values:
    //
    //   * e has been decremented until it points at the first fractional 1.
    //   * m has been shifted left and masked accordingly.

    // ———————————————————————————————————————————————————————————————————

    // First, mask away the bits of m we don't need: the ones "to the left
    // of the decimal point". Which ones are that? Suppose (e - 127) is 4,
    // for example. Then the formula says:
    //
    //     float = (−1)^s * 2^4 * 0b1.mmmmmmmmmmmmmmmmmmmmmmm
    //           = (−1)^s * 0b1mmmm.mmmmmmmmmmmmmmmmmmm
    //
    // So we want to mask away the top 4 bits of m, which we'll do by AND-ing
    // m with 0b00001111111111111111111.  That mask is 0x7f_ffff >> 4.
    // (Don't worry about the leading 1 bit: we'll deal with that next.)
    //
    // In general, we'll AND m by 0x7f_ffff >> (e - 127). Unfortunately,
    // "x >> y" can act weird ("underflow") when y >= bit_width(x) == 32.
    // (Some processors will take y modulo 32 in such a SHR instruction.)
    // So we limit (e - 127) to the size of the mantissa.
    m &= 0x7f_ffff >> (e - 127).min(23);

    // At this point, m might be zero, which means the input was an integer.
    // (It only had 1-bits to the left of the decimal point, and we cleared
    // them all.) In that case, we return 0 (which is the floating point
    // representation of 0.0f32).
    if m == 0 { return 0; }

    // Now we want to shift the highest fractional bit in m into the place
    // of the implicit "1." bit in our floating-point formula. Accordingly
    // we adjust the value of e.
    //
    //     Before:  (−1)^s * 2^(e−127)   * 1.00000001000111000000000
    //                                       ~~~~~~~^~~~~~~~~~~~~~~~
    //                                                  m
    //
    //     After:   (−1)^s * 2^(e−127−8) * 1.00011100000000000000000
    //                                     ^ ~~~~~~~~~~~~~~~~~~~~~~~
    //                                               (m << 8)
    //
    // We shift by as many places as there are leading zeroes in m,
    // plus one. (We discount the leading zeros that merely arise from
    // storing m, a "23-bit" number, in a u32 variable.)
    let t = m.leading_zeros() + 1 - (32 - 23);
    m <<= t; e -= t;

    // Make sure to actually clear the bit we "pushed out" of m. (In other
    // words, reinforce its 23-bit-ness.)
    m &= 0x7f_ffff;

    // OK, seeeeeeeemmmmmmmmmmmmmmmmmmmmmmms legit. Reassemble the value:
    s | (e << 23) | m

    // ———————————————————————————————————————————————————————————————————

    // To recap, here's the computation visualized for
    //
    //     frac(101010.001000111) = 0.001000111.
    //
    // First, there's the masking step.
    //
    //     2^5  * 1.01010001000111000000000    Initial value = 101010.001000111
    //     2^5  * 1.00000001000111000000000    "Left of ." bits of `m` masked out
    //
    // Then we slide m into the "1." bit's place, and adjust e accordingly.
    //
    //     2^5  * 1.00000001000111000000000    7 leading zeros
    //              ~~~~~~~^                   So: m <<= 8 and e -= 8.
    //     2^−3 * 1.00011100000000000000000    Final value = 0.001000111
}

fn main() {
    for &test in &[
        123.45,
        0.0001,
        -65.5,
        4509835.001f32,
        333333333333333.001f32,
        0.0,
        1.0,
        std::f32::NAN,
        std::f32::INFINITY,
        std::f32::NEG_INFINITY,
    ] {
        let mine = frac(test);
        let reference = test % 1.0;
        assert!(mine.is_nan() && reference.is_nan() || mine == reference);
        println!("OK: {}", test);
    }
}

This is a really fun implementation! I've always wanted to do bit shuffling on floats and this is a really nice and clean example. All of your comments are very useful + ALSO theyre cute, the:

    let mut m = i & 0x7f_ffff;  // I keep thinking "mantissa" would
                                // make for a cute given name.

was a surprise and it made me very happy. My other favorite comments are the s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm, the decorative ~* ruler, and the recap comment at the end working through a full example.

I've got a few suggestions I want to share because they use Rust function that I like and also cool floating-point facts.

First, you don't need to be unsafe when you're turning your float to bits, because there are two methods for taking floats to and from bits. You can write the frac function as:

pub fn frac(f: f32) -> f32 {
    f32::from_bits(_frac(f.to_bits()))
}

The next one is a NaN thing. I was reading IEEE-754 for SIGBOVIK prep and I was surprised how much structure they leave in NaNs. One, you probably want to handle signaling and quiet NaNs, and always return a quiet NaN. (I think?) You also want to preserve the "payload", which is all of the extra bits in the NaN.

An operation that propagates a NaN operand to its result and has a single NaN as an input should produce a NaN with the payload of the input NaN if representable in the destination format.

Did you know that IEEE-754 suggests that implementations might display the payload if it's not canonical? I remembered seeing examples of values like -sNaNx09af3 or whatever but I can't find them in the standard when I go back to check… (This is just a recommendation and is defined by languages anyway.)

The top bit of the mantissa is usually used as the "is-quiet" flag, and is the recommended bit to set when making the mantissa non-zero. This won't be completely portable because it's implementation-defined what is signaling and what is quiet, but this bit as an "is-quiet" flag is true on all of the rust Tier-1 platforms, and I think true on all supported platforms in fact.

It's unspecified whether operations should preserve the sign of NaNs, but if you want to, you can just stick s | back in there…

    if e == 0xff {
        return (0xff << 23) | (1 << 22) | m;
    }

Rust also makes it really easy to count leading/trailing zeros, so you can use the .leading_zeros() method instead of a loop when shifting the bits out without a bunch of fuss. Instead of the loop:

    while m & 0x80_0000 == 0 { m <<= 1; e -= 1; }

you can have:

    let shift = m.leading_zeros() - 8;
    m <<= shift;
    e -= shift;

Thank you so much for making a version of this with such great explanations, I really enjoyed it. I hope my comments are interesting & useful!! 💖