## 2008-04-11

### Fixnum Figurations

The polytopic numbers P.k.n are a kind of figurate numbers giving the number of lattice points on a k-dimensional simplex (a right triangle, tetrahedron, …) x0 + x1 + … + xk-1n of integer side n. For k = 2 we have the Pythagorean triangular numbers 1, 3, 6, 10…, the last one called the tetractys. In general, the way to calculate P.k.n is to stack n simplices of dimension k - 1 with decreasing sides:

```P.0.n = 1
P.k.n = ⟨∑ i : 0 ≤ i < n : P.(k-1).i⟩
```

This can be put in closed form by taking antidifferences:

```P.k.n = n↑k / k!
```

where nk = n⋅(n + 1)⋅…⋅(n + k - 1) is the rising factorial power. For positive, integer k, nk = (n + k - 1)↓k, which is the falling factorial power, defined analogously. Since by definition the binomial coefficient C.n.k = nk / k!, we have that:

```P.k.n = C.(n + k - 1).k
```

and the problem of calculating polytopic number can be reduced to the calculation of binomial coefficients. Mark-Jason Dominus has posted a while ago about the best method to compute binomial coefficients of machine-size precision without integer overflow. I'll weave his insight in the formal derivation of a program satisfying these requirements.

What is the largest binomial coefficient we can compute? Let B be the bit width of a positive machine integer (in OCaml, B = 30). For a given n, the largest binomial coefficient is C.n.⌊k/2⌋, so that:

```
C.n.k < 2B
⇐
C.(2⋅k).k < 2B
≡
(2⋅k)↓k / k! < 2B
≡
⟨∏ i : 0 ≤ i < k : 2⋅k - i⟩/⟨∏ i : 1 ≤ i ≤ k : i⟩ < 2B
≡ { Changing the index of summation }
⟨∏ i : 1 ≤ i ≤ k : k + i⟩/⟨∏ i : 1 ≤ i ≤ k : i⟩ < 2B
≡ { Algebra }
⟨∏ i : 1 ≤ i ≤ k : (k + i)/i⟩ < 2B
⇐ { Taking logarithms }
⟨∑ i : 1 ≤ i ≤ k : lg.((k + i)/i)⟩ < B
```

Define S.k = ⟨∑ i : 1 ≤ ik : lg.((k + i)/i)⟩. Then S.0 = 0, and S.(k + 1) - S.k = 1 + lg.(2⋅k + 1) - lg.(k + 1) (the calculations are tedious but straightforward.) Hence, the following program is immediate:

```let maxcomb b =
let rec iter s k =
let d = log (float (4*k + 2) /. float (k + 1)) in
if s < d then k else iter (s -. d) (k + 1)
in iter (log 2. *. float b) 0
```

With this, we can determine that for B = 30, k must be at most 16. This is a bit too strong, since C.33.15 < 230 < C.33.16. I'm now set to derive the program to compute binomial coefficients. But first, I'll need this later (foreshadowing!):

```let rem m n =
let r = m mod n in
if r < 0 then r + abs n else r

let rec gcd m n =
if n = 0 then m else gcd n (rem m n)
```

(cf Daan Leijen, Division and Modulus for Computer Scientists.) The program must fulfill the following specification:

```let binomial n k =
(* 0 ≤ k ≤ n ≤ 32 *)
let r = ref ... in
(* Binomial *)
!r
```

However, with the stated conditions, C.n.k = C.n.(n - k), and I can take advantage of that to minimize work by forcing kn - k:

```let rec binomial n k =
if k > n - k then binomial n (n - k) else
(* 0 ≤ k ≤ n - k < n ≤ 32 *)
let r = ref ... in
(* Binomial *)
!r
```

From the definition of binomial coefficient:

```P ≡ r = ⟨∏ i : 0 ≤ i < k : n - i⟩/⟨∏ i : 0 ≤ i < k : i + 1⟩
```

Changing the constant k by a variable h:

```P.h ≡ r = ⟨∏ i : 0 ≤ i < h : n - i⟩/⟨∏  i : 0 ≤ i < h : i + 1⟩
```

with P.k ≡ P. The initialization:

```  let r, h = ref 1, ref 0 in ...
```

establishes P.0. We have:

```let rec binomial n k =
if k > n - k then binomial n (n - k) else
(* 0 ≤ k ≤ n - k < n ≤ 32 *)
let r, h = ref 1, ref 0 in
(* P.h *)
while !h != k do
(* increment h under invariance of P.h *)
done;
(* P.k *)
!r
```

Now:

```   P.(h+1)
≡
r = ⟨∏ i : 0 ≤ i < h+1 : n - i⟩/⟨∏ i : 0 ≤ i < h+1 : i + 1⟩
≡ { Splitting the range }
r = (⟨∏ i : 0 ≤ i < h : n - i⟩⋅(n - h)) / (⟨∏ i : 0 ≤ i < h : i + 1>⋅(h+1))
```

which is restored with:

```  r := !r * (n - !h) / (!h + 1);
h := !h + 1
```

The assignment to r might unfortunately cause overflow. But here's the trick: if (n-h) / (h+1) = a/b in lowest terms, then r⋅(n-h) / (h+1) = (r/g)⋅a / (b/g), where g = (r, b). This is the most simplifying we can do while computing the binomial coefficient stepwise. Hence:

```let rec binomial n k =
if k > n - k then binomial n (n - k) else
(* 0 ≤ k ≤ n - k < n ≤ 32 *)
let r, h = ref 1, ref 0 in
(* P.h *)
while !h != k do
let f = gcd (n - !h) (!h + 1) in
let a, b = (n - !h) / f, (!h + 1) / f in
let g = gcd !r b in
r := (!r / g) * a / (b / g);
incr h
(* P.(h+1) *)
done;
(* P.k *)
!r
```

is the desired program. With this, computing polytopic numbers is just a matter of defining:

```let polytopic k n =
if n = 0 then 0 else
binomial (n + k - 1) k
```