如何用原始类型进行模乘法

Question

有没有办法构建例如(853467 * 21660421200929) % 100000000000007没有BigInteger库(请注意,每个数字都适合64位整数,但乘法结果不适合)？

这个解决方案效率低下:

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;
}

Answer 1

你应该使用俄罗斯农民的乘法.它采用加倍重复计算所有的值(b*2^i)%m,而如果将其添加在i个位的位a设置.

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;
}

它改进了算法,因为它需要O(log(a))时间,而不是O(a)时间.

注意事项:无符号,仅当m63位或更少时才有效.

Answer 2

Keith Randall的答案很好,但正如他所说,一个警告是它只有在m63位或更少的情况下才有效.

这是一个有两个优点的修改:

1. 它即使m是64位也可以工作.
2. 它不需要使用模运算,这在某些处理器上可能很昂贵.

(注意,res -= m和temp_b -= m行依赖于64位无符号整数溢出,以便给出预期的结果.这应该没问题,因为无符号整数溢出在C和C++中是很好定义的.因此,使用无符号整数类型很重要. )

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;
}

Answer 3

两种方法都对我有用。第一个与你的相同，但我将你的数字更改为显式 ULL。第二个使用汇编符号，它应该运行得更快。还有密码学中使用的算法（我猜主要是 RSA 和基于 RSA 的密码学），就像已经提到的蒙哥马利约简一样，但我认为实现它们需要时间。

#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;
}

Answer 4

对重复加倍算法的一项改进是检查一次可以计算多少位而不发生溢出。可以对两个参数进行早期退出检查——加速 N 不是素数的（不太可能？）事件。

例如 100000000000007 == 0x00005af3107a4007，它允许每次迭代计算 16（或 17）位。在本例中，实际迭代次数为 3。

// 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;
}