为什么这个模块化 GCD 在大输入时会失败？

Question

这个问题实际上来自一个名为 codechef 的竞争性编程网站。

问题如下。

给定整数 A、B 和 N，你应该计算 A^N+B^N 的 GCD 和 |A-B|。 （假设任何正整数 a 的 GCD(0,a)=a）。从此 数字可能非常大，以 1000000007 (109+7) 为模计算。

完整的问题是here

值的限制如下：

1

1

B

当 A、B、N 的值非常大时，我的解决方案通过了第一个子任务但第二个失败。

这是我的解决方案：

#include <bitset>
#include <iostream>

using std::bitset;
using std::cin;
using std::cout;

typedef unsigned long long ull;
constexpr size_t bit_size = sizeof(ull) * 8;

ull modular_exponetiation(ull a, bitset<bit_size> n, ull mod) {
  ull c = 0;
  ull d = 1;
  for (int t = n.size() - 1; t >= 0; t--) {
    c *= 2;
    d = ((d % mod) * (d % mod)) % mod;  //(d*d)%mod
    if (n[t] == 1) {
      c++;
      d = ((d % mod) * (a % mod)) % mod;  //(d*a)%mod
    }
  }
  return d;
}

ull euclid_gcd(ull a, ull b) {
  if (b == 0)
    return a;
  else
    return euclid_gcd(b, a % b);
}

int main() {
  int test;
  cin >> test;
  while (test--) {
    ull a, b, n;
    cin >> a >> b >> n;
    ull modder = a - b;
    if (modder != 0) {
      ull out_res = 0;
      bitset<bit_size> bin_rip(n);
      ull first_mod_res = (modular_exponetiation(a, bin_rip, modder) +
                           modular_exponetiation(b, bin_rip, modder)) %
                          modder;
      if (first_mod_res == 0)
        out_res = modder;
      else
        out_res = euclid_gcd(modder, first_mod_res);
      cout << out_res % 1000000007 << std::endl;
    } else {
      // mod by 0 is not defined using the problem defined result.
      // GCD(0,a) == a;
      ull modder = 1000000007;
      bitset<bit_size> bin_rip(n);
      ull res = (modular_exponetiation(a, bin_rip, modder) +
                 modular_exponetiation(b, bin_rip, modder)) %
                modder;
      cout << res << std::endl;
    }
  }
  return 0;
}

请求

这不是作业，我也不需要精确的答案或代码更正。我明白所有这些，但不明白为什么它会在更大的值上失败？

任何方向或提示都会很有用。

Answer 1

如果modder = 1e12 那么你的模数不起作用。导致 1e12 * 1e12。会有溢出。

查看this 取模而不溢出。

你可以试试这个。这里的乘法是求和。

long long multiply(long long a,long long b,long long m){
if(b == 0){
    return 0;
}
if(b==1){
    return a%m;
}
if(b&1){
    return ((a%m)+multiply(a,b-1,m))%m;
}
long long x = multiply(a,b>>1,m);
return multiply(x,2,m);
}

long long bigmod(long long a,long long b, long long m){
if(b == 0){
    return 1;
}
if(b == 1){
    return a%m;
}
if(b & 1){
    return multiply(a%m,bigmod(a,b-1,m),m);
}
long long x = bigmod(a,b>>1,m);
return multiply(x,x,m);
}

Answer 2

特别感谢 SAJIB 指出问题。

这里我会自己回答，以便更好地理解。

d = ((d % mod) * (d % mod)) % mod; 
d = ((d % mod) * (a % mod)) % mod;

是中断并导致溢出的两行。所有功劳归于@sajib。

他还指出了解决这个问题的正确方法（加法乘法）。我很高兴理解他的代码的作用。

---

所以我在这里详细解释一下。

(a * b) % m

如果 a 和 b 都是非常大的 long long 值，则会导致溢出。

一种简单的方法是使用任意精度的数据类型，例如 python 中的 int 或 Java 中的 Biginteger 类。但是这种方法不会有成效，因为内部将string转换为int再进行运算会导致二进制数系统中加法和乘法的计算变慢。

有效的解决方案：由于a和b可能是非常大的数，如果直接相乘肯定会溢出。因此，我们使用乘法的基本方法，即，

a * b = a + a + … + a (b times)

所以我们可以很容易地计算加法的值（在模 m 下），而不需要任何 计算中溢出。但是如果我们尝试将 a 的值重复加到 b 次，那么对于较大的 b 值肯定会超时，因为这种方法的时间复杂度会变成 O(b)。

所以我们将a的上述重复步骤以更简单的方式划分，即，

If b is even then 
  a * b = 2 * a * (b / 2), 
otherwise 
  a * b = a + a * (b - 1)

实现如下：

// Returns (a * b) % mod
long long moduloMultiplication(long long a,
                               long long b,
                               long long mod)
{
    long long res = 0;  // Initialize result

    // Update a if it is more than
    // or equal to mod
    a %= mod;

    while (b)
    {
        // If b is odd, add a with result
        if (b & 1)
            res = (res + a) % mod;

        // Here we assume that doing 2*a
        // doesn't cause overflow
        a = (2 * a) % mod;

        b >>= 1;  // b = b / 2
    }

    return res;
}