Skip to content

Instantly share code, notes, and snippets.

@primus-lab
Created January 3, 2020 09:23
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save primus-lab/707ef2b98f74d665f6ca503a62be2c11 to your computer and use it in GitHub Desktop.
Save primus-lab/707ef2b98f74d665f6ca503a62be2c11 to your computer and use it in GitHub Desktop.
Compositeness test for numbers of the form N=k*b^n+c
# Author: Pedja
from mpmath import *
from sympy import *
print(" ***** LUX *****\n\n\n")
while True:
k=int(input("Enter the coefficient : "))
b=int(input("Enter the base : "))
n=int(input("Enter the exponent : "))
c=int(input("Enter the constant : "))
def polynomial(m,x):
if m==0:
return 2
elif 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 test(k,b,n,c):
N=k*b**n+c
d=3
while not(jacobi_symbol(d-2,N)==c//abs(c) and jacobi_symbol(d+2,N)==1):
d=d+1
s=polynomial(k*b,d)%N
ctr=1
while ctr<=n-1:
s=polynomial(b,s)%N
ctr=ctr+1
if int(s)==polynomial(abs(c)-1,d):
return "probably prime"
else:
return "composite"
if k<1:
print("Coefficient must be natural number greater than zero")
elif b<2:
print("Base must be natural number greater than one")
elif n<2:
print("Exponent must be natural number greater than one")
elif c==0:
print("Constant must be nonzero integer")
elif c>0:
print(str(k)+"*"+str(b)+"^"+str(n)+"+"+str(c)+" is "+test(k,b,n,c))
else:
print(str(k)+"*"+str(b)+"^"+str(n)+str(c)+" is "+test(k,b,n,c))
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