1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
import random
import math
import sys
def even(n):
return n % 2 == 0
# Use the Miller-Rabin algorithm to test primality of x
def is_prime(x, k):
# Base cases
if even(x):
return False
if x in [2, 3]:
return True
# Determine s and t
s = 0
t = x - 1
while even(t):
s += 1
t //= 2
# Gather witnesses
witnesses = [random.randrange(2, x - 1) for _ in range(k)]
for w in witnesses:
if power_mod(w, x - 1, x) != 1: # Basic Fermat test
return False
for i in range(1, s + 1): # Advanced test if Fermat fails
if power_mod(w, 2 ** i * t, x) == 1 and power_mod(w, 2 ** (i - 1) * t, x) not in [1, x - 1]:
return False
return True # All witnesses must agree to get here
# Use repeated-squaring method to determine: base^power % mod
def power_mod(base, power, mod):
binary = bin(power)
squares = [base % mod] # Start with base case
result = squares[0] if binary[-1] == "1" else 1
for i in range(1, math.floor(math.log2(power)) + 1):
squares.append(squares[i - 1] ** 2 % mod)
if binary[-i - 1] == "1":
result *= squares[i]
return result % mod
# Generate a twin prime of bit length m. Only tests values of 6x±1.
def generate_twin_prime(m):
twin = False
while not twin:
start, end = bit_range(m)
n = random.randrange(start + 6 - (start % 6), end, 6) # Only test m bit ints that are multiples of 6
num_witnesses = max(4, math.ceil(math.log(n, 4))) # Based on probability 4^{-k} that a witness will be wrong.
twin = is_prime(n - 1, num_witnesses) and is_prime(n + 1, num_witnesses) # Check twin on either side of 6x
return n - 1, n + 1
# Tuple containing range for bit length m
def bit_range(m):
return 2 ** (m - 1), 2 ** m
# Write a twin in the correct format to file
def write_twin(twin):
with open("output_twin_prime.txt", "w") as file:
file.write(f"{twin[0]}\n{twin[1]}")
def main():
assert len(sys.argv) >= 2
twin = generate_twin_prime(int(sys.argv[1]))
write_twin(twin)
print(twin)
if __name__ == "__main__":
main()
|
