Commit | Line | Data |
---|---|---|
6775a491 MW |
1 | #! /usr/local/bin/sage |
2 | ### -*-python-*- | |
3 | ||
4 | from itertools import count, izip | |
5 | from random import randrange | |
6 | ||
7 | ## Constants. | |
8 | EC_YBIT = EC_XONLY = 1 | |
9 | EC_CMPR = EC_LSB = 2 | |
10 | EC_EXPLY = 4 | |
11 | EC_SORT = 8 | |
12 | ||
13 | ## Conversion primitives. | |
14 | def poly_to_int(x, p): | |
15 | return ZZ(sum([ZZ(a)*p**i for a, i in izip(x.list(), count())])) | |
16 | def int_to_poly(n, p, t): | |
17 | return sum([a*t**i for a, i in izip(n.digits(p), count())]) | |
18 | ||
19 | def fe2ip(x): | |
20 | return poly_to_int(x.polynomial(), x.parent().characteristic()) | |
21 | def i2fep(n, k): | |
22 | return int_to_poly(n, k.characteristic(), k.gen()) | |
23 | ||
24 | def fe2osp(x): | |
25 | k = x.parent() | |
26 | n = fe2ip(x) | |
27 | nb = ceil((k.cardinality() - 1).nbits()/8) | |
28 | v = n.digits(256, padto = nb); v.reverse(); return v | |
29 | ||
30 | def os2sp(v): | |
31 | return ''.join('%02x' % i for i in v) | |
32 | ||
33 | def fe2sp(x): | |
34 | if x == 0 or x == 1: | |
35 | return '%d' % x | |
36 | if x.parent().degree() > 1: | |
37 | return '0x%x' % fe2ip(x) | |
38 | else: | |
39 | n = ZZ(x + 99) | |
40 | if n < 199: return '%d' % (n - 99) | |
41 | else: return '0x%x' % x | |
42 | ||
43 | def os_to_hex(v): | |
44 | return ''.join('%02x' % b for b in v) | |
45 | ||
46 | ## Other utilities. | |
47 | def find_2torsion_point(E, n): | |
48 | while is_even(n): n //= 2 | |
49 | while True: | |
50 | Q = E.random_point() | |
51 | R = n*Q | |
52 | if R: break | |
53 | while True: | |
54 | S = 2*R | |
55 | if not S: break | |
56 | R = S | |
57 | return R | |
58 | ||
59 | def pick_at_random(choices): | |
60 | return choices[randrange(len(choices))] | |
61 | ||
62 | ## Output an elliptic curve description. | |
63 | def ecdescr(E): | |
64 | k = E.base_ring() | |
65 | m = k.degree() | |
66 | if m == 1: | |
67 | fdesc = '%s: %d' % (pick_at_random(['prime', 'niceprime']), | |
68 | k.characteristic()) | |
69 | cdesc = '%s: %s, %s' % (pick_at_random(['prime', 'primeproj']), | |
70 | fe2sp(E.a4()), fe2sp(E.a6())) | |
71 | elif k.characteristic() == 2: | |
72 | fdesc = '%s: 0x%x' % ('binpoly', poly_to_int(k.polynomial(), 2)) | |
73 | cdesc = '%s: %s, %s' % (pick_at_random(['bin', 'binproj']), | |
74 | fe2sp(E.a2()), fe2sp(E.a6())) | |
75 | else: | |
76 | raise ValueError, 'only prime or binary fields supported' | |
77 | return '"%s; %s"' % (fdesc, cdesc) | |
78 | ||
79 | ## Output a point description. | |
80 | def ecpt(Q): | |
81 | if Q is None: return 'FAIL' | |
82 | elif not Q: return 'inf' | |
83 | else: return '"%s, %s"' % (fe2sp(Q[0]), fe2sp(Q[1])) | |
84 | ||
85 | ## Compress a point. | |
86 | def ptcompr_lsb(Q): | |
87 | x, y = Q[0], Q[1] | |
88 | if Q.curve().base_ring().characteristic() != 2: | |
89 | return is_odd(fe2ip(y)) and EC_YBIT or 0 | |
90 | elif x == 0: | |
91 | return 0 | |
92 | else: | |
93 | return (y/x).polynomial()[0] == 1 and EC_YBIT or 0 | |
94 | ||
95 | def ptcompr_sort(Q): | |
96 | y = fe2ip(Q[1]); yy = fe2ip((-Q)[1]) | |
97 | return y > yy and EC_YBIT or 0 | |
98 | ||
99 | ## Some elliptic curves. Chosen to have 2-torsion, so as to test edge cases | |
100 | ## where Q = -Q; in particular, 2-torsion points can't have their | |
101 | ## y-coordinates compressed to 1. | |
102 | p = previous_prime(2^192); k_p = GF(p) | |
103 | E_p = EllipticCurve(k_p, [-3, 6]) | |
104 | n_p = 6277101735386680763835789423293484020607766684840576538738 | |
105 | T_p = find_2torsion_point(E_p, n_p) | |
106 | ||
107 | k_2.<t> = GF(2^191) | |
108 | E_2 = EllipticCurve(k_2, [1, 1, 0, 0, 1]) | |
109 | n_2 = 3138550867693340381917894711593325610042432240066118385966 | |
110 | T_2 = find_2torsion_point(E_2, n_2) | |
111 | ||
112 | def test_ec2osp(E, Q): | |
113 | ec = ecdescr(E) | |
114 | pt = ecpt(Q) | |
115 | xs, ys = fe2osp(Q[0]), fe2osp(Q[1]) | |
116 | ybit_lsb = ptcompr_lsb(Q) | |
117 | ybit_sort = ptcompr_sort(Q) | |
118 | ||
119 | def out(f, s): | |
120 | print ' %s\n %d %s\n %s;' % (ec, f, pt, os_to_hex(s)) | |
121 | ||
122 | out(EC_XONLY, [EC_XONLY] + xs) | |
123 | out(EC_CMPR, [EC_CMPR | ybit_lsb] + xs) | |
124 | out(EC_CMPR | EC_SORT, [EC_CMPR | EC_SORT | ybit_sort] + xs) | |
125 | out(EC_EXPLY, [EC_EXPLY] + xs + ys) | |
126 | out(EC_CMPR | EC_EXPLY, [EC_CMPR | EC_EXPLY | ybit_lsb] + xs + ys) | |
127 | out(EC_CMPR | EC_SORT | EC_EXPLY, | |
128 | [EC_CMPR | EC_SORT | EC_EXPLY | ybit_sort] + xs + ys) | |
129 | ||
130 | print 'ec2osp {' | |
131 | for i in xrange(20): test_ec2osp(E_p, E_p.random_point()) | |
132 | test_ec2osp(E_p, T_p) | |
133 | for i in xrange(20): test_ec2osp(E_2, E_2.random_point()) | |
134 | test_ec2osp(E_2, T_2) | |
135 | print '}' | |
136 | ||
137 | def test_os2ecp(E, Q): | |
138 | ec = ecdescr(E) | |
139 | k = E.base_ring() | |
140 | m = k.degree() | |
141 | x, y = Q[0], Q[1] | |
142 | xs, ys = fe2osp(Q[0]), fe2osp(Q[1]) | |
143 | ybit_lsb = ptcompr_lsb(Q) | |
144 | ybit_sort = ptcompr_sort(Q) | |
145 | ||
146 | def out(f, s, Q): | |
147 | sufflen = randrange(16) | |
148 | suff = [randrange(256) for i in xrange(sufflen)] | |
149 | print ' %s\n %d\n %s\n %s\n %d;' % ( | |
150 | ec, f, os_to_hex(s + suff), ecpt(Q), sufflen) | |
151 | ||
152 | ## This is the algorithm from `gfreduce_quadsolve'. | |
153 | B = x^3 + E.a2()*x^2 + E.a4()*x + E.a6() | |
154 | A = E.a1()*x + E.a3() | |
155 | if A == 0: | |
156 | yy = sqrt(B) | |
157 | else: | |
158 | xx = B/A^2 | |
159 | while True: | |
160 | rho = k.random_element() | |
161 | if rho.trace() == 1: break | |
162 | z = sum(rho^(2^i)*sum(xx^(2^j) for j in xrange(i)) for i in xrange(m)) | |
163 | if z.polynomial()[0]: z -= 1 | |
164 | yy = A*z | |
165 | ||
166 | out(EC_XONLY, [EC_XONLY] + xs, E(x, yy)) | |
167 | out(EC_EXPLY, [EC_EXPLY] + xs + ys, Q) | |
168 | ||
169 | out(EC_LSB, [EC_CMPR | ybit_lsb] + xs, Q) | |
170 | out(EC_LSB | EC_EXPLY, [EC_EXPLY | EC_CMPR | ybit_lsb] + xs + ys, Q) | |
171 | out(EC_LSB | EC_EXPLY, | |
172 | [EC_EXPLY | EC_CMPR | (not ybit_lsb)] + xs + ys, None) | |
173 | ||
174 | out(EC_SORT, [EC_CMPR | EC_SORT | ybit_sort] + xs, Q) | |
175 | out(EC_SORT | EC_EXPLY, | |
176 | [EC_EXPLY | EC_CMPR | EC_SORT | ybit_sort] + xs + ys, Q) | |
177 | out(EC_SORT | EC_EXPLY, | |
178 | [EC_EXPLY | EC_CMPR | EC_SORT | (not ybit_sort)] + xs + ys, None) | |
179 | ||
180 | def test_os2ecp_2torsion(E, Q): | |
181 | ec = ecdescr(E) | |
182 | k = E.base_ring() | |
183 | m = k.degree() | |
184 | x, y = Q[0], Q[1] | |
185 | xs, ys = fe2osp(Q[0]), fe2osp(Q[1]) | |
186 | ||
187 | def out(f, s, Q): | |
188 | sufflen = randrange(16) | |
189 | suff = [randrange(256) for i in xrange(sufflen)] | |
190 | print ' %s\n %d\n %s\n %s\n %d;' % ( | |
191 | ec, f, os_to_hex(s + suff), ecpt(Q), sufflen) | |
192 | ||
193 | out(EC_LSB, [EC_CMPR] + xs, Q) | |
194 | out(EC_LSB, [EC_CMPR | EC_YBIT] + xs, None) | |
195 | out(EC_SORT, [EC_SORT | EC_CMPR] + xs, Q) | |
196 | out(EC_SORT, [EC_SORT | EC_CMPR | EC_YBIT] + xs, None) | |
197 | ||
198 | def test_os2ecp_notpoints(E): | |
199 | ec = ecdescr(E) | |
200 | k = E.base_ring() | |
201 | ||
202 | def out(f, s): | |
203 | print ' %s\n %d\n %s\n FAIL\n 0;' % (ec, f, os_to_hex(s)) | |
204 | ||
205 | while True: | |
206 | x = k.random_element() | |
207 | if not E.is_x_coord(x): break | |
208 | xs = fe2osp(x) | |
209 | out(EC_XONLY, [EC_XONLY] + xs) | |
210 | out(EC_LSB, [EC_CMPR] + xs) | |
211 | out(EC_SORT, [EC_CMPR | EC_SORT] + xs) | |
212 | ||
213 | print 'os2ecp {' | |
214 | for i in xrange(20): test_os2ecp(E_p, E_p.random_point()) | |
215 | test_os2ecp_2torsion(E_p, T_p) | |
216 | test_os2ecp_notpoints(E_p) | |
217 | for i in xrange(20): test_os2ecp(E_2, E_2.random_point()) | |
218 | test_os2ecp_2torsion(E_2, T_2) | |
219 | test_os2ecp_notpoints(E_2) | |
220 | print '}' |