Skip to content

Instantly share code, notes, and snippets.

@primus-lab
Created October 30, 2019 08:07
Show Gist options
  • Save primus-lab/10b0425b4ed01e149db28ca964830499 to your computer and use it in GitHub Desktop.
Save primus-lab/10b0425b4ed01e149db28ca964830499 to your computer and use it in GitHub Desktop.
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