-
Notifications
You must be signed in to change notification settings - Fork 0
/
generate_pseudoprimes.sage
executable file
·154 lines (130 loc) · 4.54 KB
/
generate_pseudoprimes.sage
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#!/usr/bin/env sage
import random
import sys
proof.arithmetic(False)
def millerrabinliar(n,a):
r = 0
m = n-1
while m % 2 == 0:
r += 1
m = m/2
if Mod(a,n)^m == Mod(1,n):
return True
for i in range(r):
if Mod(a,n)^(m*2^i) == Mod(-1,n):
return True
return False
# returns three primes whose product is a pseudoprime
def arnault3(bases, k2, k3):
try:
(ks_ok, list_of_intersections) = check_ks(bases, k2, k3)
except:
return None
if not ks_ok:
return None
while True:
p1 = construct_p1(list_of_intersections, k2, k3)
if not is_prime(p1):
continue
print(".", file=sys.stderr, end="", flush=True)
p2 = k2*(p1-1)+1
if not is_prime(p2):
continue
print(";", file=sys.stderr, end="", flush=True)
p3 = k3*(p1-1)+1
if is_prime(p3):
print("#", file=sys.stderr, end="", flush=True)
break
# try to factor, p1-1, p2-1, p3-1, if any has a factor of correct size, then output
n = p1*p2*p3
assert ((n-1)%(p1-1) == 0 and (n-1)%(p2-1) == 0 and (n-1)%(p3-1) == 0 and is_prime(p1) and is_prime(p2) and is_prime(p3))
return n, p1, p2, p3
def compute_Sa_intersection(a, k2, k3):
Sa = []
for x in range(4*a):
if kronecker(a,x) == -1 and gcd(x,4*a) == 1 and x%4 == 3:
Sa.append(x)
k2i = k2.inverse_mod(4*a)
Sa2 = [k2i*(s+k2-1)%(4*a) for s in Sa]
k3i = k3.inverse_mod(4*a)
Sa3 = [k3i*(s+k3-1)%(4*a) for s in Sa]
return (Set(Sa).intersection(Set(Sa2))).intersection(Set(Sa3)).list()
def check_ks(bases, k2, k3):
if k2 == k3:
return False, []
list_of_intersections = []
for a in bases:
intersection = compute_Sa_intersection(a, k2, k3)
if intersection == []:
return False, []
list_of_intersections.append((a, intersection))
return True, list_of_intersections
def construct_p1(list_of_intersections, k2, k3):
p1 = crt(-1* k3.inverse_mod(k2), -1 * k2.inverse_mod(k3), k2, k3)
modulus = k2*k3
p1 = crt(p1, 3, modulus, 8)
modulus = lcm(4, modulus)
for a, intersection in list_of_intersections:
if a % 2 == 0:
continue
z = random.choice(intersection)
p1 = crt(p1, z, modulus, 4 * a)
modulus = modulus * a
return p1
base_k2 = 10
base_k3 = 10
#heuristically try specific baseBound, k2, k3 based on target bitsize
def generate_pseudoprime_near_given_size(bitsize, optimize=False):
global base_k2, base_k3
assert bitsize in [160, 192, 224, 256, 384, 512, 521, 1024, 2048, 4096]
if bitsize == 160:
baseBound, k2, k3 = 11, 73, 101
elif bitsize == 192:
baseBound, k2, k3 = 13, 61, 101
elif bitsize == 224:
baseBound, k2, k3 = 14, 197, 257
elif bitsize == 256:
baseBound, k2, k3 = 16, 233, 101
elif bitsize == 384:
baseBound, k2, k3 = 23, 137, 157
elif bitsize == 512:
baseBound, k2, k3 = 30, 137, 157
elif bitsize == 521:
baseBound, k2, k3 = 30, 137, 157
elif bitsize == 1024:
baseBound, k2, k3 = 52, 241, 281
elif bitsize == 2048:
baseBound, k2, k3 = 92, 673, 733
elif bitsize == 4096:
baseBound, k2, k3 = 164, 977, 997
if optimize:
k2, k3 = base_k2, base_k3
base_k3 += 1
if base_k3 >= 500:
base_k2 += 1
base_k3 = base_k2
print("Doing", baseBound, k2, k3, file=sys.stderr)
bases = primes_first_n(baseBound)
result = arnault3(bases, k2, k3)
return result, baseBound, k2, k3
def generate_pseudoprime_of_exactly_given_size(bitsize, optimize=False):
while True:
arnault, baseBound, k2, k3 = generate_pseudoprime_near_given_size(bitsize, optimize)
if arnault == None:
continue
print(arnault[0].nbits(), arnault[1].nbits(), arnault[2].nbits(), arnault[3].nbits(), file=sys.stderr)
if arnault[0].nbits() == bitsize:
break
return arnault, baseBound, k2, k3
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser()
parser.add_argument("bits", type=int)
parser.add_argument("count", type=int)
parser.add_argument("--optimize", action="store_true", default=False)
args = parser.parse_args()
#write pseudoprime data ppCount times
for i in range(args.count):
arnault, baseBound, k2, k3 = generate_pseudoprime_of_exactly_given_size(args.bits, args.optimize)
pp, p1, p2, p3 = arnault
print(",".join([str(i), str(pp), str(p1)]), p2, p3)