diff --git a/optmal-test-order-math/index.html b/optmal-test-order-math/index.html
index ed4356618cb685f4b2bc0ef060265bf2e5566a2e..962ca5f25a795d0552662a6db1006721fe5b23c1 100644
--- a/optmal-test-order-math/index.html
+++ b/optmal-test-order-math/index.html
@@ -61,35 +61,58 @@ The quantities in [square brackets] in the last line can be precomputed for all
 
 <p><strong>Algorithm 1</strong> <code>best_order</code> should compute the best order to the tests, using the assumptions above.
 
-<div><pre><code class="language-python">import numpy.typing as npt
+<div><pre><code class="language-python">from typing import List
+import numpy.typing as npt
 import numpy as np
 
 def reverse(array: npt.NDArray[np.float64]) -&gt; npt.NDArray[np.float64]:
     """Reverse a numpy array"""
     return array[::-1]
 
-def cumprod_noninclusive(array: npt.NDArray[np.float64]) -&gt; npt.NDArray[np.float64]:
-    ret = np.ones_like(array)
-    for i in range(1, len(array)):
-        ret[i] = ret[i-1] * array[i-1]
+def cumsum(
+        array: npt.NDArray[np.float64],
+        inclusive: bool = True,
+        one_larger_output: bool = False,
+) -&gt; npt.NDArray[np.float64]:
+    if one_larger_output:
+        assert not inclusive
+    ret = np.zeros(len(array) + int(one_larger_output))
+    for i in range(int(not inclusive), len(ret)):
+        ret[i] = ret[i-1] + array[i - int(not inclusive)]
+    return ret
+
+def cumprod(
+        array: npt.NDArray[np.float64],
+        inclusive: bool = True,
+        one_larger_output: bool = False,
+) -&gt; npt.NDArray[np.float64]:
+    if one_larger_output:
+        assert not inclusive
+    ret = np.ones(len(array) + int(one_larger_output))
+    for i in range(int(not inclusive), len(ret)):
+        ret[i] = ret[i-1] * array[i - int(not inclusive)]
     return ret
 
-def next_best(nt: int, d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]):
+def next_best(d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]) -&gt; List[int]:
     """
     Returns the index of the best test to run.
-
-    nt, d, and p should already include a fictitious initial and fictitious final test.
     """
 
-    # Compute D and P for the unmodified order.
-    # np.cumsum is inclusive, D[k] = d[0] + d[1] + ... + d[k]
-    # cumprod_noninclusive is not, P[k] = p[0] * p[1] * ... * p[k-1] * (1 - p[k])
-    D = np.cumsum(d)
-    P = cumprod_noninclusive(p) * (1 - p)
-
-    # Precompute intermediate values.
-    A = np.cumsum(D * P) - d[0] * p[0]
-    B = np.cumsum(P) - p[0]
+    if len(d) == 1:
+        return 0
+
+    assert len(d) == len(p)
+
+    # D[ell] := d[0] + d[1] + ... + d[ell]
+    D = cumsum(d)
+    # P[ell] := p[0] * p[1] * ... * p[ell-1] * (1 - p[ell])
+    P = cumprod(p, inclusive=False) * (1 - p)
+    # A[ell] := D[1] * P[1] + D[2] * P[2] + ... + D[ell - 1] * P[ell - 1] for ell &gt; 0
+    A = cumsum(D * P, inclusive=False) - D[0] * P[0]
+    # B[ell] := P[1] + P[2] + ... + P[ell - 1] for ell &gt; 0
+    B = cumsum(P, inclusive=False) - P[0]
+    # C[ell] := D[ell+1] * P[ell+1] + D[ell+2] * P[ell+2] + ... + D[len(D) - 1] * P[len(D) - 1]
+    # TODO: this
     C = np.cumsum(reverse(D * P))
 
     # Compute E of the swapped order 0 &lt;-&gt; k for all k.
@@ -102,42 +125,33 @@ def next_best(nt: int, d: npt.NDArray[np.float64], p: npt.NDArray[np.float64]):
         + C
     )
 
-    # Find min, excluding the trailing fictitious test, which would otherwise be optimal.
-    return np.argmin(E_of_swap[:-1])
+    return np.argmin(E_of_swap)
 
 def best_order(
-        nt: int,
         d: npt.NDArray[np.float64],
         p: npt.NDArray[np.float64],
-) -&gt; npt.NDArray[np.float64]:
+) -&gt; npt.NDArray[int]:
     """
     Returns a list of indices which consist the best order to run.
-
-    nt, d, and p should not contain fictitious tests.
     """
 
-    order = np.zeros(nt)
+    assert len(d) == len(p)
 
-    # Add the fictitious test to these.
-    nt = nt + 2
-    d = np.concatenate([0], d, [0])
-    p = np.concatenate([1], p, [0])
+    tests = np.arange(len(d), dtype=int)
 
     # Loop through each test and find the best order.
-    for k in range(nt - 1):
+    for k in range(len(d)):
         # The first k have already been assigned.
-        best = next_best(nt - k, d[k:], p[k:])
-        order[k] = best
+        best = next_best(d[k:], p[k:])
 
-        # Swap best with k + 1, the first real test in the next round.
-        np.swap(p, k + 1, best)
-        np.swap(d, k + 1, best)
+        # Swap best with k, so that arr[k+1:] are the remaining tests
+        if k != best:
+            p[[k, best]] = p[[best, k]]
+            d[[k, best]] = d[[best, k]]
+            tests[[k, best]] = tests[[best, k]]
 
-        # Set the first test in the next round to be fictitious.
-        # This deletes the test we select for this round.
-        d[k + 1] = 0
-        p[k + 1] = 1
+    assert list(sorted(tests)) == list(range(len(d)))
 
-    return order
+    return tests
 </code></pre></div></p>
 </body></html>
\ No newline at end of file
diff --git a/optmal-test-order-math/template.html b/optmal-test-order-math/template.html
index e62a6f39f259a97acb8e711aaa4967db9fead562..45d4ee6226763d586276ec3bcbca8f768c67e483 100644
--- a/optmal-test-order-math/template.html
+++ b/optmal-test-order-math/template.html
@@ -1,15 +1,19 @@
+â—Šbegin[
+  (require
+   txexpr
+   "lib/txexpr.rkt")
+  (define-values (rest template-data) (splitf-txexpr* doc 'template-data))
+  (define title (get-child template-data 'title #:default '("Working Title")))
+  (define head (get-child template-data 'head #:default '()))
+  (define body (get-elements rest))
+]
 <!DOCTYPE html>
-<html lang="en">
-<head>
-  <meta charset="utf-8" />
-  <script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
-  <script
-	id="MathJax-script"
-	async
-	src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"
-	></script>
-</head>
-<body>
-  â—Š(->html doc #:splice? #t)
-</body>
-</html>
+â—Š(->html
+  `(html ((lang "en"))
+    (head
+      (meta ((charset "utf-8")))
+      (title ,@title)
+      ,@head)
+    (body
+     (h1 ,@title)
+     ,@body)))