Calculate (a^b)%c where 0<=a,b,c<=10^18

Question

How can I calculate (a ^ b) % c, where 0 <= a, b, c <= 10^18. Here, (a ^ b) means a to the power b, not a xor b.

My current code for the problem is:

unsigned long long bigMod(unsigned long long b,
                          unsigned long long p,
                          unsigned long long m){
    if(b == 1) return b;
    if(p == 0) return 1;
    if(p == 1) return b;

    if(p % 2 == 0){  
        unsigned long long temp = bigMod(b, p / 2ll, m);
        return ((temp) * (temp) )% m;
    }else return (b  *  bigMod(b, p-1, m)) % m;
}

For this input:

a = 12345 b = 123456789 and c = 123456789012345

the expected output should be:

59212459031520

[Exponentiation by squaring](https://en.wikipedia.org/wiki/Exponentiation_by_squaring) might be a useful starting point. It's not too hard to add modulus to the formula, if you observe that `(a**b)%c == ((a%c)**b)%c)`. — Kevin, Sep 09 '15 at 17:27
Possible duplicate of [calculating a^b mod c](http://math.stackexchange.com/questions/26722/calculating-ab-mod-c). — Andrew Morton, Sep 09 '15 at 17:30
I was trying ((a%m)*(b%m))%m but may be my approach is not right. — Abu Hanifa, Sep 09 '15 at 17:31
@AndrewMorton may be I tried that approach isn't correct for this limit — Abu Hanifa, Sep 09 '15 at 17:33
possible duplicate of [Calculating pow(a,b) mod n](http://stackoverflow.com/questions/8496182/calculating-powa-b-mod-n) — phuclv, Sep 09 '15 at 17:37
[Calculating (a^b)%MOD](http://stackoverflow.com/q/11272437/995714). For a and b that are 64 bits, use [`__int128`](https://gcc.gnu.org/onlinedocs/gcc/_005f_005fint128.html) for the intermediate multiplication. If it's not available, write your own 128-bit multiplication. Reference: [Fastest way to calculate a 128-bit integer modulo a 64-bit integer](http://stackoverflow.com/q/2566010/995714), [Most accurate way to do a combined multiply-and-divide operation in 64-bit?](http://stackoverflow.com/q/8733178/995714) — phuclv, Sep 09 '15 at 17:42
@AbuHanifa The approach ((a%m)^(b%m)%m fails because the exponent must be modulo'ed by phi(m), where phi is Euler's totient function. If m is prime, then phi(m)=m-1 and you could reduce as ((a%m)^(b%(m-1)))%m. — Nico Schertler, Sep 09 '15 at 19:41

score 2 · Answer 1 · answered Sep 09 '15 at 18:59

You have a problem with temp*temp (long long overflow). You can omit this problem using algorithm of fast mod power to multiply them mod m. Here You have working code:

unsigned long long bigMultiply(unsigned long long b,unsigned long long p, unsigned long long m)
{
  if(p == 0 )return b;
  if(p%2 == 0)
  {  
     unsigned long long temp = bigMultiply(b,p/2ll,m);
     return ((temp)+(temp))%m;
  }
  else
    return (b  +  bigMultiply(b,p-1,m))%m;
}    
unsigned long long bigMod(unsigned long long b,unsigned long long p, unsigned long long m)
{

  if(b == 1)
    return b;
  if(p == 0 )return 1;
  if( p == 1)return b;
  if(p%2 == 0)
  {  
     unsigned ll temp = bigMod(b,p/2ll,m);
     return bigMultiply(temp,temp,m);
  }
  else
    return (b  *  bigMod(b,p-1,m))%m;
}

Trouble with some corner cases: `bigMod(1, p, 1)` is 1 and should be 0. Suggest `if(b == 1) return b%m;` Similar for `if(p == 0 )return 1;` --> `if(p == 0 )return 1%m;` to handle `m==1` and `if( p == 1)return b;` --> `if( p == 1)return b%m;` to handle `b >= m`. — chux - Reinstate Monica, Sep 20 '17 at 11:38
I am quite certain `b * bigMod(b,p-1,m)` can readily overflow. That aside, this code is very nice yet has some limitations: `(temp)+(temp)` can still overflow when `m > ULLONG_MAX/2`, yet that is beyond OP's `c <= 10^18`. — chux - Reinstate Monica, Sep 20 '17 at 11:48

score 0 · Answer 2 · edited Sep 09 '15 at 17:44

0

I use this code in c++:

long long power(long long a, long long b, long long c)
{
    if (b==0)
    {
        return 1;
    }

    if (b % 2 == 0)
    {
        long long w = power(a, b/2, c);
        return (w*w) % c;
    }
    else
    {
        int w = power(a, b-1, c);
        return (a*w) % c;
    }
}

It has logarithmic complexity.

edited Sep 09 '15 at 17:44

phuclv

37,963
15
156
475

answered Sep 09 '15 at 17:38

Kamil

170
3
9

2

it doesn't work if a and b are both 64 bits, as `w*w` needs a 128-bit type to store – phuclv Sep 09 '15 at 17:39

Calculate (a^b)%c where 0<=a,b,c<=10^18

2 Answers2

Linked