Carlson Bimbuh
Carlson Bimbuh

Reputation: 126

GCD Exponentiation

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

  1. (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)

  2. 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

Upvotes: 0

Views: 311

Answers (2)

Rory Daulton
Rory Daulton

Reputation: 22564

I see two ways to speed up your code.

First, use the mathematical fact that

gcd(r, s) = gcd(r % s, s)

(if s is not zero). So you do not need to compute a**n + b**n in full, you just need it modulo a - b. And you can do that by finding (a**n) % (a-b) and (b**n) % (a-b) then adding those modulo a - b.

Now, find a**n by the exponentiation by squaring method. That involves a loop that executes log2(n) times. In each pass through the loop, take the remainder mod a - b to keep your numbers low and speed up the calculations.

So there is your algorithm. Find (a**n) % (a-b) and (b**n) % (a-b) by exponentiation by squaring and modulus at each step. Then add them and take the modulus again. Finally, find the GCD of that value with a - b.


In some situations, such as a - b prime, I can see some shortcuts. The modulus of the powers of a number do repeat, as you note. However, finding out just when they repeat is a non-trivial problem for large values of a - b, especially if a - b is composite and hard to factor. Unless you have some additional information on the values of a - b and your other parameters, I suggest you not use the repetitions. If the values of a and b are small and known in advance (as in your example of a = 100 and b = 4, the repetitions are more attractive and you could pre-compute the values of powers modulus 96.


Rather than use this code, you should probably use Python's built-in pow function. See here for documentation. Hat tip to @DSM.

By request, here is my routine for exponentiation by squaring modulo a given number. There are some variations that could be done, of course. This version does no error checking on the parameters an does a little bit-twiddling for minor efficiencies.

def exp_by_squaring_mod(x, n, mod):
    """Calculate x**n % mod by squaring. This assumes x and n are non-
    negative integers and mod is a positive integer. This returns 1 for
    0**0.
    """
    result = 1
    x %= mod
    # Reduce n and keep constant the value of result * x**n % mod
    while n:  # while n is not zero
        if n & 1:  # n is odd
            result = result * x % mod
        x = x * x % mod
        n >>= 1  # integer divide by 2
    return result

Upvotes: 2

hdu
hdu

Reputation: 31

Try to optimize your exponentiation using this method:

def powerBySquaring(x,n):
    if n < 0:
        x = 1 / x
        n = -n
    if n == 0: return 1
    y = 1
    while n > 1:
        if n % 2 == 0:
            x = x * x
            n = n / 2
        else:
            y = x * y
            x = x * x
            n = (n-1)/2
    return x * y

Upvotes: 0

Related Questions