Archive for January, 2024

On the N Days Of Christmas, part 4

January 19, 2024

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.

FunctionTime (seconds)
#1: single loop with sum function0.0005
#2: two loops24
#3: two sum functions0.0004
#4: dot product0.0007
#5: dot product refinement0.0007
#6: parallelised version of #213
#7: remembering all previous day totals 0.0005
#8: summation over triangular numbers0.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.

Posted in Mathematics, Programming | 1 Comment »

On the N Days of Christmas, part 3

January 14, 2024

man in grey sweater holding yellow sticky note
Photo by hitesh choudhary on Pexels.com

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 day total 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 n and 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:

source: http://realmwalkofsoul.weebly.com/developer-diary/development-update-2-optimization-the-root-of-all-evil

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

FunctionTime (seconds)
#10.0005
#224
#30.0004
#40.0007
#50.0007
Run time of daysofxmas(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.

Posted in Mathematics, Programming, Programming Quotes | 1 Comment »