Commit | Line | Data |
---|---|---|
ee39a683 MW |
1 | #! /usr/local/bin/sage |
2 | ### -*- mode: python; coding: utf-8 -*- | |
3 | ||
4 | ###-------------------------------------------------------------------------- | |
fc2d44af MW |
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 | ###-------------------------------------------------------------------------- | |
14 | ### Define the curve. | |
ee39a683 MW |
15 | |
16 | p = 2^255 - 19; k = GF(p) | |
fc2d44af MW |
17 | A = k(486662); A0 = (A - 2)/4 |
18 | E = EllipticCurve(k, [0, A, 0, 1, 0]); P = E.lift_x(9) | |
19 | l = 2^252 + 27742317777372353535851937790883648493 | |
20 | ||
21 | assert is_prime(l) | |
22 | assert (l*P).is_zero() | |
23 | assert (p + 1 - 8*l)^2 <= 4*p | |
24 | ||
25 | ###-------------------------------------------------------------------------- | |
26 | ### Example points from `Cryptography in NaCl'. | |
27 | ||
28 | x = ld(map(chr, [0x70,0x07,0x6d,0x0a,0x73,0x18,0xa5,0x7d | |
29 | ,0x3c,0x16,0xc1,0x72,0x51,0xb2,0x66,0x45 | |
30 | ,0xdf,0x4c,0x2f,0x87,0xeb,0xc0,0x99,0x2a | |
31 | ,0xb1,0x77,0xfb,0xa5,0x1d,0xb9,0x2c,0x6a])) | |
32 | y = ld(map(chr, [0x58,0xab,0x08,0x7e,0x62,0x4a,0x8a,0x4b | |
33 | ,0x79,0xe1,0x7f,0x8b,0x83,0x80,0x0e,0xe6 | |
34 | ,0x6f,0x3b,0xb1,0x29,0x26,0x18,0xb6,0xfd | |
35 | ,0x1c,0x2f,0x8b,0x27,0xff,0x88,0xe0,0x6b])) | |
36 | X = x*P | |
37 | Y = y*P | |
38 | Z = x*Y | |
39 | assert Z == y*X | |
ee39a683 MW |
40 | |
41 | ###-------------------------------------------------------------------------- | |
42 | ### Arithmetic implementation. | |
43 | ||
44 | def sqrn(x, n): | |
45 | for i in xrange(n): x = x*x | |
46 | return x | |
47 | ||
48 | def inv(x): | |
49 | t2 = sqrn(x, 1) # 1 | 2 | |
50 | u = sqrn(t2, 2) # 3 | 8 | |
51 | t = u*x # 4 | 9 | |
52 | t11 = t*t2 # 5 | 11 | |
53 | u = sqrn(t11, 1) # 6 | 22 | |
54 | t = u*t # 7 | 2^5 - 1 = 31 | |
55 | u = sqrn(t, 5) # 12 | 2^10 - 2^5 | |
56 | t2p10m1 = u*t # 13 | 2^10 - 1 | |
57 | u = sqrn(t2p10m1, 10) # 23 | 2^20 - 2^10 | |
58 | t = u*t2p10m1 # 24 | 2^20 - 1 | |
59 | u = sqrn(t, 20) # 44 | 2^40 - 2^20 | |
60 | t = u*t # 45 | 2^40 - 1 | |
61 | u = sqrn(t, 10) # 55 | 2^50 - 2^10 | |
62 | t2p50m1 = u*t2p10m1 # 56 | 2^50 - 1 | |
63 | u = sqrn(t2p50m1, 50) # 106 | 2^100 - 2^50 | |
64 | t = u*t2p50m1 # 107 | 2^100 - 1 | |
65 | u = sqrn(t, 100) # 207 | 2^200 - 2^100 | |
66 | t = u*t # 208 | 2^200 - 1 | |
67 | u = sqrn(t, 50) # 258 | 2^250 - 2^50 | |
68 | t = u*t2p50m1 # 259 | 2^250 - 1 | |
69 | u = sqrn(t, 5) # 264 | 2^255 - 2^5 | |
70 | t = u*t11 # 265 | 2^255 - 21 | |
71 | return t | |
72 | ||
73 | assert inv(k(9))*9 == 1 | |
74 | ||
fc2d44af MW |
75 | ###-------------------------------------------------------------------------- |
76 | ### The Montgomery ladder. | |
77 | ||
78 | A0 = (A - 2)/4 | |
79 | ||
80 | def x25519(n, x1): | |
81 | ||
82 | ## Let Q = (x_1 : y_1 : 1) be an input point. We calculate | |
83 | ## n Q = (x_n : y_n : z_n), returning x_n/z_n (unless z_n = 0, | |
84 | ## in which case we return zero). | |
85 | ## | |
86 | ## We're given that n = 2^254 + n'_254, where 0 <= n'_254 < 2^254. | |
87 | bb = n.bits() | |
88 | x, z = 1, 0 | |
89 | u, w = x1, 1 | |
90 | ||
91 | ## Initially, let i = 255. | |
92 | for i in xrange(len(bb) - 1, -1, -1): | |
93 | ||
94 | ## Split n = n_i 2^i + n'_i, where 0 <= n'_i < 2^i, so n_0 = n. | |
95 | ## We have x, z = x_{n_{i+1}}, z_{n_{i+1}}, and | |
96 | ## u, w = x_{n_{i+1}+1}, z_{n_{i+1}+1}. | |
97 | ## Now either n_i = 2 n_{i+1} or n_i = 2 n_{i+1} + 1, depending | |
98 | ## on bit i of n. | |
99 | ||
100 | ## Swap (x : z) and (u : w) if bit i of n is set. | |
101 | if bb[i]: x, z, u, w = u, w, x, z | |
102 | ||
103 | ## Do the ladder step. | |
104 | xmz, xpz = x - z, x + z | |
105 | umw, upw = u - w, u + w | |
106 | xmz2, xpz2 = xmz^2, xpz^2 | |
107 | xpz2mxmz2 = xpz2 - xmz2 | |
108 | xmzupw, xpzumw = xmz*upw, xpz*umw | |
109 | x, z = xmz2*xpz2, xpz2mxmz2*(xpz2 + A0*xpz2mxmz2) | |
110 | u, w = (xmzupw + xpzumw)^2, x1*(xmzupw - xpzumw)^2 | |
111 | ||
112 | ## Finally, unswap. | |
113 | if bb[i]: x, z, u, w = u, w, x, z | |
114 | ||
115 | ## Almost done. | |
116 | return x*inv(z) | |
117 | ||
118 | assert x25519(y, k(9)) == Y[0] | |
119 | assert x25519(x, Y[0]) == x25519(y, X[0]) == Z[0] | |
120 | ||
ee39a683 | 121 | ###----- That's all, folks -------------------------------------------------- |