Conway's FRACTRAN

Question

I've implemented fractran and playing with sequences iteratively in C.

http://mathworld.wolfram.com/FRACTRAN.html

Sloan: http://oeis.org/A007542

I can generate correct results, but I would like to so do in time measured in hours not weeks... I have a 2 cpu Ultrasparc III machine (zone). pthreads improves throughput slightly.

I have to iterate though each "fraction" to determine the first successful operation. code snippet -

for(j=0; j< 14 ;j++)
{
  mpz_mul(tmp, elem, num[j]);
  mpz_tdiv_qr(q, r, tmp, denom[j]);
  if ( mpz_cmp(r, zero) == 0)
  {
         mpz_set(elem, q);
         break;
      }
}

One gmp multiplication, one division, one compare per "fraction". Yes, this is brute force. Suggestions?

Answer 1

If you've profiled it, you'll find that the problem is division. That's the really slow operation that you'd like to avoid, but avoiding multiplication wouldn't hurt either.

Real FRACTRAN programs (yes, I'm going to use the phrase "real FRACTRAN programs" as if there is such a thing) tend to work at a deep level in terms of prime factorization. For example, consider this program (written by me; you're welcome) which implements GCD:

2/7 3/11 13/17 19/23 51/65 1/13 46/95 1/19 5/6 91/2 209/3

You give it 2^p 3^q and it returns 5^gcd(p,q).

This suggests that a better way of implementing FRACTRAN is to work with completely factorized numbers. Rather than 46/95 (an instruction in the above program), you would represent this as an array of integers indexed by prime numbers:

#  2,  3,  5,  7, 11, 13, 17, 19, 23
{  1,  0, -1,  0,  0,  0,  0, -1,  1 }
# Verify that 2^1 5^-1 19^-1 23^1 = 46/95

The machine state would similarly be an array indexed by prime numbers. To execute an instruction, you make sure that performing the subtractions would not make the state underflow, and if that checks out, simply add the two arrays together elementwise.

There are two outstanding issues to consider.

Firstly, how long should the arrays be? That is, how many primes do you need? Thankfully that can be determined statically from the instruction stream. If there are primes in the input which are not mentioned in the program, they can just be copied to the output. So pick the array size for the program, and you should be fine there.

The other problem is how large the capacity for each array entry should be. For the instructions, once again, that can be determined statically. For the machine state, however, it's not clear. Of course, 64 bits ought to be enough for anybody (famous last words), but if you're trying to implement FRACTRAN in FRACTRAN, it may not be. Complicating things is that real FRACTRAN programs tend to use a lot of flags (i.e. primes whose power are always 0 or 1 in the machine state) as pseudo program counters, so using an MPZ for each entry would probably be wasteful.

So for the machine state, I'd consider a hybrid. If a prime power will fit in a machine word, use a machine word. Otherwise, use a MPZ or (since you only need addition and subtraction) a home-grown equivalent. Bonus points if you can use SIMD instructions.

Good luck, and please keep me posted. I'd love to see this on github.