Given the natural numbers, we can do basic arithmetic. For example, we can add up all the natural numbers less than n to get
add that to a second copy of itself, but with the terms in reverse order, and divide by two; this may seem complicated, but it clearly gets us the same answer:
now interleave the first half with the second:
but now each pair of terms adds one drawn from the increasing sequence 0, 1, … to one drawn from the decreasing sequence, n−1, n−2, …, hence each pair adds up simply to n−1; and there are n such pairs, so the total comes to:
and, since either n is even (so n/2 is whole) or n is odd (so (n−1)/2 is whole), this always is a whole number, despite dividing by 2. We can confirm this result inductively. First notice that sum(0) = sum({}) has no values to add together, so it's 0, which is indeed 0.(0−1)/2, in so far as we can make any sense of 0−1 (which isn't a natural thing to do to the natural numbers – but we can do the arithmetic in the integers, which subsume the naturals); in case that leaves anyone feeling nervous, we can separately look at sum(1) = sum({0}) = 0 because that's the single term being summed; and 1.(1−1)/2 is indeed 0. Now suppose we have some n (e.g. 0 or 1) for which sum(n) is n.(n−1)/2; then sum(1+n) = n+sum(n) = n +n.(n−1)/2 = n.(2 +n−1)/2 = n.(n+1)/2 = (n+1).((n+1)−1)/2, so our formula works for n+1 whenever it works for n; but, given that it works for 0 (and 1) this implies, inductively, that it works for any natural n.
Well, that was nice and easy; so how about summing the squares of the
naturals less than n, for given n ? Well, I can't think of a trick like
the duplicate in reverse, add pairwise and halve
trick that got us the
right formula for the sum of the naturals themselves, so let's see if we can
come up with an inductive approach. In fact, it turns out there's a slightly
more complex-looking problem that, with a little ingenuity, is easy to solve,
and from which we can infer the sum of squares; best of all, it generalizes
nicely to give us the sums of any natural power. (As before for 0−1,
we'll have some negative integers show up; you can chose to do arithmetic in the
integers to resolve these, or you can note that they always arise in
terms that have at least one factor of 0, meaning that we could have shuffled
the sets over which we're summing to simply omit the terms that have a 0 factor,
hence also all the terms that involve an unnatural term in them.) So, instead
of sum(: i.i ←i :n), let's look at:
because the extra factor is just a long-winded way of writing 3/3, but we can now re-arrange this to:
Now, the j = 0 term in the second sum is zero, so we can discard it; thus each j for which the second sum contributes a non-zero value is 1+i for some i – and look what happens if I write j = 1+i:
the two terms are summing the same values ! Of course, the second is summing over a different set of values for i, but we can at least cancel all the terms from values of i that show up for both sums; that is, every i in n for which 1+i is also in n; which is, in fact, every i in n except i = n−1, so this difference of sums just reduces to:
and we can get sum(: i.i ←i :n) as sum(: i.(i−1) +i ←i :n) = n.(n−1).(n−2)/3 +n.(n−1)/2 = n.(n−1).(2.n−1)/6. Note that, once again, we divide n.(n−1).(n−2) by 3 but necessarily get a whole number because n, n−1 and n−2 are three consecutive integers, so one of them is a multiple of three.
Well, I said we could generalize that to get the sum of any power; you might have guessed already that I'll do it by actually summing products of successive integers, instead of powers; once we know the sums for powers up to k, we can compute the sum of products of k+1 terms, i.(i−1)…(i−k), and deduce the sum of (k+1)-th powers from this, using our earlier knowledge of the sums of powers up to k. Now, because we're dealing with products of consecutive naturals, it actually simplifies matters to introduce a function (the number of subsets having m members in a set with n members; but I'll here use the transpose of the usual chose(n, m) = C(m, n)) related to Pascal's triangle:
which lets us re-write the previous result as sum(: C(2) :n) = C(3, n) and (since C(1) is simply the identity on {naturals}) the initial sum(n) result as sum(: C(1) :n) = C(2, n); indeed I'll now show that, in general, sum(: C(k) :n) = C(1+k, n). So, let's see:
and, as before, C(k+1, 0) is 0, so we can ignore the i=0 term in the second sum and replace it with sum(: C(k+1, j+1) ←j+1 :n), which cancels every term in the first sum that has i+1 in n, leaving only the i = n−1 term:
As promised, we can now use this to obtain the sums of powers of naturals:
which gets rather fiddly to work out, so naturally I wrote some code to do the job for me; here are the next few (asuming I've finally fixed all the bugs in the code, of course):
Using these, we can of course also sum any polynomial over any arithmetic sequence (that is, any sequence in which the difference between successive entries is constant). Note that (although the way I've displayed the sum polynomials above doesn't make it obvious) the leading order term in sum(: power(k) :n) is power(k+1, n)/(k+1); this can fairly easilly be derived, independently, by observing that the difference between its values at n and 1+n is just power(n, k), the single term added to the sum.
We can use the results above to obtain the mean and variance for uniform discrete random variates with values ranging from 1 to n (i.e. fair n-sided die-rolls), as functions of natural n. The mean is sum(1+n)/n = (1+n)/2 and the expected value of the square is sum(: i.i ←i :1+n)/n = (n+1).(2.n+1)/6, giving a variance of
yielding, for the dice commonly used in rôle-playing games:
| die | d4 | d6 | d8 | d10 | d12 | d20 |
|---|---|---|---|---|---|---|
| mean | 2.5 | 3.5 | 4.5 | 5.5 | 6.5 | 10.5 |
| variance | 5/4 | 35/12 | 21/4 | 33/4 | 143/12 | 133/4 |
| standard deviation | 1.12 | 1.71 | 2.29 | 2.87 | 3.45 | 5.77 |
Note that, when you add several die rolls, the sum's mean and
variance are the sums of the means and variances of the individual rolls added;
the standard deviation of the sum is then the square root of its variance and
the typical spread is one or two standard deviations either side of the
mean. Thus 3d6 has mean 3×3.5 = 10.5, variance 3×35/12 = 35/4 =
8.75 and standard deviation 2.96, so nearly always falls in the range from 5 to
16 (for 3 and 18 the probabilities are 1/216; for 4 and 17 they are 1/72; for a
total of 8/216 = 1/27; so nearly always
, here, is over 96% of the time)
and usually falls in the range from 8 to 13. For 3d4+3, we get mean 10.5 and
variance 3.75, hence standard deviation 1.94, so results are usually in the
range 9 through 12 and nearly always in the range 7 through 14 (with probability
1/64 for each of 6 and 15, so nearly always
is almost 97% in this
case).
Written by Eddy.