Recently I’ve been learning to use Mathematica, a piece of software I was curious about for a long time. Luckily my university has a licence for all staff, so I snagged a key as soon as I learned this and have been mucking around with some serious stuff — namely it’s surprisingly good neural network features — and some less serious stuff to get my head around the Wolfram Language.
This weekend I’ve been messing around with FRACTRAN, a fascinating esoteric programming language and model of computation developed by John H Conway, who’s perhaps most famous amongst computer types for inventing the Game of Life cellular automaton. FRACTRAN is a language in which programs consist entirely of lists of fractions, like so:
primeProg = {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};
Now that may seem nonsensical, but bear with me here. FRACTRAN works with really simple rules; given a list of fractions and an initial input N:
- Find the first fraction F in the program listing which becomes an integer when multiplied by N, then replace N by N*F
- Repeat until N doesn’t produce any integers when multiplied by any fraction in the program, then halt.
Easy peasy! So we can write a FRACTRAN interpreter in Mathematica quite easily. This one outputs each stage of program execution in order to a list called outputFrac, to allow us to manipulate the results later if need be:
fracRunList[fracProg_, input_, steplimit_] := Module[{j, state},
j = 0;
state = input;
outputFrac = {};
While[j <= steplimit, newProg = state *fracProg;
integerList = IntegerQ[#] & /@ newProg;
intSpots = Position[integerList, True];
AppendTo[outputFrac, state];
If[Length[intSpots] == 0, Break[]];
state = newProg[[intSpots[[1, 1]]]]; j++]];
So, when you call this function with fracRunList[{list of fractions}, N, timestep limit], it multiplies N through the list, checks that new list for integer values, appends that value to the list outputFrac, then starts again. The function will halt either when it reaches the timestep limit you specified, or when no more integers result from multiplying N through the list.
When we run the program above — suspiciously called ‘primeProg’ — with an initial N=2 for 50 steps, we get this:
{2, 15, 825, 725, 1925, 2275, 425, 390, 330, 290, 770, 910, 170, 156,
132, 116, 308, 364, 68, 4, 30, 225, 12375, 10875, 28875, 25375,
67375, 79625, 14875, 13650, 2550, 2340, 1980, 1740, 4620, 4060,
10780, 12740, 2380, 2184, 408, 152, 92, 380, 230, 950, 575, 2375,
9625, 11375, 2125, 1950, 1650, 1450, 3850, 4550, 850, 780, 660, 580,
1540, 1820, 340, 312, 264, 232, 616, 728, 136, 8, 60}
That may look like nonsense, but note that scattered through that list of numbers we have 2, 4, and 8 — which are respectively 2^1, 2^2, and 2^3. So what PrimeProg does is actually output all the prime numbers, in the form of prime exponents of 2!
We can see this easily if we run a simple filter on the list outputFrac after running the program for 50,000 steps:
findPrimes2[list_] :=
Log[2, Select[list, IntegerQ[Log[2, #]] &]];
fracRunList[primeProg, 2, 50000]
findPrimes2[outputFrac]
Output: {1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31}
After 50,000 steps, our clever little list of fractions has produced the first 13 prime numbers! Theoretically we can run this forever, and produce every prime number. Although each one takes longer to come out than the last, unsurprisingly, so I don’t recommend it.
Another variation of the prime-finder program I found from the Esolang wiki is more efficient, using only 9 fractions to output prime exponents of 10. We’ll test it out below, this time filtering the resulting output list for prime exponents of 10:
prime10Short = {3/11, 847/45, 143/6, 7/3, 10/91, 3/7, 36/325, 1/2,
36/5};
findPrimes10[list_] :=
Log[10, Select[list, IntegerQ[Log[10, #]] &]];
fracRunList[prime10Short, 10, 50000]
findPrimes10[outputFrac]
Output: {1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}
This prime-finder has managed to dig up 15 primes in 50,000 steps, rather than 13 like the original.
The really remarkable thing about FRACTRAN, though, is that it’s actually Turing-complete — it can in principle calculate anything calculable by any other programming language. A simple example of a multiplication programme shows off how this works:
multiFrac = {455/33, 11/13, 1/11, 3/7, 11/2, 1/3};
fracRunList[multiFrac, 72, 50]
Output: {72, 396, 5460, 4620, 63700, 53900, 4900, 2100, 900, 4950, 68250,
57750, 796250, 673750, 61250, 26250, 11250, 61875, 853125, 721875,
9953125, 8421875, 765625, 328125, 140625, 46875, 15625}
What’s happened here is that we gave our FRACTRAN program a single number that actually represents our two input numbers — in this case 3 and 2 — as a product of prime numbers raised to appropriate powers — 2^3 * 3^2 = 8 * 9 = 72. FRACTRAN then outputs the result as the power of a different prime — 15,625 = 5^6 = 5^(2 * 3). A fairly roundabout way to get multiplication done, but it works!
What this shows is that it’s possible to have FRACTRAN programs operate on multiple inputs, so long as those inputs are encoded as products of prime powers. In fact, we can assign specific primes to be our data registers and use carefully-constructed fractions to operate on those registers, and even construct complicated programs with loops! As this StackOverflow answer shows, it turns out FRACTRAN is exactly equivalent to a Minsky Register Machine, which have been proven to be Turing-equivalent — hence confirming that FRACTRAN is actually a Turing-complete language.
As a consequence some intrepid folk have built some impressive constructs in FRACTRAN. One of my favourites is a FRACTRAN interpreter which is itself written in FRACTRAN! Using just 48 fractions, this program takes as input an encoded FRACTRAN program and initial state, and correctly interprets the program and outputs the result. Here’s everything you need to try it:
fracInFrac = {5/19, 1558654261983398483185/122130132904968017083,
185/1822837804551761449, 4996917562403854655/41, 272365/67, 43/5,
43/71, 125173/47, 145915005923554298917151/952809757913927,
950886101246622507133/41426511213649, 160585150715989139597/13,
8752951/23, 17/43, 17/29, 6409/47, 5/17, 31/53, 17042839/7,
1829/41, 59/73, 331639/23, 4307/41, 89/59, 3713/31, 79/83,
268837/23, 31/79, 8633/7, 101/97, 68579/11, 9797/13, 9797/47,
35/101, 9167/13, 103/107, 1774381/47, 109/103, 109/113, 578899/23,
11227/13, 127/109, 127/131, 16637/47, 16637/11, 1114679/61, 2/127,
5/2, 3/37};
(* Encode a FRACTRAN program as a base-11 number for input into
fracInFrac interpreter *)
base2 = 11;
pad2[f_] :=
Block[{n = IntegerDigits[Numerator[f], base2 - 1],
d = IntegerDigits[Denominator[f], base2 - 1], len},
len = Max[Length[n], Length[d]]; n = PadLeft[n, len];
d = PadLeft[d, len]; Flatten[{0, Riffle[n, d], base2 - 1}]];
digits2[progList_] := Join[Flatten[pad2 /@ progList], {base2 - 1}];
encode2[progList_] := FromDigits[Reverse[digits2[progList]], base2];
(* FracInFrac input state encoder function *)
fracInput[fracProg_, init_] := 5*7^init*67^(encode2[fracProg]);
In order to encode our FRACTRAN program into a format the interpreter can understand, we first need to encode the program as a single number. In this case the encode2 function reverses the list of fractions in the program, then encodes the digits as base-10 numbers within a base-11 number. Then we need to combine that with our initial state into a single number that we can pass to the interpreter. We do this using the fracInput function, which gives us a ridiculous huge number that consists of 5 * 7^(initial state) * 67^(encoded program). In fact, the resulting numbers are far too huge to print here as examples even for the simple adding program (which is just {3/2}!), and Mathematica can’t even cope with encoding larger programs and simply spits out an error. Changing the encoding to use a smaller prime for the program is possible, but I leave that as an exercise for the reader.
Another intrepid StackOverflow commenter produced a FRACTRAN interpreter for FRACTRAN in 84 fractions, which has a slightly less ginormous program encoding:
fracInFrac3 = {197*103/(2^11*101), 101/103, 103*127/(2*101), 101/103,
109/101, 2*23/(197*109), 109/23, 29/109, 197*41*47/(31*59),
11^10*53/(127*197), 197/53, 37/197, 7^10*43/(11^10*37), 37/43,
59/(37*47), 59/47, 41*61/59, 31*67/(41*61), 61/67, 7*67/(127*61),
61/67, 101/71, 73/(127^9*29), 79/(127^2*73), 83/(127*73),
89/(2*29), 163/29, 127^11*89/79, 337/83, 2*59/89, 71/61,
7*173/(127*163), 163/173, 337*167/163, 347/(31*337), 337/347,
151/337, 1/71, 19*179/(3*7*193), 193/179, 157/(7*193), 17*181/193,
7*211/(19*181), 181/211, 193/181, 157/193, 223/(7*157), 157/223,
281*283/239, 3*257*269/(7*241), 241/269, 263/241, 7*271/(257*263),
263/271, 281/263, 241/(17*281), 1/281, 307/(7*283), 283/307,
293/283, 71*131/107, 193/(131*151), 227/(19*157), 71*311/227,
233/(151*167*311), 151*311/229, 7*317/(19*229), 229/317,
239*331/217, 71*313/157, 239*251/(151*167*313), 239*251/(151*313),
149/(251*293), 107/(293*331), 137/199,
2^100*13^100*353/(5^100*137), 2*13*353/(5*137), 137/353, 349/137,
107/349, 5^100*359/(13^100*149), 5*359/(13*149), 149/359, 199/149};
To encode the program, we reverse the order of the fractions in the program, put ’10’ between each numerator and denominator, and encode the whole list as a base-11 number:
(* Encode a FRACTRAN program in reversed base-11 digits *)
baseConvert[frac_] :=
{Numerator[frac], 10, Denominator[frac], 10};
baseDigits[fracList_] :=
Most[Flatten[Reverse[Join[baseConvert /@ fracList]]]];
baseEncode[fracCode_] := FromDigits[baseDigits[fracCode], 11];
(* Sample encoding: Simple adding program {3/2} *)
encodedAdder = baseEncode[addProg]
Output: 475
That seems manageable, right? But then we have to encode the initial state 72, as well, and this program uses the format (3^(state) * 5^(encoded program) * 199)…
(* Complete encoded initial state for FRACTRAN interpreter *)
(* 3^initial_state*5^encoded_program*199 *)
encodeFracProgState[fracList_, init_] := (3^init)*(5^baseEncode[fracList])*199;
(* Sample complete encoding: Adder, initial state 72 *)
encodeFracProgState[addProg, 72]
Output: 4595528627302514457847822534456305637274485006848124607416562426715142
2223431883566346800946031635723355137742896740401730524435178459281271
4997262424758884123382086634982376119146502282302917806400755768226744
1829230852502546067970809133366591936071662069537535458381098407899190
9230824160175285935800563682805991749525746897352823649995912091981153
9351940155029296875
…Yikes. Well, because I love you guys so much, I tested this out. I ran the encoded program through the 84-fraction interpreter, and piped the output to a text file which rapidly blew up to 4.8MB of numbers. The correct answer pops out after 6,030 lines of gibberish and looks like this:
1331690264838856002293720794993649380710435784303038760424759391768076 4462808226237499620786133558911763356798565240461593352120399425752403 5496560518676882088942416619737846495343367625845408913052935028877605 1678370577471585522503962029887317190284834394621170248161337171206903 0172712146158755426661553931885800712601024170227364957005027071924897 2557777777294510307702117543602421742701115006472526481890849420208766 6735678794793784618377685546875
Which is actually the expanded form of:
3^243 * 5^475 * 149
And you’ll note that 3^243, which is also 3^(3^(3+2)), giving us the answer to our original request to add 3 and 2, buried up in the second layer of exponents. Whew!
Anyway, as you can see it’s easy to get lost in FRACTRAN despite its apparent simplicity. It’s really interesting to play with, though, and it’s an odd moment when you realise that these simple lists of fractions are actually capable of some remarkable things. The weekend is nearly over now and I have to prepare for actual work, but perhaps next weekend I’ll return to this and construct a compiler for FRACTRAN, making it possible to write programs in a higher-level language and squish them into FRACTRAN form.
Or I might start fooling around with different nerd stuff instead, who knows?
Dr Eric Silverman 