I wrote the code below, to get the Lucas-Lehmer series up to p, for p the exponent of a Mersenne number. After checking it I found that it doesn't work for some primes p such as 11, 23, 29 etc. Any help will be very valuable!
Here's the code:
def ll_series (p):
ll_list=[4]
print 4
for i in range(1, p+1):
ll_list.append((ll_list[i-1]**2 - 2) % (2**p-1))
print(ll_list[i])
return ll_list