lab/ds/25-1/1e/secp256k1.py

# secp256k1_pure_python.py
# Pure Python scalar->public (affine x,y) for secp256k1
# No external libs required.

# Curve parameters for secp256k1
P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F
A = 0  # a = 0 for secp256k1
B = 7
# Base point G
Gx = 0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798
Gy = 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8
# Order of base point
N = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141

# ---------- Field helpers ----------
def inv_mod(x, p=P):
    """Multiplicative inverse modulo p (Python built-in pow with -1 not supported)"""
    return pow(x, p - 2, p)

# ---------- Jacobian coordinates helpers ----------
# Represent point as tuple (X, Y, Z) meaning affine (X/Z^2, Y/Z^3).
INFINITY = (0, 1, 0)  # Z == 0 means point at infinity in these helpers

def is_infinity(Pj):
    return Pj[2] == 0

def to_jacobian(x, y):
    return (x % P, y % P, 1)

def from_jacobian(Pj):
    X, Y, Z = Pj
    if Z == 0:
        return None  # infinity
    Z_inv = inv_mod(Z)
    Z_inv2 = (Z_inv * Z_inv) % P
    Z_inv3 = (Z_inv2 * Z_inv) % P
    x = (X * Z_inv2) % P
    y = (Y * Z_inv3) % P
    return (x, y)

# Point doubling in Jacobian coords
def jacobian_double(Pj):
    X1, Y1, Z1 = Pj
    if Z1 == 0 or Y1 == 0:
        return INFINITY
    # Formula for a = 0
    S = (4 * X1 * pow(Y1, 2, P)) % P
    M = (3 * pow(X1, 2, P)) % P  # since a=0
    X3 = (pow(M, 2, P) - 2 * S) % P
    Y3 = (M * (S - X3) - 8 * pow(Y1, 4, P)) % P
    Z3 = (2 * Y1 * Z1) % P
    return (X3, Y3, Z3)

# Mixed addition: add Jacobian Pj and affine Q=(x2,y2)
def jacobian_add_affine(Pj, Qx, Qy):
    X1, Y1, Z1 = Pj
    if Z1 == 0:
        return to_jacobian(Qx, Qy)
    if Y1 == 0 and X1 == 0 and Z1 == 0:
        return to_jacobian(Qx, Qy)

    Z1z = (Z1 * Z1) % P
    U2 = (Qx * Z1z) % P
    S2 = (Qy * Z1z * Z1) % P  # Qy * Z1^3

    H = (U2 - X1) % P
    r = (S2 - Y1) % P

    if H == 0:
        if r == 0:
            return jacobian_double(Pj)
        else:
            return INFINITY

    HH = (H * H) % P
    HHH = (H * HH) % P
    V = (X1 * HH) % P

    X3 = (r * r - HHH - 2 * V) % P
    Y3 = (r * (V - X3) - Y1 * HHH) % P
    Z3 = (Z1 * HH) % P
    return (X3, Y3, Z3)

# ---------- Scalar multiplication ----------
def scalar_mult(k, Px=Gx, Py=Gy):
    """
    Multiply base point (Px,Py) by scalar k.
    Returns affine (x,y) or None for infinity.
    k should be integer >=0.
    """
    if k % N == 0:
        return None
    if k < 0:
        # use -P
        return scalar_mult(-k, Px, (-Py) % P)

    # Ensure scalar reduced mod N and >0
    k = k % N
    if k == 0:
        return None

    # Initialize result as infinity (Jacobian)
    R = INFINITY
    # Take base point as affine (we'll use mixed add)
    Px = Px % P
    Py = Py % P

    # Left-to-right binary method
    bits = k.bit_length()
    for i in range(bits - 1, -1, -1):
        # R = 2*R
        R = jacobian_double(R)
        if (k >> i) & 1:
            # R = R + P
            R = jacobian_add_affine(R, Px, Py)
    # Convert R to affine
    return from_jacobian(R)

# ---------- Utility outputs ----------
def pubkey_bytes_compressed(x, y):
    prefix = b'\x02' if (y % 2 == 0) else b'\x03'
    xb = x.to_bytes(32, 'big')
    return prefix + xb

def pubkey_bytes_uncompressed(x, y):
    return b'\x04' + x.to_bytes(32, 'big') + y.to_bytes(32, 'big')

# ---------- Public function as user asked ----------
def priv_to_pub(priv_int):
    """
    Input: priv_int (Python int) — private key (assumed 0 < priv < N but function reduces)
    Output: (x, y) tuple of Python ints (affine coordinates). Returns None if point at infinity (shouldn't for valid keys).
    """
    if not isinstance(priv_int, int):
        raise TypeError("priv must be int")
    k = priv_int % N
    if k == 0:
        raise ValueError("private key is zero modulo curve order")
    pt = scalar_mult(k, Gx, Gy)
    if pt is None:
        raise RuntimeError("Result is point at infinity (shouldn't happen for valid priv)")
    return pt

# ---------- Example / quick test ----------
if __name__ == "__main__":
    # Test vector: priv = 1 -> pub = G
    priv = 1
    x, y = priv_to_pub(priv)
    assert x == Gx and y == Gy
    print("priv=1 -> pub == G : OK")

    # Test with a random private (example)
    import secrets
    priv = secrets.randbelow(N - 1) + 1
    x, y = priv_to_pub(priv)
    print("priv (hex):", hex(priv))
    print("pub x:", hex(x))
    print("pub y:", hex(y))
    print("compressed pubkey:", pubkey_bytes_compressed(x, y).hex())