GCD of numbers of the form 2^i-1

Question

How do I get GCD(2^a[i]-1,2^a[j]-1) with 1<=a[x]<=100

from fractions import gcd
powj=pow(2,n[j])-1
powk=pow(2,n[k])-1
gcdjk=gcd(powj,powk)

leads to trouble for large numbers and gives a run-time error.
I can't see a pattern in the 2^i-1 values except for primes having no other factors other than 1 and themselves.

i  2^i -1
--------------
1  1 = 1
2  3 = 1,3
3  7 = 1,7
4  15 = 1,3,5,15
5  31 = 1,31
6  63 = 1,3,7,9,21,63
7  127= 1,127
8  255= 1,3,5,15,17,51,85,255

Edit: Need to solve this for numbers of the form 2^i-1 ONLY. Following is the code:

import sys
import math
from fractions import gcd

t=int(input())
for i in range(0,t):
    door=0
    c=int(input())
    n = map(int,sys.stdin.readline().split(' '))
    for j in range(0,c-1):
        for k in range(j+1,c):
            if( gcd(n[j],n[k]) == n[k]):
                powj=pow(2,n[j])-1
                powk=pow(2,n[k])-1
                gcdjk=gcd(powj,powk)
                if(gcdjk==powk):
                    door = door+1
                else:
                    door = door-gcdjk
    print (door)

Input Sample:

2
3
10 2 3
2
3 5

Constraints:

1<=T<=20
1<=ArraySize<=10^5
1<=a[i]<=100

Answer 1

Consider the binary GCD algorithm. If both operands are of the form 2ⁱ-1, it can be substantially simplified.

To start with, there are obviously no zeroes at the end in the first step, so you go straight to the loop.

In the loop, in the subtraction, you have two numbers that are of the form 2ⁱ-1, and the left hand side is bigger than the right hand side, so the subtraction just resets as many low bits in y as there are bits set in x, that is, the subtraction is equivalent to y &= ~x. The subtraction is immediately followed by shifting y right by the number of trailing zeroes in it, so you have a number of the form 2ⁱ-1 again, but popcnt(x) shorter.

From this it should be obvious that only the lengths (ie exponents) ever matter, and the identity
gcd(2^a-1, 2^b-1) = 2^{gcd(a, b)}-1 follows from it.

Answer 2

These numbers are pretty small. With Python's built-in bignum handling, they're well within the reach of the Euclidean algorithm fractions.gcd uses:

>>> fractions.gcd(2**50-1, 2**100-1)
1125899906842623L

Your error is coming from somewhere else. You may even just be timing out when you try to iterate over all pairs of numbers in a 10000-element list. There are almost 50 million such pairs. Depending on how much time you get, your algorithm may simply be too slow.

Answer 3

Here is simple way u can use euclid's algorithm to solve for powers of two without actually evaluating them :-

we need to find a%b to solve using euclids algorithm for GCD :-

a = 2^x-1 b = 2^y-1

and a>b

we need to express a = k*b + m where m < b then a%b = m

suppose k = 2^(x-y)

2^x - 1 = 2^(x-y)*(2^y-1) + m , m = 2^(x-y)-1

hence

a%b = m = 2^(x-y) -1

hence m is again in similar power of two minus 1 form hence we can apply euclids algorithm on it.

Further Analysis :-

a = 2^x-1
b = 2^y-1 

GCD(a,b) = F(x,y)

where 

F(x,y) = x         if x==y
F(x,y) = F(x-y,y)  if x > y
F(x,y) = F(x,y-x)  if y < x

From further analysis F(x,y) = GCD(x,y)

Reference :- GCD