Start here

The number of possible bingo cards

During a recent holiday I was sitting through the bingo session that some younger relatives were playing.  My mind wandered onto the question of how many possible cards exist.

A bingo card (as played in the UK, maybe elsewhere too?) has 27 boxes arranged as 3 rows by 9 columns.  Each row has 5 numbers and 4 blanks.  No column may be all blank.  The first column contains numbers between 1-9, the second contains numbers between 10-19, and so on.  The last 2 columns contain numbers between 70-79 and 80-90 respectively.  Note the uneven distribution of possible numbers in the columns.  Numbers are in increasing order downwards if more than one appears in a column.  I don’t believe it’s a requirement for there to be (or not to be) a column with 3 numbers in, but both cases are computed below for completeness.

So moving straight to the answer: 3,669,688,706,217,187,500

If exactly one column must have 3 numbers: 1,957,437,874,940,625,000

If at least one column must have 3 numbers: 2,494,937,545,340,625,000

The numbers were computed with the Python script below:


# Look up table of how many ways there are to fill the nth column with numbers
# from the appropriate range, based on how many non-blank v blank entries there
# are
perms_per_col =      [( 9,  9* 8//2,  9* 8*7//6)] # 1-9, 9 numbers
perms_per_col += 7 * [(10, 10* 9//2, 10* 9*8//6)] # 2nd-8th cols, 10 numbers
perms_per_col +=     [(11, 11*10//2, 11*10*9//6)] # 80-90, 11 numbers

# Build list of 4-tuples containing the possible indices (0...8) of 4 blanks
# within a row of 9 cells
patterns = []
for x0 in range(0,9):
    for x1 in range(x0+1,9):
        for x2 in range(x1+1,9):
            for x3 in range(x2+1,9):
                patterns += [(x0, x1, x2, x3)]

total = 0
total_one_three = 0
total_any_threes = 0

# 3 nested loops, over the 3 rows in the card
for i0, p0 in enumerate(patterns):
    # 'comfort' message to show progress
    print ("%d" % (i0,))
    for p1 in patterns:
        for p2 in patterns:
            # How many blank boxes there are in each column?
            blanks = [0] * 9
            for i in p0: blanks[i] += 1
            for i in p1: blanks[i] += 1
            for i in p2: blanks[i] += 1
            threes = 0 # Number of columns with 3 non-blank cells
            cards_with_this_layout = 1
            for col in range(0,9):
                # No column is allowed to have 3 blanks.  If one does, bail out of this pattern
                if blanks[col] == 3:
                    cards_with_this_layout = 0
                if blanks[col] == 0:
                    threes += 1
                table_index = 2 - blanks[col]
                cards_with_this_layout *= perms_per_col[col][table_index]
            total += cards_with_this_layout
            if threes == 1:
                total_one_three += cards_with_this_layout
            if threes > 0:
                total_any_threes += cards_with_this_layout

print ("There are {:,d} cards possible".format(total))
print ("There are {:,d} cards possible with exactly one column having three numbers".format(total_one_three))
print ("There are {:,d} cards possible with at least one column having three numbers".format(total_any_threes))

Approximating the Gudermannian function

The what?  It’s the formula for inverting the Mercator projection in the latitude direction.  I’ve recently added to LogMyGsm a display of the latitude at the point in the centre of the map display.  The internal coordinates are in Web Mercator space to ease the map tile handling.  To get the latitude for display, the gd(Y) operation has to be applied.

The standard formula to return degrees latitude from the Y coordinate of the projection is gd(Y), coded thus:

static inline double exact_lat(double y)
    double t = (1.0 - 2.0*y)*M_PI;
    double M = (2.0* atan(exp(t))) - 0.5*M_PI;
    return M*(180.0/M_PI);

Given that this has to be run for every update as the map is dragged (so tens of updates a second possibly), the atan and exp might be a little wasteful. And optimising is fun…
The following is an approximation

#define P1 179.9989063857
#define P3 507.2276380744
#define P5 176.2675623673
#define Q0 1.0000000000
#define Q2 4.4623636863
#define Q4 4.2727924855
#define Q6 0.4175728442

static inline double approx_lat(double y)
    double z = 1.0 - 2.0*y;
    double z2 = z*z;
    double z4 = z2*z2;
    double t, b0, b4, b, t2;
    b4 = Q4 + Q6*z2;
    b0 = Q0 + Q2*z2;
    t2 = P3 + P5*z2;
    b = b0 + b4*z4;
    t = P1 + t2*z2;
    return z*t/b;

This looks like more code.  However, there are no calls to the transcendental libm functions exp and atan here.  On an x86-64 machine, this is more than 10x faster than the exact version.  I’m hoping the conclusion would be similar on the Cortex-A9 and ARM11 in my 2 smartphones.  To be confirmed!

As for the accuracy: it’s within 5e-5 degrees to beyond 71 degrees N/S, and within 7e-4 degrees between 71 and 85 degrees.  So the accuracy up to 71 degrees is enough to give +/-1 in  the 4th decimal place.  There isn’t much habitable landmass beyond 71 degrees, so the accuracy there was compromised to improve the accuracy below 71 degrees.


The degrees error versus latitude is shown in the above.  Notice the nice equiripple error behaviour up to 71 degrees.

At some point I will need to write up the derivation properly.  In outline : a pade approximation was computed by solving a set of linear equations at a set of Y values in [0,1].  The initial set of Y values were Chebyshev nodes.  After each iteration, the Y values were adjusted to focus attention around the place with the greatest remaining error, and the procedure repeated.  Luckily, this converged to the minimax solution shown.








Updating Arch Linux, self-inflicted disaster, and a fingers-on-the-cliff-edge recovery

Yesterday evening I finally got around to updating the Arch Linux installation on my main machine.  I tend to put this task off – after a few recent updates I’ve been stuck with an unbootable machine and had to reach for a rescue CD.  Once it was the switch to systemd, another time it was a change in how LVM started.  Last night was to be no exception, only this time the disaster was for once entirely of my own creation, although it stemmed from the move of all the binaries out of /bin, /sbin and so on into /usr/bin.

The initial pacman -Syu failed because of those binaries moving. The fix was clearly explained here. Of course the first stage took ages; downloading and installing nearly 500 packages is not a quick job.  So I went off to do something else.  By the time this was done, it was already 1am and time to think about some sleep.  By this time I’d forgotten there were two more (perhaps important!) steps in the recipe.  But I thought I’d better check the machine would still boot, just so I knew what I’d be waking up to in the morning.  One systemctl reboot later, and I was staring at a brick.  Booting is hard with no /bin/mount, /bin/fsck and so on.  It suddenly dawned on me what I had done.

So SystemRescueCD was booted, the method of choice in the past.  I was hopeful of running pacman with a suitable set of switches to fix up the mess.  Unfortunately, the systemrescue CD I had in the cupboard was only a 32-bit system, and the Arch system is 64-bit.  Oh dear.

Next, I remembered there was an ancient Ubuntu 10.04 still present on another partition.  I got that booted, fortunately remembered the passwords, and was in.  By suitable use of setting LD_LIBRARY_PATH and using /mnt/tmp/usr/lib64/ld-linux*.so to run pacman, I even got pacman to run queries.  But it wouldn’t quite do the rest of the installation – I think it was still looking in the wrong place for something.  By this time it was very late and I gave up for the night.

Sleep works wonders for the creativity.  This morning, I decided to try manually creating symlinks in the Arch filesystem for /bin/mount and friends from within Ubuntu, to see if that was enough to get Arch bootable again.  Viola!  Well almost.  It booted, but seemed to be lacking whatever file tells systemd to start some gettys on the consoles.  But luckily, sshd was running.  I used connectbot on my phone to get into the machine, remove the temporary symlinks from /bin, and complete the install by running the 2nd and 3rd steps of the recipe linked to earlier.

Now there is just the usual pile of .pacnew files to deal with.  Or ignore.

This has to be the closest I’ve come to needing a clean install in quite a while…

Minimax approximation to arctan / atan / atan2

Even during 10 years in a past life working in a group that designed FIR filters, I never grappled with the Remez exchange algorithm hands-on.  This week I finally found the excuse I needed.

A while ago, I needed a crude but very fast approximation to arctan.  In fact it was atan2, but I’ll start with the easier problem first.  With a release of the application imminent, I wanted to credit the original source I got the formula from, but I couldn’t find it any longer.  So I set about recomputing the formula for myself.

The version I found a while ago looked like this:

atan(x) ≈ (0.97179803008 * x) - (0.19065470515 * x**3)

for x in the range [-1,+1], i.e. the result is in the range [-pi/4,+pi/4].  The maximum error in the result is just less than the equivalent of 0.3 degrees.

[If the person behind these coefficients gets in touch, I’ll be more than happy to give credit and link to the original work here.]

Having now got the Remez exchange algorithm working, I found a slightly improved version is

atan(x) ≈ (0.97239 * x) - (0.19195 * x**3)

which has a worst-case error of 0.2837 degrees.  Even better, a plot shows the 6 extrema all have this absolute error, as expected.

Moving to the original atan2 problem, the requirement was rather unusual.  Starting with the x and y arguments as integers, the answer is to be expressed as a 32-bit integer, where 2**32 represents 360 degrees.  This allows subsequent rotations to be done using integer add and subtract, with the convenient property that multiples of 360 degrees are lost by the modulo wraparound at 2**32.  The resulting function looks like this, with the coefficients as shown earlier:

#define K1 (0.97239)
#define K2 (0.19195)
#define K3 ((float)((double)(1<<29) * 4.0 * K1 / M_PI))
#define K4 ((float)((double)(1<<29) * 4.0 * K2 / M_PI))

int32_t atan2i_approx2(int32_t y, int32_t x)
  int32_t u, v, s;
  uint32_t w;
  float f, f2, g;
  int32_t r;
  float fx, fy;
  float tt, bb;

  fx = (float) x;
  fy = (float) y;
  u = x + y;
  v = x - y;
  w = (u ^ v) >> 31;
  tt = (w) ? -fx : fy;
  bb = (w) ? fy : fx;
  f = tt / bb;
  f2 = f * f;
  g = f * (K3 - K4*f2);
  s = (((u >> 30) & 3) ^ ((v >> 31) & 1)) << 30;
  r = s + (int32_t) g;
  return r;

LogMyGSM explanation

The picture below tries to explain the on-screen display of LogMyGSM.


The source code is on github.  If you’ve got an Android development environment, you can download and build it for yourself.  I haven’t released it on the market yet.


Climbing to the top of the Mendips

Earlier on, I was at the top of the Mendips for the 2nd time – this time with my daughter (I was on my own last time in Feb 2007 when this picture was taken).  Sadly the weather did little to add to the occasion – as seems to be happening to our walks all too often at the moment.  Let’s just say by the time we reached Beacon Batch we were in a howling gale, driving drizzle, and my phone was wrapped in a plastic bag.  I thought it best not to attempt a Foursquare check-in or a “tweet from the top” as I generally like to do! So I’m left with this screen shot off my phone after we got back to the car as a momento.


We were parked to the right of “Mendip Lodge Wood”, at the car park next to the Burrington Inn.  We did the route clockwise.  The trail going off the left of the screen is the “long” route that we’d hoped to do, but we decided to abort half-way due to the conditions.