Demonstrates why you can get a malleable ECDSA (secp256k1) signature in Python

secp256k1_malleable.py

"""
This code is adapted from one of Vitalik's earliest implementations of Ethereum,
written in python surprisingly. It helps demonstrate the process behind ECDSA
"""
import numpy as np
import hashlib, hmac
import sys
if sys.version[0] == '2':
    safe_ord = ord
else:
    safe_ord = lambda x: x
np.set_printoptions(formatter={'object':hex})

# Elliptic curve parameters (secp256k1)
# from https://www.secg.org/sec2-v2.pdf
P = 2**256 - 2**32 - 977
# order of finite cyclic field
N = 115792089237316195423570985008687907852837564279074904382605163141518161494337
A = 0
B = 7
Gx = 55066263022277343669578718895168534326250603453777594175500187360389116729240
Gy = 32670510020758816978083085130507043184471273380659243275938904335757337482424
G = np.array([Gx, Gy])

def bytes_to_int(x):
    o = 0
    for b in x:
        o = (o << 8) + safe_ord(b)
    return o

# Extended Euclidean Algorithm
def inv(a, n):
    if a == 0:
        return 0
    lm, hm = 1, 0
    low, high = a % n, n
    while low > 1:
        r = high//low
        nm, new = hm-lm*r, high-low*r
        lm, low, hm, high = nm, new, lm, low
    return lm % n

def to_jacobian(p):
    o = (p[0], p[1], 1)
    return o

def jacobian_double(p):
    if not p[1]:
        return (0, 0, 0)
    ysq = (p[1] ** 2) % P
    S = (4 * p[0] * ysq) % P
    M = (3 * p[0] ** 2 + A * p[2] ** 4) % P
    nx = (M**2 - 2 * S) % P
    ny = (M * (S - nx) - 8 * ysq ** 2) % P
    nz = (2 * p[1] * p[2]) % P
    return (nx, ny, nz)

def jacobian_add(p, q):
    if not p[1]:
        return q
    if not q[1]:
        return p
    U1 = (p[0] * q[2] ** 2) % P
    U2 = (q[0] * p[2] ** 2) % P
    S1 = (p[1] * q[2] ** 3) % P
    S2 = (q[1] * p[2] ** 3) % P
    if U1 == U2:
        if S1 != S2:
            return (0, 0, 1)
        return jacobian_double(p)
    H = U2 - U1
    R = S2 - S1
    H2 = (H * H) % P
    H3 = (H * H2) % P
    U1H2 = (U1 * H2) % P
    nx = (R ** 2 - H3 - 2 * U1H2) % P
    ny = (R * (U1H2 - nx) - S1 * H3) % P
    nz = (H * p[2] * q[2]) % P
    return (nx, ny, nz)

def from_jacobian(p):
    z = inv(p[2], P)
    return ((p[0] * z**2) % P, (p[1] * z**3) % P)

def jacobian_multiply(a, n):
    if a[1] == 0 or n == 0:
        return (0, 0, 1)
    if n == 1:
        return a
    if n < 0 or n >= N:
        return jacobian_multiply(a, n % N)
    if (n % 2) == 0:
        return jacobian_double(jacobian_multiply(a, n//2))
    if (n % 2) == 1:
        return jacobian_add(jacobian_double(jacobian_multiply(a, n//2)), a)

def multiply(a, n):
    return from_jacobian(jacobian_multiply(to_jacobian(a), n))

def add(a, b):
    return from_jacobian(jacobian_add(to_jacobian(a), to_jacobian(b)))

def privtopub(privkey):
    return multiply(G, privkey)

pubkey = "0xf39Fd6e51aad88F6F4ce6aB8827279cffFb92266"
prvkey  = "0xac0974bec39a17e36ba4a6b4d238ff944bacb478cbed5efcae784d7bf4f2ff80"

prvkey2 = hex(N - int(prvkey, 16))

pair1 = np.array(privtopub(int(prvkey,16)))
pair2 = np.array(privtopub(int(prvkey2,16)))

assert pair1[0] == pair2[0]