Modular exponentiation fails for large mod in C++

Question

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.

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