chiark - git - mdw - catacomb/blob - utils/ec-compr-test.sage

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