Modular exponentiation fails for large mod in C++
Asked Answered
A

4

4

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.

Allseed answered 22/10, 2015 at 5:33 Comment(1)
Use a multi-precision library like GMP.Davis
H
6

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;
}
Hom answered 23/10, 2015 at 8:28 Comment(1)
Amazing. Kudos to your approach.Allseed
B
1

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.

Blunderbuss answered 22/10, 2015 at 5:52 Comment(5)
my n and p, are both less than 2*10^6. So, the product should never be a problem I guess. Am I wrong somewhere?Allseed
ans = (ans*n)%mod; and n = (n*n)%mod;, n and ans may be larger than 2*10^6Blunderbuss
Does this mean that Project Euler #48 on HackerRank (hackerrank.com/contests/projecteuler/challenges/euler048) cannot be solved in C++?Allseed
@saint1729 You can implement BigInteger yourself in C++, but I suggest to use pythonBlunderbuss
@saint1729 You don't need to implement BigInteger for this specific scenario. See my answer.Hom
D
0

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.

Donohoe answered 22/10, 2015 at 5:39 Comment(1)
I there. I tried it. In my case a*b <= 4*10^12 < 2^64-1 and so that is not giving me any advantage.Allseed
R
0

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.

Reposeful answered 22/10, 2015 at 6:30 Comment(0)

© 2022 - 2025 — McMap. All rights reserved.