libm: implement accelerated computation of (x << e) % y

[quaternic]

Builds on #1011, both were split from #1002

The computation is implemented in linear_mul_reduction.

[quaternic]

I've now rebased and applied the previous feedback.

@tgross35

Could you add a module-level doc comment with some of the common names used here? E.g. R, RR (if that's not R*R) r, m, q, xq. I'm unfortunately a bit lost :) (but I don't need to understand it in detail)

To make the module easier to follow, I moved the public interface, fn linear_mul_reduction, to the start of the file and added more explanatory comments. That should be understandable on its own. The implementation is essentially unchanged from before, other than one minor tweak.

[quaternic]

That minor tweak:
Previously the right shift was done first, and then the checked_sub was with the original y. Now that's done first with the left-shifted y, and the shift is last.

[quaternic]

I've gone through and rewritten the comments for the implementation of Reducer, which should actually be understandable now. (At least this time I didn't get too much of a headache from it.)

[tgross35]

The comments are super thorough and I was able to (to a reasonable level) understand what's going on, thank you for all the updates!