Here's a more compact and mathematical description of what is going on. Let $a$ and $b$ be the input, already reduced modulo $m$, so $a < m$ and $b < m$. (Code-wise, this means after the b %= m
line.) We want to calculate $ab \mod m$, which is to say, setting $x=ab$, we want to find $r_x$ such that $x = q_xm + r_x$ for some $q_x$.
The Non-Overflowing Case
In the non-overflowing case, we could just calculate the modulus and be done with it. Instead, the code calculates: $$y = \left\lfloor\frac{x}{m}\right\rfloor = \left\lfloor q_x + \frac{r_x}{m}\right\rfloor = q_x + \left\lfloor\frac{r_x}{m}\right\rfloor = q_x$$
Now $x - ym = q_xm + r_x - q_xm = r_x$. I elided a concern here though. It could be that $x$ doesn't fit in the mantissa of our floating point variable. As long as $m$ fits, though, this will not be a problem as $x = ab < m^2$. We'll end up with $x' = x + e$ where $e$ is some round-off error whose magnitude is less than $m$. We proceed as before: $$y = \left\lfloor\frac{x'}{m}\right\rfloor = \left\lfloor q_x + \frac{r_x+e}{m}\right\rfloor = q_x + \left\lfloor\frac{r_x+e}{m}\right\rfloor$$
In this case, we can't eliminate the $\left\lfloor\frac{r_x+e}{m}\right\rfloor$ because, while $r_x < m$ and $e < m$, it may be the case that $r_x+e > m$ or $r_x+e < 0$, though it will certainly be (strictly) between $-m$ and $2m$, so $\left\lfloor\frac{r_x+e}{m}\right\rfloor$ is either $0$ or $\pm1$. Now $$x - ym = q_x m + r_x - q_x m - \left\lfloor\frac{r_x+e}{m}\right\rfloor m = r_x - \left\lfloor\frac{r_x+e}{m}\right\rfloor m$$
Performing the modulus operation now will get rid of that extra term, however, due to the convention C chooses, (-1)%m = -1
. To get a convention where we always return a positive number, we can add $m$ to the result if negative.
The Overflowing Case
Let's assume that we are doing the integer arithmetic mod $N$, e.g. $2^{64}$, and let's assume $ab > N$ which means $m^2 > N > m$. Now the multiplication is going to wrap. We'll write $x = ab = z + kN$ for an integer $k < m$. That means the computation a*b
will give the number $z$. The floating point computation is as before assuming, again, that the $m$ fits in the mantissa (which may require 80-bit extended precision floats for larger $m$). For the modulus, define $z = q_z m + r_z$ and $kN = q_N m + r_N$ so $x = q_z m + q_N m + r_z + r_N$. As before storing $x$ as a floating point variable may produce some round-off error $|e| < m$ so again have $x' = x + e$. If we do the same calculation as before it looks like: $$\begin{align}
z - \left\lfloor\frac{x'}{m}\right\rfloor m
& = q_zm + r_z - q_z m - q_N m - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\
& = r_z - q_N m - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\
& = r_z + r_N - \left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor m \\
\end{align}$$
I've pulled a rabbit out of the hat here. Here, we're doing addition mod $N$, and mod $N$ $q_N m + r_N = kN = 0$ so $q_N m = -r_N$. As before, we mod by $m$ to get $r_z + r_N = x \mod m$, and, as before, the C convention may lead to this being negative which will need correcting.
There's one potential remaining issue. $\left\lfloor\frac{r_z+r_N+e}{m}\right\rfloor$ could be $2$. This causes no problem unless $2m > N$ in which case you'd get the wrong answer. For $N = 2^{64}$ this means $m > 2^{63}$. (As an aside, if $2m = N$ it causes no issue since we'll end up with $0$.) If we configure the floating point hardware to round down so that $e \leq 0$, this case will not come up. (Though don't forget my restriction that the mantissa can hold $m$!).
Connecting this to most/least significant bits
To specifically address the part that you quoted, consider an integer, $B$, greater than 2 which we'll think of as a base, as in "base $10$". In the usual $B=10$ case, the number $21$ is represented as $21 = 2B+1$. In a floating point representation we'd write this as $2.1\times B^1$. Now let's say we wanted to multiply two (positive) single base-$B$ digit numbers, the result would require at most two digits, say $cB+d$ where (as required by a [standard] base-$B$ representation) $0 \leq c < B$ and $0\leq d < B$. If we're restricted to only being able to store one base-$B$ digit of the result, there are two obvious choices, either $c$ or $d$.
Choosing $d$ is equivalent to working mod $B$ as $cB + d = d\mod B$, this is what happens with integer arithmetic. (Incidentally, at the assembly level, integer multiplication often does produce both of these digits.)
Floating point arithmetic, on the other hand, effectively corresponds to choosing $c$, but compensating that by incrementing the exponent. In effect, we represent the result as $c.d\times B^1$ but since we can store only one base-$B$ digit, this becomes just $c\times B^1$. (In practice, we'll consider numbers as multi-digit numbers in a small base (i.e. 2), rather than 1- or 2-digit numbers in a large base. This allows us to save some of the higher digits of $d$ if they aren't needed to store $c$, but in the worst-case scenario all of $d$ is lost. None of $c$ is lost until we start running out of room in the exponent. For the code above, this is not an issue.)
As long as $m$ can be represented faithfully in the floating point format, the expression $\left\lfloor\frac{ab}{m}\right\rfloor$ can be viewed as extracting that upper digit in base-$m$. You can view the code and math above as the interplay between base-$N$ and base-$m$ representations of a number.
Practicalities
Based on section 5.2.4.2.2 of this draft, the C11 standard appears to only require long double
to have a mantissa roughly 33 bits in length. (In particular, it appears to only specify the minimum number of decimal digits that can be faithfully represented.) In practice, most C compilers when targeting general purpose CPUs and particularly x86-family CPUs, will use IEEE754 types. In this case double
will effectively have a 53-bit mantissa. x86-family CPUs support an 80-bit format with a mantissa with effectively 64-bits, and several but not all compilers will have long double
indicate that when targeting x86. The range of validity of the code depends on these implementation details.