Why this modular GCD fails on large inputs?

Question

This question actually comes from a Competitive Programming site called codechef.

The question is as follow.

Given integers A, B and N, you should calculate the GCD of A^N+B^N and |A−B|. (Assume that GCD(0,a)=a for any positive integer a). Since this number could be very large, compute it modulo 1000000007 (109+7).

The complete question is here

Constraints on values are as follow :

1 <= A,B,N <= 1e9 (Sub task #1)

1 <= A,B,N <= 1e12 (Sub task #2)

B <= A always

My Solution passes the first Subtask but Fails on Second when Values of A, B, N are very large.

Here is my solution :

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

Request

It is not a Homework nor I want a precise answer or code correction. I do understand all these but cannot understand why would it fail on larger values?

Any direction or hint would be useful.

Answer 1

SPECIAL THANKS TO SAJIB for Pointing out the Issue.

Here I will answer it myself for better understanding.

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

Are the two lines that breaks and causes overflow. All credits to @sajib.

He also pointed out the correct way to solve this issue (Multiplication with summations). I had great time understanding what his code does.

---

So here I explain it in detail.

(a * b) % m

will cause overflow if both a and b are very large long long values.

A naive approach is to use arbitrary precision data type such as int in python or Biginteger class in Java. But that approach will not be fruitful because internal conversion of string to int and then perform operation will lead to slow down the calculations of addition and multiplications in binary number system.

Efficient solution: Since a and b may be very large numbers, if we try to multiply directly then it will definitely overflow. Therefore we use the basic approach of multiplication i.e.,

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

So we can easily compute the value of addition (under modulo m) without any overflow in the calculation. But if we try to add the value of a repeatedly up to b times then it will definitely timeout for the large value of b, since the time complexity of this approach would become O(b).

So we divide the above-repeated steps of a in simpler way i.e.,

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

Implementation is as follow :

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

Answer 2

If modder = 1e12 then your modulo not work. cause 1e12 * 1e12. there will be overflow.

Look at this to modulo without overflow.

You can try with this. Here multiplication done by summation.

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