mdw@git.distorted.org.uk Git - catacomb/blob - utils/curve25519.sage

   1 #! /usr/local/bin/sage
   2 ### -*- mode: python; coding: utf-8 -*-
   3
   4 ###--------------------------------------------------------------------------
   5 ### Some general utilities.
   6
   7 def ld(v):
   8   return 0 + sum(ord(v[i]) << 8*i for i in xrange(len(v)))
   9
  10 def st(x, n):
  11   return ''.join(chr((x >> 8*i)&0xff) for i in xrange(n))
  12
  13 def piece_widths_offsets(wd, n):
  14   o = [ceil(wd*i/n) for i in xrange(n + 1)]
  15   w = [o[i + 1] - o[i] for i in xrange(n)]
  16   return w, o
  17
  18 def pieces(x, wd, n, bias = 0):
  19
  20   ## Figure out widths and offsets.
  21   w, o = piece_widths_offsets(wd, n)
  22
  23   ## First, normalize |n| < bias/2.
  24   if bias and n >= bias/2: n -= bias
  25
  26   ## First, collect the bits.
  27   nn = []
  28   for i in xrange(n - 1):
  29     m = (1 << w[i]) - 1
  30     nn.append(x&m)
  31     x >>= w[i]
  32   nn.append(x)
  33
  34   ## Now normalize them to the appropriate interval.
  35   c = 0
  36   for i in xrange(n - 1):
  37     b = 1 << (w[i] - 1)
  38     if nn[i] >= b:
  39       nn[i] -= 2*b
  40       nn[i + 1] += 1
  41
  42   ## And we're done.
  43   return nn
  44
  45 def combine(v, wd, n):
  46   w, o = piece_widths_offsets(wd, n)
  47   return sum(v[i] << o[i] for i in xrange(n))
  48
  49 ###--------------------------------------------------------------------------
  50 ### Define the curve.
  51
  52 p = 2^255 - 19; k = GF(p)
  53 A = k(486662); A0 = (A - 2)/4
  54 E = EllipticCurve(k, [0, A, 0, 1, 0]); P = E.lift_x(9)
  55 l = 2^252 + 27742317777372353535851937790883648493
  56
  57 assert is_prime(l)
  58 assert (l*P).is_zero()
  59 assert (p + 1 - 8*l)^2 <= 4*p
  60
  61 ###--------------------------------------------------------------------------
  62 ### Example points from `Cryptography in NaCl'.
  63
  64 x = ld(map(chr, [0x70,0x07,0x6d,0x0a,0x73,0x18,0xa5,0x7d
  65 ,0x3c,0x16,0xc1,0x72,0x51,0xb2,0x66,0x45
  66 ,0xdf,0x4c,0x2f,0x87,0xeb,0xc0,0x99,0x2a
  67 ,0xb1,0x77,0xfb,0xa5,0x1d,0xb9,0x2c,0x6a]))
  68 y = ld(map(chr, [0x58,0xab,0x08,0x7e,0x62,0x4a,0x8a,0x4b
  69 ,0x79,0xe1,0x7f,0x8b,0x83,0x80,0x0e,0xe6
  70 ,0x6f,0x3b,0xb1,0x29,0x26,0x18,0xb6,0xfd
  71 ,0x1c,0x2f,0x8b,0x27,0xff,0x88,0xe0,0x6b]))
  72 X = x*P
  73 Y = y*P
  74 Z = x*Y
  75 assert Z == y*X
  76
  77 ###--------------------------------------------------------------------------
  78 ### Arithmetic implementation.
  79
  80 def sqrn(x, n):
  81   for i in xrange(n): x = x*x
  82   return x
  83
  84 sqrtm1 = sqrt(k(-1))
  85
  86 def inv(x):
  87   t2 = sqrn(x, 1)        #   1 | 2
  88   u = sqrn(t2, 2)        #   3 | 8
  89   t = u*x                #   4 | 9
  90   t11 = t*t2             #   5 | 11
  91   u = sqrn(t11, 1)       #   6 | 22
  92   t = u*t                #   7 | 2^5 - 1 = 31
  93   u = sqrn(t, 5)         #  12 | 2^10 - 2^5
  94   t2p10m1 = u*t          #  13 | 2^10 - 1
  95   u = sqrn(t2p10m1, 10)  #  23 | 2^20 - 2^10
  96   t = u*t2p10m1          #  24 | 2^20 - 1
  97   u = sqrn(t, 20)        #  44 | 2^40 - 2^20
  98   t = u*t                #  45 | 2^40 - 1
  99   u = sqrn(t, 10)        #  55 | 2^50 - 2^10
 100   t2p50m1 = u*t2p10m1    #  56 | 2^50 - 1
 101   u = sqrn(t2p50m1, 50)  # 106 | 2^100 - 2^50
 102   t = u*t2p50m1          # 107 | 2^100 - 1
 103   u = sqrn(t, 100)       # 207 | 2^200 - 2^100
 104   t = u*t                # 208 | 2^200 - 1
 105   u = sqrn(t, 50)        # 258 | 2^250 - 2^50
 106   t = u*t2p50m1          # 259 | 2^250 - 1
 107   u = sqrn(t, 5)         # 264 | 2^255 - 2^5
 108   t = u*t11              # 265 | 2^255 - 21
 109   return t
 110
 111 def quosqrt(x, y):
 112
 113   ## First, some preliminary values.
 114   y2 = sqrn(y, 1)        #   1 | 0, 2
 115   y3 = y2*y              #   2 | 0, 3
 116   xy3 = x*y3             #   3 | 1, 3
 117   y4 = sqrn(y2, 1)       #   4 | 0, 4
 118   w = xy3*y4             #   5 | 1, 7
 119
 120   ## Now calculate w^(p - 5)/8.  Notice that (p - 5)/8 =
 121   ## (2^255 - 24)/8 = 2^252 - 3.
 122   u = sqrn(w, 1)         #   6 | 2
 123   t = u*w                #   7 | 3
 124   u = sqrn(t, 1)         #   8 | 6
 125   t = u*w                #   9 | 7
 126   u = sqrn(t, 3)         #  12 | 56
 127   t = u*t                #  13 | 63 = 2^6 - 1
 128   u = sqrn(t, 6)         #  19 | 2^12 - 2^6
 129   t = u*t                #  20 | 2^12 - 1
 130   u = sqrn(t, 12)        #  32 | 2^24 - 2^12
 131   t = u*t                #  33 | 2^24 - 1
 132   u = sqrn(t, 1)         #  34 | 2^25 - 2
 133   t = u*w                #  35 | 2^25 - 1
 134   u = sqrn(t, 25)        #  60 | 2^50 - 2^25
 135   t2p50m1 = u*t          #  61 | 2^50 - 1
 136   u = sqrn(t2p50m1, 50)  # 111 | 2^100 - 2^50
 137   t = u*t2p50m1          # 112 | 2^100 - 1
 138   u = sqrn(t, 100)       # 212 | 2^200 - 2^100
 139   t = u*t                # 213 | 2^200 - 1
 140   u = sqrn(t, 50)        # 263 | 2^250 - 2^50
 141   t = u*t2p50m1          # 264 | 2^250 - 1
 142   u = sqrn(t, 2)         # 266 | 2^252 - 4
 143   t = u*w                # 267 | 2^252 - 3
 144   beta = t*xy3           # 268 |
 145
 146   ## Now we have beta = (x y^3) (x y^7)^((p - 5)/8) =
 147   ## x^((p + 3)/8) y^((7 p - 11)/8) = (x/y)^((p + 3)/8).
 148   ## Suppose alpha^2 = x/y.  Then beta^4 = (x/y)^((p + 3)/2) =
 149   ## alpha^(p + 3) = alpha^4 = (x/y)^2, so beta^2 = ±x/y.  If
 150   ## y beta^2 = x then alpha = beta and we're done; if
 151   ## y beta^2 = -x, then alpha = beta sqrt(-1); otherwise x/y
 152   ## wasn't actually a square after all.
 153   t = y*beta^2
 154   if t == x: return beta
 155   elif t == -x: return beta*sqrtm1
 156   else: raise ValueError, 'not a square'
 157
 158 assert inv(k(9))*9 == 1
 159 assert 5*quosqrt(k(4), k(5))^2 == 4
 160
 161 ###--------------------------------------------------------------------------
 162 ### The Montgomery ladder.
 163
 164 A0 = (A - 2)/4
 165
 166 def x25519(n, x1):
 167
 168   ## Let Q = (x_1 : y_1 : 1) be an input point.  We calculate
 169   ## n Q = (x_n : y_n : z_n), returning x_n/z_n (unless z_n = 0,
 170   ## in which case we return zero).
 171   ##
 172   ## We're given that n = 2^254 + n'_254, where 0 <= n'_254 < 2^254.
 173   bb = n.bits()
 174   x, z = 1, 0
 175   u, w = x1, 1
 176
 177   ## Initially, let i = 255.
 178   for i in xrange(len(bb) - 1, -1, -1):
 179
 180     ## Split n = n_i 2^i + n'_i, where 0 <= n'_i < 2^i, so n_0 = n.
 181     ## We have x, z = x_{n_{i+1}}, z_{n_{i+1}}, and
 182     ## u, w = x_{n_{i+1}+1}, z_{n_{i+1}+1}.
 183     ## Now either n_i = 2 n_{i+1} or n_i = 2 n_{i+1} + 1, depending
 184     ## on bit i of n.
 185
 186     ## Swap (x : z) and (u : w) if bit i of n is set.
 187     if bb[i]: x, z, u, w = u, w, x, z
 188
 189     ## Do the ladder step.
 190     xmz, xpz = x - z, x + z
 191     umw, upw = u - w, u + w
 192     xmz2, xpz2 = xmz^2, xpz^2
 193     xpz2mxmz2 = xpz2 - xmz2
 194     xmzupw, xpzumw = xmz*upw, xpz*umw
 195     x, z = xmz2*xpz2, xpz2mxmz2*(xpz2 + A0*xpz2mxmz2)
 196     u, w = (xmzupw + xpzumw)^2, x1*(xmzupw - xpzumw)^2
 197
 198     ## Finally, unswap.
 199     if bb[i]: x, z, u, w = u, w, x, z
 200
 201   ## Almost done.
 202   return x*inv(z)
 203
 204 assert x25519(y, k(9)) == Y[0]
 205 assert x25519(x, Y[0]) == x25519(y, X[0]) == Z[0]
 206
 207 ###----- That's all, folks --------------------------------------------------