[TẶNG BẠN] TRỌN BỘ Bí kíp học tốt 08 môn
Chắc suất Đại học top - Giữ chỗ ngay!!
ĐĂNG BÀI NGAY để cùng trao đổi với các thành viên siêu nhiệt tình & dễ thương trên diễn đàn.
Cách đây n tháng mình có làm code kiểm tra tính nguyên tố của một số bằng phương pháp Baillie–PSW (Baillie–PSW Primality Test). Cái này mình làm khá lâu rồi nên cũng chả nhớ mình code có đúng không:>
Mọi người có thể chạy thử nếu muốn:v (Lưu ý đây là Python nhé)
P/S: Mình comment khá nhiều trong code để tránh bị loạn:v
Mọi người có thể chạy thử nếu muốn:v (Lưu ý đây là Python nhé)
P/S: Mình comment khá nhiều trong code để tránh bị loạn:v
Mã:
from math import *
from numpy import *
def baillie_psw(n: int): # Baillie–PSW primality test
"""
The Baillie–PSW primality test is a probabilistic primality testing algorithm that determines
whether a number is composite or is a probable prime.
The Baillie-PSW test is a combination of
a strong Fermat probable prime test to base 2 and a strong Lucas probable prime test.
Let n be the odd positive integer that we wish to test for primality.
1. Optionally, perform trial division to check if n is divisible by a small prime number less than some convenient limit.
2. Perform a base 2 strong probable prime test. If n is not a strong probable prime base 2, then n is composite; quit.
3. Find the first D in the sequence 5, −7, 9, −11, 13, −15, ... for which the Jacobi symbol (D/n) is −1.
Set P = 1 and Q = (1 − D) / 4.
4. Perform a strong Lucas probable prime test on n using parameters D, P, and Q.
If n is not a strong Lucas probable prime, then n is composite. Otherwise, n is almost certainly prime.
"""
# 1. Optionally, perform trial division
# to check if n is divisible by a small prime number less than some convenient limit before using this test.
if not quick_trial_division(n):
return False
# 2. Perform a base 2 strong probable prime test. (This is actually Miller-Rabin primality test using base 2)
# If n is not a strong probable prime base 2, then n is composite; quit.
if not miller_rabin_base_2(n):
return False
# 3. Find the first D in the sequence 5, −7, 9, −11, 13, −15, ... for which the Jacobi symbol (D/n) is −1.
# Set P = 1 and Q = (1 − D) / 4.
# If n is a perfect square, then step 3 will never yield a D with (D/n) = −1.
# So we need to check whether n is a perfect square or not before starting step 3
sqrt_of_n = int(sqrt(n))
if sqrt_of_n * sqrt_of_n == n:
return False
D, P, Q = D_finder(n)
# 4. Perform a strong Lucas probable prime test on n using parameters D, P, and Q.
# If n is not a strong Lucas probable prime, then n is composite.
# Otherwise, n is almost certainly prime.
if not strong_lucas_probable_prime(n, D, P, Q):
return False
return True
def quick_trial_division(n):
"""
This one is for Baillie–PSW primality test only
This will only check if n is divisible by a small prime number (less than 100),
so don't use this as the main test.
"""
if n % 2 == 0:
return False
f = 3
while f * f <= n:
if n % f == 0:
return False
elif f > 997:
break
else:
f += 2
return True
def miller_rabin_base_2(n):
"""
Just the same as original Miller-Rabin Primality Test,
except for the value of a will always be 2.
This is used for Baillie–PSW primality test only
"""
if n == 2 or n == 3:
return True
if n % 2 == 0:
return False
s = n - 1
d = 0
r = 0
while True:
if s % 2 == 0:
r += 1
s //= 2
else:
d = s
break
x = (2**d) % n
if x == 1 or x == n-1:
return True
for _ in range(r-1): # Repeat r-1 times
t = 0
x = (x * x) % n
if x == n-1:
return True
return False
def D_finder(n):
"""
Find the first D in the sequence 5, −7, 9, −11, 13, −15, ... for which the Jacobi symbol (D/n) is −1.
Then let P = 1 and Q = (1 − D) / 4.
"""
D = 5
while True:
if jacobi_symbol(D, n) == -1:
P = 1
Q = (1-D)/4
return D, P, Q
else:
if D > 0:
D = (-1)*D - 2
else:
D = (-1)*D + 2
def strong_lucas_probable_prime(n, D, P, Q):
"""
Perform a strong Lucas probable prime test on n using parameters D, P, and Q.
"""
delta_n = n + 1 # It's actually n - (D/n), but since (D/n) = −1, we can just assign delta_n to n+1
# Now, factor delta_n into the form d * 2^s, where is odd
s = 0
while True:
if delta_n % 2 == 0:
s += 1
delta_n //= 2
else:
d = delta_n
break
# A strong Lucas pseudoprime for a given (P, Q) pair is an odd composite number n with GCD(n, D) = 1,
# satisfying one of the conditions
# U_d ≡ 0 (mod n) or V_{d*2^r} ≡ 0 (mod n)
# for 0 <= r < s
U_d = U_V_finder(n, d)
if U_d % n == 0:
return True
for r in range(s):
delta_prime = (1 << r) * d
V_delta = U_V_finder(n, delta_prime)
if V_delta % n == 0:
return True
return False
def U_V_finder(n, t, D, P, Q):
"""
Find U_{n+1} and V_{n+1}
Use for strong lucas probable prime test only
"""
k = 0
# Let U_{2k}, V_{2k}, U_{2k+1} and V_{2k+1} be an array in the form of [a, b]
# where a denotes the subscript and b is their value
U_even = [2*k, 0] # U_{2k}
V_even = [2*k, 0] # V_{2k}
U_temp = [] # These are made for calculation of U_{2k} and V_{2k}
V_temp = [] # since they depend on U_{k} and V_{k}
U_odd = [2*k+1, 1] # U_{2k+1}
V_odd = [2*k+1, 1] # V_{2k+1}
t_range = (t // 2) + 1
for k in range(1, t_range):
if k == 1:
U_even[1] = (U_odd[1] * V_odd[1]) % n
V_even[1] = (((V_odd[1] * V_odd[1]) % n) - ((2*(Q**k)) % n)) % n
else:
U_even[1] = (U_temp[k-2] * V_temp[k-2]) % n
V_even[1] = (((V_temp[k-2] * V_temp[k-2]) % n) - ((2*(Q**k)) % n)) % n
if (P*U_even[1] + V_even[1]) % n == 0:
U_odd[1] = (P*U_even[1] + V_even[1])/2
else:
U_odd[1] = (P*U_even[1] + V_even[1] + n)/2
V_odd[1] = (D*U_even[1] + P*V_even[1])/2
U_temp.append(U_even[1])
U_temp.append(U_odd[1])
V_temp.append(V_even[1])
V_temp.append(V_odd[1])
if t % 2 == 1:
return U_odd
else:
return V_even
def jacobi_symbol(k, n): #Jacobi Symbol
"""
Jacobi Symbol (k/n)
Return -1, 0 or -1
"""
if n % 2 == 0:
raise ValueError("'n' must be odd.")
k %= n
result = 1
while k != 0:
while k % 2 == 0:
k /= 2
n_mod_8 = n % 8
if n_mod_8 in (3, 5):
result = -result
k, n = n, k
if k % 4 == 3 and k % 4 == 3:
result = -result
k %= n
if n == 1:
return result
else:
return 0