大数乘法取模

Question

我知道这个问题可能看起来很幼稚，但我已经在谷歌和这个网站上搜索了很多次，但都没有找到令人满意的答案。 我只是想计算(A*B)%MOD，假设a、b和MOD都是long long类型的。 假设MOD比A和B都大，以至于A%MOD = A且B%MOD = B，但A*B大于64位。如何正确地计算(A*B)%MOD的值？

- unrealsoul007

---

马特：我不同意复制。unrealsoul: A*B = (A-X)B + XB，你可以总是用这种方式把A拆分成较小的数字。例如，将X设置为floor(A/2)。然后如果子结果仍然太大，可以应用相同的过程。 - Cheers and hth. - Alf

这是被采纳答案的人的评论：“如果long的最大值为2^63 - 1，则只要x < 290331368171，1768431 * x就不会溢出。”但这正是我的疑问，如果A*B溢出了怎么办。 - unrealsoul007

分而治之 - Cheers and hth. - Alf

你可以使用GMP库。这将允许使用大于64位的整数。 - SethB

这篇与模1000000007问题不同，正如OP所指出的。其他问题的答案都不完全适用，因为它们都建议仅使用GMP或执行某些缓慢、错误或x86特定的操作。大家应该让这个问题保持开放状态。 - tmyklebu

Answer 1

这里的基本思路是首先定义一个不会溢出的addmod函数，它利用了负数在算术中的特性。然后使用位运算来定义timesmod函数。时间复杂度为O(N)，其中N是使用的位数（在本例中为64位）。

#include <iostream>
using namespace std;

typedef long long BigInt; // must be signed, to detect overflow

BigInt A = 0x7fffffffffffff01;
BigInt B = 0x7fffffffffffff02;
BigInt M = 0x7fffffffffffff03;

// For simplicity it is assumed x, y, and m are all positive.
BigInt addmod( BigInt x, BigInt y, BigInt m )
{
  x %= m;
  y %= m;
  BigInt sum = x-m+y; // -m <= sum < m-1
  return sum < 0 ? sum + m : sum;
}

BigInt timesmod( BigInt x, BigInt y, BigInt m )
{
  x %= m;
  y %= m;
  BigInt a = x < y ? x : y; // min
  BigInt b = x < y ? y : x; // max
  BigInt product = 0;
  for (; a != 0; a >>= 1, b = addmod(b,b,m) )
    if (a&1) product = addmod(product,b,m);
  return product;
}

int main()
{
  cout << "A = " << A << endl;
  cout << "B = " << B << endl;
  cout << "M = " << M << endl;
  cout << "A*B mod M = " << timesmod(A,B,M) << endl;
  return 0;
}

输出：

A = 9223372036854775553
B = 9223372036854775554
M = 9223372036854775555
A*B mod M = 2

这很容易确认，因为A=-2且B=-1 mod M。

注意：此代码未经过优化。

- Matt

---

问题在于您正在使用有符号类型的 BigInt，因此无法可靠地检测溢出（它会导致未定义的行为）。您需要使用无符号类型才能获得可靠的溢出行为。 - Chris Dodd

@ChrisDodd，你说的完全正确；感谢你指出这一点。我已经更新了addMod函数以解决这个问题。 - Matt

你的方法几乎与我编写的最终方法完全相同，如下所示：long long int multiply(long long a, long long b, long long mod) { long long result; if(b==0) return 0LL; result = multiply(a,b>>1,mod); result = (result+result)%mod; if(b&1) result = (result+a)%mod; return result; } - unrealsoul007

@unrealsoul007 我得到的结果是multiply(0x7fffffffffffff01,0x7fffffffffffff02,0x7fffffffffffff03)返回了-9223372036854711038。你得到什么结果？（我认为答案应该是 2。） - Matt

@Matt 只需稍微修改代码，将 long long 替换为 unsigned long long，你就能得到答案 2。 - unrealsoul007

Answer 2

我认为你可以将128位的乘积分成两部分（高64位和低64位），并对每个部分取模p。假设p大约是4^k，那么你可以通过将hi64 / (p>>k)相除来大致计算出这个数字中有多少个p；这应该会给你大约k-1位正确答案。从整个数字中减去这么多个p，现在hi64就少了大约k-1位。再次进行此操作，但计算(hi64 << k-1) / (p >> k)。然后再次进行此操作，计算(hi64 << k+k-2) / (p >> k)。

另一个帖子建议使用Schrage的技巧，听起来更好，但我不理解它。希望那个帖子的作者回来并完成他的回答！

- tmyklebu