I've been studying equational unification. I decided to test my
understanding of it by implementing unification and matching in
Abelian groups. I am quite surprised by how little code it takes.
Let me share it with you.
John
Test cases:
2x+y=3z
2x=x+y
64x-41y=a
Code:
> -- Unification and matching in Abelian groups
> -- John D. Ramsdell -- August 2009
>
> module Main (main, test) where
>
> import Data.Char (isSpace, isAlpha, isAlphaNum, isDigit)
> import Data.List (sort)
> import System.IO (isEOF)
>
> -- Chapter 8, Section 5 of the Handbook of Automated Reasoning by
> -- Franz Baader and Wayne Snyder describes unification and matching in
> -- communtative/monoidal theories. This module refines the described
> -- algorithms for the special case of Abelian groups.
>
> -- In this module, an Abelian group is a free algebra over a signature
> -- with three function symbols,
> --
> -- * the binary symbol +, the group operator,
> -- * a constant 0, the identity element, and
> -- * the unary symbol -, the inverse operator.
> --
> -- The algebra is generated by a set of variables. Syntactically, a
> -- variable is an identifer such as x and y.
>
> -- The axioms associated with the algebra are:
> --
> -- * x + y = y + x Commutativity
> -- * (x + y) + z = x + (y + z) Associativity
> -- * x + 0 = x Group identity
> -- * x + -x = 0 Cancellation
>
> -- A substitution maps variables to terms. A substitution s is
> -- extended to a term as follows.
> --
> -- s(0) = 0
> -- s(-t) = -s(t)
> -- s(t + t') = s(t) + s(t')
>
> -- The unification problem is given the problem statement t =? t',
> -- find a substitution s such that s(t) = s(t') modulo the axioms of
> -- the algebra. The matching problem is to find substitution s such
> -- that s(t) = t' modulo the axioms.
>
> -- A term is represented as the sum of factors, and a factor is the
> -- product of an integer coeficient and a variable or the group
> -- identity, zero. In this representation, every coeficient is
> -- non-zero, and no variable occurs twice.
>
> -- A term can be represented by a finite map from variables to
> -- non-negative integers. To make the code easier to understand,
> -- association lists are used instead of Data.Map.
>
> newtype Lin = Lin [(String, Int)]
>
> -- Constructors
>
> -- Identity element (zero)
> ide :: Lin
> ide = Lin []
>
> -- Factors
> var :: Int -> String -> Lin
> var 0 _ = Lin []
> var c x = Lin [(x, c)]
>
> -- Invert by negating coefficients.
> neg :: Lin -> Lin
> neg (Lin t) =
> Lin $ map (\(x, c) -> (x, negate c)) t
>
> -- Join terms ensuring that coefficients are non-zero, and no variable
> -- occurs twice.
> add :: Lin -> Lin -> Lin
> add (Lin t) (Lin t') =
> Lin $ foldr f t' t
> where
> f (x, c) t =
> case lookup x t of
> Just c' | c + c' == 0 -> remove x t
> | otherwise -> (x, c + c') : remove x t
> Nothing -> (x, c) : t
>
> -- Remove the first pair in an association list that matches the key.
> remove :: Eq a => a -> [(a, b)] -> [(a, b)]
> remove _ [] = []
> remove x (y@(z, _) : ys)
> | x == z = ys
> | otherwise = y : remove x ys
>
> canonicalize :: Lin -> Lin
> canonicalize (Lin t) =
> Lin (sort t)
>
> -- Convert a linearized term into an association list.
> assocs :: Lin -> [(String, Int)]
> assocs (Lin t) = t
>
> term :: [(String, Int)] -> Lin
> term assoc =
> foldr f ide assoc
> where
> f (x, c) t = add t $ var c x
>
> -- Unification and Matching
>
> newtype Equation = Equation (Lin, Lin)
>
> newtype Maplet = Maplet (String, Lin)
>
> -- Unification is the same as matching when there are no constants
> unify :: Monad m => Equation -> m [Maplet]
> unify (Equation (t0, t1)) =
> match $ Equation (add t0 (neg t1), ide)
>
> -- Matching in Abelian groups is performed by finding integer
> -- solutions to linear equations, and then using the solutions to
> -- construct a most general unifier.
> match :: Monad m => Equation -> m [Maplet]
> match (Equation (t0, t1)) =
> case (assocs t0, assocs t1) of
> ([], []) -> return []
> ([], _) -> fail "no solution"
> (t0, t1) ->
> do
> subst <- intLinEq (map snd t0) (map snd t1)
> return $ mgu (map fst t0) (map fst t1) subst
>
> -- Construct a most general unifier from a solution to a linear
> -- equation. The function adds the variables back into terms, and
> -- generates fresh variables as needed.
> mgu :: [String] -> [String] -> Subst -> [Maplet]
> mgu vars syms subst =
> foldr f [] (zip vars [0..])
> where
> f (x, n) maplets =
> case lookup n subst of
> Just (factors, consts) ->
> Maplet (x, g factors consts) : maplets
> Nothing ->
> Maplet (x, var 1 $ genSym n) : maplets
> g factors consts =
> term (zip genSyms factors ++ zip syms consts)
> genSyms = map genSym [0..]
>
> -- Generated variables start with this character.
> genChar :: Char
> genChar = 'g'
>
> genSym :: Int -> String
> genSym i = genChar : show i
>
> -- So why solve linear equations? Consider the matching problem
> --
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] =?
> -- d[0]*a[0] + d[1]*a[1] + ... + d[m-1]*a[m-1]
> --
> -- with n variables and m constants. We seek a most general unifier s
> -- such that
> --
> -- s(c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1]) =
> -- d[0]*a[0] + d[1]*a[1] + ... + d[m-1]*a[m-1]
> --
> -- which is the same as
> --
> -- c[0]*s(x[0]) + c[1]*s(x[1]) + ... + c[n-1]*s(x[n-1]) =
> -- d[0]*a[0] + d[1]*a[1] + ... + d[m-1]*a[m-1]
> --
> -- Notice that the number of occurrences of constant a[0] in s(x[0])
> -- plus s(x[1]) ... s(x[n-1]) must equal d[0]. Thus the mappings of
> -- the unifier that involve constant a[0] respect integer solutions of
> -- the following linear equation.
> --
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = d[0]
> --
> -- To compute a most general unifier, a most general integer solution
> -- to a linear equation must be found.
>
> -- Integer Solutions of Linear Inhomogeneous Equations
>
> type LinEq = ([Int], [Int])
>
> -- A linear equation with integer coefficients is represented as a
> -- pair of lists of integers, the coefficients and the constants. If
> -- there are no constants, the linear equation represented by (c, [])
> -- is the homogeneous equation:
> --
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = 0
> --
> -- where n is the length of c. Otherwise, (c, d) represents a
> -- sequence of inhomogeneous linear equations with the same
> -- left-hand-side:
> --
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = d[0]
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = d[1]
> -- ...
> -- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = d[m-1]
> --
> -- where m is the length of d.
>
> type Subst = [(Int, LinEq)]
>
> -- A solution is a partial map from variables to terms, and a term is
> -- a pair of lists of integers, the variable part of the term followed
> -- by the constant part. The variable part may specify variables not
> -- in the input. For example, the solution of
> --
> -- 64x = 41y + 1
> --
> -- is x = -41z - 16 and y = -64z - 25. The computed solution is read
> -- off the list returned as an answer.
> --
> -- intLinEq [64,-41] [1] =
> -- [(0,([0,0,0,0,0,0,-41],[-16])),
> -- (1,([0,0,0,0,0,0,-64],[-25]))]
>
> -- Find integer solutions to linear equations
> intLinEq :: Monad m => [Int] -> [Int] -> m Subst
> intLinEq coefficients constants =
> intLinEqLoop (length coefficients) (coefficients, constants) []
>
> -- The algorithm used to find solutions is described in Vol. 2 of The
> -- Art of Computer Programming / Seminumerical Alorithms, 2nd Ed.,
> -- 1981, by Donald E. Knuth, pg. 327.
>
> -- On input, n is the number of variables in the original problem, c
> -- is the coefficients, d is the constants, and subst is a list of
> -- eliminated variables.
> intLinEqLoop :: Monad m => Int -> LinEq -> Subst -> m Subst
> intLinEqLoop n (c, d) subst =
> -- Find the smallest coefficient in absolute value
> let (i, ci) = smallest c in
> case () of
> _ | ci < 0 -> intLinEqLoop n (invert c, invert d) subst
> -- Ensure the smallest coefficient is positive
> | ci == 0 -> fail "bad problem"
> -- Lack of non-zero coefficients is an error
> | ci == 1 ->
> -- A general solution of the following form has been found:
> -- x[i] = sum[j] -c'[j]*x[j] + d[k] for all k
> -- where c' is c with c'[i] = 0.
> return $ eliminate n (i, (invert (zero i c), d)) subst
> | divisible ci c ->
> -- If all the coefficients are divisible by c[i], a solution is
> -- immediate if all the constants are divisible by c[i],
> -- otherwise there is no solution.
> if divisible ci d then
> let c' = divide ci c
> d' = divide ci d in
> return $ eliminate n (i, (invert (zero i c'), d')) subst
> else
> fail "no solution"
> | otherwise ->
> -- Eliminate x[i] in favor of freshly created variable x[n],
> -- where n is the length of c.
> -- x[n] = sum[j] (c[j] div c[i] * x[j])
> -- The new equation to be solved is:
> -- c[i]*x[n] + sum[j] c'[j]*x[j] = d[k] for all k
> -- where c'[j] = c[j] mod c[i] for j /= i and c'[i] = 0.
> let c' = map (\x -> mod x ci) (zero i c)
> c'' = divide ci (zero i c)
> subst' = eliminate n (i, (invert c'' ++ [1], [])) subst in
> intLinEqLoop n (c' ++ [ci], d) subst'
>
> -- Find the smallest coefficient in absolute value
> smallest :: [Int] -> (Int, Int)
> smallest xs =
> foldl f (-1, 0) (zip [0..] xs)
> where
> f (i, n) (j, x)
> | n == 0 = (j, x)
> | x == 0 || abs n <= abs x = (i, n)
> | otherwise = (j, x)
>
> invert :: [Int] -> [Int]
> invert t = map negate t
>
> -- Zero the ith position in a list
> zero :: Int -> [Int] -> [Int]
> zero _ [] = []
> zero 0 (_:xs) = 0 : xs
> zero i (x:xs) = x : zero (i - 1) xs
>
> -- Eliminate a variable from the existing substitution. If the
> -- variable is in the original problem, add it to the substitution.
> eliminate :: Int -> (Int, LinEq) -> Subst -> Subst
> eliminate n m@(i, (c, d)) subst =
> if i < n then
> m : map f subst
> else
> map f subst
> where
> f m'@(i', (c', d')) = -- Eliminate i in c' if it occurs in c'
> case get i c' of
> 0 -> m' -- i is not in c'
> ci -> (i', (addmul ci (zero i c') c, addmul ci d' d))
> -- Find ith coefficient
> get _ [] = 0
> get 0 (x:_) = x
> get i (_:xs) = get (i - 1) xs
> -- addnum n xs ys sums xs and ys after multiplying ys by n
> addmul 1 [] ys = ys
> addmul n [] ys = map (* n) ys
> addmul _ xs [] = xs
> addmul n (x:xs) (y:ys) = (x + n * y) : addmul n xs ys
>
> divisible :: Int -> [Int] -> Bool
> divisible small t =
> all (\x -> mod x small == 0) t
>
> divide :: Int -> [Int] -> [Int]
> divide small t =
> map (\x -> div x small) t
>
> -- Input and Output
>
> instance Show Lin where
> showsPrec _ (Lin []) =
> showString "0"
> showsPrec _ x =
> showFactor t . showl ts
> where
> Lin (t:ts) = canonicalize x
> showFactor (x, 1) = showString x
> showFactor (x, -1) = showChar '-' . showString x
> showFactor (x, c) = shows c . showString x
> showl [] = id
> showl ((s,n):ts)
> | n < 0 =
> showString " - " . showFactor (s, negate n) . showl ts
> showl (t:ts) = showString " + " . showFactor t . showl ts
>
> instance Read Lin where
> readsPrec _ s0 =
> [ (t1, s2) | (t0, s1) <- readFactor s0,
> (t1, s2) <- readRest t0 s1 ]
> where
> readPrimary s0 =
> [ (t0, s1) | (x, s1) <- scan s0, isVar x,
> let t0 = var 1 x ] ++
> [ (t0, s1) | ("0", s1) <- scan s0,
> (s, _) <- scan s1, not (isVar s),
> let t0 = ide ] ++
> [ (t0, s2) | (n, s1) <- scan s0, isNum n,
> (x, s2) <- scan s1, isVar x,
> let t0 = var (read n) x ] ++
> [ (t0, s3) | ("(", s1) <- scan s0,
> (t0, s2) <- reads s1,
> (")", s3) <- scan s2 ]
> readFactor s0 =
> [ (t1, s2) | ("-", s1) <- scan s0,
> (t0, s2) <- readPrimary s1,
> let t1 = neg t0 ] ++
> [ (t0, s1) | (s, _) <- scan s0, s /= "-",
> (t0, s1) <- readPrimary s0 ]
> readRest t0 s0 =
> [ (t2, s3) | ("+", s1) <- scan s0,
> (t1, s2) <- readFactor s1,
> (t2, s3) <- readRest (add t0 t1) s2 ] ++
> [ (t2, s3) | ("-", s1) <- scan s0,
> (t1, s2) <- readPrimary s1,
> (t2, s3) <- readRest (add t0 (neg t1)) s2 ] ++
> [ (t0, s0) | (s, _) <- scan s0, s /= "+" && s /= "-" ]
>
> isNum :: String -> Bool
> isNum (c:_) = isDigit c
> isNum _ = False
>
> isVar :: String -> Bool
> isVar (c:_) = isAlpha c && c /= genChar
> isVar _ = False
>
> scan :: ReadS String
> scan "" = [("", "")]
> scan (c:s)
> | isSpace c = scan s
> | isAlpha c = [ (c:part, t) | (part,t) <- [span isAlphaNum s] ]
> | isDigit c = [ (c:part, t) | (part,t) <- [span isDigit s] ]
> | otherwise = [([c], s)]
>
> instance Show Equation where
> showsPrec _ (Equation (t0, t1)) =
> shows t0 . showString " = " . shows t1
>
> instance Read Equation where
> readsPrec _ s0 =
> [ (Equation (t0, t1), s3) | (t0, s1) <- reads s0,
> ("=", s2) <- scan s1,
> (t1, s3) <- reads s2 ]
>
> instance Show Maplet where
> showsPrec _ (Maplet (x, t)) =
> showString x . showString " -> " . shows t
>
> -- Test Routine
>
> -- Given an equation, display a unifier and a matcher.
> test :: String -> IO ()
> test prob =
> case readM prob of
> Err err -> putStrLn err
> Ans (Equation (t0, t1)) ->
> do
> putStr "Problem: "
> print $ Equation (canonicalize t0, canonicalize t1)
> subst <- unify $ Equation (t0, t1)
> putStr "Unifier: "
> print subst
> putStr "Matcher: "
> case match $ Equation (t0, t1) of
> Err err -> putStrLn err
> Ans subst -> print subst
> putStrLn ""
>
> readM :: (Read a, Monad m) => String -> m a
> readM s =
> case [ x | (x, t) <- reads s, ("", "") <- lex t ] of
> [x] -> return x
> [] -> fail "no parse"
> _ -> fail "ambiguous parse"
>
> data AnsErr a
> = Ans a
> | Err String
>
> instance Monad AnsErr where
> (Ans x) >>= k = k x
> (Err s) >>= _ = Err s
> return = Ans
> fail = Err
>
> main :: IO ()
> main =
> do
> done <- isEOF
> case done of
> True -> return ()
> False ->
> do
> prob <- getLine
> test prob
> main
