Specific modular multiplication algorithm [duplicate]

Question

I have 3 large 64 bit numbers: A, B and C. I want to compute:

(A x B) mod C

considering my registers are 64 bits, i.e. writing a * b actually yields (A x B) mod 2⁶⁴.

What is the best way to do it? I am coding in C, but don't think the language is relevant in this case.

---

After getting upvotes on the comment pointing to this solution:

(a * b) % c == ((a % c) * (b % c)) % c

let me be specific: this isn't a solution, because ((a % c) * (b % c)) may still be bigger than 2⁶⁴, and the register would still overflow and give me the wrong answer. I would have:

(((A mod C) x (B mod C)) mod 2⁶⁴) mod C

Mantinea

(A*B)%C it works ... or not ? – Iniquity

@GrijeshChauhan A*B overflows in C (N-bit number multiplied by N-bit number should ideally give you 2N-bit number, but that's not how arithmetic works in C), so you can't just use this truncated product. – Greathouse

added C tag , I Belive C people can help you better – Iniquity

#5916308 see the 2nd answer (by KennyTM) – Scandent

@AlexeyFrunze so may times I got very good answers from you. Thanks for this too :) – Iniquity

@AlexeyFrunze won't using usigned long long actually solve the issue? – Colettacolette

Too lazy to derive and check but wouldn't (A*B) mod C == ((A mod C) * (B mod C)) mod C make it easier? – Inefficiency

@Thrustmaster does not help at all. C may still be up to 2^64. – Colettacolette

@Scandent (a * b) % c == ((a % c) * (b % c)) % c doesn't help at all: ((a % c) * (b % c)) may still overflow the 64 bit register. – Mantinea

The Chinese remainder theorem comes to mind: en.wikipedia.org/wiki/Chinese_remainder_theorem - but this is pretty involved and doesn't really guarantee you'll get the exact solution. Also, I'm not entirely sure it even avoids the overflow issue while implementing it - just thought I'd share it. – Joggle

en.wikipedia.org/wiki/Karatsuba_algorithm - this involves splitting long integers into two shorter ones, turning 64-bit operands to a set of 32-bit operands – Soph

@harold How 128 bit output is stored in x86-64? – Mantinea

@harold You need mov rax, rdx before ret. – Greathouse

@harold Btw, there may be a division overflow (say, a=b=0xffffffffffffffff (their product will be 0xfffffffffffffffe0000000000000001) and c=1). – Greathouse

@harold Make c bigger, say 0xfffffffffffffffe, and with the same a=b=0xffffffffffffffff you will still have a division overflow. – Greathouse

Answer 1

As I have pointed in comment, Karatsuba's algorithm might help. But there's still a problem, which requires a separate solution.

Assume

A = (A1 << 32) + A2

B = (B1 << 32) + B2.

When we multiply those we get:

A * B = ((A1 * B1) << 64) + ((A1 * B2 + A2 * B1) << 32) + A2 * B2.

So we have 3 numbers we want to sum and one of this is definitely larger than 2^64 and another could be.

But it could be solved!

Instead of shifting by 64 bits once we can split it into smaller shifts and do modulo operation each time we shift. The result will be the same.

This will still be a problem if C itself is larger than 2^63, but I think it could be solved even in that case.

Soph

how ? A1 * B1 can be up to 2^64 again, so even if you shift it by 1 you'll have data loss. – Tribalism

They can't, they are both < 2^32. – Soph

So what ? They are < 2^32, their product is < 2^64, and it can easily overflow if you shift it. – Tribalism

Correct. But we can check that in advance and if we overflow we know that it was exactly 1 bit. Assuming we're shifting by 1 bit each time. 1 bit overflow is the known value and we even know how much it is mod C. – Soph

@Soph I'm probably missing it, but what about division/modulo, how do you suggest to do it once the product of A and B is computed? – Greathouse

Once we have calculated A1*B1, A1*B2, A2*B1 and A2*B2 we can start combining it applying modulo on each step to avoid overflow. There is no way we can "over-apply" it - we simply keep all our computations in modulo C. – Soph

Thank you very much for your explanation, it was really helpful. But I did not understand how we can avoid the overflow issue when C is 64 bits – Jacky