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