modular multiplication

Question

I was reading the Modular Multiplication page on wikipedia...and could not understand the algorithm to compute $a \cdot b \pmod{m}$.

uint64_t mul_mod(uint64_t a, uint64_t b, uint64_t m)
{
  long double x;
  uint64_t c;
  int64_t r;
  if (a >= m) a %= m;
  if (b >= m) b %= m;
  x = a;
  c = x * b / m;
  r = (int64_t)(a * b - c * m) % (int64_t)m;
  return r < 0 ? r + m : r;
}

I get that:

$$ a \cdot b \pmod{m} \equiv ((a \mod{m}) \cdot (b \mod{m})) \mod{m}$$

the step I didn't get is c=x*b/mand later doing r = (int64_t)(a * b - c * m) % (int64_t)m; what is that for?

The article says:

by employing the trick that, by hardware, floating-point multiplication results in the most significant bits of the product kept, while integer multiplication results in the least significant bits kept

I can't see how this is related !!!

Answer 1

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.

Answer 2

@D.W I noticed that for the case where e is negative, after the long double division, the value is being roundup, if that's the least significant bit, making it look like it represents exact value negating the rounding error. See this example. (2^63-2)(2^63-2) = 85070591730234615828950163710522949636. but this is represented as 85070591730234615828950163710522949632 in long double. This is a case where e is negative. when u divide it by (2^63-1), the answer should be 9223372036854775804.99999999999999999967. but this is being rounded up to 9223372036854775805. so we really are not getting the case of (r+m) as the quotient we getting is the actual q but not q-1(which we are expecting since e is negative). To avoid rounding up, we need more precision to hold the decimal. The max error u can get when m is of order 2^63 is 2^60 when ab is in the range of (2^125-2^126). Which requires at least 3 bits to hold the decimal value(e/m), but we are using 63 bits already to hold the mantissa of (a*b)/m. If we use lower order m, the error will be lesser too, even though we use lesser bits to hold the mantissa, we would need more bits to hold the decimal value. for example, if m is in the order of 2^58, the max error between 2^115 and 2^116 is 2^50. Now the decimal value will be of order e/m (1/2^8) which requires 8 bits. So it looks like the total bits necessary (58+8 = 66) is always 66 bits but we only have 63 bits and hence its rounding up always when e is negative.