随便写了一点代码

主要是练手的..

# SAGE
def SModQ(n,q,T):
	if n.nbits()>=q:
		n-=T
	return n

def SModMers(n,q,T):
	if n.nbits()==1:
		n-=2
		return T+n
	n=n*n
	n=(n>>q)+(n&T)-2
	return SModQ(n,q,T)

def Mersenne_Test(n):
	a=ZZ(4)
	T=1<<n
	T-=1
	for i in xrange(0,n-2):
		a=SModMers(a,n,T)
	return a==0

for i in primes(10000):
	if Mersenne_Test(i):
		print "2^%d-1 is a prime"%i