diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,24 @@ +2009-09-17 John D. Ramsdell + + * agum.cabal (Version): Released as version 2.2. + +2009-09-14 John D. Ramsdell + + * src/Algebra/AbelianGroup/IntLinEq.hs: Integer solutions to + linear equation solver was placed in its own module. + +2009-09-13 John D. Ramsdell + + * src/Algebra/AbelianGroup/UnificationMatching.hs (unify): Changed + the result to be a substitution since unification always succeeds. + 2009-09-05 John D. Ramsdell * src/Algebra/AbelianGroup/UnificationMatching.hs: Added - reference to Andrew Kennedy's Ph.D. thesis as in contains a proof + reference to Andrew Kennedy's Ph.D. thesis as it contains a proof of correctness of the implemented matching algorithm. + + * agum.cabal (Version): Released as version 2.1. 2009-08-29 John D. Ramsdell diff --git a/Makefile b/Makefile --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ # Requires GNU Make # The all target creates a default configuration if need be. -PACKAGE = agum +PACKAGE := $(wildcard *.cabal) CONFIG = dist/setup-config SETUP = runhaskell Setup.hs @@ -12,10 +12,10 @@ Makefile: @echo make $@ -$(PACKAGE).cabal: +$(PACKAGE): @echo make $@ -$(CONFIG): $(PACKAGE).cabal +$(CONFIG): $(PACKAGE) $(SETUP) configure --ghc --user --prefix="${HOME}" %: force diff --git a/agum.cabal b/agum.cabal --- a/agum.cabal +++ b/agum.cabal @@ -1,5 +1,5 @@ Name: agum -Version: 2.1 +Version: 2.2 Maintainer: ramsdell@mitre.org Cabal-Version: >= 1.2 License: GPL @@ -21,6 +21,7 @@ Library Build-Depends: base <= 4.1.0.0, containers Exposed-Modules: Algebra.AbelianGroup.UnificationMatching + Algebra.AbelianGroup.IntLinEq Hs-Source-Dirs: src GHC-Options: -Wall -fno-warn-name-shadowing -fwarn-unused-imports @@ -29,6 +30,7 @@ Main-Is: Algebra/AbelianGroup/Main.hs Build-Depends: base <= 4.1.0.0, containers Other-Modules: Algebra.AbelianGroup.UnificationMatching + Algebra.AbelianGroup.IntLinEq Hs-Source-Dirs: src GHC-Options: -Wall -fno-warn-name-shadowing -fwarn-unused-imports diff --git a/readme.txt b/readme.txt --- a/readme.txt +++ b/readme.txt @@ -1,5 +1,5 @@ This package contains a library for unification and matching in -Abelian groups and a program that exercises the library. +an Abelian group and a program that exercises the library. $ agum Abelian group unification and matching -- :? for help @@ -15,7 +15,7 @@ agum> 64x-41y=a Problem: 64x - 41y = a -Unifier: [x : g0,y : g1,a : 64g0 - 41g1] -Matcher: [x : -16a - 41g6,y : -25] +Unifier: [a : 64g1 - 41g2,x : g1,y : g2] +Matcher: [x : -16a - 41g6,y : -25a - 64g6] agum> :quit diff --git a/src/Algebra/AbelianGroup/IntLinEq.hs b/src/Algebra/AbelianGroup/IntLinEq.hs new file mode 100644 --- /dev/null +++ b/src/Algebra/AbelianGroup/IntLinEq.hs @@ -0,0 +1,206 @@ +-- Integer Solutions of Linear Inhomogeneous Equations +-- +-- Copyright (C) 2009 John D. Ramsdell +-- +-- This program is free software: you can redistribute it and/or modify +-- it under the terms of the GNU General Public License as published by +-- the Free Software Foundation, either version 3 of the License, or +-- (at your option) any later version. + +-- This program is distributed in the hope that it will be useful, +-- but WITHOUT ANY WARRANTY; without even the implied warranty of +-- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +-- GNU General Public License for more details. + +-- You should have received a copy of the GNU General Public License +-- along with this program. If not, see . + +-- | +-- Module : Algebra.AbelianGroup.IntLinEq +-- Copyright : (C) 2009 John D. Ramsdell +-- License : GPL +-- +-- Integer Solutions of Linear Inhomogeneous Equations +-- +-- A linear equation with integer coefficients is represented as a +-- pair of lists of non-zero 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 the +-- inhomogeneous equation: +-- +-- > c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = g +-- +-- where g = gcd(d[0], d[1], ..., d[m-1]), and m is the length of d. +-- Thus g is the greatest common denominator of the elements of d. +-- +-- 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. In other words, the length of the coefficents in the +-- answer may exceed the length of the coefficients 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]))] +-- +-- 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. To show sums, we write +-- +-- > sum[i] c[i]*x[i] for c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1]. +-- +-- The algorithm's initial values are the linear equation (c,d) and an +-- empty substitution s. +-- +-- 1. Let c[i] be the smallest non-zero coefficient in absolute value. +-- +-- 2. If c[i] < 0, multiply c and d by -1 and goto step 1. +-- +-- 3. If c[i] = 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. Use the equation to eliminate x[i] +-- from the range of the current substitution s. If variable x[i] is +-- in the original equation, add the mapping to substitution s. +-- +-- 4. If c[i] divides every coefficient in c, +-- +-- * if c[i] divides every constant in d, divide c and d by c[i] +-- and goto step 3, +-- +-- * otherwise fail because there is no solution. +-- +-- 5. Otherwise, eliminate x[i] as above 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]) +-- +-- Goto step 1 and solve the equation: +-- +-- > c[i]*x[n] + sum[j] (c[j] mod c[i])*x[j] = d[k] for all k + +module Algebra.AbelianGroup.IntLinEq + (LinEq, Subst, intLinEq) where + +-- | A linear equation with integer coefficients is represented as a +-- pair of lists of non-zero integers, the coefficients and the +-- constants. +type LinEq = ([Int], [Int]) + +-- | A solution to a linear equation 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. In other words, the length +-- of the coefficents in the answer may exceed the length of the +-- coefficients in the input. +type Subst = [(Int, LinEq)] + +-- | Find integer solutions to a linear equation or fail when there +-- are no solutions. +intLinEq :: Monad m => LinEq -> 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 non-zero 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] mod c[i])*x[j] = d[k] for all k + intLinEqLoop n (map (\x -> mod x ci) c ++ [ci], d) subst' + where + subst' = eliminate n (i, (invert c' ++ [1], [])) subst + c' = divide ci (zero i c) + +-- Find the smallest non-zero 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 diff --git a/src/Algebra/AbelianGroup/Main.hs b/src/Algebra/AbelianGroup/Main.hs --- a/src/Algebra/AbelianGroup/Main.hs +++ b/src/Algebra/AbelianGroup/Main.hs @@ -31,7 +31,7 @@ do putStr "Problem: " print $ Equation (t0, t1) - subst <- unify $ Equation (t0, t1) + let subst = unify $ Equation (t0, t1) putStr "Unifier: " print subst putStr "Matcher: " diff --git a/src/Algebra/AbelianGroup/UnificationMatching.hs b/src/Algebra/AbelianGroup/UnificationMatching.hs --- a/src/Algebra/AbelianGroup/UnificationMatching.hs +++ b/src/Algebra/AbelianGroup/UnificationMatching.hs @@ -1,4 +1,4 @@ --- Unification and matching in Abelian groups +-- Unification and matching in an Abelian group -- -- Copyright (C) 2009 John D. Ramsdell -- @@ -72,6 +72,7 @@ import Data.Char (isSpace, isAlpha, isAlphaNum, isDigit) import Data.Map (Map) import qualified Data.Map as Map +import Algebra.AbelianGroup.IntLinEq -- Chapter 8, Section 5 of the Handbook of Automated Reasoning by -- Franz Baader and Wayne Snyder describes unification and matching in @@ -113,7 +114,7 @@ -- factors. A factor is the product of a non-zero integer coefficient -- and a variable. In this representation, no variable occurs twice. -- Thus a term is represented by a finite map from variables to --- non-negative integers. +-- non-zero integers. -- | A term in an Abelian group is represented by the group identity -- element, or as the sum of factors. A factor is the product of a @@ -213,10 +214,12 @@ -- | Given 'Equation' (t0, t1), return a most general substitution s -- such that s(t0) = s(t1) modulo the equational axioms of an Abelian --- group. -unify :: Monad m => Equation -> m Substitution +-- group. Unification always succeeds. +unify :: Equation -> Substitution unify (Equation (t0, t1)) = - match $ Equation (add t0 (neg t1), ide) + case match $ Equation (add t0 (neg t1), ide) of + Nothing -> error "Internal error--unification failed" + Just s -> s -- Matching in Abelian groups is performed by finding integer -- solutions to linear equations, and then using the solutions to @@ -231,7 +234,7 @@ ([], _) -> fail "no solution" (t0, t1) -> do - subst <- intLinEq (map snd t0) (map snd t1) + 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 @@ -295,139 +298,8 @@ -- 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 the --- inhomogeneous equation: --- --- c[0]*x[0] + c[1]*x[1] + ... + c[n-1]*x[n-1] = g --- --- where g = gcd(d[0], d[1], ..., d[m-1]), and m is the length of d. --- Thus g is the greatest common denominator of the elements 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 non-zero 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] mod c[i])*x[j] = d[k] for all k - intLinEqLoop n (map (\x -> mod x ci) c ++ [ci], d) subst' - where - subst' = eliminate n (i, (invert c' ++ [1], [])) subst - c' = divide ci (zero i c) - --- Find the smallest non-zero 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 +-- to a linear equation must be found. See module +-- Algebra.AbelianGroup.IntLinEq. -- Elementary Abelian group matching is equivalent to unification with -- constants. A proof of correctness of this algorithm, cast as