Forget π (pi), τ (tau) is where it's at

Below is Vi Hart's video explaining why τ (= 2π) is the more natural "circle constant". Everything just works out more neatly with τ and thinking about a full turn rather than a half turn just makes sense. And yes,

e^iτ = 1

is how it should be. The Tau Manifesto has more details but suffice to say, while Mar 14th never did much for me, on Jun 28th I might just make an effort this year.

No this is not entirely a waste of time, using τ with beginners would actually eliminate quite a bit of confusion.

Bit of a maths puzzle.

Ugh! This comes out nicely on my blog but looks like crap in Google reader, presumably because the javascript can't run.
Define a function:
\[ f(x) = \frac{2x}{1 - x^2} \] In case that maths layout thing isn't working, that's f(x) = 2x/(1-x^2)
Now, \[ f(0) = 0 \] so that repeats. Also \[ f(\sqrt 3) = -\sqrt 3 \] and \[ f(\sqrt -3) = \sqrt 3 \] so that repeats too.
The question is, what values repeat if you keep applying f to the result? I thought I knew the answer but now I think maybe it's more tricky.

Why you can't have infinite data compression.

I was leaving a comment on this blog post and it got kind of long so I thought I'd make it a post (and got a whole lot longer!).

Even if you can squeeze a quart into a pint pot...

Compression methods make lots of files smaller. Zip can shrink text because they have lots of redundancy. In fact most things created by humans have lots of redundancy in them. Repeating something in a different way it helps us grasp understand more easily. Redundancy also defends against small errors. So if, in my first sentence, I wrote "compresion", with just one "s" you would still know that I probably meant "compression". Maybe I meant "comprehension" but from the title and the rest of the sentence you could rule that out.

It seems like maybe, if you were clever enough, you could somehow compress a file, then compress the compressed file and then do that again and again and end using lots less disk space (at the expense of doing lots of compression and decompression all the time).

... you can't squeeze 511 quarts into 255 pint pots.

To see why this cannot be done, no matter how clever you are, imagine you had a compressor that always made things smaller. Let's consider all possible strings with 8 or fewer bits. There are 511 (2^9 - 1) of them. If they all get smaller then they must all turn into strings of 7 or fewer bits. There are 255 (2^8 - 1) of those.

Now take all of these 7 bit strings and uncompress them. You have 255 inputs so you get 255 outputs, but that's not enough. You started off by compressing 511 things and you only got half of them back. Some of your input strings must have compressed down to the same output string. That's no good, you need to always get the original back. A compressor that doesn't give the original back is called "lossy" and while that can be great for pictures and movies where the human eye discards lots of the fine detail anyway, it's not the kind of compression that you need for legal documents or computer software.

What if I only compress the ones that get smaller?

It seems like you can still do better by only compressing the ones that get smaller and leaving the others uncompressed. The problem here is that while the strings get smaller on average, you now have to remember which ones you compressed and which you didn't. So now you have to add an extra bit to signify whether it was compressed. There are 256 exactly-8-bit strings and 255 7-or-fewer-bit strings which means at least one of the 8 bit is left as an 8-bit string. Now add the "was I compressed" bit and you get a 9-bit string. It got bigger! Actually it's a little bit more subtle than this because you might say that for an 8-bit, uncompressed string you don't need the extra bit, it obviously wasn't compressed but you must remember that some of the 7-bit strings weren't compressed either and so became 8-bit strings. You have to be able to tell them apart from an 8-bit string that wasn't compressed. Either way, you have 511 inputs and you don't want to mix any up so you need 511 outputs and in the best case, that requires strings up to 8 bits long.

Working at it from the ground up.

Another way to come at this is by thinking about what such a miracle compressor would do to short strings. First, consider the 1-bit strings. There are 2 of them, "0" and "1". Since this compressor must make everything shorter it's stuck already. So let's make an exception for 1-bit strings and say that they always compress to 1-bit strings. So what about 2-bit strings? Well, if it turned any of those into 1-bit strings we'd have a problem because the 1-bit outputs are already taken (for the 1-bit strings). So let's make an exception for the 2-bit strings as well and let them compress to 2-bit outputs. What about the 3-bit strings? Well, all the 2-bit outputs are already taken, so we'll have to make an exception for the 3-bit strings too... Eventually we'll make an exception for every input, no exceptions!

Back in the real world.

There appears to be a loop-hole in real-world computing. It seems you can always get a win-win solution by only compressing the files that compress well and applying the right method. So Zip your text files, PNG your images and don't compress anything that doesn't compress well. You always save space or break even. This seems to contradict what I wrote above. That's because it ignores the fact that filenames take up space. It only works because your computer has already allocated some number of bytes for the name you're going to give the file.. So even if you give it a really short name, the extra space is used anyway. Also, if you already had a really long name, adding .zip on the end may not be possible as it would exceed the size-limit on filenames. Some filesystems don't have a limit and only use as much space as the length of the name but then adding .zip takes up 4 more bytes.

Just being silly now.

If you use a filesystem that does not have a limit on the length of a file name then the following compression method "works". Convert the file you want to compress to a number (that's essentially what your computer does when you save a file anyway - the numbers are big, this blog post has about 6,000 letters in it so would be saved as a number with about 14,000 digits). Now compress this number by subtracting 1 from it and sticking a .z on the end of the file name. Now do that again and again and again. Eventually the number will be 0 and the file data will take no space on disk. When you want to read your file you just keep uncompressing it until all the .zs are gone from the file name. Great, now if only you could find a way to compress the filenames you'd be set...

The truth.

Do you promise to tell the truth, the whole truth and nothing but the truth?

I do... hold on, do you want it depth first or breadth first?

Computing combinatorials

I was commenting on this fun post about efficiently computing n-choose-c but my attempt to post some code seems to have been eaten by the comment system (I suspect use of square and angle brackets upset it), so I'll post something here instead.

Here's why c always takes integer values and more generally why (k + 1) * (k + 2) * ... * (k + n) is always divisible by factorial(n). There are three ways to show this.

Apologies for the crappy maths typesetting, I should really figure out how to do this nicely.

Proof by hand-waving.

First off, something hand-wavy. One of (k + 1)*(k + 2) is even so 2 divides into that product. One of (k + 1)*(k + 2)*(k + 3) is threeven (not a real word!) so 3 divides into that product. 4 is not so simple. One of (k + 1)*(k + 2)*(k + 3)*(k + 4) is fourven so 4 divides into that product - however maybe the fourven number was also the even number we used to cancel the 2 earlier on so it only has one 2 left for division but if that's the case then since one of the remaining 3 numbers is also even and you can get a 2 from that to make up the 4. 5 works just like 2 and 3. 6 contends with both 2 and 3 but you can resolve that too. This gets messy quickly but essentially, for any prime p, every time you come to the next multiple of p in the divisor you will have already added at least one multiple of p in the k+i product and by the time you hit a multiple of p² you will have already hit a corresponding one in the k+i product too, so you'll still have enough ps in the product overall to cancel all the ps in the divisor. Still a bit hand-wavy but it could be made rigorous, however there is another proof that's easily rigorous but less insightful.

Proof by counting prime factors.

There is a well known formula for finding the power of p in factorial(n). It's

sum(i=1..infinity, [n/pⁱ])

where [] means take the integer part, throwing away any fractions. So if p = 5, and n=30 then this is

[30/5] + [30/25} + [30/125] + [30/625] + ... = 6 + 1 + 0 + 0 + ...

It's all zeros after that, so the end result is 7. So 5⁷ divides factorial(30) and 7 is the highest power of 5 that divides it.

So rewriting the k + i product as factorial(k + n)/factorial(k) we now know that for any prime p, it contains exactly

sum(i=1..infinity, [(n+k)/pⁱ]) - sum(i=1..infinity, [(k)/pⁱ])

powers of p. If this is greater than

sum(i=1..infinity, [(k)/pⁱ])

then we know that factorial(n) divides factorial(k + n)/factorial(k).

Well, [x + y] >= [x] + [y] so

[(n+k)/pⁱ] >= [(n)/pⁱ] + [(k)/pⁱ]

now bring one term over to the left hand side to get

[(n+k)/pⁱ] - [(n)/pⁱ] >= [(k)/pⁱ]

Now sum over i=1..infinity and you get exactly what we needed to show. So for any prime p, there are at least as many powers of p in the k+i product as there are in factorial(n), all the ps in the divisor cancel out and the divisor divides in evenly.

Proof by "it just is".

The final and least insightful way is to note that factorial(n+k)/(factorial(k) * factorial(n)) is n+k choose n and so must be an integer

Code

Here's a python snippet that checks that the division is always even with an assert.

#! /usr/bin/python2.4 import sys def comb(n, k): # k = max(k, n - k) # would make it go faster c = 1 for i in range(k + 1, n + 1): c *= i assert c % (i - k) == 0 c /= i - k return c n = int(sys.argv[1]) for i in range(0, n + 1): print "comb(%d, %d) = %d" % (n, i, comb(n, i))

Give it an n and it will print out comb(n, i) for all i. The assert never fires for me.

As it loops, c gets the values comb(k + 1, 1), comb(k + 2, 2), ..., comb(n, n - k), all of which are integers (and remember that comb(n, n - k) = comb(n, k) so c ends up with the desired value.)

If your lanugage does tail-call optimisation, you could define

comb(n, 0) = 1

comb(n, i) = comb(n - 1, i - 1) * n / (n - k)

and be done. With memoisation you could avoid lots of repeated calculations.

Independent nonsense: negligent doctors

As usual, the Indo provides some useless statistics, this time on malpractice cases. What they tell us is:

A senior house officer -- a grade of junior doctor -- was most likely to be involved (74pc) as against 14pc for registrars and 8pc for consultants.

What they don't tell us is what percentage of patients are dealt with by each type of doctor. So if only 50% of patients were dealt with by junior doctors but they were involved in 74% of malpractice cases then that's a black mark against junior doctors and would be something worth knowing. Elsewhere in the article it's stated that they deal with "the majority" of cases. Without these other figures, the statistics are not just useless, they're possibly quite misleading.

in reference to: http://www.independent.ie/national-news/ae-blunders-responsible-for-one-in-seven-medical-claims-1951185.html (view on Google Sidewiki)

Hyperbinary numbers

Update: The proof is over here.

I've been following along the a discussion on The Math Less Travelled about hyperbinary numbers (here, here and here). I'm posting this here because it's a bit too big to fit in the margins of the other blog (although I have found no remarkable proof). Since I have very little time to work on this (and have already given myself a few nights with insufficient sleep!) I thought I'd do a brain-dump of what I did in case someone else can finish it. Apologies for the awful formatting, I don't post much maths so have no nice styles for it and no time to find them.

Quick summary

h(n) is the number of ways of writing n as a sum of powers of 2 where each power of 2 can only occur 0, 1 or 2 times. See the articles for more details but the key facts are (I'll refer to them by these numbers):

1. h(2n) = h(n) + h(n - 1)
2. h(2n + 1) = h(n)

From those you can get several relationships, like h(k * 2^n - 1) = h(k-1). Other interesting stuff is in those posts. The main question is how to find the inverse of h(n). That is, given k, what values of n have h(n) = k. It turns out that each even solution to that leads to an infinite number of odd solutions and all odd solutions can be traced back to an even solution. So the question then is what are the even solutions - the "primary occurrences"? My contribution so far has been to notice that there seem to be φ(k) even solutions for a given k although I am still stumped on how to prove that.

Another useful things is that h(2n + 1) = h(n) means that all the odd positions are just repeats of earlier values. Applying fact 1 to the right hand side of fact 2 gets you h(2n) = h(2n + 1) + h(2n - 1). So all the even positions are just the sum of the values that come immediately before and after.

Finally, h(2^n - 1) = 1 and h(2^n) = n + 1 and from that you can get h(2^n + 1) = n

New stuff

The last bits got me thinking about h(2^n + k). It's fairly easy to figure that out for more and more values of k. For h(2^n + 2) just apply fact 1 to get that as the sum of 2 things we already know:

h(2^n + 2) = h(2^(n-1) + 1) + h(2^(n-1)) = (n-1) + (n-1) + 1 # note this line for discussion below. = 2n - 1.

Looking at the table of values primary occurrences (here) you can see that all of the odd entries do indeed have 2^n + 2. In fact after a bit, it's the lowest value for each odd entry.

What about 1? 2*1 - 1 = 1 but if n = 1 the line commented above is not valid because h(2^n + 1) = n is only valid for n > 1.

To get a bit further, consider h(2n) = h(2n + 1) + h(2n - 1) from above. Moving things around a little that becomes h(2n + 1) = h(2n) - h(2n - 1). So we can get the next odd entry by subtracting the previous odd entry from the previous even entry. So h(2^n + 3) = n - 1. Since 2^n + 3 is odd this doesn't tell us anything about primary occurrences but it's still useful as all odd entries are duplicated later on.

Now we can get h(2^n + 4) = h(2^(n-1) + 2) + h(2^(n-1) + 1) and again we know what they are = (2n - 1) + n = 3n - 1. And low and behold, all the entries (except 2) in the primary occurrences table of the form 3n - 1, include have 2^n + 4 (5 -> 12, 8 -> 20, 11 -> 36, ...). What about 2? 2 = 3*1 -1 but if n=1 then the relation above breaks because the 2n - 1 rule is only valid for n > 1.

Continuing on in the same vein you get the first 100:

h(2^n + 0) = 1n + 1 h(2^n + 2) = 2n - 1 h(2^n + 4) = 3n - 4 h(2^n + 6) = 3n - 5 h(2^n + 8) = 4n - 9 h(2^n + 10) = 5n - 12 h(2^n + 12) = 5n - 13 h(2^n + 14) = 4n - 11 h(2^n + 16) = 5n - 16 h(2^n + 18) = 7n - 23 h(2^n + 20) = 8n - 27 h(2^n + 22) = 7n - 24 h(2^n + 24) = 7n - 25 h(2^n + 26) = 8n - 29 h(2^n + 28) = 7n - 26 h(2^n + 30) = 5n - 19 h(2^n + 32) = 6n - 25 h(2^n + 34) = 9n - 38 h(2^n + 36) = 11n - 47 h(2^n + 38) = 10n - 43 h(2^n + 40) = 11n - 48 h(2^n + 42) = 13n - 57 h(2^n + 44) = 12n - 53 h(2^n + 46) = 9n - 40 h(2^n + 48) = 9n - 41 h(2^n + 50) = 12n - 55 h(2^n + 52) = 13n - 60 h(2^n + 54) = 11n - 51 h(2^n + 56) = 10n - 47 h(2^n + 58) = 11n - 52 h(2^n + 60) = 9n - 43 h(2^n + 62) = 6n - 29 h(2^n + 64) = 7n - 36 h(2^n + 66) = 11n - 57 h(2^n + 68) = 14n - 73 h(2^n + 70) = 13n - 68 h(2^n + 72) = 15n - 79 h(2^n + 74) = 18n - 95 h(2^n + 76) = 17n - 90 h(2^n + 78) = 13n - 69 h(2^n + 80) = 14n - 75 h(2^n + 82) = 19n - 102 h(2^n + 84) = 21n - 113 h(2^n + 86) = 18n - 97 h(2^n + 88) = 17n - 92 h(2^n + 90) = 19n - 103 h(2^n + 92) = 16n - 87 h(2^n + 94) = 11n - 60 h(2^n + 96) = 11n - 61 h(2^n + 98) = 16n - 89 h(2^n + 100) = 19n - 106

With the 2^n + k only being valid for n such that 2^n > k.

Obviously calculating these by hand gets a bit boring. To write some code to do it I had to play with f(n, m) = h(2^n + m) and it breaks down into 2 cases for m even and odd.

For the odd case, let m = 2k + 1. So

f(n, m) = h(2^n + m) = h(2^n + 2k + 1) = h(2n + 1) = h(2^n + 2k) - h(2^n + 2k - 1)

For the even case, let m = 2k. So

f(n, m) = h(2^n + m) = h(2^n + 2k) = h(2(2^(n-1) + k)) = h(2^(n-1) + k) + h(2^(n-1) + k -1) (applying fact 1) = h(2^(n-1) + m/2) + h(2^(n-1) + m/2 -1) = f(n - 1, m/2) + f(n -1, m/2 - 1)

Turning that into python and adding a couple of base cases gives the following code (it tries to also cope with 2^n - m but there's a bug somewhere). #! /usr/bin/python2.4 import sys class ANPlusB(object): def __init__(self, a, b): self.a = a self.b = b def Sub(self, n_plus): return ANPlusB(self.a, self.b + self.a * n_plus) def Eval(self, n): return self.a * n + self.b def __str__(self): if self.b < 0: s = "-" b = - self.b else: s = "+" b = self.b return "% 4sn %s % 4s" % (self.a, s, b) def __add__(self, other): return ANPlusB(self.a + other.a, self.b + other.b) def __sub__(self, other): return ANPlusB(self.a - other.a, self.b - other.b) TABLE = {-2: ANPlusB(1, -1), -1: ANPlusB(0, 1), 0: ANPlusB(1, 1), 1: ANPlusB(1, 0), } biggest = 1 smallest = -1 def F(n_plus, k): """Find the a*n +b for 2^(n + n_plus) + k .""" # Note something is broken for negaitve n_plus, going below -10 # causes infinite recursion. global smallest, biggest if k not in TABLE: if k < smallest: val = _FSmaller(k) smallest = k else: val = _FBigger(k) biggest = k TABLE[k] = val val = TABLE[k] return val.Sub(n_plus) def _FBigger(k): if k % 2: return F(0, k - 1) - F(0, k - 2) else: return F(-1, k / 2) + F(-1, (k / 2) - 1) def _FSmaller(k): if k % 2: return F(0, k + 2) - F(0, k + 1) else: return F(-1, k / 2) + F(-1, (k / 2) - 1) def main(): min, max = [int(x) for x in sys.argv[1:]] for i in range(min, max+1, 2): print "h(2^n + % 4d) = %s" % (i, F(0, i)) if __name__ == "__main__": main()

You can run it as ./hyperbinary.py 0 100 or to get the first 50 even entries. Looking at the output it seems that h(2^n + m) = h(m) * n + c. I haven't proved it formally but it should be easy enough, just show that the coefficient of n obeys the fact 1 and fact 2. I didn't find a good relation between m or h(m) and the c above (the best I can spot is that int(c/h(m)) is always log_2(m) and when m is 2^k c=k*h(m) + 1). Nor did I relate any of this to φ although at least now there appears to be some modulo arithmetic going on.

Answer to 1 of the 2 puzzles and a new puzzle.

I take Nicole's comment to mean she gives up, although Nicole, it turns out that your comment has a deep relation to this problem and raises a new question.

To find out on average how many items are left unmoved by a random permutation, you can consider every permutation and how many are left unmoved by each, add them all up and divide by the total number of permutations. To avoid mathsy notations I'm going to write a tiny program that calculates this total.

total = 0 # Check what every permutation does to every item, if it leaves it # in the same place then increase the total by 1 for p in perms: for i in 1..N: if p(i) == i: total = total + 1

Now, here's the key to this problem - swap the loops. The final answer stays the same we just end up adding things up in a different order.

total = 0 # Check what every permutation does to every item, if it leaves it # in the same place then increase the total by 1 for i in 1..N: for p in perms: if p(i) == i: total = total + 1

So now what we're doing is looping over the items and counting how many permutations leave that item in place and that's something we can work out. The number of permutations that fix a given item is (N-1)! (where ! is factorial) because when you fix 1 item, you're left with N-1 other items which you move around however you like so it's just the number of permutations of N-1 items which is (N-1)!. So now the code can be just

total = 0 # Check what every permutation does to every item, if it leaves it # in the same place then increase the total by 1 for i in 1..N: total = total + (N-1)!

which is just adding (N-1)! N times. That's N x (N-1)! which is just N!. Now we have to divide that by the total number of permutations which is also N! to get an average of 1.

You can also get to this answer by drawing a grid with all the permutations across the top and all the items down the side. You put a tick in a cell in the grid if the permutation on the top leaves the item on the left unmoved. Now count the ticks. Doing it column by column is hard but doing it row by row is easy because there will be (N-1)! ticks in each row.

So the conclusion is that if you do a perfect shuffle of a deck of cards, on average 1 card will end up in the same places as before.

So where does Nicole's comment come in? I knew it was Ralph Wiggum who said it but I couldn't remember the context turns out. I found it here. It's

  Lisa: Hey Ralph, want to come with me and Alison to play "Anagrams"?
Alison: We take proper names and rearrange the letters to form a
        description of that person.
 Ralph: My cat's breath smells like cat food.

and what we've just shown is that given a random anagram, you'd expect 1 letter on average to end up in the same place - well kind of, that's only true if there are no repeated letters. I hadn't thought of this in terms of anagrams before and that adds a nice twist in terms of repeated letters. Thanks Nicole!

So the new puzzle is - given a string written in scrabble tiles (it can have repeated letters), if you mess up the tiles and make a new string at random, how many letters do you expect to be in the same place as before?

2 fun maths puzzles

I guess they're fun because I figured them out! If there's anyone actually reading this, post your solutions in the comments. If I get any responses I'll post solutions and if I don't, I might not.

1. On average, given a random permutation of N items, how many elements do we expect to stay put? E.g. if N=3 there are only 6 perumtations, 1 leaves all 3 in place, 3 leave 1 in place and 2 shift all of them (left or right). So that's (1*3 + 3*1) / 6 = 1.
2. On average, given a random function from 1..N->1..N which may or may not be injective or surjective (that is there may be X and Y such that f(X) = F(Y) and there maybe be a Z such that f(X) != Z for any X), how many elements do we expect to stay put - f(X)=X? E.g. if N=2 there are 4 possible functions, f(X) = X, f(X) = 1, f(X) = 2 and f(X) = 3-X. These fix 2, 1, 1 and 0 elements respectively. So on average that's (2+1+1)/4 = 1 but what about larger N?

Proof by taxation

A little story

One day your boss says, "Great job, we want to give you a bonus, were you planning a holiday this year?". You say "well actually, I wanted to go here", you show him a website, "but I can't afford it". "How much is it?". "1000 euro", you say. "Great, we'll pay for that". "Thanks boss!". Happy days.

Later you realise that if they just give you 1000e you'll have to pay tax on that. The tax rate is 40% so that's a 400e tax bill, you still can't afford the holiday.

You head up to your boss's office and explain the problem. He's still in a great mood and says, "Don't worry, we'll pay your tax bill too". "Wow but there's another bit of a problem, I'll have to pay tax on that extra 400 too. 40% of 400e is 160e, I still can't afford it". You boss is looking less happy now and gets out a piece of paper. "Right, we'll pay all your tax, no matter what" and starts writing down some figures 1000 400 160 64 25.60 10.24 4.096 1.6384 0.65536 you hear him mumbling as he starts adding them all up. Suddenly Joan, his secretary, who's been quiet all this time blurts out "1666e and 66c". You both stare at her in amazement then after a lot more mumbling the boss says "my total is 1666e and .23c but I suppose if I added a few more lines to the sum it'd probably be 1666.66. How did you get it Joan?"

"Well, tax is 40% so he gets to keep 60% of anything you pay him. So the final amount in his pocket is the what you give him times 0.6 . So you're looking for a number that when you multiply it by 0.6 gives you 1000. So 1000 ÷ .6 is the number you're looking for because when you multiply that by 0.6 the two .6s cancel out and only the 1000 is left and 1000 ÷ 0.6 is 1666.666666666...", says Joan.

So what we have proved is that 1000 + 1000 x 0.4 + 1000 x 0.4^2 + 1000 x 0.4^3 + ... = 1000 ÷ 0.6 . There was nothing special about our choice of 40%, so replacing 0.4 with r (and 0.6 by 1-r) and dividing out the 1000 on both sides gives 1 + r + r^2 + r^3 + ... = 1 ÷ (1 - r) Of course there are other ways to prove this but, I like my proof by taxation because it feels like it explains why they are equal.

Listing the rational numbers

Here's a nice series of posts about how to list every rational number in its lowest terms, with no repeats (ie. only include 2/5 not 4/10). The key is to make a tree with 1/1 as the root and for each node i/j, put i+j/j as the left child and i/i+j as the right child. It's all laid out very clearly in the series of post which is based on a terse maths paper.

In the comments I found a link to another paper which comes at this using matrices and relates this tree with another previously known tree that also lists the rationals. It's a fairly readable paper too but a bit more technical.