use sage to calculate the a^p-1 mod p
power_mod(a, p-1, p)
eg:
result = power_mod(5, 322, 323)
print(result)
result1 = power_mod(2, 2008, 2009)
print(result1)
system congruences:
x = 0
while(1):
x = x+1
if( mod(x,3)==1 and mod(x,4)==2 and mod(x,5)==3):
print(x)
break
import random
def primality_test(n, k):
prime_count = 0
for _ in range(k):
a = random.randint(1, n - 1)
if gcd(a, n) == 1 and power_mod(a, n - 1, n) == 1:
prime_count += 1
probability = prime_count / k
return probability
# Number to be tested for primality
n = 323 # Replace with your desired number
# Number of repetitions
k = 2009 # Adjust as needed
probability = primality_test(n, k)
print(f"Probability that {n} is prime after {k} repetitions: {probability}")
Lecture 5:
#Find a Rabin-Miller witness a > 8 for n = 1729
n=1729
fs=factor(n-1)
fs
print(fs)
#write n-1=2^k*q
k=fs[0][1]
q=(n-1)/2^k
print(k,q)
#use two Boolean variables
#to check the two conditions
#of Rabin-Miller Theorem
for a in [9..(n-1)]:
lst=[]
check1=mod(a^q,n)==1
lst.append(check1)
#print(check1)
check2=0
for i in range(k):
x=(mod(a^(q*2^i),n)==-1)
lst.append(x)
check2=check2+x
#print(lst)
if check1+check2==0:
print(a, "is a RM witness")
Lecture 6:
Find all primes p < 100 for which 3 is a primitive root modulo p.
list = prime_range(1,100)
for p in prime_range(1,100):
for k in range(1,p-1):
if(mod(3^k,p)==1):
list.remove(p)
break
print(f"3 is a primitive root modulo {list}")
2
print("k g^k (mod p)")
for k in range(1,53):
residue = mod(2^k,53)
print(f"{k} {residue}")
Lab7:
##Problem 7: Does n = 409537 pass the Solovay-Strassen test to base a = 345678?
n = 409537
#print(is_prime(n))
a = 345678
if(gcd(a,n)==1) & is_odd(n):
if(mod(jacobi_symbol(a,n),n)==power_mod(a,(n-1)/2,n)):
print(f'{n} pass the Solovay-Strassen test to base a ={a}')
else:
print(f'{n} does not pass the Solovay-Strassen test to base a ={a}')
##Is a = 1234567345679 an Euler witness for n = 10714934881993?
a = 1234567345679
n = 10714934881993
if(gcd(a,n)==1) & is_odd(n):
if(mod(jacobi_symbol(a,n),n)==power_mod(a,(n-1)/2,n)):
print(f'{a} not an Euler witness for {n}')
else:
print(f'{a} an Euler witness for {n}')
lecture 10:
Fermat
Problem 1
n = 4717
s=ceil(sqrt(n))
print(n,s,s^2)
t2=s^2-n
t=sqrt(t2)
print(t)
while t2!=(floor(t))^2:
s=s+1
t2=s^2-n
t=sqrt(t2)
print('s=',s,'t=',t)
if s^2>100*n: #note that for odd n, the algorithm works.
break
print('a=',s+t,'b=',s-t,n==(s+t)*(s-t))
#for n=10, it does not work.
#check, if for some even numbers n, Fermat succeeds.
n = 123123
s = ceil(sqrt(n))
t2 = s^2 - n
t = sqrt(t2)
while t2!=(floor(t))^2:
s = s+1
t2 = s^2 - n
t = sqrt(t2)
print('a=',s+t,'b=',s-t,n==(s+t)*(s-t))
2021 first question
list = prime_range(3,100)
for i in list:
k = (i-1)/2
print(i,(-1)^k)
print(jacobi_symbol(i-1,i))