ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 14 Nov 2018 10:50:59 +0100Implementing the AKS primality testhttps://ask.sagemath.org/question/44275/implementing-the-aks-primality-test/I am implementing the infamous AKS deterministic primality test using SageMath.
My code is as follows:
# AKS Primality Test
def fast_exponentation(a, n):
"""
Computes a^n fast
"""
ans = 1
while n:
if n & 1 : ans = ans * a
a = a * a
n >>= 1
return ans
# Determines whether n is a power a ^ b, b > 1
def is_perfect_power(n):
lgn = 1 + ( len( bin ( abs ( n ) ) ) - 2)
for b in range(2,lgn):
# we use binary search to check if n root b is an integer
lowa = 0
higha = n
while lowa < higha - 1:
mida = (lowa + higha) // 2
ab = fast_exponentation(mida,b)
if ab > n: higha = mida
elif ab < n: lowa = mida
else: return True # mida ^ b
return False
def AKS_primality_test(n):
if n == 1: return None
if n == 2: return True
# if n is of the form a^b for a > 1 and b > 1 return False
if is_perfect_power(n): return False
print("Passed perfect power test")
# find the smallest integer r such that
# a. gcd(n,r) > 1 or
# b. [n]_r has multiplicative order > 4 log2(n)^2
for r in xrange(2,n):
# if r shares a non trivial factor with n, n is composite
if gcd(r,n) > 1: return False
# compute the multiplicative order of n in Z mod r
ord_r = mod(n,r).multiplicative_order()
if ord_r > 4*log(n,2)^2: break
print("min r found! r={}".format(r))
if r == n: return True
# We get the ring where we will check for primality
ZnX = IntegerModRing(n)[x]
ZnX_div_f = ZnX.quo(x^r-1)
for j in range(1, 2*log(n,2)*sqrt(r)+2):
if ZnX_div_f((x+j)^n) != ZnX_div_f(x^n+j):
print("The freshman dream test failed for j={}".format(j))
return False
print("Passed final test")
return True
print("271 is prime = " + str(AKS_primality_test(271)))
For the implementation I am following the scheme in page 550 of *A Computational Introduction to Number Theory
and Algebra*, by Victor Shoup ([available online in his webpage](https://shoup.net/ntb/ntb-v2.pdf)).
![image description](/upfiles/1542188648972937.png)
However, the code above, which should test whether 271 is prime, gives us the following output:
Passed perfect power test
min r found! r=269
The freshman dream test failed for j=1
271 is prime = False
This is wrong - 271 is in fact prime. So I assume I am misunderstanding how to use quotient rings in SageMath.
Why does the equality test fail?
EDIT: as pointed out by @rburing, the line `if ZnX_div_f((x+j)^n) != ZnX_div_f(x+j^n): ` should read `if ZnX_div_f((x+j)^n) != ZnX_div_f(x^n+j): ` I have edited the code above to reflect this, but the behavior of the program is the same nonetheless.JsevillamolWed, 14 Nov 2018 10:50:59 +0100https://ask.sagemath.org/question/44275/