Fractran in Haskell

One of the blogs I tend to follow is Good Math, Bad Math which had a post on one of those oddball things that can be lots of fun to play with, in this case Fractran, a Turing Complete language by John Horton Conway. I implemented this myself in Haskell back in (hmmm, not at all sure) 1997(?) in
Haskell. The quality of the Haskell novice I was then (and remain) shows though my Haskell has improved since then. Read that post for better information on what Fractran is and how it works.

This implementation factors the numbers used and just keeps track of the factors in the numerator and denominator of the fractions involved. I probably reinvented the wheel to get the primes and factors and all. I tested it in the primes program (named primegame) and
the addition program (named addergame), but not on much else as actually writing Fractran programs was not something I tried very hard to master.

import Ratio
--
-- run takes an integer i
-- and a list of fractions (the program)
-- it returns a result list of the integers generated
--

run :: Integer -> [Rational] -> [Integer]
run i p =      takeWhile (>0) res
        where
      res = i:(map runstep res)
      runstep j        = runl j p

runl :: Integer -> [Rational] -> Integer
runl    j []     = 0
runl    j (f:fs) = let
                     ifv = (fromInteger j)*f
                  in  if   isInteger1 ifv
                      then (numerator ifv)  -- essentially a toInteger
                      else runl j fs

isInteger i = (i == (truncate i))
isInteger1 i = ((denominator i) == 1) -- works in this context

--
-- the numbers themselves dont show much, so write them as products
-- of powers of primes
--
primes      :: Integral a => [a]
primes       = map head (iterate sieve [2..])
sieve (p:xs) = [ x | x<-xs, x `rem` p /= 0 ] 
powersOfTwo = 2:(map (2*) powersOfTwo)  
--
-- returns 0 if not a power of two
-- else returns the power 
-- 
whichPowerOfTwo x =     l                     
 where                  
       (a,b) = span (< x) powersOfTwo
       l = x == (head b) then 1 + (length a) else 0

--
-- brute force factorization
--

factor  x   = factor1 x primes 
factor1 1 _      = []
factor1 x (p:ps) = let (c,q) = fp x p
                  in  if c > 0
                      then (c,p):(factor1 q ps)
                      else factor1 q ps

--
-- multiply two lists of prime, power pairs
--

mult [] l = l
mult l [] = l
mult l@((p,pow):ps) l'@((p',pow'):ps')
 |  p == p'              = (p,pow+pow'):(mult ps ps')
 |  p < p'             = (p,pow):(mult ps l')
 |  p >  p'              = (p',pow):(mult l ps')

--
-- divide one list of prime,power pairs by another
--

divvy []  l = []
divvy l  [] = l
divvy l@((p,pow):ps) l'@((p',pow'):ps')
| p == p' && pow == pow'       = (divvy ps ps')
| p == p'                      = (p, pow - pow'):(divvy ps ps')
| p <  p'                      = (p,pow):(divvy ps l')
| p >  p'                      = (p, -pow'):(divvy ps ps')

--
-- isIntegerL returns true if the list represents an integer
-- when all the powers are >= 0
--

isIntegerL l = and (map ((>0).snd) l)

--
-- fp takes two integers and returns the largest
-- power of the second that evenly divides the first
-- and the quotient thus determined
--

fp :: Integer -> Integer -> (Integer, Integer)
fp x p     = if   x `mod` p == 0
               then let (c,q) = fp (x `div` p) p
                    in  (c+1,q)
               else (0, x)

primegame = [17%91, 78%85, 19%51, 23%38, 29%33, 77%29, 95%23, 77%19, 1%17, 11%13, 13%11, 15%2, 1%7,5
5%1 ]

--
-- start with 2^a * 3^b
-- ends with 2^a+b
-- so started with (8=2^3) * (243 = 3^5) = 1944 should result in 256 = 2^8
-- that is, run 1944 addergame => 256
--
addergame = [2%3]

p1 = filter (/= 0) (map whichPowerOfTwo (run 2 primegame))

-- main = putStr (show (take 20 p1))
main = putStr (unlines (map show (map factor (take 10000 (run 2 primegame)))))