Stirling Numbers – a Case Study
Vector 8:2, page 44
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 44
Theme¶
Several different APL functions, for the evaluation of Stirling numbers of the second kind, are presented. Their advantages and disadvantages are discussed.
Exposition¶
A colleague doing some theoretical statistical work on the analysis of an experiment which had resulted in counts from the truncated Poisson distribution being observed, recently asked me if I had an APL function to evaluate the Stirling numbers of the second kind. I did not have such a function, but he told me they were defined for positive integers \(n\) and \(t\) by
a formula obtainable from Abramowitz and Stegun (1964). With APL it is easy to turn this expression into executable form, for example:
However this function can be improved: it ignores the fact that the first term is identically zero, it does not make it clear that the result is always positive, and APL has a more concise way of summing alternating series. A better function is therefore given by
It may or may not be numerically preferable to factorise out the \(n!\) and use
Of course APL’s power comes from its ability to handle several different values simultaneously. We could easily have a vector for T.
However, it is likely that the values of T we would require, at least to form a table of the Stirling numbers (which might be useful in the long run to save calculation) would be all the positive integers up to T. This leads to
and we could have, had we wished, developed STIR2 to this form instead of STIR3. If we wish to produce and save a table of the numbers we may build STIR5 into a loop:
It is now time to look at the table of numbers this function produces; this is shown in Table 1 for N=6 and T=6. It becomes very clear that a function which wastes its time calculating values for N>T is unsatisfactory. This, of course, should have become evident if a little more algebra had been done before turning to programming, or if STIR2, say, had been properly tested, but the functions are so simple it hardly seemed necessary to test them all individually!
Table 1
| n \ t | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 0 | 1 | 3 | 7 | 15 | 31 |
| 3 | 0 | 0 | 1 | 6 | 25 | 90 |
| 4 | 0 | 0 | 0 | 1 | 10 | 61 |
| 5 | 0 | 0 | 0 | 0 | 1 | 15 |
| 6 | 0 | 0 | 0 | 0 | 0 | 1 |
Development¶
It is important to define more coherently the problem we are tackling. We now want to produce a table of Stirling numbers in which the element at row \(n\) and column \(p\) contains \(\sigma^n_{n+p}\), and to solve the practical problem underlying the computing problem we would like \(n\) to take values up to at least 100 and \(p\) to reach at least \(2n\) for as many values of \(n\) as possible.
This means rewriting STIR5 as
Running this in a loop produces Table 2 for N=8 and P=6. The magnitude of the values appearing in the table gives rise to concern for the function: when will it overflow? This is certainly a matter which will have to be investigated, but there are other interesting aspects to the table.
Table 2
| n \ p | 1 | 2 | 3 | 4 | 5 | 6 |
|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 3 | 7 | 15 | 31 | 63 | 127 |
| 3 | 6 | 25 | 90 | 301 | 966 | 3025 |
| 4 | 10 | 65 | 350 | 1701 | 7770 | 34105 |
| 5 | 15 | 140 | 1050 | 6951 | 42525 | 246738 |
| 6 | 21 | 266 | 2646 | 22827 | 179487 | 1323652 |
| 7 | 28 | 462 | 5880 | 63987 | 627396 | 5715424 |
| 8 | 36 | 750 | 11880 | 159027 | 1899612 | 20912320 |
From row 2 it is clear that \(\sigma^2_{t+2} = 2\sigma^2_{t+1} + 1\), and from row 3 it is not too difficult to deduce that \(\sigma^3_{t+3} = 3\sigma^3_{t+2} + \sigma^2_{t+3}\), and hence to discover the recurrence relation
This enables us to write a function without calculating the factorials or binomial coefficients in our original definition.
∇ Z←N STIR8 P;I;J
[1] Z←(N,P+1)⍴1
[2] I←2
[3] NEWROW:J←2
[4] NEWVAL:Z[I;J]←Z[I-1;J]+I×Z[I;J-1]
[5] →((P+1)≥J←J+1)/NEWVAL
[6] →(N≥I←I+1)/NEWROW
[7] Z← 0 1 ↓Z
∇
But what is this? An APL function with a double loop? Surely that is unnecessary. Instead of working a row at a time we must evaluate the table by columns.
∇ Z←N STIR9 P;J
[1] Z←(N,P+1)⍴1
[2] J←2
[3] NEWCOL:Z[;J]←+\Z[;J-1]×⍳N
[4] →((P+1)≥J←J+1)/NEWCOL
[5] Z← 0 1 ↓Z
∇
Now we have a much faster function (taking about a quarter of the time STIR8 does for N=80 and P=80), which is approaching its final stage, but let us look more closely at the table and the function we have written. Column 1 is just +\⍳N and column 2 is therefore +\(⍳N)×+\⍳N. Indeed column P is ⍎1↓(7×P)⍴'×+\(⍳N)' an expression so concise it is hardly worth writing a function for. Notice that ⍎1↓(9×P)⍴'×⎕←+\(⍳N)' prints out the entire table given by STIR9, although it transposes it and prints each line to a different format. This is avoidable, for small values of N and P, by ⍎2↓((12×P)⍴'×⍎⎕←10 0⍕+\K'),'←⍳N' but we wish to store the table, not print it, a task which is done by our final version
which is slightly faster even than STIR9.
Recapitulation¶
We have seen how it is worthwhile to look at the problem from a theoretical point of view so that we may obtain different algorithms for calculation. Even with an apparently obvious formula, simply translated into APL, it may be more appropriate to use a recurrence relation. This is not only for reasons of time, which are not of vital importance if the values obtained are to be stored in any case, but for reasons of numerical accuracy. In passing we suggested that STIR2, which involves the calculation of a factorial and a set of binomial coefficients, might have different numerical properties from STIR3 which involves two sets of factorials.
To investigate this suppose N=50 and T=60, the sort of numbers which might indeed arise as the result of a small experiment. STIR1 and STIR2 give \(\sigma^{50}_{60} = 2.02\times10^{25}\) while STIR3 and STIR4 give \(5.06\times10^{25}\). The sensitivity of each function to the order of summing the series can be looked at by changing -/ in STIR2 and STIR3 to -/⌽; the results change to \(1.71\times10^{26}\) and \(7.09\times10^{25}\) respectively. With STIR4 we can even alter the calculated value depending on which values other than 60 we include in vector T. STIR (and STIR9 which uses the same arithmetic) gives \(9.55\times10^{24}\) for \(\sigma^{50}_{60}\), which I believe to be correct. Not only is the algorithm contained in ⍎1↓(7×P)⍴'×+\(⍳N)' remarkably succinct, it is numerically preferable to the others. STIR works whenever 350>N+2×P and occasionally outside that range. When the function fails it is sometimes because the workspace becomes full, and sometimes because a domain error occurs.
The figure shows an isometric plot of the common logarithm of the Stirling numbers of the second kind, for values of \(n\) and \(t-n\) up to 80. The values given above, except for that of \(9.55\times10^{24}\), are system-dependent; the ones quoted were obtained from IBM APL on a Personal Computer, the same phenomenon occurs with different values using TRYAPL2. Other systems give a domain error instead of inaccurate results.
Figure: the Stirling Numbers of the Second Kind
Coda¶
The style of programming illustrated in STIR can be utilised in many instances of recursion, though the Stirling numbers of the second kind may be the most elegant. Readers might like to deduce what the following expressions evaluate.
⍎'+',((9×N-2)⍴'Z←Z,+/¯2↑'),'Z←1 1'
⍎1↓((16×P)⍴'×+\(-(⍳N)+K←K+1)'),'+0×K←¯1'
⍎'+',((31×P)⍴'Z←Z+N⍴+⍀((⌈N÷K),K)⍴((K←K+1)⍴0),'),'Z←N⍴K←1'
Reference¶
[1] Abramowitz, M. & Stegun, IA (eds). Handbook of mathematical functions with formulas, graphs, and mathematical tables. Washington: US Govt. Print. Off. 1964.