### -*- mode: python; coding: utf-8 -*-
###--------------------------------------------------------------------------
-### Define the field.
+### Some general utilities.
+
+def ld(v):
+ return 0 + sum(ord(v[i]) << 8*i for i in xrange(len(v)))
+
+def st(x, n):
+ return ''.join(chr((x >> 8*i)&0xff) for i in xrange(n))
+
+def piece_widths_offsets(wd, n):
+ o = [ceil(wd*i/n) for i in xrange(n + 1)]
+ w = [o[i + 1] - o[i] for i in xrange(n)]
+ return w, o
+
+def pieces(x, wd, n, bias = 0):
+
+ ## Figure out widths and offsets.
+ w, o = piece_widths_offsets(wd, n)
+
+ ## First, normalize |n| < bias/2.
+ if bias and n >= bias/2: n -= bias
+
+ ## First, collect the bits.
+ nn = []
+ for i in xrange(n - 1):
+ m = (1 << w[i]) - 1
+ nn.append(x&m)
+ x >>= w[i]
+ nn.append(x)
+
+ ## Now normalize them to the appropriate interval.
+ c = 0
+ for i in xrange(n - 1):
+ b = 1 << (w[i] - 1)
+ if nn[i] >= b:
+ nn[i] -= 2*b
+ nn[i + 1] += 1
+
+ ## And we're done.
+ return nn
+
+def combine(v, wd, n):
+ w, o = piece_widths_offsets(wd, n)
+ return sum(v[i] << o[i] for i in xrange(n))
+
+###--------------------------------------------------------------------------
+### Define the curve.
p = 2^255 - 19; k = GF(p)
+A = k(486662); A0 = (A - 2)/4
+E = EllipticCurve(k, [0, A, 0, 1, 0]); P = E.lift_x(9)
+l = 2^252 + 27742317777372353535851937790883648493
+
+assert is_prime(l)
+assert (l*P).is_zero()
+assert (p + 1 - 8*l)^2 <= 4*p
+
+###--------------------------------------------------------------------------
+### Example points from `Cryptography in NaCl'.
+
+x = ld(map(chr, [0x70,0x07,0x6d,0x0a,0x73,0x18,0xa5,0x7d
+,0x3c,0x16,0xc1,0x72,0x51,0xb2,0x66,0x45
+,0xdf,0x4c,0x2f,0x87,0xeb,0xc0,0x99,0x2a
+,0xb1,0x77,0xfb,0xa5,0x1d,0xb9,0x2c,0x6a]))
+y = ld(map(chr, [0x58,0xab,0x08,0x7e,0x62,0x4a,0x8a,0x4b
+,0x79,0xe1,0x7f,0x8b,0x83,0x80,0x0e,0xe6
+,0x6f,0x3b,0xb1,0x29,0x26,0x18,0xb6,0xfd
+,0x1c,0x2f,0x8b,0x27,0xff,0x88,0xe0,0x6b]))
+X = x*P
+Y = y*P
+Z = x*Y
+assert Z == y*X
###--------------------------------------------------------------------------
### Arithmetic implementation.
for i in xrange(n): x = x*x
return x
+sqrtm1 = sqrt(k(-1))
+
def inv(x):
t2 = sqrn(x, 1) # 1 | 2
u = sqrn(t2, 2) # 3 | 8
t = u*t11 # 265 | 2^255 - 21
return t
+def quosqrt(x, y):
+
+ ## First, some preliminary values.
+ y2 = sqrn(y, 1) # 1 | 0, 2
+ y3 = y2*y # 2 | 0, 3
+ xy3 = x*y3 # 3 | 1, 3
+ y4 = sqrn(y2, 1) # 4 | 0, 4
+ w = xy3*y4 # 5 | 1, 7
+
+ ## Now calculate w^(p - 5)/8. Notice that (p - 5)/8 =
+ ## (2^255 - 24)/8 = 2^252 - 3.
+ u = sqrn(w, 1) # 6 | 2
+ t = u*w # 7 | 3
+ u = sqrn(t, 1) # 8 | 6
+ t = u*w # 9 | 7
+ u = sqrn(t, 3) # 12 | 56
+ t = u*t # 13 | 63 = 2^6 - 1
+ u = sqrn(t, 6) # 19 | 2^12 - 2^6
+ t = u*t # 20 | 2^12 - 1
+ u = sqrn(t, 12) # 32 | 2^24 - 2^12
+ t = u*t # 33 | 2^24 - 1
+ u = sqrn(t, 1) # 34 | 2^25 - 2
+ t = u*w # 35 | 2^25 - 1
+ u = sqrn(t, 25) # 60 | 2^50 - 2^25
+ t2p50m1 = u*t # 61 | 2^50 - 1
+ u = sqrn(t2p50m1, 50) # 111 | 2^100 - 2^50
+ t = u*t2p50m1 # 112 | 2^100 - 1
+ u = sqrn(t, 100) # 212 | 2^200 - 2^100
+ t = u*t # 213 | 2^200 - 1
+ u = sqrn(t, 50) # 263 | 2^250 - 2^50
+ t = u*t2p50m1 # 264 | 2^250 - 1
+ u = sqrn(t, 2) # 266 | 2^252 - 4
+ t = u*w # 267 | 2^252 - 3
+ beta = t*xy3 # 268 |
+
+ ## Now we have beta = (x y^3) (x y^7)^((p - 5)/8) =
+ ## x^((p + 3)/8) y^((7 p - 11)/8) = (x/y)^((p + 3)/8).
+ ## Suppose alpha^2 = x/y. Then beta^4 = (x/y)^((p + 3)/2) =
+ ## alpha^(p + 3) = alpha^4 = (x/y)^2, so beta^2 = ±x/y. If
+ ## y beta^2 = x then alpha = beta and we're done; if
+ ## y beta^2 = -x, then alpha = beta sqrt(-1); otherwise x/y
+ ## wasn't actually a square after all.
+ t = y*beta^2
+ if t == x: return beta
+ elif t == -x: return beta*sqrtm1
+ else: raise ValueError, 'not a square'
+
assert inv(k(9))*9 == 1
+assert 5*quosqrt(k(4), k(5))^2 == 4
+
+###--------------------------------------------------------------------------
+### The Montgomery ladder.
+
+A0 = (A - 2)/4
+
+def x25519(n, x1):
+
+ ## Let Q = (x_1 : y_1 : 1) be an input point. We calculate
+ ## n Q = (x_n : y_n : z_n), returning x_n/z_n (unless z_n = 0,
+ ## in which case we return zero).
+ ##
+ ## We're given that n = 2^254 + n'_254, where 0 <= n'_254 < 2^254.
+ bb = n.bits()
+ x, z = 1, 0
+ u, w = x1, 1
+
+ ## Initially, let i = 255.
+ for i in xrange(len(bb) - 1, -1, -1):
+
+ ## Split n = n_i 2^i + n'_i, where 0 <= n'_i < 2^i, so n_0 = n.
+ ## We have x, z = x_{n_{i+1}}, z_{n_{i+1}}, and
+ ## u, w = x_{n_{i+1}+1}, z_{n_{i+1}+1}.
+ ## Now either n_i = 2 n_{i+1} or n_i = 2 n_{i+1} + 1, depending
+ ## on bit i of n.
+
+ ## Swap (x : z) and (u : w) if bit i of n is set.
+ if bb[i]: x, z, u, w = u, w, x, z
+
+ ## Do the ladder step.
+ xmz, xpz = x - z, x + z
+ umw, upw = u - w, u + w
+ xmz2, xpz2 = xmz^2, xpz^2
+ xpz2mxmz2 = xpz2 - xmz2
+ xmzupw, xpzumw = xmz*upw, xpz*umw
+ x, z = xmz2*xpz2, xpz2mxmz2*(xpz2 + A0*xpz2mxmz2)
+ u, w = (xmzupw + xpzumw)^2, x1*(xmzupw - xpzumw)^2
+
+ ## Finally, unswap.
+ if bb[i]: x, z, u, w = u, w, x, z
+
+ ## Almost done.
+ return x*inv(z)
+
+assert x25519(y, k(9)) == Y[0]
+assert x25519(x, Y[0]) == x25519(y, X[0]) == Z[0]
###----- That's all, folks --------------------------------------------------
~mdw / catacomb / blobdiff