Here's a simpler, but more technical, attempt at an explanation.

---

    while (true) {

Treat your initial 64 bit random value as a fixed point fraction in the range from 0 to just under 1.

        u64 fract = random64();

Convert it to a 128 bit type, still treating it as a fraction, and multiply it by another 64-bit nonzero integer, S. The result is a value between 0 and S (excluding S), with the fractional part in the lower 64 bits and the integer part in the upper 64 bits.

        u128 upto_s = (u128)fract * (u128)s;
        u64  fract_part = (u64)upto_s;
        u64  int_part   = (u64)(upto_s >> 64);

Because we multiplied by S, the possible resulting fraction bits increase by S between each sequential random value. This affects the number of possible fractions for each integer part. For example, with S=3, you could pack like

    . X . . X . . X . . X . . X . .

or

    X . . X . . X . . X . . X . . X

This is NOT the case if the number of possible fractions (here 2⁶⁴) is a multiple of S. We can therefore reject some of our possible fractions to make it a multiple of S. This can be done by rejecting the first 2⁶⁴ % s fractional values, and thus only accepting fractional values at least 2⁶⁴ % s.

        # Faster early test to avoid modulo
        # s > (2⁶⁴) % s
        if (fract_part >= s) {
            return int_part;
        }

        # 2⁶⁴ % s == (2⁶⁴ - s) % s == (-s) % s
        u64 min_fract = (-s) % s;
        if (fract_part >= min_fract) {
            return int_part;
        }
    }