Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Primality test for generalized Fermat numbers
'''
This script was written in python 3.x.
In order to run this script, please make sure your python version is 3.x or above.
How to run:
python fermat.py
or if it doesn't work use this one:
python3 fermat.py
Author: Pedja Terzic <pedja.terzic@hotmail.com>
'''
from mpmath import *
print(" ***** FERMAT *****\n\n\n")
while True:
b=int(input("Enter the base : "))
n=int(input("Enter the exponent : "))
F=b**2**n+1
def polynomial(m,x):
if m==1:
return x
else:
p0=2
p1=x
l=2
while l<=m:
p=x*p1-p0
p0=p1
p1=p
l=l+1
return p
def jacobi(a,q):
j=1
while a != 0:
while a%2==0:
a=a/2
if q%8==3 or q%8==5:
j=-j
#interchange(a,q)
c=a
a=q
q=c
if a%4==3 and q%4==3:
j=-j
a=fmod(a,q)
if q==1:
return j
else:
return 0
if n < 2:
print("Exponent must be greater than one")
elif b%2==1:
print("Base must be an even number")
else:
d=3
while not(jacobi(d-2,F)==-1 and jacobi(d+2,F)==-1):
d=d+1
s=polynomial(b/2,polynomial(b/2,d))%F
ctr=1
while ctr<=2**n-2:
s=polynomial(b,s)%F
ctr=ctr+1
if s==0:
print(str(b) + "^2^" + str(n) + "+1 is prime")
else:
print(str(b) + "^2^" + str(n) + "+1 is composite")
try_again = ""
# Loop until users opts to go again or quit
while not(try_again == "1") and not(try_again == "0"):
try_again = input("Press 1 to try again, 0 to exit. ")
if try_again in ["1", "0"]:
continue # a valid entry found
else:
print("Invalid input- Press 1 to try again, 0 to exit.")
# at this point, try_again must be "0" or "1"
if try_again == "0":
break
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.