Let’s see a quick implementation of the DUAL_EC_DRBG algorithm.
In [1]:
# Constants p = 0xffffffff00000001000000000000000000000000ffffffffffffffffffffffff a = 0xffffffff00000001000000000000000000000000fffffffffffffffffffffffc # -3 b = 0x5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b Px = 0x6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296 Py = 0x4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5 Qx = 0xc97445f45cdef9f0d3e05e1e585fc297235b82b5be8ff3efca67c59852018192 Qy = 0xb28ef557ba31dfcbdd21ac46e2a91e3c304f44cb87058ada2cb815151e610046 n = 0xffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551
In [2]:
# Curve equation: y^2 = x^3 + ax + b is_on_curve = lambda x, y: (y**2) % p == (x**3 + a*x + b) % p print("P is on curve: ", is_on_curve(Px, Py)) print("Q is on curve: ", is_on_curve(Qx, Qy))
Out:
P is on curve: True Q is on curve: True
Elliptic curve arithmetic
In [3]:
def extended_euclidian_algorithm(a: int, b: int) -> tuple[int, int, int]: old_r, r = a, b old_s, s = 1, 0 old_t, t = 0, 1 while r != 0: quotient = old_r // r old_r, r = r, old_r - quotient * r old_s, s = s, old_s - quotient * s old_t, t = t, old_t - quotient * t return old_r, old_s, old_t
In [4]:
def modular_multiplicative_inverse(a: int, n: int) -> int: g, x, _ = extended_euclidian_algorithm(a, n) assert g == 1, f"{a} has no multiplicative inverse modulo {n}" return x % n
In [5]:
# The point at infinity INF = -1, -1
In [6]:
def ec_add(x1: int, y1: int, x2: int, y2: int) -> tuple[int, int]: if (x1, y1) == INF: return x2, y2 if (x2, y2) == INF: return x1, y1 if x1 == x2 and y1 != y2: return INF if x1 == x2: m = (3 * x1**2 + a) * modular_multiplicative_inverse(2 * y1, p) else: m = (y1 - y2) * modular_multiplicative_inverse(x1 - x2, p) x3 = (m**2 - x1 - x2) % p y3 = (m * (x1 - x3) - y1) % p return x3, y3
In [7]:
def ec_point_int_mul(x: int, y: int, n: int) -> tuple[int, int]: res = INF addend = x, y while n: if n & 1: res = ec_add(*res, *addend) addend = ec_add(*addend, *addend) n >>= 1 return res
In [8]:
seed = 0xd530b913e6f2ef88b21616fd34a603f203d0578c outsize = p.bit_length() - 16 # Chop off first 16 bits outmask = (1 << outsize) - 1 output = "" for it in range(25): seed, _ = ec_point_int_mul(Px, Py, seed) if it == 0: continue # Do not output first iteration r, _ = ec_point_int_mul(Qx, Qy, seed) r &= outmask hex_out = hex(r)[2:].zfill(outsize // 4) output += hex_out + "\n" print(output.strip())
Out:
be58a87ed729a0585d9af5e845e604c7ec2783f2b40b9dbba8cc36d9e3f0 230ce69cca1336e5d70dfca682665d0c2176d040ec693d8c6c936bf7a546 94ffaad3700ac9651c933280b9f46abf15f6b2402aec8b9c634a750add48 e206ecd890a0f00a067ad7626fd9b352561fc8a5550fccb42d80f9032fde 96c844ba91a7a8b67b5ef4b3c7920fcd62f5bd4ce0f850a2d1ff65070328 36d4f3f653b6515e3c58f51863189b92eca2bbd4321acb21ea5c448d784c aff823a29659e6086e667dbcb045e9fdf499f64ad6f45bbb6ead7ecc20e9 c9ce423afbac6b7396e93f776a3675cd807b0b403d146966b6de4e157f09 10969f99625df3aa2720c57d74056be5abe148db737ef0b56d99de62b3c4 7b2ddc14165ebcc34183e5f183b5d26b97c35abe087ee2b720b3ee5da9e8 a694e58543452b0b24eb09dd0fc8f254908055774997f503d727c04951ad 4d6f14bcbd77162289be8903405f74dac4e0782d7924565f84b0ab615133 14ff391b7da36f7fb11957b9c89f626c6e3c15b1c256c67b7a26e34b3efd 2a21fdc8243b0f77d1bfc2750f321b535d1c9c8219b7ab05ed64eab15d69 896449de8a415895193684d94cde16b2c78e476afe95b6bd2c443a5396aa c0dee1488e178e2126df8233be0a3d149b5ec44edbc6fb2601cc019057e4 08dc47177f6f9d5c45611d53e3415f4ffea3f5856c0a0753b17fcd363bec 225c3f62f15534aa0867d16c30a834b05a1bb6ca3355c8182a454ce154e6 73d12602be56fdea575215e2c50aace0cb6c803ee4ddd9f930ee4f80cff6 2838a063167235fe4762c91a810be06b0d8daddcf7653608eb924de636b0 e8ea631bb338a590a541a781fa734a3e3293b8ef6541680ff5a5584dc49b 8b7baefb785845f13c74275566eb919ba7f667a6dca96374195a5ef277bb 4dc9b8d717de71e8b4be821eb242dc065fe1ba6fc4385d92ce3bb392a3e1 7b4a93bcf41d93038bc3fd39c8485b75c4bf9e39ef4528b1b5c5c2235646