| 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 '}' |