Daily Archives: 10/27/2006

Prime Number Pathology Courtesy of Good Math, Bad Math

I look forward to reading the Good Math, Bad Math blog, especially on Fridays. Often there is a pointer to an interesting computer language that somehow stretches your understanding of programming languages and computer science. Today was no exception, as he introduced Fractran, and interesting “language” that does computation by multiplying fractions.

It turns out that it isn’t the first time I’d seen this: Conway mentioned it in his Book of Numbers, a terrific little book full of gems of knowledge like this. He presents it just as an aside, by giving you the following list of fractions, labelled A-N,

17/91, 78/85, 19/51, 23/38, 29/33, 77/29, 95/23, 
77/19, 1/17, 11/13, 13/11, 15/14, 15/2, 55/1

You begin by initializing an accumulator to 2. You then proceed by trying to multiply your accumulator by each number. If it results in a whole number, you
update your accumulator, and start again. Every once in a while, this computation generates a result which is an exact power of two. The sequence of
these numbers forms a monotonically increasing series, with only the prime powers of two. I wrote a little Python program to demonstrate this.

#!/usr/bin/python

program = [(17,91), (78,85), (19,51), (23,38), (29,33), (77,29),
           (95,23), (77,19), (1,17), (11,13), (13,11), (15,14),
           (15,2), (55,1) ]

acc = 2

def poweroftwo(acc):
        while acc > 1 and acc % 2 == 0:
                acc = acc // 2
        return acc == 1

def log2(acc):
        cnt = 0
        while acc>1:
                acc = acc // 2
                cnt = cnt + 1
        return cnt

while True:
        for n, d in program:
                if acc % d == 0:
                        break
        acc = n * acc / d
        if poweroftwo(acc):
                print "2^%d = %d" % (log2(acc), acc)

When you run this, you get output like….

2^2 = 4
2^3 = 8
2^5 = 32
2^7 = 128
2^11 = 2048
2^13 = 8192
2^17 = 131072
2^19 = 524288
2^23 = 8388608
2^29 = 536870912
2^31 = 2147483648
2^37 = 137438953472
2^41 = 2199023255552
2^43 = 8796093022208
...

When I first saw this, I must admit that I was mystified, but it’s really not that hard (although far from obvious) to figure out what is going on. Your accumulator is literally encoding the state of a computer. In the case of the calculation above, the “machine” has 10 registers, each one labelled with a prime (basically the
primes less than thirty). If you wish to increment the value in a particular register, you multiply the accumulator by the corresponding prime. If you wish to decrement a counter, you divide. You aren’t allowed to divide if the results aren’t integers, which basically means that each register contains a positive non-negative value. The state of the machine is just the product of all the machine values. Pretty neat!

Using this idea, you can convert the list of fractions given by Conway into a pretty neat little .signature program to compute primes (and runs a lot faster):

int j=1,h,c,o,n,w,a,y,m,v;main(){for(;;)o&&w?a++,o--,w--:(c&&a?j++,h++,w++,c--,
a--:(h&&a?y++,h--,a--:(j&&y?m++,j--,y--:(h&&n?v++,h--,n--:(v?o++,n++,v--:(m?c++
,y++,m--:(y?o++,n++,y--:(a?a--:(w?n++,w--:(n?w++,n--:(j&&o?h++,c++,j--,o--:(j?h
++,c++,j--:(c++,n++))))))))))))),h+c+o+n+w+a+y+m+v||printf("%d\\n",j);}

Not sure about the character mangling, if you have difficulty cutting and pasting this, try clicking this link to get the text.

How’s that for obfuscation? Thanks to Tom for helping me compress this so it will fit down into 4 lines.

Addendum: Here’s an even shorter version. Tom wrote a slightly shorter one by using recursion and avoiding the main loop, but I think that’s vaguely cheating, as C compilers don’t necessarily implement tail recursion as simple loops.

Technorati Tags: , ,