forked from sethtroisi/gmp-ecm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
check.sage
118 lines (109 loc) · 3.51 KB
/
check.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
def FindGroupOrder(p,s):
K = GF(p)
v = K(4*s)
u = K(s^2-5)
x = u^3
b = 4*x*v
a = (v-u)^3*(3*u+v)
A = a/b-2
x = x/v^3
b = x^3 + A*x^2 + x
E = EllipticCurve(K,[0,b*A,0,b^2,0])
return factor(E.cardinality())
def FindGroupOrderA(p,A):
K = GF(p)
d = K((A+2)/4)
a = K(4*d-2)
b = K(16*d+2)
E = EllipticCurve(K,[0,a/b,0,1/b^2,0])
return factor(E.cardinality())
# for parameter sigma = 1:s
def FindGroupOrderParam1(p,s):
return FindGroupOrderA (p, 4*s^2/2^64-2)
# for parameter sigma = 2:s
def FindGroupOrderParam2(p,s):
K = GF(p)
E = EllipticCurve(K,[0,36])
P = s*E(-3,3)
x,y = P.xy()
x3 = (3*x+y+6)/(2*(y-3))
A = -(3*x3^4+6*x3^2-1)/(4*x3^3)
d = K((A+2)/4)
a = K(4*d-2)
b = K(16*d+2)
E = EllipticCurve(K,[0,a/b,0,1/b^2,0])
return factor(E.cardinality())
# for parameter sigma = 3:s
def FindGroupOrderParam3(p,s):
return FindGroupOrderA (p, 4*s/2^32-2)
def FindGroupOrderParam (p, sigma, param):
if param == 0:
return FindGroupOrder (p, sigma)
elif param == 1:
return FindGroupOrderParam1 (p, sigma)
elif param == 2:
return FindGroupOrderParam2 (p, sigma)
elif param == 3:
return FindGroupOrderParam3 (p, sigma)
else:
print "Invalid parametrization: ", param
raise ValueError
# check if the prime p is found with B1,B2,param,sigma, or raises an error
# check_found_aux ("./ecm", 31622776601683800097, 11000, 1873422, 1, 800667805)
# check_found_aux ("./ecm", 31622776601683800097, 11000, 1873422, 1, 800667806)
def check_found_aux (ecm, p, B1, B2, param, sigma):
f = open("/tmp/inxyz", "w")
f.write(str(p) + "\n")
f.close()
f = open("/tmp/doitxyz", "w")
f.write(ecm + " -param " + str(param) + " -sigma " + str(sigma) + " " + str(B1) + " " + str(B2) + " < /tmp/inxyz > /tmp/outxyz\n")
f.close()
os.system("chmod +x /tmp/doitxyz")
os.system("/tmp/doitxyz")
f = open("/tmp/outxyz", "r")
l = f.readlines()
f.close()
n = len(l)
if l[n-1] <> 'Found input number N\n':
print "prime p=", p, "not found with B1=", B1, "B2=", B2, "param=", param, "sigma=", sigma
raise ValueError
def is_found(l, B1, B2):
n = len(l)
if l[n-1][0] > B2:
return False
for i in range(n-2,-1,-1):
if l[i][0]^l[i][1] > B1:
return False
return True
# check if a prime p is found with bounds B1 and B2,
# for parameter 'param' and sigma in [sigma_min,sigma_max-1]
# check_found ("./ecm", 31622776601683800097, 11000, 1873422, 0, 1000)
# check_found ("./ecm", 31622776601683800097, 11000, 1873422, 1, 1000)
# check_found ("./ecm", 31622776601683800097, 11000, 1873422, 2, 1000)
# check_found ("./ecm", 31622776601683800097, 11000, 1873422, 3, 1000)
def check_found (ecm, p, B1, B2, param, sigma_max):
assert (is_prime (p))
e2 = 0
e3 = 0
tested = 0
tries = 0
for sigma in range(sigma_max):
try:
l = FindGroupOrderParam (p, sigma, param)
except ArithmeticError:
continue
tries += 1
assert (l[0][0] == 2)
e2 += l[0][1]
if l[1][0] == 3:
e3 += l[1][1]
if is_found (l, B1, B2):
tested += 1
# check the factor is really found
check_found_aux (ecm, p, B1, B2, param, sigma)
print tries, tested, 1.0*e2/tries, 1.0*e3/tries, 2.0^(e2/tries)*3.0^(e3/tries)
# check all parametrizations 0, 1, 2, 3
#check_found_all ("./ecm", 31622776601683800097, 11000, 1873422, 200)
def check_found_all (ecm, p, B1, B2, sigma_max):
for param in range(4):
check_found (ecm, p, B1, B2, param, sigma_max)