Elliptic curves

An elliptic curve $E$ is a set of points satisfying an equation of form $$ y^2 = x^3 + ax + b, $$ where $4a^3 + 27b^2 \ne 0$, together with a "point at infinity", denoted by $\mathcal{O}$.

Let us first draw a graph of a cubic polynomial, followed by the corresponding elliptic curve.

In [ ]:
@interact
def _(a=slider(-5, 5, default=-2),
      b=slider(-5, 5, default=-2),
      xr=range_slider(-10, 10, default=(-3, 3), label="x range"),
      yr=range_slider(-10, 10, default=(-5, 5), label="y range")):
    x, y = var('x y')
    E = x^3 + a*x + b
    xlim = (x, *xr)
    ylim = (y, *yr)
    show(plot(E, xlim))
    show(implicit_plot(y^2 == E, xlim, ylim, frame=False, axes=True))

Let us now demonstrate the group addition operation on the elliptic curve.

In [ ]:
@interact
def _(a=slider(-5, 5, default=-2),
      b=slider(-5, 5, default=-2),
      Px=slider(-10, 10, default=2),
      Ps=slider(["-", "+"], default="+", label="Py sign"),
      Qx=slider(-10, 10, default=3),
      Qs=slider(["-", "+"], default="-", label="Qy sign"),
      xr=range_slider(-10, 10, default=(-3, 3), label="x range"),
      yr=range_slider(-10, 10, default=(-5, 5), label="y range")):
    x, y = var('x y')
    s = {"+": 1, "-": -1}
    E = x^3 + a*x + b
    pts = []
    Py = s[Ps] * sqrt(E.subs(x == Px))
    if Py.is_real():
        pts.append((Px, Py))
    Qy = s[Qs] * sqrt(E.subs(x == Qx))
    if Qy.is_real():
        pts.append((Qx, Qy))
    xlim = (x, *xr)
    ylim = (y, *yr)
    p = implicit_plot(y^2 == E, xlim, ylim, frame=False, axes=True)
    if pts:
        p += point(pts, color="red", size=50)
    if len(pts) > 1:
        if pts[0] == pts[1] and Py != 0:
            lm = (3*Px^2 + a)/(2*Py)
        elif Px != Qx:
            lm = (Py - Qy)/(Px - Qx)
        else:
            pts += [(Px, yr[0]), (Px, yr[1])]
        if len(pts) == 2:
            Rx = lm^2 - Px - Qx
            Ry = lm * (Px - Rx) - Py
            R = [(Rx, Ry), (Rx, -Ry)]
            pts.append(R[1])
            p += point(R[0], color="magenta", size=50)
            p += point(R[1], color="brown", size=50)
            if Ry != 0:
                p += line2d(R, color="cyan")
        p += line2d(pts, color="green")
    show(p, xmin=xr[0], xmax=xr[1], ymin=yr[0], ymax=yr[1])