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, …) `x _{0}` +

`x`+ … +

_{1}`x`≤

_{k-1}`n`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 `n`↑`k` = `n`⋅(`n` + 1)⋅…⋅(`n` + `k` - 1) is the rising factorial power. For positive, integer `k`, `n`↑`k` = (`n` + `k` - 1)↓`k`, which is the falling factorial power, defined analogously. Since by definition the binomial coefficient C.`n`.`k` = `n`↓`k` / `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 ≤ `i` ≤ `k` : 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 < 2^{30} < 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 `k` ≤ `n` - `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 (* incrementhunder 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