{-
goodstein.hs
Ed Wells
Version 3
July 10, 2006
This module contains the code necessary to produce the Goodstein's sequence
starting with a natural number. For more information please see the URL:
http://astronomy.swin.edu.au/~pbourke/fun/goodstein/
The main routine in this module is capable of producing the two sequences
listed at that URL, which start with the numbers 25 and 4.
The method used is to create a data structure representing the series of
power expressions for each term in the Goodstein sequence. A simple data
structure, Expr, is used to contain Numbers, Addition of elements in a term,
and the Power elements for each term.
To reduce the complexity of the data structure and a dramatic increase in the
number of rules necessary for pattern matching, general multiplication and
subtraction are disallowed.
In the Goodstein sequence we only ever subtract by one, and only when generating
the nex term in the sequence. This simplified subtraction is handled by the
function sub and called in the function nextTerm.
Additionally, the only form of multiplication in the Goodstein sequence is the
coefficient of a power element. For this reason, the power elements contain the
coefficient as part of its constructor.
The process of subtraction to obtain the next term is handled by a recurssive match
against Num, Add, and Pow elements. Overall, the data structure forms a tree structure
in which Num elements are leaf elements, and Add and Pow elements are branching nodes.
Pattern matching in the sub function will perform a depth-first walk of the tree when
an Add element is found.
We are careful to always construct the tree such that the lowest-powered element is
encountered first by depth-first traversal. This element may be a scalar Num or
a more complex Pow element. The recursive rules for handling subtraction from a
Pow element are chosen to be minimal, with the exception of Rule 2, which is so common
in practice to warrant expansion from Rule 3.
Running the main function produces a simple Read-Evaluate-Print-Loop interpreter
in which the user can input a number and produce subsequent terms in the number's
Goodstein sequence.
The function integerLogBase is originally defined in the Numerics section
of the Haskell 98 Library Report. It is viewable for the following URL.
http://www.informatik.uni-freiburg.de/~thiemann/haskell/haskell98-library-html/numeric.html#sect4
The function formatExpr is inspired by a similar function in online course
"Programming in Haskell" at Chalmers, viewable at the following URL.
http://www.cs.chalmers.se/Cs/Grundutb/Kurser/d1pt/d1pta/Defining%20New%20Types/sld017.htm
All other code is authored by myself.
Compiles using GHC 6.4.2
GHCi works.
To use as a module, comment the function main and uncomment the module statement.
-}
{-
module Goodstein (
Expr, thrd, eval, add, pow, mutate, sub,
nextTerm, getTerm, getTermF, putExpr,
intToExpr2, exprToInteger, formatExpr,
integerLogBase, readInteger, main)
where
-}
data Expr = Num Integer
| Add Arg2 -- Function and two expressions.
| Pow Arg3
deriving Show
type Arg2 = (Expr, Expr)
type Arg3 = (Expr, Expr, Expr)
-- Simple function to extract third element from a triple.
thrd :: (a, b, c) -> c
thrd (_, _, c) = c
evalArgs :: Arg2 -> Arg2
evalArgs (e1,e2) = (eval e1, eval e2)
evalArgs3 :: Arg3 -> Arg3
evalArgs3 (e1, e2, e3) = (eval e1, eval e2, eval e3)
-- Walking the data structure.
eval :: Expr -> Expr
eval (Num a) = Num a
eval (Add as) = add as
eval (Pow as) = pow as
{-
The following arithmetic operations are carefully chosen so that only the minimal
number of operations are possible when constructing and mutating our expressions.
For instance, a term like "a-b" is not allowed because we only subtract with
each mutation and that can be handled using "sub a", where a :: Expr.
We are careful to construct the expression such that our power terms are nested from left to right
(outer to inner) as least to greatest.
Also, in general, we do not want Add (Num a, Num b) -> Num (a+b). The only place we have such
a pattern should be in the exponents of the powers, which we do NOT want to add together.
-}
-- Simplify and evaluate additions.
add :: Arg2 -> Expr
add args = let (e1,e2) = evalArgs args
in case (e1,e2) of
( a, Num 0) -> a
(Num 0, b) -> b
(Num a, _) -> Add (e1, e2)
( _, Num b) -> Add (e2, e1) -- Push numbers to left.
_ -> Add (e1, e2)
-- Simplify and evaluate powers.
pow :: Arg3 -> Expr
pow args = let eArgs = evalArgs3 args -- (coef, base, exp)
in case eArgs of
(Num 0, _, _) -> Num 0
( _, Num 0, _) -> Num 0
( coef, _, Num 0) -> coef
-- ( coef, Pow (Num 1, base, base'), exp)
-- -> let upperPow = pow (Num 1, base', exp)
-- in Pow (coef, base, upperPow) -- Note: recursive evocation of pow.
( coef, base, Num 1) -> Pow eArgs
_ -> Pow eArgs
mutate :: Expr -> Expr
mutate (Num a) = Num a
mutate (Add (a, b)) = add (mutate a, mutate b)
mutate (Pow (coef, Num b, exp))
= let b' = b + 1
coef' = mutate' (Num b') coef
expr' = pow (coef', Num b, exp)
in mutate' (Num b') expr'
mutate (Pow (coef, Pow (c, Num b, e), exp))
= let b' = b + 1
coef' = mutate' (Num b') coef
c' = mutate' (Num b') c
base' = mutate' (Num b') (pow (c', Num b, e))
expr' = pow (coef', base', exp)
in mutate' (Num b') expr'
mutate expr = expr -- An error term could go here.
mutate' :: Expr -> Expr -> Expr
mutate' new (Add (e1, e2)) = add (mutate' new e1, mutate' new e2)
mutate' new (Pow (coef, old, Num 1)) = Pow (coef, new, Num 1)
mutate' new (Pow (coef, old, Num e)) = Pow (coef, new, Num e)
mutate' new (Pow (coef, old, Add (Num e, expr))) = let newExpr = mutate' new expr
in Pow (coef, new, Add (Num e, newExpr))
mutate' new (Pow (coef, old, expr)) = let newExpr = mutate' new expr
in Pow (coef, new, newExpr)
mutate' new expr = expr -- An error term could go here.
-- Subtract one from an expression. This is only used when mutating
-- an expression to the next term in the Goodstein sequence.
sub :: Expr -> Expr
sub arg = case arg of
-- Rule 0. Pull subtractions into first element of
-- of addition and evaluate.
Add as -> let expr = sub (fst as)
in add (expr, snd as)
-- Rule 1. Subtract one from a number.
Num a -> Num (a-1)
-- Rule 2. Optimiztion to Rule 3.
-- B^1 - 1 => [ B - 1 ]
Pow (Num 1, base, Num 1)
-> sub base
-- Rule 3. Power to a number.
-- B^e - 1 => [ B^(e-1) - 1 ] + [ B - 1 ] * B^(e-1)
Pow (Num 1, base, Num e)
-> let e' = (sub (Num e))
coefPow = sub base
rightPow = pow (coefPow, base, e')
leftPow = sub (pow (Num 1, base, e'))
in add (leftPow, rightPow)
-- Rule 4. Power to number plus Expression.
-- B^(e+E) - 1 => [B^((e-1)+E) - 1] + [B - 1] * B^((e-1)+E)
Pow (Num 1, base, Add (Num e, expr))
-> let e' = sub (Num e)
coefPow = sub base
addPow = add (e', expr)
rightPow = pow (coefPow, base, addPow)
leftPow = sub (pow (Num 1, base, addPow))
in add (leftPow, rightPow)
-- Rule 5. Power to Power Expression.
-- B^E - 1 => Recurssive call to sub based on Expression.
Pow (Num 1, base, Pow (Num 1, base', exp))
-> let base'' = Pow (Num 1, base, base')
subPow = sub (Pow (Num 1, base'', exp))
in eval subPow
{- -> let base'' = Pow (Num 1, base, base') -- base^base' -- Do not evaluate with pow.
in case exp of
-- Expression is Power to number.
-- b^b'^ne => (b^b')^ne
ne@(Num e)
-> pow (sub (Pow (Num 1, base'', ne)))
-- Expression is Power to number plus Expression.
-- b^b'^an => (b^b')^an
an@(Add (Num e, expr))
-> pow (sub (Pow (Num 1, base'', an)))
-- Expression is Power to Expression.
_ -> pow (sub (Pow (Num 1, base'', exp)))
-}
-- Rule 6. Coefficient times Power to Expression.
-- c*B^E - 1 => [B^E - 1] + (c-1)*B^E
Pow (coef, base, exp)
-> let coef' = sub coef
rightPow = pow (coef', base, exp)
leftPow = sub (pow (Num 1, base, exp))
in add (leftPow, rightPow)
nextTerm = eval.sub.mutate
putExpr :: String -> Expr -> IO ()
putExpr name expr = do
putStrLn (concat ["\n", name, " is: \n"])
putStrLn (formatExpr expr)
putStrLn (concat ["\n", "Mutate ", name, ": \n"])
putStrLn (formatExpr (nextTerm expr))
putStrLn "\n"
getTerm :: Integer -> Integer -> Expr
getTerm i n = getTerm' n (intToExpr2 i)
getTermF i n = formatExpr (getTerm i n)
getTerm' :: Integer -> Expr -> Expr
getTerm' n expr | n == 0 = expr
| n > 0 = getTerm' (n-1) (nextTerm expr)
formatExpr :: Expr -> String
formatExpr (Num a)
= show a
formatExpr (Add (e1, e2))
= concat [formatExpr e1, " + ", formatExpr e2]
formatExpr (Pow (Num 1, b, Num e))
= concat [formatExpr b, "^", show e]
formatExpr (Pow (Num 1, b, e))
= concat [formatExpr b, "^(", formatExpr e, ")"]
formatExpr (Pow (Num c, b, Num e))
= concat ["(", show c, "*", formatExpr b, "^", show e, ")"]
formatExpr (Pow (c, b, Num e))
= concat ["(", "(", formatExpr c, ")", "*", formatExpr b, "^", show e, ")"]
formatExpr (Pow (Num c, b, e))
= concat ["(", show c, "*", formatExpr b, "^(", formatExpr e, ")", ")"]
formatExpr (Pow (c, b, e))
= concat ["(", "(", formatExpr c, ")", "*", formatExpr b, "^(", formatExpr e, ")", ")"]
-- The following definition for integerLogBase
-- and all notes are from the Numerics section
-- of the Haskell 98 Libraries Report.
--
-- Compute the (floor of the) log of i in base b.
-- Simplest way would be just divide i by b until it's smaller then b,
-- but that would be very slow! We are just slightly more clever.
integerLogBase :: Integer -> Integer -> Integer
integerLogBase b i =
if i < b then
0
else
-- Try squaring the base first to cut down the number of divisions.
let l = 2 * integerLogBase (b*b) i
doDiv :: Integer -> Integer -> Integer
doDiv i l = if i < b then l else doDiv (i `div` b) (l+1)
in doDiv (i `div` (b^l)) l
-- Transform an Integer into an Expression of base 2.
intToExpr2 :: Integer -> Expr
intToExpr2 0 = Num 0
intToExpr2 1 = Num 1
intToExpr2 i = let b = 2
nextPow = integerLogBase b i
diff = i - (b^nextPow)
nextExpr = intToExpr2 nextPow
diffExpr = intToExpr2 diff
baseExpr = pow (Num 1, Num b, nextExpr)
-- eval gets rid of lone Num 0 and Pow (_,_,Num 0) expressions.
in add (diffExpr, baseExpr)
-- Transform an Expression into an Integer.
exprToInteger :: Expr -> Integer
exprToInteger expr = exprToInteger' expr 0
exprToInteger' :: Expr -> Integer -> Integer
exprToInteger' (Num a) acc = acc + a
exprToInteger' (Add (a,b)) acc = let left = exprToInteger' a 0
right = exprToInteger' b 0
in acc + left + right
exprToInteger' (Pow (c, b, e)) acc = let coef = exprToInteger' c 0
base = exprToInteger' b 0
exp = exprToInteger' e 0
in acc + (coef * (base^exp) )
-- Convert a Char into a Maybe Integer, using Nothing to signal non-numeric input.
char2Integer :: Char -> Maybe Integer
char2Integer '0' = Just 0
char2Integer '1' = Just 1
char2Integer '2' = Just 2
char2Integer '3' = Just 3
char2Integer '4' = Just 4
char2Integer '5' = Just 5
char2Integer '6' = Just 6
char2Integer '7' = Just 7
char2Integer '8' = Just 8
char2Integer '9' = Just 9
char2Integer _ = Nothing
-- As simple function to convert our input string into a number.
-- We use char2Integer to get a list of Maybe Integer types.
-- If we find a Nothing value, then we denote the condition by returning the Integer 0.
-- This is not elegant, but it works for our simply code.
-- Also, notice we pass the helper function the list of Maybe Integers in reverse order
-- to make multiplying by the base number easier to compute.
readInteger :: String -> Integer
readInteger [] = 0
readInteger str = let intStr = reverse (map (char2Integer) str)
in readInteger' intStr 1 0
readInteger' :: [Maybe Integer] -> Integer -> Integer -> Integer
readInteger' [] base acc = acc
readInteger' (Nothing:ss) base acc = 0
readInteger' ((Just s):ss) base acc = let base' = 10 * base
acc' = acc + (s * base)
in readInteger' ss base' acc'
{- ---------------------------------
Main Functions.
------------------------------------ -}
main = main0
-- User inputs a number.
main0 = do
putStrLn ""
putStrLn "Input a number > "
line <- getLine
main1 line
-- After a line is input, try to get a number from it.
main1 [] = main0
main1 line = let num = readInteger line
expr = intToExpr2 num
in
do
putStrLn (concat ["The start of the Goodstein sequence is, ", (show num), ": \n"])
putStrLn (formatExpr expr)
main2 expr
-- Loop over choice of next term in sequence.
main2 expr =
do
putStrLn ""
putStrLn "Would you like to: "
putStrLn " (e)valuate this term "
putStrLn " (s)how this term "
putStrLn " (g)et the next term "
putStrLn " (l)eap to later term "
putStrLn " (i)nput a new number "
putStrLn " (q)uit? \n"
choice <- getLine
let makeChoice [] = ' '
makeChoice (c:cs) = c
case (makeChoice choice) of
'e' -> let num = exprToInteger expr
in do
putStrLn ""
putStrLn "This term evaluates to the Integer: \n"
putStrLn (show num)
putStrLn ""
main2 expr
's' -> do
putStrLn (show expr)
main2 expr
'g' -> let expr' = nextTerm expr
in do
putStrLn (formatExpr expr')
main2 expr'
'l' ->do
putStrLn "Input the number of terms to leap over: \n"
line <- getLine
putStrLn "The term evaluates to the expression: \n"
let expr' = getTerm' (readInteger line) expr
putStrLn (formatExpr expr')
main2 expr'
'i' -> main0
'q' -> return ()
_ -> do
putStrLn "Choice not recognized, please try again! \n"
main2 expr
{- ---------------------------------
End of Code: Update Log
------------------------------------
2006-07-10 -- Fixed:
- Rule 4 in function sub
from: B^(e+E) - 1 => [B^(e-1) - 1] + [B - 1] * B^((e-1)+E)
to: B^(e+E) - 1 => [B^((e-1)+E) - 1] + [B - 1] * B^((e-1)+E)
This correction fixes all known outstanding bugs.
- Simplified Rule 5 from a case expression to a simple recursion on sub.
The original code is left commented out for comparative purposes.
- Added functions getTerm, which takes a number to start a Goodstein
sequence and the number of terms to evaluate.
- Added function getTerm', a helper for getTerm. This function takes
an expression and the number of terms to evaluate.
- Added function getTermF, which provides formated output by calling
formatExpr on the output from getTerm.
- Added option in main2 to show the Expr data structure of a term.
- Fixed mutate and mutate' so that Pow terms with Pow in the base,
Add in coefficient, and Pow in the coefficient will properly mutate.
This can be exhibited by (getTermF 256 2).
- Fixed formatExpr so that there are proper parenthesis around coefficients
of power elements.
- Added choices to the REPL interpreter to leap ahead using getTerm and
to improve error handling.
2006-07-07 -- Original posting, work in progress.
------------------------------------ -}