20

Is there a way to build e.g. (853467 * 21660421200929) % 100000000000007 without BigInteger libraries (note that each number fits into a 64 bit integer but the multiplication result does not)?

This solution seems inefficient:

int64_t mulmod(int64_t a, int64_t b, int64_t m) {
    if (b < a)
        std::swap(a, b);
    int64_t res = 0;
    for (int64_t i = 0; i < a; i++) {
        res += b;
        res %= m;
    }
    return res;
}
8
  • It looks like in this case you could cheat because you don't care about anything greater than parameter __int64 m (or uint64_t for those that favor it) hence you could deal only with 64-bit types. Commented Aug 28, 2012 at 22:29
  • 1
    Have you read about the Montgomery reduction algorithm? Commented Aug 28, 2012 at 22:31
  • 1
    @ildjarn, Montgomery reduction is unlikely to be a win in this case, since the product needs to be converted to Montgomery representation and back. Commented Aug 28, 2012 at 22:52
  • 2
    Funny, this would be trivial in x64 assembly. Commented Aug 28, 2012 at 23:01
  • 1
    GCC, Clang, and Intel ICC support __int128. I thought I remembered a similar name for MSVC, but struggle to find it now. Commented May 3, 2021 at 20:09

8 Answers 8

28

You should use Russian Peasant multiplication. It uses repeated doubling to compute all the values (b*2^i)%m, and adds them in if the ith bit of a is set.

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    int64_t res = 0;
    while (a != 0) {
        if (a & 1) res = (res + b) % m;
        a >>= 1;
        b = (b << 1) % m;
    }
    return res;
}

It improves upon your algorithm because it takes O(log(a)) time, not O(a) time.

Caveats: unsigned, and works only if m is 63 bits or less.

Sign up to request clarification or add additional context in comments.

1 Comment

Should res be declared uint64_t?
24

Keith Randall's answer is good, but as he said, a caveat is that it works only if m is 63 bits or less.

Here is a modification which has two advantages:

  1. It works even if m is 64 bits.
  2. It doesn't need to use the modulo operation, which can be expensive on some processors.

(Note that the res -= m and temp_b -= m lines rely on 64-bit unsigned integer overflow in order to give the expected results. This should be fine since unsigned integer overflow is well-defined in C and C++. For this reason it's important to use unsigned integer types.)

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t res = 0;
    uint64_t temp_b;

    /* Only needed if b may be >= m */
    if (b >= m) {
        if (m > UINT64_MAX / 2u)
            b -= m;
        else
            b %= m;
    }

    while (a != 0) {
        if (a & 1) {
            /* Add b to res, modulo m, without overflow */
            if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
                res -= m;
            res += b;
        }
        a >>= 1;

        /* Double b, modulo m */
        temp_b = b;
        if (b >= m - b)       /* Equiv to if (2 * b >= m), without overflow */
            temp_b -= m;
        b += temp_b;
    }
    return res;
}

4 Comments

I like this, since it handles full 64 bit values. If you test whether b < a at the top, and if so, swap a and b, it can significantly speed up the time since it's more likely the while loop can exit early.
Why can't you do b %= m if m>UINT64_MAX / 2u? Does modulo operation magically becomes unstable?
You can definitely do b %= m. But, modulo operations can be slow (depending on the processor) so are worth avoiding if possible. So, if (m > UINT64_MAX / 2u) b -= m; is a possible optimisation to avoid a modulo operation in the case that m is large, so the modulo can be simplified down to a simple subtraction.
Might be a tad late with this comment, but the res += b and b += temp_b do also overflow, which is not mentioned in your answer, even though the -= operations are specifically mentioned. Not from a C++ background, so I'm not sure if this is standard behaviour, but maybe add that in your answer?
5

Both methods work for me. The first one is the same as yours, but I changed your numbers to excplicit ULL. Second one uses assembler notation, which should work faster. There are also algorithms used in cryptography (RSA and RSA based cryptography mostly I guess), like already mentioned Montgomery reduction as well, but I think it will take time to implement them.

#include <algorithm>
#include <iostream>

__uint64_t mulmod1(__uint64_t a, __uint64_t b, __uint64_t m) {
  if (b < a)
    std::swap(a, b);
  __uint64_t res = 0;
  for (__uint64_t i = 0; i < a; i++) {
    res += b;
    res %= m;
  }
  return res;
}

__uint64_t mulmod2(__uint64_t a, __uint64_t b, __uint64_t m) {
  __uint64_t r;
  __asm__
  ( "mulq %2\n\t"
      "divq %3"
      : "=&d" (r), "+%a" (a)
      : "rm" (b), "rm" (m)
      : "cc"
  );
  return r;
}

int main() {
  using namespace std;
  __uint64_t a = 853467ULL;
  __uint64_t b = 21660421200929ULL;
  __uint64_t c = 100000000000007ULL;

  cout << mulmod1(a, b, c) << endl;
  cout << mulmod2(a, b, c) << endl;
  return 0;
}

8 Comments

I don't know about inline assembler, does it also use a loop?
@ChristianAmmer no and it doesn't need one. It uses a double-width multiplication and division is always double-width. It's only in high level languages that the high part of the multiplication suddenly gets lost.
This is OK for the example, but if (a * b > (2^64 - 1) * c) it will fail. But I assume the OP means that the implied quotient is a 64-bit value as well.
About loop: I don't know how assembler calculates the multiplication,I mean under the hoos, but on C++ there is no need for a loop, because due to %uint64 we know, that the result is 64bits at most. @Brett 3 numers are 64bit.
@Benjamin - I'm just pointing out that the assembly implementation isn't general. Try: a=8534670000000000000, b=216604212009290. Both are 64-bit, but divq will result in an exception.
|
5

An improvement to the repeating doubling algorithm is to check how many bits at once can be calculated without an overflow. An early exit check can be done for both arguments -- speeding up the (unlikely?) event of N not being prime.

e.g. 100000000000007 == 0x00005af3107a4007, which allows 16 (or 17) bits to be calculated per each iteration. The actual number of iterations will be 3 with the example.

// just a conceptual routine
int get_leading_zeroes(uint64_t n)
{
   int a=0;
   while ((n & 0x8000000000000000) == 0) { a++; n<<=1; }
   return a;
}

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n)
{
     uint64_t result = 0;
     int N = get_leading_zeroes(n);
     uint64_t mask = (1<<N) - 1;
     a %= n;
     b %= n;  // Make sure all values are originally in the proper range?
     // n is not necessarily a prime -- so both a & b can end up being zero
     while (a>0 && b>0)
     {
         result = (result + (b & mask) * a) % n;  // no overflow
         b>>=N;
         a = (a << N) % n;
     }
     return result;
}

2 Comments

+1, Nice speed improvement, but would recommend a check for n==0 to prevent an infinite loop in get_leading_zeroes().
I'd check if n <= 4^16 - when n == 0 return nan or some error, when n == 1 always return 0, otherwise if 1 < n <= 4^16, then (a % n) * (b % n) is, by definition, guaranteed to fit within uint64_t, so just directly return (a % n) * (b % n) % n instead of wasting time looping. Another trick I like using is when a == b and they're both a < n < (a + a) (i.e. over the half-way mark), then I'd square-mod the smaller half of things by first setting a = b = n - a.
3

You could try something that breaks the multiplication up into additions:

// compute (a * b) % m:

unsigned int multmod(unsigned int a, unsigned int b, unsigned int m)
{
    unsigned int result = 0;

    a %= m;
    b %= m;

    while (b)
    {
        if (b % 2 != 0)
        {
            result = (result + a) % m;
        }

        a = (a * 2) % m;
        b /= 2;
    }

    return result;
}

4 Comments

+1 for a working solution, I have to think about it to fully understand, but it works because (a * b) == (a * 2) * (b / 2), right?
It will actually fail for some inputs. a * 2 can overflow and be reduced incorrectly if m is bigger than 1 << 63 (or 1 << 31 if ints are 32 bit).
You can in fact reduce by (~0ULL/m) in each step. Eg. for 100000000000007, you can use 131072 (1<<17) instead of 2. That also explains harold's comment; for such large m the stepsize becomes 1 and you make no progress.
@harold: You're right: The first factor must not have its top bit set. I believe that's the only limitation on the function argument values of the present algorithm, though.
3

a * b % m equals a * b - (a * b / m) * m

Use floating point arithmetic to approximate a * b / m. The approximation leaves a value small enough for normal 64 bit integer operations, for m up to 63 bits.

This method is limited by the significand of a double, which is usually 52 bits.

uint64_t mod_mul_52(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t c = (double)a * b / m - 1;
    uint64_t d = a * b - c * m;

    return d % m;
}

This method is limited by the significand of a long double, which is usually 64 bits or larger. The integer arithmetic is limited to 63 bits.

uint64_t mod_mul_63(uint64_t a, uint64_t b, uint64_t m) {
    uint64_t c = (long double)a * b / m - 1;
    uint64_t d = a * b - c * m;

    return d % m;
}

These methods require that a and b be less than m. To handle arbitrary a and b, add these lines before c is computed.

a = a % m;
b = b % m;

In both methods, the final % operation could be made conditional.

return d >= m ? d % m : d;

Comments

1

I can suggest an improvement for your algorithm.

You actually calculate a * b iteratively by adding each time b, doing modulo after each iteration. It's better to add each time b * x, whereas x is determined so that b * x won't overflow.

int64_t mulmod(int64_t a, int64_t b, int64_t m)
{
    a %= m;
    b %= m;

    int64_t x = 1;
    int64_t bx = b;

    while (x < a)
    {
        int64_t bb = bx * 2;
        if (bb <= bx)
            break; // overflow

        x *= 2;
        bx = bb;
    }

    int64_t ans = 0;

    for (; x < a; a -= x)
        ans = (ans + bx) % m;

    return (ans + a*b) % m;
}

4 Comments

Can't you use x=(1<<63-m)/b ? That's rounded down, so b*x <= 1<<63 - m, and it doesn't take a loop to calculate that. Doesn't change big-O, as the number of for-loop iterations decreases by <50%.
@MSalters: doesn't that introduce a division, which is potentially expensive?
@CraigMcQueen: Yes, but only one, and there's already a modulo in a loop.
That's not improvement -- I don't know the O() for the last loop, but taking just 3 "random" numbers the loop ran over 300000 iterations.
0

Just so this point is not lost: some compilers like GCC and Clang support __int128 without external libraries. So as long as your inputs are 64-bit you could really do the multiplication in one x86-64 instruction, then a 128-bit mod.

Note: I have not rigorously tested this for negative inputs, 0 mod, etc.

int64_t modmul(int64_t a, int64_t b, int64_t m) 
{
    return ((__int128)a * b) % m;
}
modmul:
        mov     rcx, rsi
        mov     rax, rdi
        mov     rsi, rdx
        sub     rsp, 8
        imul    rcx       ; the actual multiply
        mov     rdi, rsi
        sar     rdi, 63
        mov     rcx, rdi
        mov     r8, rax
        mov     rax, rdx
        mov     rdx, rsi
        mov     rdi, r8
        mov     rsi, rax
        call    __modti3 ; GCC hands off mod to a call
        add     rsp, 8
        ret

https://godbolt.org/z/d7E8az3Mr

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.