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-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 < 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 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
No comments:
Post a Comment