Chapter 31: Elliptic Curves and the Group Law (SageMath)#

Python Version Available

This notebook contains the SageMath implementation of Chapter 31, using native EllipticCurve, PolynomialRing, and built-in algebraic geometry. A pure Python/NumPy version is available in Chapter 31: Elliptic Curves and the Group Law.

31.1 Historical Background#

The application of elliptic curves to cryptography emerged from two independent proposals in the mid-1980s:

  • In 1985, Victor Miller presented the idea of using the group of rational points on an elliptic curve over a finite field as the basis for a public-key cryptosystem at a crypto conference, with the paper appearing in 1986.

  • In 1987, Neal Koblitz independently proposed essentially the same idea, publishing in Mathematics of Computation.

Both recognized a key insight: the discrete logarithm problem on an elliptic curve group appeared to be significantly harder than the corresponding problem in the multiplicative group of a finite field. This meant that the same level of security could be achieved with dramatically smaller key sizes. For example, a 256-bit elliptic curve key provides roughly the same security as a 3072-bit RSA key.

Elliptic curves had been studied as abstract mathematical objects for over a century – by Weierstrass, Mordell, Weil, and many others – but their application to practical cryptography was revolutionary. Today, Elliptic Curve Cryptography (ECC) underpins TLS 1.3, Bitcoin, Signal, and numerous other protocols.

Note

This chapter focuses on the mathematics of elliptic curves: defining the group law, implementing arithmetic over finite fields, and verifying the group axioms computationally. Cryptographic protocols built on these foundations (ECDH, ECDSA, etc.) are treated in later chapters.

31.2 The Weierstrass Equation#

The discriminant condition \(\Delta \neq 0\) ensures the cubic \(x^3 + ax + b\) has no repeated roots, which is equivalent to the curve having no cusps or self-intersections. When \(\Delta = 0\), the curve degenerates into a singular curve that does not carry a group structure in the usual sense.

Over the real numbers \(\mathbb{R}\), elliptic curves appear as smooth plane curves. Over a finite field \(\mathbb{F}_p\), the “curve” is a discrete set of points. In both cases, the same algebraic formulas define the group operation.

# Plot y^2 = x^3 + ax + b over the reals for several parameter choices
# SageMath's EllipticCurve provides a native .plot() method

curves_params = [(-1, 0), (0, 1), (-2, 2)]
titles = [r"$y^2 = x^3 - x$", r"$y^2 = x^3 + 1$", r"$y^2 = x^3 - 2x + 2$"]

plots = []
for (a_coeff, b_coeff), title in zip(curves_params, titles):
    E = EllipticCurve(QQ, [a_coeff, b_coeff])
    disc = E.discriminant()
    p = E.plot(xmin=-3, xmax=3, ymin=-5, ymax=5)
    p += text(title, (0, 4.5), fontsize=12, color='black')
    p += text(r"$\Delta = %s$" % disc, (-1.5, -4), fontsize=10, color='brown')
    plots.append(p)

show(graphics_array(plots, nrows=1, ncols=3), figsize=[14, 4.5])
../_images/93f1ab2a5bc68ca3a7977d51320688aedac93e9ff01f505d6d87529e703f222e.png

31.3 The Group Law: Point Addition#

The remarkable fact about elliptic curves is that their points form an abelian group under a geometric operation called “point addition.”

# Geometric illustration of point addition on y^2 = x^3 - 2x + 2
E = EllipticCurve(QQ, [-2, 2])

# Choose two points on the curve (work numerically for the plot)
a_coeff, b_coeff = -2, 2
x1_val = -1.0
y1_val = sqrt(RR(x1_val^3 + a_coeff*x1_val + b_coeff))
x2_val = 1.5
y2_val = sqrt(RR(x2_val^3 + a_coeff*x2_val + b_coeff))

# Compute the addition geometrically
lam = (y2_val - y1_val) / (x2_val - x1_val)
x3_val = lam^2 - x1_val - x2_val
y3_intersect = lam*(x1_val - x3_val) - y1_val  # P+Q y-coordinate
y3_third = -y3_intersect  # third intersection point (before reflection)

# Build the plot
p = E.plot(xmin=-2.0, xmax=3.5, ymin=-5, ymax=5, color='blue', thickness=2)

# Secant line through P and Q
t = var('t')
p += parametric_plot((t, lam*(t - x1_val) + y1_val), (t, -2, 3.5),
                      color='red', linestyle='dashed', thickness=1.2)

# Vertical reflection line
p += line([(x3_val, y3_third), (x3_val, y3_intersect)],
          color='green', linestyle='dashed', thickness=1.2)

# Plot points P, Q, R', P+Q
p += point([(x1_val, y1_val)], color='red', size=60, zorder=5)
p += text("$P$", (x1_val - 0.3, y1_val + 0.4), fontsize=14, color='red')

p += point([(x2_val, y2_val)], color='green', size=60, zorder=5)
p += text("$Q$", (x2_val + 0.25, y2_val + 0.3), fontsize=14, color='green')

p += point([(x3_val, y3_third)], color='black', size=50, zorder=5)
p += text("$R'$", (x3_val + 0.3, y3_third - 0.4), fontsize=13, color='black')

p += point([(x3_val, y3_intersect)], color='purple', size=80, zorder=5, marker='^')
p += text("$P+Q$", (x3_val + 0.35, y3_intersect + 0.4), fontsize=14, color='purple')

p.axes_labels(['$x$', '$y$'])
show(p, figsize=[7, 7], title="Point Addition on an Elliptic Curve")

print(f"P = ({float(x1_val):.4f}, {float(y1_val):.4f})")
print(f"Q = ({float(x2_val):.4f}, {float(y2_val):.4f})")
print(f"P + Q = ({float(x3_val):.4f}, {float(y3_intersect):.4f})")
../_images/52b82d8542e8d2b298f7efb986160665265d84d5fd78815b41ff7915bebe3330.png
P = (-1.0000, 1.7321)
Q = (1.5000, 1.5411)
P + Q = (-0.4942, -1.6934)

31.4 The Group Structure#

Closure and commutativity are relatively straightforward from the addition formulas. The identity and inverse properties follow from the definitions. Associativity is the hardest to prove – it requires a careful algebraic verification (or an argument via the Riemann-Roch theorem). We will verify it computationally in Section 31.11.

31.5 Symbolic Addition Formulas with SageMath#

One of the most elegant features of SageMath is its ability to derive the addition formulas symbolically. Rather than implementing the chord-and-tangent formulas by hand, we define a generic elliptic curve over a quotient ring and let SageMath compute the formulas for us.

Note

This approach uses PolynomialRing, ideal(), and quotient_ring() to work in the ring where \(y_1^2 = x_1^3 + ax_1 + b\) and \(y_2^2 = x_2^3 + ax_2 + b\) are identities. The resulting addition expressions are the general formulas valid for any elliptic curve in short Weierstrass form.

#some technical setup
R.<x1,y1,x2,y2,a,b> = PolynomialRing(QQ, 'x1,y1,x2,y2,a,b')
Ideal1 = ideal(y1^2-(x1^3+a*x1+b),y2^2-(x2^3+a*x2+b))
Rquot.<X1,Y1,X2,Y2,A,B>=R.quotient_ring(Ideal1)
Fquot=Rquot.fraction_field()
E=EllipticCurve(Fquot,[A,B]);
#two abstract points on the curve E
P=E.point([X1,Y1])
Q=E.point([X2,Y2])
P
(X1 : Y1 : 1)
Q
(X2 : Y2 : 1)
#addition of two distinct points
pretty_print(P+Q)
\(\displaystyle \left(\frac{X_{1}^{2} X_{2} + X_{1} X_{2}^{2} - 2 Y_{1} Y_{2} + X_{1} A + X_{2} A + 2 B}{X_{1}^{2} - 2 X_{1} X_{2} + X_{2}^{2}} : \frac{-3 X_{1} Y_{1} X_{2}^{2} + 3 X_{1}^{2} X_{2} Y_{2} + Y_{1}^{2} Y_{2} - Y_{1} Y_{2}^{2} - X_{1} Y_{1} A - 2 Y_{1} X_{2} A + 2 X_{1} Y_{2} A + X_{2} Y_{2} A - 3 Y_{1} B + 3 Y_{2} B}{-3 X_{1}^{2} X_{2} + 3 X_{1} X_{2}^{2} + Y_{1}^{2} - Y_{2}^{2} - X_{1} A + X_{2} A} : 1\right)\)
#doubling formula for a point P
pretty_print(2*P)
\(\displaystyle \left(\frac{X_{1} Y_{1}^{2} - 3 X_{1}^{2} A + A^{2} - 9 X_{1} B}{4 Y_{1}^{2}} : \frac{Y_{1}^{4} + 3 X_{1} Y_{1}^{2} A - 9 X_{1}^{2} A^{2} - A^{3} + 18 Y_{1}^{2} B - 27 X_{1} A B - 27 B^{2}}{8 Y_{1}^{3}} : 1\right)\)
#opposite of a point
pretty_print(-P)
\(\displaystyle \left(X_{1} : -Y_{1} : 1\right)\)
#triplication formula of a point P
pretty_print(3*P)
\(\displaystyle \left(\frac{256 Y_{1}^{14} - 3840 X_{1} Y_{1}^{12} A + 14592 X_{1}^{2} Y_{1}^{10} A^{2} - 1792 Y_{1}^{10} A^{3} + 4096 X_{1} Y_{1}^{8} A^{4} - 25344 Y_{1}^{12} B + 50688 X_{1} Y_{1}^{10} A B - 20736 X_{1}^{2} Y_{1}^{8} A^{2} B + 3840 Y_{1}^{8} A^{3} B + 62208 Y_{1}^{10} B^{2} - 34560 X_{1} Y_{1}^{8} A B^{2} - 20736 Y_{1}^{8} B^{3}}{2304 X_{1}^{2} Y_{1}^{12} + 4608 Y_{1}^{12} A - 3840 X_{1} Y_{1}^{10} A^{2} - 3840 X_{1}^{2} Y_{1}^{8} A^{3} + 13824 X_{1}^{2} Y_{1}^{10} B + 256 Y_{1}^{8} A^{4} + 9216 Y_{1}^{10} A B - 20736 X_{1} Y_{1}^{8} A^{2} B + 20736 X_{1}^{2} Y_{1}^{8} B^{2} - 13824 Y_{1}^{8} A B^{2}} : \frac{-65536 Y_{1}^{25} - 1179648 X_{1} Y_{1}^{23} A + 14745600 X_{1}^{2} Y_{1}^{21} A^{2} - 19660800 Y_{1}^{21} A^{3} + 37945344 X_{1} Y_{1}^{19} A^{4} - 12386304 X_{1}^{2} Y_{1}^{17} A^{5} - 14155776 Y_{1}^{23} B + 81395712 X_{1} Y_{1}^{21} A B - 159252480 X_{1}^{2} Y_{1}^{19} A^{2} B + 196608 Y_{1}^{17} A^{6} + 110886912 Y_{1}^{19} A^{3} B - 100859904 X_{1} Y_{1}^{17} A^{4} B + 159252480 Y_{1}^{21} B^{2} - 329121792 X_{1} Y_{1}^{19} A B^{2} + 238878720 X_{1}^{2} Y_{1}^{17} A^{2} B^{2} - 84934656 Y_{1}^{17} A^{3} B^{2} - 254803968 Y_{1}^{19} B^{3} + 286654464 X_{1} Y_{1}^{17} A B^{3} + 143327232 Y_{1}^{17} B^{4}}{-1769472 Y_{1}^{24} - 3538944 X_{1} Y_{1}^{22} A + 1769472 X_{1}^{2} Y_{1}^{20} A^{2} + 7077888 Y_{1}^{20} A^{3} - 4128768 X_{1} Y_{1}^{18} A^{4} - 4128768 X_{1}^{2} Y_{1}^{16} A^{5} - 14155776 Y_{1}^{22} B - 10616832 X_{1} Y_{1}^{20} A B + 31850496 X_{1}^{2} Y_{1}^{18} A^{2} B + 65536 Y_{1}^{16} A^{6} + 21233664 Y_{1}^{18} A^{3} B - 33619968 X_{1} Y_{1}^{16} A^{4} B - 31850496 Y_{1}^{20} B^{2} + 31850496 X_{1} Y_{1}^{18} A B^{2} + 79626240 X_{1}^{2} Y_{1}^{16} A^{2} B^{2} - 28311552 Y_{1}^{16} A^{3} B^{2} + 95551488 X_{1} Y_{1}^{16} A B^{3} + 47775744 Y_{1}^{16} B^{4}} : 1\right)\)

The Power of Symbolic Computation

Notice that SageMath derives the full addition, doubling, negation, and triplication formulas automatically from the curve equation. In the Python version, these formulas had to be implemented manually. SageMath’s EllipticCurve handles all the algebraic complexity internally, allowing us to focus on the mathematics rather than the implementation.

31.6 Elliptic Curves over Finite Fields#

When we work over a finite field \(\mathbb{F}_p\) (where \(p > 3\) is prime), all arithmetic is performed modulo \(p\). The curve equation becomes:

\[ E(\mathbb{F}_p): \quad y^2 \equiv x^3 + ax + b \pmod{p},\]

and the set of points \(E(\mathbb{F}_p)\) is finite. The following theorem gives a tight bound on how many points the curve can have.

In other words, the number of points is approximately \(p + 1\), with an error of at most \(2\sqrt{p}\). The quantity \(t = p + 1 - \#E(\mathbb{F}_p)\) is called the trace of Frobenius.

Tip

Hasse’s theorem tells us that for large \(p\), the group \(E(\mathbb{F}_p)\) has roughly \(p\) elements. This is crucial for cryptographic applications: if we choose \(p \approx 2^{256}\), the group has approximately \(2^{256}\) elements, making brute-force attacks infeasible.

31.7 Creating Elliptic Curves in SageMath#

SageMath’s EllipticCurve constructor works seamlessly over any field. We simply specify the base field and the coefficients [a, b] for the short Weierstrass form \(y^2 = x^3 + ax + b\). This replaces the manual Python class from the Python version with a single line of code.

# Define an elliptic curve over F_97
F = GF(97)
E = EllipticCurve(F, [2, 3])
print(E)
print(f"Discriminant: {E.discriminant()}")
print(f"j-invariant: {E.j_invariant()}")
print(f"#E(F_97) = {E.cardinality()}")
Elliptic Curve defined by y^2 = x^3 + 2*x + 3 over Finite Field of size 97
Discriminant: 62
j-invariant: 36
#E(F_97) = 100

31.8 Verifying Curve Membership#

SageMath can enumerate all rational points on a curve over a finite field. Let us list them and verify membership.

# Enumerate all points on E(F_97)
E = EllipticCurve(GF(97), [2, 3])
pts = E.points()

print(f"All {len(pts)} points satisfy the curve equation (by construction).")
print(f"\nFirst 10 points:")
for P in pts[:10]:
    print(f"  {P}")

# Verify a few points manually
print(f"\nSample verifications on E: y^2 = x^3 + 2x + 3 (mod 97):")
for P in pts[1:8]:
    x, y = P.xy()
    lhs = y^2 % 97
    rhs = (x^3 + 2*x + 3) % 97
    print(f"  P = ({int(x):2d}, {int(y):2d})  |  y^2 mod 97 = {int(lhs):2d}  |  "
          f"x^3+2x+3 mod 97 = {int(rhs):2d}  |  match: {lhs == rhs}")
All 100 points satisfy the curve equation (by construction).

First 10 points:
  (0 : 1 : 0)
  (0 : 10 : 1)
  (0 : 87 : 1)
  (1 : 43 : 1)
  (1 : 54 : 1)
  (3 : 6 : 1)
  (3 : 91 : 1)
  (4 : 47 : 1)
  (4 : 50 : 1)
  (10 : 21 : 1)

Sample verifications on E: y^2 = x^3 + 2x + 3 (mod 97):
  P = ( 0, 10)  |  y^2 mod 97 =  3  |  x^3+2x+3 mod 97 =  3  |  match: True
  P = ( 0, 87)  |  y^2 mod 97 =  3  |  x^3+2x+3 mod 97 =  3  |  match: True
  P = ( 1, 43)  |  y^2 mod 97 =  6  |  x^3+2x+3 mod 97 =  6  |  match: True
  P = ( 1, 54)  |  y^2 mod 97 =  6  |  x^3+2x+3 mod 97 =  6  |  match: True
  P = ( 3,  6)  |  y^2 mod 97 = 36  |  x^3+2x+3 mod 97 = 36  |  match: True
  P = ( 3, 91)  |  y^2 mod 97 = 36  |  x^3+2x+3 mod 97 = 36  |  match: True
  P = ( 4, 47)  |  y^2 mod 97 = 75  |  x^3+2x+3 mod 97 = 75  |  match: True

31.9 Plotting the Curve over \(\mathbb{F}_p\) and \(\mathbb{R}\)#

Over the reals, \(E\) is a smooth plane curve. Over \(\mathbb{F}_p\), it becomes a discrete scatter of points. Let us compare the two views for the same curve parameters using SageMath’s point() and EllipticCurve.plot().

# Plot E(F_97) as a scatter alongside E(R) using SageMath plotting

# Left: E(F_97)
E_Fp = EllipticCurve(GF(97), [2, 3])
pts = E_Fp.points()
affine_pts = [(ZZ(P.xy()[0]), ZZ(P.xy()[1])) for P in pts if P != E_Fp(0)]
p_scatter = point(affine_pts, color='blue', size=12, alpha=0.7)
p_scatter += text(
    r"$E(\mathbb{F}_{97}): y^2 = x^3 + 2x + 3$" + f"\n|E| = {len(pts)} points",
    (50, 90), fontsize=11, color='black'
)

# Right: E(R)
E_R = EllipticCurve(QQ, [2, 3])
p_real = E_R.plot(xmin=-2, xmax=5, ymin=-12, ymax=12, color='blue', thickness=2)
p_real += text(r"$E(\mathbb{R}): y^2 = x^3 + 2x + 3$", (1.5, 11), fontsize=11, color='black')

show(graphics_array([p_scatter, p_real], nrows=1, ncols=2), figsize=[14, 5.5])
../_images/370028e9c3cb2447a06f5fafea1c253818e33568b8025cfd28a6f6668b0fca25.png

31.10 Step-by-Step Point Addition over \(\mathbb{F}_p\)#

In SageMath, point addition is as simple as using the + operator. Let us pick two points on \(E(\mathbb{F}_{97})\) and demonstrate the elegance of SageMath’s built-in arithmetic, while also tracing the steps manually for comparison.

# Point addition in SageMath -- clean and simple
E = EllipticCurve(GF(97), [2, 3])
pts = E.points()

# Pick two distinct affine points
P = pts[1]
Q = pts[5]
print(f"Curve: {E}")
print(f"P = {P}")
print(f"Q = {Q}")

# SageMath addition -- just use +
R = P + Q
print(f"\nP + Q = {R}")
print(f"On curve? {R in E}")

# Manual step-by-step for comparison
x1, y1 = P.xy()
x2, y2 = Q.xy()
F = GF(97)

dx = F(x2 - x1)
dy = F(y2 - y1)
lam = dy / dx
x3 = lam^2 - x1 - x2
y3 = lam*(x1 - x3) - y1

print(f"\nManual step-by-step:")
print(f"  dx = {dx}")
print(f"  dy = {dy}")
print(f"  lambda = dy/dx = {lam}")
print(f"  x3 = lambda^2 - x1 - x2 = {x3}")
print(f"  y3 = lambda*(x1 - x3) - y1 = {y3}")
print(f"  Manual result: ({x3}, {y3})")
print(f"  Matches SageMath: {E((x3, y3)) == R}")
Curve: Elliptic Curve defined by y^2 = x^3 + 2*x + 3 over Finite Field of size 97
P = (0 : 10 : 1)
Q = (3 : 6 : 1)

P + Q = (85 : 71 : 1)
On curve? True

Manual step-by-step:
  dx = 3
  dy = 93
  lambda = dy/dx = 31
  x3 = lambda^2 - x1 - x2 = 85
  y3 = lambda*(x1 - x3) - y1 = 71
  Manual result: (85, 71)
  Matches SageMath: True

31.11 Verifying the Group Axioms#

We now computationally verify all five group axioms on a small curve \(E(\mathbb{F}_p)\). SageMath’s native point arithmetic via the + operator makes this verification concise and readable. For a small prime, we can exhaustively check associativity for all triples – a luxury we will not have for cryptographic-size primes.

# Use a smaller prime for exhaustive verification
E = EllipticCurve(GF(11), [1, 6])
pts = E.points()
n = len(pts)
print(f"Curve: {E}")
print(f"#E(F_11) = {n}")
print(f"Points: {pts}")

O = E(0)  # point at infinity

# 1. CLOSURE: P + Q is on the curve for all P, Q
closure_ok = all((P + Q) in E for P in pts for Q in pts)
print(f"\n1. Closure: {'PASS' if closure_ok else 'FAIL'}")

# 2. IDENTITY: P + O = P for all P
identity_ok = all(P + O == P and O + P == P for P in pts)
print(f"2. Identity (O): {'PASS' if identity_ok else 'FAIL'}")

# 3. INVERSES: P + (-P) = O for all P
inverse_ok = all(P + (-P) == O for P in pts)
print(f"3. Inverses: {'PASS' if inverse_ok else 'FAIL'}")

# 4. COMMUTATIVITY: P + Q = Q + P for all P, Q
comm_ok = all(P + Q == Q + P for P in pts for Q in pts)
print(f"4. Commutativity: {'PASS' if comm_ok else 'FAIL'}")

# 5. ASSOCIATIVITY: (P + Q) + R = P + (Q + R) for all P, Q, R
assoc_ok = True
count = 0
for P in pts:
    for Q in pts:
        for R in pts:
            if (P + Q) + R != P + (Q + R):
                assoc_ok = False
                print(f"  ASSOC FAIL: P={P}, Q={Q}, R={R}")
            count += 1
print(f"5. Associativity: {'PASS' if assoc_ok else 'FAIL'}  ({count} triples checked)")
Curve: Elliptic Curve defined by y^2 = x^3 + x + 6 over Finite Field of size 11
#E(F_11) = 13
Points: [(0 : 1 : 0), (2 : 4 : 1), (2 : 7 : 1), (3 : 5 : 1), (3 : 6 : 1), (5 : 2 : 1), (5 : 9 : 1), (7 : 2 : 1), (7 : 9 : 1), (8 : 3 : 1), (8 : 8 : 1), (10 : 2 : 1), (10 : 9 : 1)]

1. Closure: PASS
2. Identity (O): PASS
3. Inverses: PASS
4. Commutativity: PASS
5. Associativity: PASS  (2197 triples checked)

31.12 Scalar Multiplication and the Double-and-Add Algorithm#

Given a point \(P\) and a positive integer \(k\), the scalar multiple \(kP\) is defined as:

\[ kP = \underbrace{P + P + \cdots + P}_{k \text{ times}}.\]

Computing this naively requires \(k - 1\) additions, which is infeasible when \(k \approx 2^{256}\). The double-and-add algorithm (analogous to square-and-multiply for modular exponentiation) computes \(kP\) in \(O(\log k)\) group operations by processing the binary representation of \(k\).

In SageMath, scalar multiplication is built in: simply write k * P.

# Scalar multiplication in SageMath -- just use k * P
E = EllipticCurve(GF(97), [2, 3])
pts = E.points()
P = pts[1]

print(f"Curve: {E}")
print(f"Base point P = {P}")
print(f"\nScalar multiples of P:")
print(f"{'k':>4s}  {'kP':>30s}  {'on curve?':>10s}")
print("-" * 50)

for k in range(1, 20):
    kP = k * P
    on_curve = kP in E
    print(f"{int(k):4d}  {str(kP):>30s}  {str(on_curve):>10s}")

# Find the order of P using SageMath's built-in method
order_P = P.order()
n_E = E.cardinality()
print(f"\nOrder of P: {order_P}")
print(f"#E(F_97) = {n_E}")
print(f"Order of P divides #E? {n_E % order_P == 0}")
Curve: Elliptic Curve defined by y^2 = x^3 + 2*x + 3 over Finite Field of size 97
Base point P = (0 : 10 : 1)

Scalar multiples of P:
   k                              kP   on curve?
--------------------------------------------------
   1                    (0 : 10 : 1)        True
   2                   (65 : 32 : 1)        True
   3                   (23 : 24 : 1)        True
   4                   (52 : 68 : 1)        True
   5                   (88 : 56 : 1)        True
   6                   (95 : 66 : 1)        True
   7                   (10 : 76 : 1)        True
   8                   (84 : 37 : 1)        True
   9                   (11 : 80 : 1)        True
  10                   (80 : 10 : 1)        True
  11                   (17 : 87 : 1)        True
  12                   (74 : 20 : 1)        True
  13                   (87 : 70 : 1)        True
  14                   (46 : 72 : 1)        True
  15                   (53 : 24 : 1)        True
  16                   (92 : 81 : 1)        True
  17                    (1 : 43 : 1)        True
  18                   (21 : 73 : 1)        True
  19                   (85 : 26 : 1)        True

Order of P: 50
#E(F_97) = 100
Order of P divides #E? True

31.13 Tracing the Double-and-Add Algorithm#

While SageMath handles scalar multiplication internally, it is instructive to trace the double-and-add algorithm step by step to understand the efficiency gain.

# Trace the right-to-left double-and-add algorithm
E = EllipticCurve(GF(97), [2, 3])
pts = E.points()
P = pts[1]

k = 53
print(f"Computing {k} * P where P = {P}")
print(f"Binary representation of {k}: {bin(k)}")
print(f"Binary digits (MSB first): {[int(b) for b in bin(k)[2:]]}")
print()

O = E(0)
result = O
addend = P
k_copy = k
step = 0

print(f"{'Step':>4s}  {'Bit':>3s}  {'Action':>15s}  {'Result':>30s}  {'Addend':>30s}")
print("-" * 90)

while k_copy > 0:
    bit = k_copy & 1
    if bit:
        result = result + addend
        action = "add"
    else:
        action = "skip"

    print(f"{int(step):4d}  {int(bit):3d}  {action:>15s}  {str(result):>30s}  {str(addend):>30s}")

    addend = 2 * addend
    k_copy >>= 1
    step += 1

print(f"\nFinal result: {k}*P = {result}")
print(f"Verification: {k * P}")
print(f"Match: {result == k * P}")
print(f"\nTotal group operations: {step} doubles + {bin(k).count('1')} adds")
print(f"Naive method would require {k - 1} additions")
Computing 53 * P where P = (0 : 10 : 1)
Binary representation of 53: 0b110101
Binary digits (MSB first): [1, 1, 0, 1, 0, 1]

Step  Bit           Action                          Result                          Addend
------------------------------------------------------------------------------------------
   0    1              add                    (0 : 10 : 1)                    (0 : 10 : 1)
   1    0             skip                    (0 : 10 : 1)                   (65 : 32 : 1)
   2    1              add                   (88 : 56 : 1)                   (52 : 68 : 1)
   3    0             skip                   (88 : 56 : 1)                   (84 : 37 : 1)
   4    1              add                   (47 : 79 : 1)                   (92 : 81 : 1)
   5    1              add                   (23 : 24 : 1)                   (21 : 24 : 1)

Final result: 53*P = (23 : 24 : 1)
Verification: (23 : 24 : 1)
Match: True

Total group operations: 6 doubles + 4 adds
Naive method would require 52 additions

31.14 The Elliptic Curve Discrete Logarithm Problem (ECDLP)#

The best known algorithms for the general ECDLP on a well-chosen curve run in \(O(\sqrt{n})\) time (e.g., Pollard’s rho algorithm). There is no known sub-exponential algorithm (unlike the ordinary DLP in \(\mathbb{F}_p^*\), which can be solved in sub-exponential time via the number field sieve). This is the fundamental reason ECC can use shorter keys than RSA or classical Diffie-Hellman.

Warning

Not all elliptic curves are suitable for cryptography. Curves susceptible to the MOV attack (embedding degree too small), anomalous curves (\(\#E(\mathbb{F}_p) = p\)), or curves with other special structure can have their ECDLP solved efficiently. Standardized curves (NIST P-256, Curve25519, etc.) are carefully chosen to avoid these pitfalls.

# ECDLP demonstration with SageMath
E = EllipticCurve(GF(97), [2, 3])
n_E = E.cardinality()

# Find a generator (point with maximal order)
G = None
for candidate in E.points():
    if candidate.order() == n_E:
        G = candidate
        break

if G is None:
    # Group may not be cyclic; find point with largest order
    G = max(E.points(), key=lambda P: P.order())

print(f"Base point G = {G}, order = {G.order()}")

# Choose a secret scalar and compute Q = kG
secret_k = 42
Q = secret_k * G
print(f"\nSecret k = {secret_k}")
print(f"Public Q = kG = {Q}")

# SageMath can solve the discrete log on small curves
# Using the built-in discrete_log method
k_found = G.discrete_log(Q)
print(f"\nSageMath discrete_log finds: k = {k_found}")
print(f"Verification: {k_found * G == Q}")

print(f"\nFor a 256-bit curve, brute force requires ~2^256 operations.")
print(f"Pollard's rho requires ~2^128 operations -- still infeasible.")
Base point G = (0 : 10 : 1), order = 50

Secret k = 42
Public Q = kG = (84 : 60 : 1)

SageMath discrete_log finds: k = 42
Verification: True

For a 256-bit curve, brute force requires ~2^256 operations.
Pollard's rho requires ~2^128 operations -- still infeasible.
/var/folders/z7/wp7m8p7x1250jzvklw5z24mm0000gn/T/ipykernel_27462/1454447375.py:26: DeprecationWarning: The syntax P.discrete_log(Q) is being replaced by Q.log(P) to make the argument ordering of logarithm methods in Sage uniform. Please update your code.
See https://github.com/sagemath/sage/issues/37150 for details.
  k_found = G.discrete_log(Q)

31.15 Verifying Hasse’s Theorem#

We now verify Hasse’s bound \(|\#E(\mathbb{F}_p) - p - 1| \leq 2\sqrt{p}\) by counting points on a fixed curve for several primes. SageMath’s cardinality() uses the efficient Schoof-Elkies-Atkin algorithm internally.

# Count points on y^2 = x^3 + x + 1 for various primes
primes_list = [7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
               71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137,
               139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199]

a_coeff, b_coeff = 1, 1
valid_primes = []
counts = []
traces_list = []

for p in primes_list:
    try:
        E_tmp = EllipticCurve(GF(p), [a_coeff, b_coeff])
        n = E_tmp.cardinality()
        valid_primes.append(p)
        counts.append(n)
        traces_list.append(p + 1 - n)
    except ArithmeticError:
        pass  # singular curve for this prime

# Plot: point count vs p+1 with Hasse band
p_points = point(list(zip(valid_primes, counts)), color='blue', size=20,
                 legend_label=r'$|E(\mathbb{F}_p)|$')
p_center = line([(valid_primes[0], valid_primes[0]+1),
                 (valid_primes[-1], valid_primes[-1]+1)],
                color='red', linestyle='dashed', thickness=1.5,
                legend_label='$p + 1$')

# Hasse band
upper = [(p, p + 1 + 2*sqrt(RR(p))) for p in valid_primes]
lower = [(p, p + 1 - 2*sqrt(RR(p))) for p in valid_primes]
p_band_u = line(upper, color='orange', linestyle='dotted', thickness=1.5,
                legend_label=r'Hasse band $p+1 \pm 2\sqrt{p}$')
p_band_l = line(lower, color='orange', linestyle='dotted', thickness=1.5)

p_plot = p_points + p_center + p_band_u + p_band_l
p_plot.axes_labels(['Prime $p$', 'Number of points'])
show(p_plot, figsize=[10, 5], title=r"Point Count $|E(\mathbb{F}_p)|$ vs Hasse Bound")

# Print table
print(f"{'p':>5s}  {'#E':>5s}  {'p+1':>5s}  {'trace':>6s}  {'2sqrt(p)':>9s}  {'Hasse OK?':>10s}")
print("-" * 48)
for p_val, c_val, t_val in zip(valid_primes, counts, traces_list):
    bound = 2 * sqrt(RR(p_val))
    ok = abs(t_val) <= bound
    print(f"{int(p_val):5d}  {int(c_val):5d}  {int(p_val+1):5d}  {int(t_val):6d}  {float(bound):9.2f}  {'YES' if ok else 'NO':>10s}")
../_images/7d892ca30e1ad7e380efcc2b0206684a2c95d53583bc38164ce11aff6092de99.png
    p     #E  p+Integer(1)   trace   2sqrt(p)   Hasse OK?
------------------------------------------------
    7      5      8       3       5.29         YES
   11     14     12      -2       6.63         YES
   13     18     14      -4       7.21         YES
   17     18     18       0       8.25         YES
   19     21     20      -1       8.72         YES
   23     28     24      -4       9.59         YES
   29     36     30      -6      10.77         YES
   37     48     38     -10      12.17         YES
   41     35     42       7      12.81         YES
   43     34     44      10      13.11         YES
   47     60     48     -12      13.71         YES
   53     58     54      -4      14.56         YES
   59     63     60      -3      15.36         YES
   61     50     62      12      15.62         YES
   67     56     68      12      16.37         YES
   71     59     72      13      16.85         YES
   73     72     74       2      17.09         YES
   79     86     80      -6      17.78         YES
   83     90     84      -6      18.22         YES
   89    100     90     -10      18.87         YES
   97     97     98       1      19.70         YES
  101    105    102      -3      20.10         YES
  103     87    104      17      20.30         YES
  107    105    108       3      20.69         YES
  109    123    110     -13      20.88         YES
  113    125    114     -11      21.26         YES
  127    126    128       2      22.54         YES
  131    128    132       4      22.89         YES
  137    126    138      12      23.41         YES
  139    126    140      14      23.58         YES
  149    136    150      14      24.41         YES
  151    154    152      -2      24.58         YES
  157    171    158     -13      25.06         YES
  163    189    164     -25      25.53         YES
  167    144    168      24      25.85         YES
  173    172    174       2      26.31         YES
  179    180    180       0      26.76         YES
  181    190    182      -8      26.91         YES
  191    217    192     -25      27.64         YES
  193    201    194      -7      27.78         YES
  197    222    198     -24      28.07         YES
  199    218    200     -18      28.21         YES

31.16 Group Structure and Point Orders#

The group \(E(\mathbb{F}_p)\) is always isomorphic to \(\mathbb{Z}/n_1\mathbb{Z} \times \mathbb{Z}/n_2\mathbb{Z}\) for some \(n_1, n_2\) with \(n_2 \mid n_1\). In many cases (especially for cryptographic curves), \(n_2 = 1\), meaning the group is cyclic. A point \(G\) whose order equals \(\#E(\mathbb{F}_p)\) is called a generator of the group.

SageMath can compute point orders directly with P.order() and the full group structure with E.abelian_group().

# Group structure analysis using SageMath
E = EllipticCurve(GF(97), [2, 3])
pts = E.points()
n_E = E.cardinality()

print(f"Curve: {E}")
print(f"#E(F_97) = {n_E}")
print(f"Group structure: {E.abelian_group()}")

# Compute order of each point
from collections import Counter

orders = {}
generators = []
for P in pts:
    ord_P = P.order()
    orders[str(P)] = ord_P
    if ord_P == n_E:
        generators.append(P)

# Distribution of orders
order_counts = Counter(orders.values())
print(f"\nOrder distribution:")
for order in sorted(order_counts.keys()):
    print(f"  Order {int(order):4d}: {int(order_counts[order]):3d} points")

print(f"\nNumber of generators (order = {n_E}): {len(generators)}")
if generators:
    print(f"Sample generators: {generators[:5]}")

# Factor the group order
print(f"\nFactorization of #E = {n_E}: {factor(n_E)}")

# Plot order distribution
order_vals = sorted(order_counts.keys())
order_freqs = [order_counts[o] for o in order_vals]
p = bar_chart(order_freqs, color='steelblue')
show(p, figsize=[8, 4],
     title=f"Distribution of Point Orders on $E(\\mathbb{{F}}_{{97}})$")
Curve: Elliptic Curve defined by y^2 = x^3 + 2*x + 3 over Finite Field of size 97
#E(F_97) = 100
Group structure: Additive abelian group isomorphic to Z/50 + Z/2 embedded in Abelian group of points on Elliptic Curve defined by y^2 = x^3 + 2*x + 3 over Finite Field of size 97

Order distribution:
  Order    1:   1 points
  Order    2:   3 points
  Order    5:   4 points
  Order   10:  12 points
  Order   25:  20 points
  Order   50:  60 points

Number of generators (order = 100): 0

Factorization of #E = 100: 2^2 * 5^2
../_images/9c681265713d7d3cb3fbcf082b5975944cabf311d2b6abf1dca608c7412d1f08.png

31.17 Singular Curves: A Cautionary Note#

When the discriminant \(\Delta = -16(4a^3 + 27b^2)\) vanishes, the curve is singular – it has a cusp or a node. Singular curves do not form a group in the usual sense, and their “discrete logarithm problem” can be solved efficiently. This is why the discriminant check is critical.

# Plot singular curves over R using SageMath's implicit_plot
var('x, y')

# Cusp: y^2 = x^3 (a=0, b=0, discriminant = 0)
p_cusp = implicit_plot(y^2 - x^3, (x, -1, 3), (y, -4, 4),
                        color='red', linewidth=2)
p_cusp += point([(0, 0)], color='black', size=60)
p_cusp += text("Cusp", (0.5, 0.8), fontsize=12, color='black')
p_cusp += text(r"$y^2 = x^3$ ($\Delta = 0$)", (1, 3.5), fontsize=11, color='red')

# Node: y^2 = (x-1)^2(x+2) has a node at x=1
p_node = implicit_plot(y^2 - (x-1)^2*(x+2), (x, -2.5, 3), (y, -4, 4),
                        color='red', linewidth=2)
p_node += point([(1, 0)], color='black', size=60)
p_node += text("Node", (1.5, 0.8), fontsize=12, color='black')
p_node += text(r"$y^2 = (x-1)^2(x+2)$ ($\Delta = 0$)", (0.5, 3.5), fontsize=11, color='red')

show(graphics_array([p_cusp, p_node], nrows=1, ncols=2), figsize=[12, 5])

# Show that SageMath correctly rejects singular curves
for a_test, b_test in [(0, 0), (-3, 2)]:
    try:
        E_bad = EllipticCurve(GF(97), [a_test, b_test])
        print(f"a={a_test}, b={b_test}: accepted (unexpected)")
    except ArithmeticError as e:
        print(f"a={a_test}, b={b_test}: correctly rejected -- {e}")
../_images/b838ad521e05f237cf3d7c83b3b4903b623d2bbf0202ec9ed8b8ea4fa70c826b.png
a=0, b=0: correctly rejected -- y^2 = x^3 defines a singular curve
a=-3, b=2: correctly rejected -- y^2 = x^3 + 94*x + 2 defines a singular curve

31.18 Comparing Curves: Point Counts and the Hasse Band#

Different curve parameters \((a, b)\) yield different group orders over the same field. Let us survey several curves over \(\mathbb{F}_{97}\) and see how their point counts distribute within the Hasse band.

# Survey point counts for various curves over F_97
p = 97
hasse_center = p + 1

print(f"p = {p}")
print(f"Hasse band: [{float(p + 1 - 2*sqrt(RR(p))):>.1f}, {float(p + 1 + 2*sqrt(RR(p))):>.1f}]")
print(f"Center: p + 1 = {hasse_center}")
print()

curve_orders = []

for a_val in range(0, 10):
    for b_val in range(1, 10):
        try:
            E_tmp = EllipticCurve(GF(p), [a_val, b_val])
            n = E_tmp.cardinality()
            curve_orders.append(n)
        except ArithmeticError:
            continue

# Plot histogram of point counts
min_o, max_o = min(curve_orders), max(curve_orders)
hist_data = [curve_orders.count(v) for v in range(min_o, max_o + 1)]
p_hist = bar_chart(hist_data, color='steelblue')
show(p_hist, figsize=[10, 5],
     title=f"Distribution of Point Counts for Various Curves over $\\mathbb{{F}}_{{{p}}}$")

print(f"Total non-singular curves tested: {len(curve_orders)}")
print(f"Min #E = {min(curve_orders)}, Max #E = {max(curve_orders)}")
print(f"All within Hasse bound: "
      f"{all(abs(n - hasse_center) <= 2*sqrt(RR(p)) for n in curve_orders)}")
p = 97
Hasse band: [78.3, 117.7]
Center: p + 1 = 98
../_images/c10e9071a7c9a7133ed982e62c131040a49394ad1575644973b6f466a812d07e.png
Total non-singular curves tested: 89
Min #E = 79, Max #E = 117
All within Hasse bound: True

31.19 Exercises#

31.20 Summary#

In this chapter we established the mathematical foundations of elliptic curve arithmetic:

  • An elliptic curve in short Weierstrass form is \(y^2 = x^3 + ax + b\) with \(4a^3 + 27b^2 \neq 0\) (non-singular).

  • The points on \(E\) over a field \(K\), together with a point at infinity \(\mathcal{O}\), form an abelian group under the chord-and-tangent addition law.

  • Over a finite field \(\mathbb{F}_p\), Hasse’s theorem constrains the group order: \(|\#E(\mathbb{F}_p) - p - 1| \leq 2\sqrt{p}\).

  • The scalar multiplication \(kP\) is computed efficiently via the double-and-add algorithm in \(O(\log k)\) group operations.

  • The Elliptic Curve Discrete Logarithm Problem (ECDLP) – given \(G\) and \(Q = kG\), find \(k\) – is believed to be computationally intractable for well-chosen curves, forming the basis of modern ECC.

SageMath’s native EllipticCurve replaces the manual Python implementation with elegant, concise operations: E = EllipticCurve(GF(p), [a, b]), P + Q, k * P, P.order(), E.cardinality(), and E.abelian_group().

Tip

Looking ahead. In subsequent chapters, we will use the group law developed here to build cryptographic protocols: ECDH (key exchange), ECDSA (digital signatures), and EdDSA (Edwards-curve signatures). We will also study specialized attacks on poorly chosen curves.

31.21 References#

  1. V. S. Miller, “Use of Elliptic Curves in Cryptography,” Advances in Cryptology – CRYPTO ‘85, LNCS 218, Springer, 1986, pp. 417–426.

  2. N. Koblitz, “Elliptic Curve Cryptosystems,” Mathematics of Computation, vol. 48, no. 177, 1987, pp. 203–209.

  3. J. H. Silverman, The Arithmetic of Elliptic Curves, 2nd edition, Graduate Texts in Mathematics 106, Springer, 2009.

  4. L. C. Washington, Elliptic Curves: Number Theory and Cryptography, 2nd edition, CRC Press, 2008.

  5. H. Hasse, “Zur Theorie der abstrakten elliptischen Funktionenkorper,” Journal fur die reine und angewandte Mathematik, 1936.

  6. D. Hankerson, A. Menezes, S. Vanstone, Guide to Elliptic Curve Cryptography, Springer, 2004.

  7. W. Stein et al., SageMath, the Sage Mathematics Software System, https://www.sagemath.org.