# Modular multiplication with machine word limitations

Imagine I have 64-bit machine and the widest integer available is 64-bit signed long. I cannot use BigInteger or similar libraries for performance reasons, and all calculations I get would me modulo $2^{64}$. I need to choose prime $p$ close to $2^{63}$ but less than $2^{63}$ (which one would it be, I think choosing $2^{61}$ would make my computations faster?) and need to implement multiplication modulo $p$. Is there any known algorithm to do so?

-
If you're working on a 64-bit architecture, chances are that an exact 64 bit multiplication with a 128 bit result is supported natively. For example the imulq and mulq instructions on x86-64 do just that and place the result in two separate 64 bit registers. Check if such a multiplication is available in a (standard) library. Once you have that you're all set (see for example the answer by Richard). – WimC Nov 4 '12 at 8:14

If you can chose a 32-bit prime as modulus, then you don't need to do anything special. If you need a prime near the max 64-bit integer, then you'll have to implement the multiplication by hand. But this is not complicated. To multiply $a$ and $b$, write $a=a_1 \cdot 2^{32} + a_0$ and $b=b_1 \cdot 2^{32} + b_0$, with $0\le a_1,a_0,b_1,b_0 \lt 2^{32}$. Expand each term in $ab$, reduce them modulo $p$, add them and reduce the result modulo $p$.

-
Yes, I can do this. But for performance reasons this is not what I would be happy with. There should be a formula to do multiplication by modulo 2^61-1 using 64 bit arithmetic. – divby0 Apr 20 '11 at 13:22

If you're doing lots of multiplications of numbers modulo your prime then it's worth implementing Montgomery Reduction. I understand you're thinking of doing calculations mod $2^{61}-1$ but I don't think you can avoid at least implicitly calculating a 128-bit product when you multiply two numbers of this size together.

-

If you are thinking of implementing a 64-bit modulo multiplication in C, this is the fastest routine I know of. This beats all algorithmic implementations due to a hardware trick: long double type holds the most significant bits after a multiplication while an integer holds the least significant ones.

if (a >= m) a %= m;
if (b >= m) b %= m;

long double x = a;

uint64_t c = x * b / m;
int64_t r = (int64_t)(a * b - c * m) % (int64_t)m;
return r < 0 ? r + m : r;


where a, b and m are 64-bit unsigned integers. (Note that m < $2^{63}$ is a must)

-