It’s my daughter’s birthday on October 25, which as I write, is today.
Happy birthday Heather!
I was idly thinking about my family’s birthdays while having a shower this morning (as you do), which either fall on the 25th (Heather) or the 26th (Karen, Nic and I) of some month.
I started thinking about the difference between the months of each birthday.
The difference between our son Nic’s birthday in February and my wife Karen’s in June is 4 months (6-2). So too between Karen’s and Heather’s (10-6).
The number of months between Karen’s and my birthday in November is 5 (11-6).
Completing the circle, between my birthday and Nic’s is 3 months, when moving from November to February (1 month each from Nov to Dec, Dec to Jan, Jan to Feb).
Curiously, within the same year, November back to February is 9 months (11-2), which is 32, 9 divided by 3 is 3, so 3 x 3 = 9, and so on. More simply, 12 -9 is 3. Three shows up in various ways.
The numbers 3, 4, 5 reminded me of the so-called “3, 4, 5” right-angled triangle, where two of the sides have length 3 and 4 with the diagonal (hypotenuse) having length 5.
My last post ended with the strong hint that I may not have been finished with this topic.
In part 4, triangular numbers were used as a first simplification to computing the total number of gifts for the 12 days of Christmas by eliminating some additions.
It turns out that the polynomial from part 5 is the equation for the class of tetrahedral numbers.
I wanted to determine the optimal solution from first principles and experimentation. The triangular numbers link was found through serendipity (listening to an unrelated podcast). I learned about tetrahedral numbers through further reading after finishing part 7.
In a Wikipedia article about the 12 days of Christmas song, Donald Knuth’s analysis of computational complexity is mentioned and then as an afterthought:
Incidentally, it is also observed that the total number of gifts after m days equals m3 / 6 + m2 / 2 + m / 3.
which is the polynomial we saw from part 5.
One thing led to another and I came across another Wikipedia page about tetrahedral numbers that are computed with this function. The name of the class of numbers arises from the fact that such numbers can be modelled by stacking spherical balls into a tetrahedron, as illustrated with this nice 3D graphic:
The number of balls in each layer is a triangular number. The number of balls in the tetrahedron is the nth tetrahedral number, and also the number of gifts for the nth day of Christmas.
gifts-for-day-of-xmasn = 1 ball in layer1 + 3 balls in layer2 + 6 balls in layer3 + … + n(n+1)/2 balls in layern
There is also a relationship to the problem of finding the optimal way to stack canon balls worked on by Johannes Kepler (one of my heroes, who realised that planetary orbits are elliptical, not circular), but that’s a different post.
Here are some other interesting things I found along the way:
The year that most of the thought went into my posts was 2024, which turns out to be the 22ndtetrahedral number.
Ian Phillipps won an Obfuscated C competition in 1988 with a C program that generates the text of the 12 days of Christmas song. The code was described by the judges as “…like what you would get by pounding on the keys of an old typewriter at random”, and as the same Wikipedia article goes on to say: “…takes advantage of the recursive structure of the song to print its lyrics with code that is shorter than the lyrics themselves.”
There’s a proof by induction for the tetrahedral numbers for all n that uses the equation for triangular numbers in the inductive step. This deserves further discussion, and perhaps I will return to this from another viewpoint in future.
The nthtetrahedral number is represented by the 3rdrising factorial of n divided by the factorial of 3. This yields a slightly more compact equation than the polynomial, but is just another way of expressing it, so is no better computationally. The equality-related expressions below partially recapitulate the explorations throughout this series of posts, with the rising factorial form at right (source: Wikipedia).
The 4th diagonal of Pascal’s triangle (built from binomial coefficients) contains the tetrahedral numbers (from left or right). The triangular numbers are in the 3rd diagonal, and the counting numbers are in the 2nd diagonal.
In this final post, I explore a curious (at least from my perspective), yet tenuous, relationship between triangular numbers and the polynomial model, from parts 4 and 5 respectively, via a brief detour through Calculus.
As a reminder, the function to compute the Nth triangular number from part 4 is shown below along with the sum of the applications of this function to the values (day number) 1 to 12:
The function to compute the number of gifts for the Nth day of Christmas from part 5 is:
and as we saw in the last post, applying this function to 12 also gives 364:
If we take the indefinite integral of the triangular numbers function, the result is as follows:
The first term of the resulting integral is the same as the polynomial.
The coefficients of the second terms differ (one half vs one quarter) between the polynomial model and the integral of the triangular numbers formula.
Ignoring the constant of integration, C (see below), this gives:
which is 40 short of what all our methods give for the number of gifts for the 12 days of Christmas. Since we are interested in the range of days from 1 to 12, the result of applying this formula to 1 should be subtracted, giving f(12) – f(1) = 324 – 0.41666 = 323.58333.
A more intuitive way to understand and compute the integral of this function is via a numerical method called the Riemann Sum, in which rectangles are used to approximate the area under the curve (which is what the integral means here):
As the number of rectangles increases and their width decreases, the Riemann Sum result becomes closer to the Actual Sum. Making a = 0 instead adds nothing since f(0) here gives 0.
Geogebra was used to create the Riemann Sum graphic above in which the actual area computed by the so-called definite integral is also shown. The definiteintegral, F, for the 12 days of Christmas does allow us to ignore C, owing to the fact that subtracting the application F to 1 from the application of F to 12 cancels C out :
As you can see, this result is the same as the Actual Area shown in the Geogebra graphic above.
So, my intuition that integrating the triangular numbers equation would give the same result as the days-of-xmas equation, and the polynomial equation in particular, turned out to be wrong. But, there’s a simple reason for that. The derivative (opposite of the integral) of the polynomial equation is not the triangular numbers equation! It is this instead:
If we apply the Riemann Sum to this, we should get the expected result:
Note that the final term (one third) does make a contribution here since f(0) is non-zero.
So ends the long strange trip down the rabbit hole this series has taken me on.
Or does it?
I’ve spent far too much time thinking about this, but it’s been a lot of fun!
Or have I?
Definitely mathematics of the recreational kind anyway!
This post summarises the results of earlier posts.
For each method, the following are given:
A worked example for n=12 days, the only one that matters in practice for The Twelve Days of Christmas as opposed to the general problem of The N Days of Christmas I’m interested in.
The mathematical notation generalising the method, sometimes shown together with the worked example.
A function in the Julia programming language showing:
how to compute the N Days of Christmas using the method;
the time taken for n=100,000 and n=1000,000,000 on my M2 Mac OS X laptop, i.e. the time taken for the function calls daysofxmas(100000) and daysofxmas(1000000000) to complete.
Where multiple functions were presented for a method in earlier posts, the most efficient and representative (closest to the mathematical notation) is given.
Subtle considerations regarding function parameter and variable types that allow larger values of n to yield correct results are omitted from the code here; assume very large (128 bit) integer types.
Headings link back to the post in which the method was first introduced.
function daysofxmas(n)
total = 0
for day in 1:n
daysum = sum(1:day)
total += daysum
end
total
end
Run-time for n=100,000: 0.0005 seconds
Run-time for n=1000,000,000: 8.5 seconds
Notes
The mathematical notation is given with application to 12 days. 12 can be replaced by n for generality.
The run-time for this method can be significantly higher. In particular, if two nested loops are used instead of a single outer loop and the efficient sum() call shown here, the time taken can be closer to 20 seconds instead of a fraction of a second. See the Row Addition post for more details.
function daysofxmas(n)
days = [1:n;]
sum(days .* reverse(days))
end
Run-time for n=100,000: 0.0007 seconds
Run-time for n=1000,000,000: 33.7 seconds
Notes
The mathematical notation is given with application to 12 days.
There’s a lot more going on in the function code than meets the eye. Reversing a large list (of 100,000 elements in this case) takes time, multiplication is time consuming, and so on. There are faster ways to do this, but few quite so concise.
The simplification shown in the worked example of doubling half the sequence of additions is not used in the function here, because as noted in the Dot Product post, it does not help the run-time and complicates the code.
Worked example for 12 days and mathematical notation
Julia function
function triangular(n)
div(n*(n+1), 2)
end
function daysofxmas(n)
total = 0
for n in 1:n
total += triangular(n)
end
total
end
Run-time for n=100,000: 0.0003 seconds
Run-time for n=1000,000,000: 6.8 seconds
Notes
See the Triangular Numbers post for the background of how we arrive at triangular numbers. Spoiler alert, it involves triangles. 🙂 A clue is that half of the width times the height of a rectangle is the area of a triangle, which if you peer closely at the triangular function above, you’ll see.
Worked example for 12 days and mathematical notation
Julia function
function daysofxmas(n)
div(n^3, 6) + div(n^2, 2) + div(n, 3)
end
Run-time for n=100,000: 0.00001 seconds
Run-time for n=1000,000,000: 0.00001 seconds
Notes
See Polynomial Model re: how the polynomial expression is arrived at.
This is the fastest of all the approaches, for any value of n.
Summary
Taking the results for n=100,000 and n=1,000,000,000 into account, we have the following ranking:
Ranking
Method
1
Polynomial Model
2
Triangular Numbers
3
Row Addition
4
Dot Product
Ranking of methods and their functions
As n gets larger, the run-time difference between the methods due to these time complexity variations becomes more obvious, as can be seen from the results.
The polynomial model method computes the N Days of Christmas in constant time vs linear or quadratic time for the other methods.
It satisfies my initial goal of a simple expression with the lowest computational time complexity and shortest run-time, for any value of n. It is more than a million times faster than its nearest competitor for n = 1 billion.
If I was looking for the fastest possible function implementations, I’d use a language like C, but Julia is better than some other scientific computing languages in this regard.
“Simplest” or “most compact” may be in the eye of the beholder when it comes to mathematical notation (they’re all fairly compact) but shortest function run-time is harder to argue with.
The final post, part 7, briefly explores a possible link between triangular numbers and the polynomial model.
Sometimes, you need to visualise data to understand it.
Here is a plot of the number of days vs the total number of gifts, ending with our favourite, traditional data point: 12 days vs 364 total gifts.
Drawing a straight line through these data points isn’t possible, so the relationship between days and total gifts is not linear. There’s clear upward curvature.
What we need to do is create a model of the data.
A first reasonable thought is that the data may represent exponential growth, something we’ve all become accustomed to in the age of COVID. That turns out not to be the case here, as a little reflection will show; there is not a pattern of doubling or tripling, except for the first few days.
There are numerous models we could use, some more complex than others, but to cut a long story short, a polynomial model turns out to fit the data best. Perfectly, in fact. In particular, the model is a so-called polynomial of degree 3.
But I’m getting slightly ahead of myself.
Spreadsheets like Excel, Numbers, and LibreOffice allow a plot and model to be created, much like the ones above, for which I used the R programming language. Here’s what this looks like in Excel:
There are a few things to notice.
The first is the equation at top right of the plot: the 3rd degree polynomial, 3rd because the largest exponent to which the variable x (day number here) is raised is 3.
The second thing to notice is R2 = 1. This R-squared value is a measure of the “goodness of fit” of the model to the data.
A polynomial of degree 2 is not a good fit, and adding more terms, e.g. to give a 4th degree polynomial, doesn’t improve the fit.
The third thing to notice is the column in the table at left titled Total (model). This is the output of applying the polynomial model to numbers in the Day column. It’s almost the same, but not quite the same as the Total column, which is the actual total number of gifts we’ve been computing through different means all along.
The R-squared goodness of fit measure is very close to 1 (perfect fit) but not exactly. It appears that it’s being rounded up. This happens no matter what spreadsheet software you use. Not only that, but the actual model equation itself differs across software used. LibreOffice gave the smallest error of the spreadsheets I tried because the model had constant multipliers in each term with more decimal places than Excel or Numbers.
So, what’s going on? The answer is that in a computer, real numbers are most often represented using some variant of so-called floating-point format. The dominant standard for floating-point numbers is IEEE-754, created in 1985. It’s very successful for general purpose computation, but errors of various sorts can creep in such that the result is an approximation. The topic of the representation of real numbers for computation would take us down another rabbit hole.
The question is: can we take this approximate model and make it exact? The answer, as Bob the Builder might say is: yes we can!
The coefficient (constant multiplier) 0.1667 in the first term is just an approximation for the rational numberone sixth, a rational number being an exact number represented as the division of two integers (whole numbers), i.e. a fraction.
The coefficient 0.5 in the second term is of course just one half, another rational number.
The coefficient 0.3333 in the 3rd term is just an approximation for the rational number one third.
The 4th term 2E-11 is engineering notation for 2-11 which is 0.00000000002, a very small number, added to the result, again, to adjust for the fact that the model coefficients are approximations of exact rational numbers.
Given this, the model simply becomes:
I’ve used n instead of x here, since this is the variable name we have used all along to denote number of days.
This turns out to give answers that are in perfect agreement with the other methods we have been using to compute the total number of gifts for a given number of days.
As the statistician Paul Box is reported to have said:
All models are wrong, but some are useful.
Sometimes they’re actually not wrong at all, as in this case. The initial polynomial model was “wrong but useful” but the simplified form is exactly right.
Even though the model was created using data for only the first 12 days, it works for larger numbers just as our other methods do, making it truly predictive, something that not all models can claim to be. Admittedly, as models go, this is a simple one.
Here is the working for f(12):
which in words is: one sixth of 12 cubed, plus one half of 12 squared, plus one third of 12, which evaluates to 288 + 72 + 4 = 364.
This polynomial of degree 3 corresponds to one sixth of a 12x12x12 cube (volume 1728) + one half of a 12×12 square (area 144) + one third of a line 12 units long.
This function, f(n), is the most compact, general expression for computing the total gifts for the Nth day of Christmas I’ve found so far.
It also turns out to be the fastest to run on a computer, since no matter how largen becomes, the function executes in constant time, O(1), meaning that computation time is not sensitive to the size of n.
The best we have been able to do before now is linear time, O(n). The last two posts have discussed computational complexity in the context of this problem, so refer to those for more detail.
The actual run-time (at least on my M2 Mac) is around one millionth of a second, for any value of n.
As in the last two posts, I wrote code in the Julia programming language to implement this function. Unlike many programming languages, Julia supports rational numbers as a first class data type. So, one way to write f(n) in Julia as daysofxmas(n) is:
function daysofxmas(n)
1//6*n^3 + 1//2*n^2 + 1//3*n
end
where 1//6 is the rational number (fraction) one sixth, for example, and n^3 is n3.
Here’s an example in which the function is applied to 12:
julia> daysofxmas(12)
364//1
The result is literally the fraction 364 over 1, simply because we were very explicit about using rational numbers as coefficients. While that’s an exact answer, we know that the answer really is just a positive whole number. We can convert the result to a suitable integer data type:
function daysofxmas(n)
Int128(1//6*n^3 + 1//2*n^2 + 1//3*n)
end
A sample run of the modified function gives the simple integer result we’re after:
julia> daysofxmas(12)
364
A reminder from previous posts that the intention is to generalise the gift total calculation to any positive number of days n, i.e. not just n = 12 days, just because we can, not because we have to! 🙂
The final twist is that integers in programming languages have limits, so that performing operations like exponentiation and multiplication on very large integers as we do here can lead to integer “overflow”, resulting in answers that seem to make no sense or give errors:
This can be solved by converting the data type of n itself (the default integer type in Julia is Int64, i.e. a 64 bit integer) to a larger integer type so that intermediate operations within the expression don’t overflow, at least where the value of n becomes too large. This also achieves the rational to integer type conversion added in the previous version of the function. The modified function is as follows:
# 9
function daysofxmas(n)
if n > 100000000
n = Int128(n)
end
1//6*n^3 + 1//2*n^2 + 1//3*n
end
giving the correct, yet ridiculous, total number of gifts for n = 1,000,000,000 (one billion) of 166,666,667,166,666,667,000,000,000 or approximately 167 million billion billion:
While writing this post, I realised that implementations of the daysofxmas function in the previous 2 posts fail to give the correct answer beyond values of n of around 100,000,000. For n = 1,000,000,000, instead of an error like the one above, the total number of gifts given is 2.4 billion billion by these functions instead of 167 million billion billion. The reason, integer overflow, is related to the error encountered above in which the use of rational numbers led to an error instead of an incorrect result due to the vagaries of the way arithmetic operations work in programming languages and the hardware they target.
The next part will summarise the different methods, results, and Julia function run-times, taking into account the correction necessary to give the correct result, even for large values of n.
Beyond this, there appears to be some kind of link between triangular numbers and the polynomial model presented here. This will be explored in the final post. First, part 6 summarises the trip we’ve been on.
I listened to an ABC Radio National Ockham’s Razor podcast episode called “The magic of storytelling… in maths?” awhile ago and the speaker reminded me about triangular numbers. I realised there was a connection between this and the current problem that leads to a more compact way to represent the first method presented in part 1.
Recall also from previous posts the intention to generalise the gift total calculation to any positive number of days n, i.e. not just n = 12 days, just because we can, not because we have to! 🙂
Imagine a triangle of height 12, the interior of which is represented by dots, and so consists of 12 rows of dots. The first row has 1 dot, the second row as 2 dots, the third row has 3 dots, the 12th row has 12 dots, and in general, the nth row has n dots.
The first so-called triangular number is 1 with each subsequent number being that row’s number of dots plus the number of dots in all the previous rows.
So for example, the second triangular number is 2+1 = 3. The 7th triangular number is 7+6+5+4+3+2+1 = 28, as shown above.
This simple function in the Julia programming language yields the nth triangular number:
function triangular(n)
num = 0
for i in 1:n
num += i
end
num
end
with the following run yielding the first 12 triangular numbers:
julia> [triangular(n) for n in 1:12]
12-element Vector{Int64}:
1
3
6
10
15
21
28
36
45
55
66
78
or even more simply:
julia> [sum(1:n) for n in 1:12]
In case you’re wondering, Int64 is just the type (64 bit integer) of the elements of the resulting list or sequence (vector) of triangular numbers.
Adding these numbers gives 364, the number of gifts in the 12 days of Christmas since the triangular numbers are just the green column of the table from part 1.
Now imagine the dots converted to green squares.
If we place another triangle, inverted on top of the green triangle with the same 12 x 12 dimension (in red, as shown), we now have a 12 x 13 rectangle, or in general, n x (n+1), or simply n(n+1).
Since the area of a triangle is half that of a rectangle, the area of the green (or red) triangle is only half of n(n+1).
Here is the mathematical notation for:
the resulting function
applying the function to 12 to get the number of gifts on day 12 (calculation detail shown), and
the resulting summation of function applications for the days 1 to 12:
Since each component of the triangle has unit size (1×1 square, or in other words, a dot), it turns out that the triangle area is equivalent to the sum of all these units, which is the same as the number of gifts for day n.
Here is some Julia code that gives the number of gifts for the 12 days of Christmas using the triangular function:
function triangular(n)
n*(n+1)/2
end
julia> sum([triangular(n) for n in 1:12])
364
We can wrap this up in the function daysofxmas, which I’ll label #8, since the last function in part 3 was #7:
# 8
function daysofxmas(n)
sum([triangular(n) for n in 1:n])
end
In the language of computational complexity and Big-O notation, the second triangular function has time complexity O(1) — or, order 1 — whereas the first implementation of triangular has time complexity O(n).
What O(1) means for n(n+1)/2 is that computing the nth triangular number takes so-called constant time: as fast the a machine can carry out the addition, multiplication, and addition operations. As n gets larger, the time to compute remains the same.
What O(n) means for the first triangular function in this post is that computing the nth triangular number has so-called linear time complexity, or time proportional to the number of days, n. So, as n gets larger, so does the time taken to complete the calculation.
For an n of 12, the time taken in the calculation for both implementations of the function triangular is very small, around a millionth of a second. For n = 1 billion, the first function (using a loop) takes about 4 seconds, whereas the second still just takes around a millionth of a second.
While the area-based triangular function is a big improvement over the brute force approach to find the number of gifts on the nth day (by replacing each row addition), it must still be coupled with a summation step to give the total number of gifts for the 12 days of Christmas, as shown in the example run above (function #8), using sum. This requires O(n) time to complete.
For n = 1 billion, that’s around 6.8 seconds to compute the total of gifts. Not bad, but we can do better, as we saw in the last post where we went down numerous implementation rabbit holes for the daysofxmas function.
Note: Before wrapping up (gift pun unintended) below, it’s worth going back to read part 3 later if you haven’t yet, to peek down those rabbit holes, but it’s quite a bit longer than this post and may be best viewed on a screen larger than a smartphone’s (or printed). While it doesn’t introduce additional mathematical approaches it’s likely to be interesting to anyone wanting to understand more about the computational approach, the tradeoffs of different algorithms, why I chose to use the Julia programming language here, optimisation, and the likelihood that your intuition regarding how to make things run faster is very likely to lead you astray. It also spells out further the relationship between the brute force approaches and the mathematics (e.g. see reference to dot product in table below).
So, let’s summarise the performance of all the functions from part 3 and part 4 now, this time with n = 100,000 — a paltry 100 thousand days of Christmas, yielding a mere 166,671,666,700,000 — about 167 million million gifts.
Function
Time (seconds)
#1: single loop with sum function
0.0005
#2: two loops
24
#3: two sum functions
0.0004
#4: dot product
0.0007
#5: dot product refinement
0.0007
#6: parallelised version of #2
13
#7: remembering all previous day totals
0.0005
#8: summation over triangular numbers
0.0003
Run time of daysofxmas(100000)for each of the functions presented in parts 3 and 4
Apart from the obvious outliers, #2 and #6, there’s not much difference for n=100,000. As n becomes larger, #2 and #6 become impractical and the difference between the others becomes more obvious. See the next post for more.
What I said in part 1 is that I want the most compact function (the simplest, with the least number of terms) possible that computes the gift total for any n. But because I also want this to be done in the shortest time possible, the function daysofxmas itself must have constant time, not linear time complexity. That means not doing any kind of summation at all!
Is there a more compact function that does not require the use of a summation step? Answering that question is the focus of part 5.
Before moving on to the next approach to the problem, this post presents code — functions — in the programming language Julia that automate the row and column methods of the first two posts.
In addition, these functions generalise to any positive number of days n. This will come in handy in the event that we are ever invaded by an alien species or artificial general intelligence is achieved, and our alien or robot overlords command that there should be say, 100 days, instead of 12 days of Christmas. 🙂
When I first wrote this post, I included examples in the R and Python languages in addition to Julia, an emerging language that’s gaining popularity, especially in scientific computing which has been dominated for decades by Fortran, MATLAB, R, Python and a handful of other languages.
In the end I decided that presenting examples in just one language was best since it avoids arbitrary differences between languages and, frankly, some languages have counter-intuitive constructs, all of which serves only to distract from the main point of the post.
For the first method in part 1, summing all the gifts received per day and adding all the the resulting values (the green column of the table in part 2), looks like this in mathematical notation:
The variable i below the right sigma (summation) symbol takes on values from 1 to day, all of which are added, leading to the day totals in the green column of part 2’s table (shown below again).
The left summation corresponds to the addition of each daytotal number in the green column to arrive at the final gift total in purple.
The table from part 2 is shown here for reference:
Here is a Julia function, daysofxmas, that is one possible implementation of this notation, taking a parameter nand returning the total gifts received in n days. Invocations of the function where n is 12 and 100 are shown. The line with “# 1” is a comment. I’ll use this to number distinct functions for reference.
The for loop corresponds to the left-most summation in the mathematical notation,
sum(1:day) corresponds to the right-most summation in which each whole number from 1 to day is added and assigned to the variable daysum, and
each day’s sum is accumulated in the variable total which is returned from the function.
# 1
function daysofxmas(n)
total = 0
for day in 1:n
daysum = sum(1:day)
total += daysum
end
total
end
julia> daysofxmas(12)
364
julia> daysofxmas(100)
171700
171,700 is obviously a lot of gifts and we wouldn’t want to have to add all those gifts per day by hand! 364 is bad enough.
Adding a formatted print call to the function shows day number, sum of gifts per day, and cumulative total, after which the final total is printed, as before. The weird-looking formatting strings (e.g. %2i) tell Julia how many spaces a number should occupy (2 or 3 here) in the output and what type of numbers to expect (integers, i.e. whole numbers).
using Printf: @printf
function daysofxmas(n)
total = 0
for day in 1:n
daysum = sum(1:day)
total += daysum
@printf("%2i: %2i => %3i\n",
day, daysum, total)
end
total
end
julia> daysofxmas(12)
1: 1 => 1
2: 3 => 4
3: 6 => 10
4: 10 => 20
5: 15 => 35
6: 21 => 56
7: 28 => 84
8: 36 => 120
9: 45 => 165
10: 55 => 220
11: 66 => 286
12: 78 => 364
364
Two for loops could have been used instead, which has the advantage of matching the notation’s summation structure more explicitly, but requires a bit more code and is not as fast as the first implementation above (which we’ll return to below).
# 2
function daysofxmas(n)
total = 0
for day in 1:n
daysum = 0
for i in 1:day
daysum += i
end
total += daysum
end
total
end
Instead, a more compact form of this function uses a language feature known as the list comprehension. Here, the inner sum yields a list of per-day gift totals, which the outer sum adds up (again, the green column of the table in part 2).
# 3
function daysofxmas(n)
sum([sum(1:day) for day in 1:n])
end
This is faster than function #2 and about the same as #1. It is also a close match to the mathematical sigma notation above, represented as Julia code.
In addition, it is more declarative, allowing us to focus more on expressing what we want to do rather than how. Another advantage is that it allows a programming language like Julia to have better control over what to optimise for best performance.
Building up the expression returned from line 3 in stages with n= 12, makes it easier to see what’s going on, if you are interested in the detail:
julia> [day for day in 1:12]
12-element Vector{Int64}:
1
2
3
4
5
6
7
8
9
10
11
12
julia> [sum(1:day) for day in 1:12]
12-element Vector{Int64}:
1
3
6
10
15
21
28
36
45
55
66
78
julia> sum([sum(1:day) for day in 1:12])
364
Now recall the second method presented in part 2: the result of summing the blue cells in the table, where each corresponding value in the sequences 1, 2, …, 12 and 12, …, 2, 1 is multiplied, with the products then being added:
In mathematics, this is known as the dot product:
Again, 12 can be generalised to n:
The following Julia code defines an implementation of the function daysofxmas that carries out this particular dot product.
Again, invocations of the function where n is 12 and 100 are shown.
#4
function daysofxmas(n)
days = [1:n;]
sum(days .* reverse(days))
end
julia> daysofxmas(12)
364
julia> daysofxmas(100)
171700
[1:n;] means: generate a sequence, a vector, of all the whole numbers from 1 to n.
Use of the reverse function here can be replaced by [n:-1:1;] which says: give me a vector of all the whole numbers from n to 1, where -1 means “step -1”, i.e. start from n and count down by 1, stopping at 1. The end result is the same.
All of the functions presented here yield their result in a tiny fraction of a second, even for a silly 1000 or 10,000 days of Christmas (166,716,670,000 or more than 166 billion gifts), both of which would require the length of the year to change as well of course!
Even for a truly stupid 100,000 days of Christmas, the time to compute the result (166,671,666,700,000 or more than 166 trillion gifts) is less than a thousandth of a second or better, faster than either Python or R on my Mac (M2 processor, 32 GB RAM).
That is, except for the function definition with the two nested loops (#2), which takes 3 seconds, more than a thousand times slower than all the other function definitions shown. There are good reasons for this that have to do with efficient use of processor resources. I’ll return to this point.
Function #3 could also be modified to implement the shortened form of the calculation, twice the first half of the series, as discussed in part 2: 2 x (12 + 22 + 30 + 36 + 40 + 42), but would need to take into account whether n is odd or even, since a sequence with an odd number of values cannot be split in half such that both halves have the same number of values and the left half is the reverse of the right half. Think of the sequences 1, 2, 3, 2, 1 and 1, 2, 3, 3, 2, 1. The second satisfies this constraint but the first doesn’t.
I was going to avoid that complication, since the code runs fast anyway compared with the effort of manual calculation. But, I couldn’t help myself. 🙂
One implementation of this modification to function #3 in Julia is shown next.
There’s no need to pay too much attention to this code unless you want to. The commentary that remains in this post is the main thing to take note of.
# 5
function daysofxmas(n)
days = [1:n;]
products = days .* reverse(days)
numproducts = length(products)
halflength = div(numproducts, 2)
total = 2 * sum(products[1:halflength])
if isodd(numproducts)
total += products[halflength+1]
end
total
end
This implementation is more complex and this greater code complexity and effort may be misplaced.
Just because one mathematical expression is shorter to write down than another, doesn’t mean that implementing it as a function will make it faster or require less code. For one thing, I’ve had to add more code to take account of the odd-length-check requirement. Most of the code after line 4 is new.
As the famous Computer Scientist, Donald Knuth said:
Well, maybe not all evil, but it’s generally a bad idea to make guesses about whether a different approach will be faster or which parts of your code need to run faster, since your intuition or common sense will very likely fail you. That’s what the discipline of computational complexity and profiling tools are for.
No, not social profiling. Code profiling.
It turns out that when you use a profiling tool on this function, most of the total runtime is spent in the dot product expression: days .* reverse(days). So, working with half of the dot product doesn’t help since it has already been computed.
Any implementation that improves the runtime of #4 or #5 would have to compute the dot product in a different way (e.g. just half of it) or compute it in parallel. Both are possible and I will briefly explore the second of these here, not with the dot multiply implementations of the function (#4 or #5), but with function #2.
Why? Interestingly, it turns out that the slowest definition of the function (#2) we’ve seen so far, the one with two loops and simple addition, can greatly benefit from a parallel implementation.
But to really stress the daysofxmas function, let’s make n = 100,000. A billion days of Christmas! This results in 166,666,667,166,666,667,000,000,000 gifts. That’s 167 million billion billion.
So, how long does this take to run for each function implemented so far? Most
Function
Time (seconds)
#1
0.0005
#2
24
#3
0.0004
#4
0.0007
#5
0.0007
Run time ofdaysofxmas(100000) for each of the 4 functions presented so far
If function #2 is modified in the following apparently minimal way…
# 6
function daysofxmas(n)
total = 0
@simd for day in 1:n
daysum = 0
@simd for i in 1:day
daysum += i
end
total += daysum
end
total
end
…the run-time for a billion days of Christmas becomes 13 seconds!
And, how has the code been changed? A Julia directive (actually a so-called macro) @simd has been added before each for loop. SIMD is an acronym for Single Instruction Multiple Data.
What this means here is that the work of the iterations (for loops) is split up, such that chunks of each vector are operated on by the same machine instructions simultaneously rather than sequentially, i.e. one after another.
This is one example of parallel execution. Note that just randomly adding such a directive to any old for loop won’t necessarily improve the run-time because of code dependencies within the loop. It may in fact make it worse. Such is the case if @simd is added before the for loop in function #1, which roughly doubles the runtime.
Exploring what actually happens when such a directive is added to a language like Julia as well as the nuances of parallel computing is, well, a whole other story.
Is there another implementation that is simpler and competes with #1 or #6 for speed? Yes there is. It turns out to be a minor variation on function #1.
# 7
function daysofxmas(n)
total = 0
daytotal = 0
for day in 1:n
daytotal += day
total += daytotal
end
total
end
Here we use a cumulative day total, remembering all the day additions previously done, only adding the current day number to the day total before updating the overall total. This takes the same amount of time to run for 100,000 days of Christmas as function #1: 0.0005 seconds.
What this last example underscores is that the algorithm matters as dictated by the discipline of computational complexity. Function #7 has time complexity of O(n) which is so-called Big-O notation meaning that there is a simple linear relationship between the number of days of Christmas, n, and the time it takes to compute the number of gifts. This is in contrast with function #2 that has a time complexity of O(n2) such that the time taken to compute the gifts for n days is the square of n, which for large n can get big very quickly, e.g. for n = 10,000, n2 is 100,000,000 instead of just n.
As an aside, #4 and #5 use more computer memory than all the other functions, because of the large vectors they hold in memory.
But enough of this!
Back to the mathematics (and a bit of code) in part 4, but still very much continuing down the path to generalisation and brevity.
The yellow row of the table above highlights an example of one of the days: the 21 gifts given on the sixth day: 6 geese a-laying, 5 golden rings, 4 calling birds, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
The right-most column (green) is the sum of each day’s gift total using the approach described in part 1.
The bottom right-most cell (purple) shows the addition of the numbers in the green cells to give the total number of gifts for the 12 days of Christmas.
The final row (blue) shows that there are:
12 lots of gift 1 = 12
11 lots of gift 2 = 22
10 lots of gift 3 = 30
9 lots of gift 4 = 36
8 lots of gift 5 = 40
7 lots of gift 6 = 42
6 lots of gift 7 = 42
5 lots of gift 8 = 40
4 lots of gift 9 = 36
3 lots of gift 10 = 30
2 lots of gift 11 = 22
1 lot of gift 12 = 12
These values also add to 364.
This amounts to multiplying each value in the sequence 12 to 1 by those numbers in reverse order, 1 to 12, then adding the resulting products:
You may notice that the first six terms in the addition above repeat, but in reverse, after the first occurrence of 42 (the meaning of life!), so the calculation can be shortened to:
This seems like a nice simplification.
But as nice as this is, it doesn’t help us with generalising to any value of n, such as 100, since even if the summation sequence can be halved for any sequence, for large values of n, that’s still a lot of numbers to multiply and add.
We’ll start exploring how to generalise to any value of n in part 3.
364 gifts, one gift for each day in a non-leap year (365 days). Every day except Christmas, you could say.
I’m not good at coming up with answers to this sort of question quickly, but I enjoy thinking about the underlying problem, different approaches, and in particular, the simplest and most general solution. This kind of recreational mathematics also often takes me down unexpected paths.
I can never remember the lyrics, so here’s a reminder of The 12 days of Christmas, a carol in which each each day, gifts are given to the recipient by a true love. The following omits a few days (verses) for brevity:
On the first day of Christmas, my true love sent to me, a partridge in a pear tree.
On the second day of Christmas, my true love sent to me, 2 turtle doves and partridge in a pear tree.
On the third day of Christmas, my true love sent to me, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
On the fourth day of Christmas, my true love sent to me, 4 calling birds, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
On the fifth day of Christmas, my true love sent to me, 5 golden rings, 4 calling birds, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
…
On the eleventh day of Christmas, my true love sent to me, 11 pipers piping, 10 lords a-leaping, 9 ladies dancing, 8 maids a-milking, 7 swans a-swimming, 6 geese a-laying, 5 golden rings, 4 calling birds, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
On the twelfth day of Christmas, my true love sent to me, 12 drummers drumming, 11 pipers piping, 10 lords a-leaping, 9 ladies dancing, 8 maids a-milking, 7 swans a-swimming, 6 geese a-laying, 5 golden rings, 4 calling birds, 3 french hens, 2 turtle doves, and a partridge in a pear tree.
For the purpose of solving the problem, we can ignore the actual gifts and just focus on the numbers, e.g. the two turtle doves becomes just 2 gifts of some kind, in this case, 2 lots of the same kind of gift (turtle dove) on day 2, or simply 2 lots of gift number 2.
A first intuition leads to adding all the gift numbers from 1 to 12 giving 78, but that’s only part of the answer.
On the first day we are given 1 gift (a partridge in a pear tree).
On the second day we are given 1 gift plus 2 more (a partridge in a pear tree, two turtle doves) giving 1+2 = 3 gifts.
On the third day we are given 1 gift plus 2 more plus 3 more (a partridge in a pear tree, two turtle doves, three french hens) giving 1+2+3 = 6 gifts.
And so on.
Of course, the ordering in the carol’s verses is the reverse of 1+2+3+…+n, given the order in which the gifts appear, but since 1+2+3+…+n = n+…+3+2+1, the summation is unaffected.
When we were all trying to come up with the right number that Christmas time, I couldn’t help but wonder whether there were simpler ways to arrive at this, more compact expressions. It turns out there are.
What interests me is this: finding the most compact mathematical expression that gives the same answer as the brute force row addition method above, whether the n in “The n days of Christmas” is 12 or 100. Imagine the time and effort required if there were 100 days of Christmas!
The second post will explore two approaches, first very briefly revisiting the row addition approach above, then moving on to column multiplication and addition of these products. It’s much easier to see this working than explain it abstractly.
The third post will generalise the approaches so far to any value of n, not just 12, via simple code in different programming languages.
The fourth post will present an approach involving the day 1 to 12 (or 1 to n) results as triangular numbers and the sum of these. I arrived at this after listening to a podcast by chance that helped me think about the problem in a different way. This approach allows the generalisation to any value of n with a mathematical expression rather than via programming.
The fifth post will present the most compact general solution for any value of n that I have found so far, arrived at by plotting and modelling the results. There will also be a short digression on rational numbers vs real numbers and problems arising from the representation of real numbers in computers.
The sixth post will summarise the approaches and their relative performance.
The seventh and final post will show that a potentially interesting link exists between the last two approaches, by employing a bit of calculus.
All of this in a single post would be a bit much. Several short posts seemed more appropriate.
Following on from my last post On errors in math and code, here’s an example (adapted from Khan Academy) of how an error can be caught with a small amount of symbolic prompting along the way towards the final numeric answer.
For the purpose of reading this post, it doesn’t matter so much whether you understand all the steps, more that taken together, the whole thing is complex enough that you can see how an error is possible. Even a simpler problem could be enough for this.
The circle of radius 4 above has the following equation:
where 16 is .
The green lines perpendicular to the x-axis each form the base of an equilateral triangle. For a given x value, the line extends the same amount, y, above and below the x-axis, so has length 2y, as shown below.
If you imagine many such triangles closely packed together over the circular base, the 3D shape at left below emerges with (a) The solid, showing a few equilateral triangles extending above the x-y plane. Note that, incidentally, in the example below, the circle’s radius is 1 not 4 as in our example (not important for how to approach a solution to the problem).
The problem: how to compute the volume of this 3 dimensional solid? The answer requires the use of integral calculus, which allows the area or volume of non-trivial functions to be computed.
To compute the volume of the 3D shape presented in the problem, we need to follow these steps:
Determine the triangle height by taking half of the yellow equilateral triangle, giving a right angle triangle, allowing us to apply the Pythagorean Theorem to determine the height. This is already shown in the yellow triangle above as the square root of 3 but we will show how to calculate it given a hypotenuse of length 2y and a base of half that (y).
Determine the area of a triangle first in terms of y, then x since we will be summing areas over a range on the x axis.
Calculate the volume of the 3D shape by summing the area of an “infinite number” of triangles by computing the definite integral over the appropriate x axis range.
Along the way I will show an error and how a sanity check could serve as a correction before proceeding to calculate the final answer.
Triangle height using the Pythagorean Theorem:
Triangle area in terms of y:
Triangle area in terms of x (from circle equation), substituting for y squared:
Volume by integration over area of triangles of width dx (delta x), where dx becomes very small, so there are “infinitely many” triangles.
This is almost right, but the bounds are wrong. Look at the circle at the start of this post. It spans -4 to 4 on the x axis, not 0 to 4. The volume should therefore be twice that calculated above. We get that from what follows, with the correct integration bounds.
But if just one half of the circle’s base is used for determining the volume, the math is simpler. Compare the last two solutions: in the first there are fewer numbers in the calculation, so less opportunity for error. If one half of the volume is computed, the whole volume is simply twice this, and that was what was initially in my mind; I simply forgot to multiply by two! Taking this into account gives:
So, a simple error (wrong bounds or not multiplying by two) was introduced in the very first step of the integration but did not become tangible until we started substituting numbers.
A graphing calculator application like Desmos first confirms the result of using the erroneous bounds and then the equivalence of the last two approaches:
After all that, the main point I’m making is that guidance during learning can be given along the way to a solution (in the textbook or online resource where the problem appears) by showing multiple choices for the equation itself, not the numeric answer, such as the first and last above along with some incorrect ones. This may be enough for someone to realise that a mistake has been made before proceeding to substitute numbers.
Online learning resources often use the multiple choice approach or the “just give the numeric answer” approach, e.g. to three decimal places or in terms of , but it seems less common to see a two stage process in which the correct equation must be selected followed by the numeric answer given. I wonder about the potential benefit of this, especially when first learning a specific, complex subject area in which mistakes are easy to make.
An argument against this is that it is better to allow failure and provide step-by-step guidance in the event of an incorrect answer, something Khan Academy excels at. Of course, this doesn’t guarantee that the same mistake won’t be made across similar problems.
.
, but it seems less common to see a two stage process in which the correct equation must be selected followed by the numeric answer given. I wonder about the potential benefit of this, especially when first learning a specific, complex subject area in which mistakes are easy to make.