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

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

while n:
if n & 1:
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)
hex_out = hex(r)[2:].zfill(outsize // 4)
output += hex_out + "\n"

print(output.strip())```
Out:
```be58a87ed729a0585d9af5e845e604c7ec2783f2b40b9dbba8cc36d9e3f0
230ce69cca1336e5d70dfca682665d0c2176d040ec693d8c6c936bf7a546
96c844ba91a7a8b67b5ef4b3c7920fcd62f5bd4ce0f850a2d1ff65070328
36d4f3f653b6515e3c58f51863189b92eca2bbd4321acb21ea5c448d784c
c9ce423afbac6b7396e93f776a3675cd807b0b403d146966b6de4e157f09
10969f99625df3aa2720c57d74056be5abe148db737ef0b56d99de62b3c4
7b2ddc14165ebcc34183e5f183b5d26b97c35abe087ee2b720b3ee5da9e8
4d6f14bcbd77162289be8903405f74dac4e0782d7924565f84b0ab615133
14ff391b7da36f7fb11957b9c89f626c6e3c15b1c256c67b7a26e34b3efd
2a21fdc8243b0f77d1bfc2750f321b535d1c9c8219b7ab05ed64eab15d69
896449de8a415895193684d94cde16b2c78e476afe95b6bd2c443a5396aa
c0dee1488e178e2126df8233be0a3d149b5ec44edbc6fb2601cc019057e4
08dc47177f6f9d5c45611d53e3415f4ffea3f5856c0a0753b17fcd363bec
225c3f62f15534aa0867d16c30a834b05a1bb6ca3355c8182a454ce154e6
73d12602be56fdea575215e2c50aace0cb6c803ee4ddd9f930ee4f80cff6