# Constants

p = 0xffffffff00000001000000000000000000000000ffffffffffffffffffffffff
a = 0xffffffff00000001000000000000000000000000fffffffffffffffffffffffc  # -3
b = 0x5ac635d8aa3a93e7b3ebbd55769886bc651d06b0cc53b0f63bce3c3e27d2604b

Px = 0x6b17d1f2e12c4247f8bce6e563a440f277037d812deb33a0f4a13945d898c296
Py = 0x4fe342e2fe1a7f9b8ee7eb4a7c0f9e162bce33576b315ececbb6406837bf51f5

Qx = 0xc97445f45cdef9f0d3e05e1e585fc297235b82b5be8ff3efca67c59852018192
Qy = 0xb28ef557ba31dfcbdd21ac46e2a91e3c304f44cb87058ada2cb815151e610046

n = 0xffffffff00000000ffffffffffffffffbce6faada7179e84f3b9cac2fc632551

# 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))

P is on curve:  True
Q is on curve:  True

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

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

# The point at infinity
INF = -1, -1

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

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

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())

be58a87ed729a0585d9af5e845e604c7ec2783f2b40b9dbba8cc36d9e3f0
230ce69cca1336e5d70dfca682665d0c2176d040ec693d8c6c936bf7a546
94ffaad3700ac9651c933280b9f46abf15f6b2402aec8b9c634a750add48
e206ecd890a0f00a067ad7626fd9b352561fc8a5550fccb42d80f9032fde
96c844ba91a7a8b67b5ef4b3c7920fcd62f5bd4ce0f850a2d1ff65070328
36d4f3f653b6515e3c58f51863189b92eca2bbd4321acb21ea5c448d784c
aff823a29659e6086e667dbcb045e9fdf499f64ad6f45bbb6ead7ecc20e9
c9ce423afbac6b7396e93f776a3675cd807b0b403d146966b6de4e157f09
10969f99625df3aa2720c57d74056be5abe148db737ef0b56d99de62b3c4
7b2ddc14165ebcc34183e5f183b5d26b97c35abe087ee2b720b3ee5da9e8
a694e58543452b0b24eb09dd0fc8f254908055774997f503d727c04951ad
4d6f14bcbd77162289be8903405f74dac4e0782d7924565f84b0ab615133
14ff391b7da36f7fb11957b9c89f626c6e3c15b1c256c67b7a26e34b3efd
2a21fdc8243b0f77d1bfc2750f321b535d1c9c8219b7ab05ed64eab15d69
896449de8a415895193684d94cde16b2c78e476afe95b6bd2c443a5396aa
c0dee1488e178e2126df8233be0a3d149b5ec44edbc6fb2601cc019057e4
08dc47177f6f9d5c45611d53e3415f4ffea3f5856c0a0753b17fcd363bec
225c3f62f15534aa0867d16c30a834b05a1bb6ca3355c8182a454ce154e6
73d12602be56fdea575215e2c50aace0cb6c803ee4ddd9f930ee4f80cff6
2838a063167235fe4762c91a810be06b0d8daddcf7653608eb924de636b0
e8ea631bb338a590a541a781fa734a3e3293b8ef6541680ff5a5584dc49b
8b7baefb785845f13c74275566eb919ba7f667a6dca96374195a5ef277bb
4dc9b8d717de71e8b4be821eb242dc065fe1ba6fc4385d92ce3bb392a3e1
7b4a93bcf41d93038bc3fd39c8485b75c4bf9e39ef4528b1b5c5c2235646