Optimize fmod with a method using integer multiplication

[quaternic]

This is kind of a retry at #898. One of the problems there was that it would have added overhead and regressed performance for typical inputs.

Unlike that PR, this doesn't aim for sub-linear scaling; the cost of evaluating fmod(x, y) is still roughly proportional to log2(|x/y|). However, the constant factor is much better. Running the random-benchmarks locally, I got walltime reductions of

fmodf16:  -56.9%
fmodf:    -85.0%
fmod:     -95.4%
fmodf128: -98.7%

New utilities in libm::support:

• trait NarrowingDiv for dividing u2N / uN when the quotient fits in uN
• a reasonable implementation of that for up to u256 / u128
• fn linear_mul_reduction<U>(x: U, mut e: u32, y: U) -> U computes (x << e) % y with the new method