I have been trying to determine a shorter method for calculating gcd((a^n + b^n), abs(a-b)). I noticed that, if I were to calculate (using the above formula) say a = 100 and b = 4, starting from 1 and ending at n (a loop), at a certain point, the answer becomes constant. For a = 100, b = 4 and n = 100, I create a loop from 1 to n, and at each point, I apply the formula, the first answer (n = 1) is 8, thereafter it is 32 till when n becomes 100. For optimization, I break out of the loop once two equal consecutive numbers are found and the latest number (32 here), becomes the answer. Does anyone know a straightforward formula for calculating gcd((a^n + b^n), a-b), or better still, my primary concern, a global formula for finding (a^n + b^n)

Note:
1. 1<=a,b,n<=10^12

2. (a^n - b^n) can be simplified as <https://math.stackexchange.com/questions/406703/how-to-simplify-an-bn>. But can't find a version for (a^n + b^n)
3. Following Rory Daulton's anwser, I have implemented Exponentiation by Squaring shown below in a function

My Python Code for the above explanation follows suite:

```
a, b, n = map(int, raw_input().split()); ans = -1

if a == b:
    ans = (a**n) + (b**n)

else:
    for j in xrange(n):
        x = gcd((a**n)+(b**n),abs(a-b))

        if x != ans:
            ans = x
        else:
            break
print ans 

```

Exponentiation by Squaring
==========================

```
def pow3(x, n):
r = 1
while n:
    if n % 2 == 1:
        r *= x
        n -= 1
    x *= x
    n /= 2
return r

```