diff --git a/examples/benchmarks/benchmark_suite.py b/examples/benchmarks/benchmark_suite.py
index 96762466..39390f54 100644
--- a/examples/benchmarks/benchmark_suite.py
+++ b/examples/benchmarks/benchmark_suite.py
@@ -20,6 +20,7 @@
 from examples.benchmarks.mario_level import MarioBenchmark
 from examples.benchmarks.lambda_calculus import LambdaCalculusBenchmark
 from examples.benchmarks.median import MedianBenchmark
+from examples.benchmarks.pell import PellsEquationBenchmark
 
 
 banknote = get_banknote()
@@ -40,6 +41,7 @@
     MarioBenchmark(),
     LambdaCalculusBenchmark(),
     MedianBenchmark(),
+    PellsEquationBenchmark(),
 ]
 
 
diff --git a/examples/benchmarks/pell.py b/examples/benchmarks/pell.py
new file mode 100644
index 00000000..bcea468e
--- /dev/null
+++ b/examples/benchmarks/pell.py
@@ -0,0 +1,57 @@
+"""A benchmark for the positive solutions of Pell's equation.
+
+The default instance is ``x² - 2y² = 1``.  Pell solutions are conventionally
+positive; allowing arbitrary signed integers would make minimizing ``x + y``
+unbounded below because the negation of every solution is also a solution.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Annotated
+
+from examples.benchmarks.benchmark import Benchmark, example_run
+from geneticengine.grammar.grammar import Grammar, extract_grammar
+from geneticengine.grammar.metahandlers.ints import IntRange
+from geneticengine.problems import Problem, SingleObjectiveProblem
+
+
+@dataclass
+class PellPair:
+    x: Annotated[int, IntRange(1, 100)]
+    y: Annotated[int, IntRange(1, 100)]
+
+    def evaluate(self, d: int) -> tuple[int, int]:
+        return self.x, self.y
+
+
+def pell_fitness(candidate: PellPair, d: int = 2) -> float:
+    """Penalize equation violations, then minimize the positive coordinate sum."""
+    x, y = candidate.evaluate(d)
+    residual = x * x - d * y * y - 1
+    return abs(residual) * 1_000_000 + x + y
+
+
+class PellsEquationBenchmark(Benchmark):
+    """Find the smallest positive solution to ``x² - 2y² = 1``."""
+
+    def __init__(self, d: int = 2) -> None:
+        if d <= 0 or int(d**0.5) ** 2 == d:
+            raise ValueError("d must be positive and nonsquare")
+        self.d = d
+        self.problem = SingleObjectiveProblem(
+            minimize=True,
+            target=5 if d == 2 else None,
+            fitness_function=lambda candidate: pell_fitness(candidate, d),
+        )
+        self.grammar = extract_grammar([PellPair], PellPair)
+
+    def get_problem(self) -> Problem:
+        return self.problem
+
+    def get_grammar(self) -> Grammar:
+        return self.grammar
+
+
+if __name__ == "__main__":
+    example_run(PellsEquationBenchmark())
diff --git a/examples/mathematical_programming/pells_equation.py b/examples/mathematical_programming/pells_equation.py
new file mode 100644
index 00000000..ab6ece7a
--- /dev/null
+++ b/examples/mathematical_programming/pells_equation.py
@@ -0,0 +1,8 @@
+"""Genetic programming example for Pell's equation."""
+
+from examples.benchmarks.benchmark import example_run
+from examples.benchmarks.pell import PellsEquationBenchmark
+
+
+if __name__ == "__main__":
+    example_run(PellsEquationBenchmark())
diff --git a/run_examples.sh b/run_examples.sh
index 6bcaf72c..cb512d4b 100755
--- a/run_examples.sh
+++ b/run_examples.sh
@@ -70,6 +70,7 @@ run_example examples/mathematical_programming/fixed_egyptian_fractions.py
 run_example examples/mathematical_programming/rational_egyptian_fractions.py
 run_example examples/mathematical_programming/reptend_fraction.py
 run_example examples/mathematical_programming/incremented_reciprocal.py
+run_example examples/mathematical_programming/pells_equation.py
 
 
 run_example examples/progsys/Number_IO.py
diff --git a/tests/benchmarks/test_mathematical.py b/tests/benchmarks/test_mathematical.py
index 3db5afdd..ff15f648 100644
--- a/tests/benchmarks/test_mathematical.py
+++ b/tests/benchmarks/test_mathematical.py
@@ -4,6 +4,9 @@
 from examples.benchmarks.mathematical import IncrementedReciprocalBenchmark
 from examples.benchmarks.mathematical import RationalEgyptianFractionBenchmark
 from examples.benchmarks.mathematical import ReptendFractionBenchmark
+from examples.benchmarks.pell import PellsEquationBenchmark
+from examples.benchmarks.pell import PellPair
+from examples.benchmarks.pell import pell_fitness
 from examples.benchmarks.mathematical import fixed_egyptian_fraction
 from examples.benchmarks.mathematical import incremented_reciprocal
 from examples.benchmarks.mathematical import reptend_fraction
@@ -42,5 +45,12 @@ def test_all_benchmark_grammars_are_constructible():
         RationalEgyptianFractionBenchmark(),
         ReptendFractionBenchmark(),
         IncrementedReciprocalBenchmark(),
+        PellsEquationBenchmark(),
     ):
         assert benchmark.get_grammar() is not None
+
+
+def test_pells_equation_minimizes_the_smallest_positive_solution():
+    assert pell_fitness(PellPair(3, 2)) == 5
+    assert pell_fitness(PellPair(1, 1)) > pell_fitness(PellPair(3, 2))
+    assert PellsEquationBenchmark().get_problem().target == [5]