Tutorial covering the concept of Elliptic Curve Cryptography and Diffie-Hellman Key Exchange, tutorial available under https://mareknarozniak.com/2020/11/30/ecdh/

ecc.py

import numpy as np


# ECC addition in R
def eccContAddition(Px, Py, Qx, Qy, a, b):
    if np.isclose([Px], [Qx])[0]:
        m = (3.*Px**2.+a)/(2.*Py)
    else:
        m = (Py-Qy)/(Px-Qx)
    Rx = m**2 - Px - Qx
    Ry = Py + m*(Rx - Px)
    return Rx, Ry, m


# EC point addition
def onCurve(X, Y, a, b, p):
    # None is the infinity point
    if X is None and Y is None:
        return True
    return np.mod(Y**2 - X**3 - a*X - b, p) == 0


def extendedEuclideanAlgorithm(x1, x2):
    # base case
    if x1 == 0:
        return x2, 0, 1
    # recursive call
    gcd, x_, y_ = extendedEuclideanAlgorithm(x2 % x1, x1)
    x = y_ - (x2 // x1) * x_
    y = x_
    return gcd, x, y


# modular multiplicative inverse
def inverse(n, p):
    if n == 0:
        raise ZeroDivisionError('Stahp dividing by zero you!!!')
    if n < 0:
        # k ** -1 = p - (-k) ** -1  (mod p)
        return p - inverse(-n, p)
    gcd, x, _ = extendedEuclideanAlgorithm(n, p)
    assert gcd == 1
    assert np.mod(n*x, p) == 1
    return np.mod(x, p)


def eccFiniteAddition(Px, Py, Qx, Qy, a, b, p):
    if None in [Px, Py]:
        return Qx, Qy
    if None in [Qx, Qy]:
        return Px, Py
    assert onCurve(Px, Py, a, b, p)
    assert onCurve(Qx, Qy, a, b, p)
    if Px == Qx and Py != Qy:
        # point1 + (-point1) = 0
        return None, None
    if Px == Qx:
        m = (3*Px**2 + a)*inverse(2*Py, p)
    else:
        m = (Py - Qy)*inverse(Px - Qx, p)
    Rx = np.mod(m**2 - Px - Qx, p)
    Ry = np.mod(-(Py + m*(Rx - Px)), p)
    assert onCurve(Rx, Ry, a, b, p)
    return Rx, Ry


def eccNegation(Px, Py, a, b, p):
    assert onCurve(Px, Py, a, b, p)
    Rx = Px
    Ry = np.mod(-Py, p)
    assert onCurve(Rx, Ry, a, b, p)
    return Rx, Ry


# EC point multiplication
def eccScalarMult(k, Px, Py, a, b, p):
    assert onCurve(Px, Py, a, b, p)
    if np.mod(k, p) == 0:
        return [], [], None, None
    if k < 0:
        # k * point = -k * (-point)
        Rx, Ry = eccNegation(Px, Py)
        return eccScalarMult(-k, Rx, Ry)
    rx, ry = Px, Py
    Rxs, Rys = [rx], [ry]
    for i in range(k):
        rx, ry = eccFiniteAddition(rx, ry, Px, Py, a, b, p)
        assert onCurve(rx, ry, a, b, p)
        Rxs.append(rx)
        Rys.append(ry)
    return Rxs, Rys, rx, ry


def bruteForceKey(Px, Py, pub_x, pub_y, a, b, p):
    for k in range(p):
        _, _, public_x, public_y = eccScalarMult(k, Px, Py, a, b, p)
        if public_x == pub_x and public_y == pub_y:
            return k

ecdh.png

ecdh_cipher.png

ecdh_cipher_alice.png

ecdh_cipher_bob.png

ecdh_curve.png

ecdh_curve_0.png

ecdh_curve_1.png

ecdh_curve_2.png

ecdh_discrete_addition.png

ecdh_finite.png

ecdh_public_key.png

ecdh_public_key_alice.png

ecdh_public_key_bob.png

ecdh_real.png

plotting.py

import numpy as np
import matplotlib.pyplot as plt


def plotR(a, b,
        line=None,
        vert=None,
        highlights=[],
        annotations=[],
        title=None,
        url=False,
        filename='curve.png'):
    # calculate the curve
    y, x = np.ogrid[-5:5:200j, -5:5:200j]
    # elliptic curve over real space
    cont = np.power(y, 2.)-np.power(x, 3.)-x*a-b
    # plot it
    fig, axs = plt.subplots(1, 1)
    if vert is not None:
        axs.axvline(x=vert, color='black', linestyle='--')
    axs.contour(x.ravel(), y.ravel(), cont, [0], colors='tab:blue')
    if line is not None:
        m, y0 = line
        Y = m*x + y0
        axs.plot(x.ravel(), Y[0], color='tab:orange')
    for xs, ys, color in highlights:
        axs.scatter(xs, ys, color=color, zorder=3)
    delta = 0.1
    for txt, Rx, Ry in annotations:
        axs.annotate(txt, (Rx+delta, Ry+delta), zorder=3)
    axs.set_title(title)
    axs.set_xlim([-5., 5.])
    axs.set_ylim([-5., 5.])
    axs.set_axisbelow(True)
    axs.grid()
    if url:
        plt.figtext(0.0, 0.01, 'https://mareknarozniak.com/2020/11/30/ecdh/', rotation=0)
    plt.tight_layout()
    fig.savefig(filename)


def plotF(sols_x, sols_y, p,
        Raxs=[],
        Rays=[],
        highlights=[],
        annotations=[],
        title=None,
        url=False,
        filename='curve_finite_field.png'):
    fig, axs = plt.subplots(1, 1)
    axs.scatter(sols_x, sols_y, zorder=1)
    delta = 0.1
    for i in range(p):
        axs.annotate(chr(i + 65), (sols_x[i] + delta, sols_y[i] + delta))
    for xs, ys, color in highlights:
        axs.scatter(xs, ys, color=color)
    for txt, Rx, Ry in annotations:
        axs.annotate(txt, (Rx - 4*delta, Ry - 6*delta))
    if None in Raxs and None in Rays:
        Raxs.remove(None)
        Rays.remove(None)
    for k in range(len(Raxs)-1):
        ox, oy = Raxs[k], Rays[k]
        dx, dy = Raxs[k+1], Rays[k+1]
        axs.arrow(ox, oy, dx-ox, dy-oy, head_width=0.2, length_includes_head=True, color='black')
    axs.set_title(title)
    axs.set_xticks(range(p))
    axs.set_yticks(range(p))
    axs.set_xlim([-1, p])
    axs.set_ylim([-1, p])
    axs.set_axisbelow(True)
    axs.grid()
    if url:
        plt.figtext(0.0, 0.01, 'https://mareknarozniak.com/2020/11/30/ecdh/', rotation=0)
    plt.tight_layout()
    fig.savefig(filename)

run.py

import matplotlib.pyplot as plt
import numpy as np

from ecc import eccContAddition, eccFiniteAddition, eccScalarMult, bruteForceKey
from plotting import plotR, plotF



curves = [
    ('secp256k1, $a=0,b=7$', 0, 7, 1.3, 2.1),
    ('$a=-3,b=1$', -3, 1, -1.5, 2.1),
    ('$a=-3,b=1$', -3, 8, -2., 2.1)]

for i, (title, a, b, Px, Qx) in enumerate(curves):
    Py = np.sqrt(np.power(Px, 3.)+Px*a+b)
    Qy = np.sqrt(np.power(Qx, 3.)+Qx*a+b)
    Rx, Ry, m = eccContAddition(Px, Py, Qx, Qy, a, b)
    Y0 = Py-m*Px
    url = i == 0
    highlights = [
        ([Px, Qx, Rx, Rx], [Py, Qy, Ry, -Ry], 'tab:orange')]
    annotations = [
        ('P', Px, Py),
        ('Q', Qx, Qy),
        ('-R', Rx, Ry),
        ('R', Rx, -Ry)]
    plotR(a, b,
        title=title,
        line=(m, Y0),
        vert=Rx,
        highlights=highlights,
        annotations=annotations,
        url=url,
        filename='ecdh_curve_{0}.png'.format(str(i)))

a = 0
b = 7
p = 17


plotR(a, b,
    title='secp256k1, $a=0,b=7$',
    url=True,
    filename='ecdh_real.png'.format(str(i)))

# secp256k1 curve over finite field
sols_x = []
sols_y = []
for cx in range(p):
    for cy in range(p):
        if np.mod(cy**2-cx**3-cx*a-b, p) == 0:
            sols_x.append(cx)
            sols_y.append(cy)

plotF(sols_x, sols_y, p,
    title='secp256k1, $a=0,b=7$',
    filename='ecdh_finite.png')

idx_p = 3
idx_q = 5

Px, Py, Qx, Qy = sols_x[idx_p], sols_y[idx_p], sols_x[idx_q], sols_y[idx_q]
Rx, Ry = eccFiniteAddition(Px, Py, Qx, Qy, a, b, p)

highlights = [
    ([Px, Qx, Rx], [Py, Qy, Ry], 'tab:orange')]
annotations = [
    ('(P)', Px, Py),
    ('(Q)', Qx, Qy),
    ('(R)', Rx, Ry)]
plotF(sols_x, sols_y, p,
    title=r'$D+F=B$',
    highlights=highlights,
    annotations=annotations,
    filename='ecdh_discrete_addition.png')

generator = 16
Px, Py = sols_x[generator], sols_y[generator]

private_keys = [4, 13]
names = ['alice', 'bob']
communicate_with = [1, 0]
public_keys = []

for i, (private_key, name) in enumerate(zip(private_keys, names)):
    Raxs, Rays, public_x, public_y = eccScalarMult(private_key, Px, Py, a, b, p)
    public_key = (public_x, public_y)
    public_keys.append(public_key)
    highlights = [
        ([Px], [Py], 'tab:green'),
        ([public_x], [public_y], 'tab:orange')]
    url = i == 0
    plotF(sols_x, sols_y, p,
        Raxs=Raxs,
        Rays=Rays,
        highlights=highlights,
        title='{0} Public Key'.format(name.capitalize()),
        url=url,
        filename='ecdh_public_key_{0}.png'.format(name))


# Diffie-Hellman Key Exchange
for private_key, name, partner in zip(private_keys, names, communicate_with):
    public_x, public_y = public_keys[partner]
    Raxs, Rays, secret_x, secret_y = eccScalarMult(private_key, public_x, public_y, a, b, p)
    highlights = [
        ([public_x], [public_y], 'tab:green'),
        ([secret_x], [secret_y], 'tab:orange')]
    url = i == 0
    plotF(sols_x, sols_y, p,
        Raxs=Raxs,
        Rays=Rays,
        highlights=highlights,
        title='{0} Cipher'.format(name.capitalize()),
        url=url,
        filename='ecdh_cipher_{0}.png'.format(name))


# break the Bob key!
Px, Py = sols_x[generator], sols_y[generator]
pub_x, pub_y = public_keys[1]
k = bruteForceKey(Px, Py, pub_x, pub_y, a, b, p)
print('Bob private key: ', k)

# break the Alice key!
pub_x, pub_y = public_keys[0]
k = bruteForceKey(Px, Py, pub_x, pub_y, a, b, p)
print('Alice private key: ', k)

# convert ecdh_curve_0.png ecdh_curve_1.png ecdh_curve_2.png +append ecdh_curve.png

# convert ecdh_public_key_alice.png ecdh_public_key_bob.png +append  ecdh_public_key.png

# convert ecdh_cipher_alice.png ecdh_cipher_bob.png +append ecdh_cipher.png

# convert ecdh_real.png ecdh_finite.png +append ecdh.png