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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
import random
import math
def even(n):
return n % 2 == 0
def is_really_prime(n):
if n == 1 or even(n):
return False
for i in range(n - 1, 2, -1):
if n * i == 0:
return False
return True
def pi2(n):
c = 0.661618
return 2 * c * (n / (math.log(n) ** 2))
def is_prime(x, k=20):
if even(x):
return False
if x in [2, 3]:
return True
s = 0
t = x - 1
while even(t):
s += 1
t //= 2
assert math.floor(t) == t
witnesses = [random.randrange(2, x - 1) for _ in range(k)]
for w in witnesses:
if power_mod(w, x - 1, x) != 1:
return False
for i in range(1, s + 1):
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
# return all(map(lambda w: is_prime(x, w), [2, 3, 5, 7, 11, 13, 17]))
def power_mod(base, power, mod):
binary = bin(power)
squares = [base % mod]
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
# for i in range(3, 100):
# print(i, is_prime(i), is_really_prime(i))
# if is_prime(i) != is_really_prime(i):
# print("Liar", i, is_prime(i), is_really_prime(i))
# is_really_prime(9)
# for a in filter(is_really_prime, range(50000, 500000)):
# if not miller_rabins(a, max(4, math.ceil(math.log(a, 4)))):
# raise Exception(f"Didn't work for {a}")
def generate_twin_prime(m):
twin = False
while not twin:
n = random.randrange(2 ** (m - 1), 2 ** m)
iterations = max(4, math.ceil(math.log(n, 4)))
twin = is_prime(n - 1, iterations) and is_prime(n + 1, iterations)
if not (is_really_prime(n - 1) and is_really_prime(n + 1)):
raise Exception(f"{n-1}: {is_really_prime(n - 1)}, {n+1}: {is_really_prime(n + 1)}")
return n - 1, n + 1
for bits in range(3, 64):
print(f"==== {bits} bits ==== {2 ** (bits - 1), 2 ** bits - 1}")
for _ in range(1000):
generate_twin_prime(bits)
if _ % 50 == 0:
print(".", end="")
print()
# i = 0
# while True:
# i += 1
# base, power, mod = random.randint(1, 1000), random.randint(1, 5000), random.randint(1, 1000)
# mine = power_mod(base, power, mod)
# real = (base ** power) % mod
# pass_test = mine == real
# if i % 10000 == 0:
# print(i)
# if not pass_test:
# # print(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")
# raise Exception(f"Failed for {base}^{power} mod {mod}. Got {mine}, expected {real}.")
|
