from random import SystemRandom # random number generator rand = SystemRandom () # strong random number generator def modadd (a, b, m): # return (a + b) mod m; a, b < m s = a + b if s > m: s = s - m return s def modsub (a, b, m): # return (a - b) mod m; a, b < m s = a - b if s < 0: s = s + m return s def modmul (a, b, m): # return (a * b) mod m; a, b < m # shift-sub (SSMM) u, v, s = a, b, 0 while v != 0: if v & 1 == 1: s = s + u if s > m: s = s - m v = v >> 1 u = u << 1 if u > m: u = u - m return s def modinv (b, a, m): # return (b * a^{-1}) mod m # proposed radix-8 modinv u, v = a, -m x, y = b, -0 while u != 1: if u & 1 == 1: tv, ty = v, y else: tv, ty = 0, 0 if v & 1 == 1: tu, tx = u, x else: tu, tx = 0, 0 tuv, txy = tu + tv, tx + ty # tuv is even if tuv & 6 == 0: # radix 8: uv = tuv // 8 if txy & 1 == 0: if txy & 2 == 0: if txy & 4 == 0: xy = txy // 8 else: xy = (txy + 4 * m) // 8 else: if txy & 4 == (m*2 & 4): xy = (txy + 6 * m) // 8 # -2m else: xy = (txy + 2 * m) // 8 else: if txy & 6 == m & 6: xy = (txy - m) // 8 else: if txy & 2 == m & 2: xy = (txy + 3 * m) // 8 else: if txy & 4 != m & 4: xy = (txy + m) // 8 else: xy = (txy + 5 * m) // 8 # -3m else: if tuv & 2 == 0: # radix 4: uv = tuv // 4 if txy & 1 == 0: if txy & 2 == 0: xy = txy // 4 else: xy = (txy + 2 * m) // 4 else: if txy & 3 == m & 3: xy = (txy - m) // 4 else: xy = (txy + m) // 4 else: # radix 2: uv = tuv // 2 if txy & 1 == 0: xy = txy // 2 else: if txy < 0: xy = (txy + m) // 2 else: xy = (txy - m) // 2 if uv < 0: v, y = uv, xy else: u, x = uv, xy if x < 0: x = x + m return x def point_addition (P, Q, m, a): # point addition R = P + Q x1, y1 = P x2, y2 = Q if x1 == -1 and y1 == -1: return Q # O + Q if x2 == -1 and y2 == -1: return P # P + O if x1 == x2: if modadd (y1, y2, m) == 0: return [-1, -1] # Point O else: return point_doubling (P, m, a) # 2P # s = ((y1 - y2) / (x1 - x2)) mod m s = modinv (modsub (y1, y2, m), modsub (x1, x2, m), m) # rx = (s * s - x1 - x2) mod m rx = modsub (modmul (s, s, m), modadd (x1, x2, m), m) # ry = (s * (x1 - rx) - y1) mod m ry = modsub (modmul (s, modsub (x1, rx, m), m), y1, m) return [int (rx), int (ry)] def point_doubling (P, m, a): # point doubling R = 2P x, y = P if y == 0: return [-1, -1] # Point O # s = ((3 * x * x + a) / (2 * y)) mod m s = modinv (modadd(a, modmul(modmul(x, x, m), 3, m), m), modadd(y, y, m), m) # rx = (s * s - 2 * x) mod m rx = modsub (modmul (s, s, m), modmul (x, 2, m), m) # ry = (s * (x - rx) - y) mod m ry = modsub (modmul (s, modsub (x, rx, m), m), y, m) return [int (rx), int (ry)] def scalar_point_multiplication (P, d, m, a): # scalar point multiplication if d == 0: return [-1, -1] # Point O k = d Q = [-1, -1] # Point O R = P while k != 0: if k & 1: Q = point_addition (Q, R, m, a) # Q = Q + R R = point_doubling (R, m, a) # R = 2R k >>= 1 return Q a = int (0x0000000000000000000000000000000000000000000000000000000000000000) b = int (0x0000000000000000000000000000000000000000000000000000000000000007) m = int (0xfffffffffffffffffffffffffffffffffffffffffffffffffffffffefffffc2f) x = int (0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798) y = int (0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8) P = [x, y] # Elliptic curve Diffie-Hellman (ECDH) key agreement: da = rand.getrandbits (256) % m # Alice's private key db = rand.getrandbits (256) % m # Bob's private key Qa = scalar_point_multiplication ( P, da, m, a) # Alice's public key Qb = scalar_point_multiplication ( P, db, m, a) # Bob's public key Qab = scalar_point_multiplication (Qb, da, m, a) # Alice calculates shared key Qba = scalar_point_multiplication (Qa, db, m, a) # Bob calculates shared key print ('da = 0x{:064x}'.format(da), end=' ') print ('db = 0x{:064x}'.format(db)) print ('Qax = 0x{:064x}'.format(Qa[0]), end=' ') print ('Qay = 0x{:064x}'.format(Qa[1])) print ('Qbx = 0x{:064x}'.format(Qb[0]), end=' ') print ('Qby = 0x{:064x}'.format(Qb[1])) print ('Qabx = 0x{:064x}'.format(Qab[0]), end=' ') print ('Qaby = 0x{:064x}'.format(Qab[1])) print ('Qbax = 0x{:064x}'.format(Qba[0]), end=' ') print ('Qbay = 0x{:064x}'.format(Qba[1])) assert (Qa [1] * Qa [1]) % m == (Qa [0] * Qa [0] * Qa [0] + a * Qa [0] + b) % m assert (Qb [1] * Qb [1]) % m == (Qb [0] * Qb [0] * Qb [0] + a * Qb [0] + b) % m assert (Qab[1] * Qab[1]) % m == (Qab[0] * Qab[0] * Qab[0] + a * Qab[0] + b) % m assert (Qba[1] * Qba[1]) % m == (Qba[0] * Qba[0] * Qba[0] + a * Qba[0] + b) % m assert Qab == Qba # verify correctness