Modular exponentiation fails for large mod in C++

问题:

This is the code I'm using for calculating (n^p)%mod. Unfortunately, it fails for large values of mod (in my case mod = 10000000000ULL) when I call it from main() method. Any idea; why?

ull powMod(ull n, ull p, ull mod) {
    ull ans = 1;
    n = n%mod;
    while(p) {
        if(p%2 == 1) {
            ans = (ans*n)%mod;
        }
        n = (n*n)%mod;
        p /= 2;
    }
    return ans;
}

Here, ull is a typedef for unsigned long long.

回答1:

Yes you can do it in C++. As others pointed out you cannot do it directly. Using a little drop of number theory it is possible to decompose the problem into two manageable sub problems.

First consider that 10^10 = 2^10 * 5^10. Both factors are coprime, so you can use the Chinese remainder theorem to find the power modulo 10^10 using the powers modulo 2^10 and modulo 5^10.

Note that in the following code the magic values u2 and u5 were found using the Extended Euclidean Algorithm. You don't need to program this algorithm yourself because these values are constants. I use maxima and its gcdex function, to compute them.

Here is the modified version:

typedef unsigned long long ull;

ull const M  = 10000000000ull;

ull pow_mod10_10(ull n, ull p) {
  ull const m2 = 1024;    // 2^10
  ull const m5 = 9765625; // 5^10
  ull const M2 = 9765625; // 5^10 = M / m2
  ull const M5 = 1024;    // 2^10 = M / m5
  ull const u2 = 841;     // u2*M2 = 1 mod m2
  ull const u5 = 1745224; // u5*M5 = 1 mod m5

  ull b2 = 1;
  ull b5 = 1;
  ull n2 = n % m2;
  ull n5 = n % m5;

  while(p) {
    if(p%2 == 1) {
      b2 = (b2*n2)%m2;
      b5 = (b5*n5)%m5;
    }
    n2 = (n2*n2)%m2;
    n5 = (n5*n5)%m5;
    p /= 2;
  }

  ull np = (((b2*u2)%M)*M2)%M;
  np    += (((b5*u5)%M)*M5)%M;
  np    %= M;
  return np;
}

回答2:

It seems that you can't avoid it.

If mod is 10000000000ULL, in (a*b)%c in your program, both a and b are smaller than mod so we treat them as 9999999999ULL, a*b will be 99999999980000000001, but unsigned long long can only express 2^64-1=18446744073709551615 < 99999999980000000001 so your method will be overflow.

回答3:

One of the possible issues here seems to be that when you do (a*b)%c, the a*b part itself could overflow, resulting in a wrong answer. One way to work around that is to use the identity that

(a*b)%c 

is equivalent to

(a%c * b%c)%c

This will prevent overflows in the intermediate multiplications as well.

回答4:

Your code line

 n = (n*n)%mod;

is repeatedly executed. As long as n is smaller than mod, this will potentially result in evaluating (mod-1)*(mod-1) at some point in time.

On input n may not be so large, but the mentioned line of code increases n in the loop.