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, …) x₀ + x₁ + … + 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³⁰ < 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
    (* 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