diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,16 @@ +2012-05-12 John D. Ramsdell + + * cmu.cabal (Version): Released as version 1.8. + + * src/Algebra/CommutativeMonoid/LinDiophEq.hs: Used the method in + the paper by Contejean and Devie for solving an inhomogeneous + equation using the algorithm for the homogeneous equation solver. + The previous version buggy. + +2012-05-02 John D. Ramsdell + + * cmu.cabal (Version): Released as version 1.7. + 2012-04-27 John D. Ramsdell * src/Algebra/CommutativeMonoid/LinDiophEq.hs: Corrected spelling diff --git a/cmu.cabal b/cmu.cabal --- a/cmu.cabal +++ b/cmu.cabal @@ -1,5 +1,5 @@ Name: cmu -Version: 1.7 +Version: 1.8 Maintainer: ramsdell@mitre.org Cabal-Version: >= 1.2 License: GPL diff --git a/src/Algebra/CommutativeMonoid/LinDiophEq.hs b/src/Algebra/CommutativeMonoid/LinDiophEq.hs --- a/src/Algebra/CommutativeMonoid/LinDiophEq.hs +++ b/src/Algebra/CommutativeMonoid/LinDiophEq.hs @@ -25,63 +25,39 @@ -- The solver uses the algorithm of Contejean and Devie as specified -- in \"An Efficient Incremental Algorithm for Solving Systems of -- Linear Diophantine Equations\", Information and Computation --- Vol. 113, pp. 143-174, 1994 after a modification explained below. +-- Vol. 113, pp. 143-174, 1994. -- -- The algorithm for systems of homogeneous linear Diophantine -- equations follows. Let e[k] be the kth basis vector for 1 <= k <= -- n. To find the minimal, non-negative solutions M to the system of --- equations sum(i=1,n,a[i]*v[i]) = 0, the modified algorithm of --- Contejean and Devie is: +-- equations sum(i=1,n,a[i]*v[i]) = 0, the algorithm of Contejean and +-- Devie is: -- -- 1. [init] A := {e[k] | 1 <= k <= n}; M := {} -- -- 2. [new minimal results] M := M + {a in A | a is a solution} -- --- 3. [breadth-first search] A := {a + e[k] | a in A, 1 <= k <= n, --- \ \< 0} --- --- 4. [unnecessary branches] A := {a in A | all m in M : some +-- 3. [unnecessary branches] A := {a in A | all m in M : some -- 1 <= k <= n : m[k] < a[k]} -- --- 5. [test] If A = {}, stop, else go to 2. +-- 4. [breadth-first search] A := {a + e[k] | a in A, 1 <= k <= n, +-- \ \< 0} -- --- The original algorithm reversed steps 3 and 4. +-- 5. [test] If A = {}, stop, else go to 2. -- -- This module provides a solver for a single linear Diophantine -- equation a*v = b, where a and v are vectors, not matrices. --- Conceptually, it uses the homogeneous solver after appending -b as --- the last element of v and by appending 1 to a at each step in the --- computation. The extra 1 is omitted when an answer is produced. -- --- Steps 3 and 4 were switched because the use of the original --- algorithm for the problem 2x + y - z = 2 produces a non-minimal --- solution. linDiophEq [2,1,-1] 2 = [[1,0,0],[0,2,0]], but the --- original algorithm produces [[1,0,0],[0,2,0],[1,1,1]]. +-- When solving an inhomogeneous equation, it uses the homogeneous +-- solver after adding -b as the first element of v and by bounding +-- the first element of a to be one at each step in the computation. +-- The first element of a solution is zero if it is a solution to the +-- associated homogeneous equation, and one if it is a solution to the +-- inhomogeneous equation. -- -- The algorithm is likely to be Fortenbacher's algorithm, the one -- generalized to systems of equations by Contejean and Devie, but I --- have not been able to verified this fact. I learned how to extend --- Contejean and Devie's results to an inhomogeneous equation by --- reading \"Effective Solutions of Linear Diophantine Equation --- Systems with an Application to Chemistry\" by David Papp and Bela --- Vizari, Rutcor Research Report RRR 28-2004, September, 2004, --- . --- --- The example that shows a problem with the original algorithm --- follows. For the problem linDiophEq [2,1,-1] 2, the value of a and --- m at the beginning of the loop is: --- --- @ --- a m --- [[0, 0, 1], [0, 1, 0], [1, 0, 0]] [] --- [[0, 1, 1], [0, 2, 0]] [[1, 0, 0]] --- [] [[1, 0, 0], [0, 2, 0]] --- @ --- --- Consider [0, 1, 1] in a. If you remove unnecessary branches first, --- the element will stay in a. After performing breadth-first search, --- a will contain [1, 1, 1], which is the unwanted, non-minimal --- solution. +-- have not been able to verified this fact. module Algebra.CommutativeMonoid.LinDiophEq (linDiophEq) where @@ -113,14 +89,32 @@ -- | The 'linDiophEq' function takes a list of integers that specifies -- the coefficients of linear Diophantine equation and a constant, --- and returns the equation's minimal, non-negative solutions. When --- solving an inhomogeneous equation, solve the related homogeneous --- equation and add in those solutions. +-- and returns the equation's minimal, non-negative solutions. +-- +-- When solving an inhomogeneous equation, the first element of a +-- solution is zero if it solves the associated homogeneous equation, +-- and one otherwise. +-- +-- Thus to solve 2x + y - z = 2, use +-- +-- @ +-- linDiophEq [2,1,-1] 2 = [[0,0,1,1],[1,1,0,0],[0,1,0,2],[1,0,2,0]] +-- @ +-- +-- The two minimal solutions to the homogeneous equation are [0,1,1] +-- and [1,0,2], so any linear combinations of these solutions +-- contributes to a solution. The solution that corresponds to +-- [1,0,0] is x = w + 1, y = v, and z = v + 2w. The solution that +-- corresponds to [0,2,0] is x = w, y = v + 2, and z = v + 2w. + linDiophEq :: [Int] -> Int -> [[Int]] linDiophEq [] _ = [] -linDiophEq v c = - newMinimalResults (vector n v) c (basis n) S.empty +linDiophEq v 0 = + newMinimalResults True (vector n v) (basis n) S.empty where n = length v +linDiophEq v c = + newMinimalResults False (vector n (negate c:v)) (basis n) S.empty + where n = 1 + length v -- Construct the basis vectors for an n-dimensional space basis :: Int -> Set (Vector Int) @@ -131,35 +125,36 @@ -- This is the main loop. -- Add elements of a that solve the equation to m and the output -newMinimalResults :: Vector Int -> Int -> Set (Vector Int) -> +-- Variable hom is true when solving a homogeneous equation +newMinimalResults :: Bool -> Vector Int -> Set (Vector Int) -> Set (Vector Int) -> [[Int]] newMinimalResults _ _ a _ | S.null a = [] -newMinimalResults v c a m = +newMinimalResults hom v a m = loop m (S.toList a) -- Test each element in a where loop m [] = -- When done, prepare for next iteration - let a' = breadthFirstSearch v c a -- Step 3 - a'' = unnecessaryBranches a' m in -- Step 4 --- The original algorithm reverses these two steps. --- let a' = unnecessaryBranches a m --- a'' = breadthFirstSearch v c a' in - newMinimalResults v c a'' m + let a' = unnecessaryBranches a m + a'' = breadthFirstSearch hom v a' in + newMinimalResults hom v a'' m loop m (x:xs) - | prod v x == c && S.notMember x m = + | prod v x == 0 && S.notMember x m = elems x:loop (S.insert x m) xs -- Answer found | otherwise = loop m xs -- Breadth-first search using the algorithm of Contejean and Devie -breadthFirstSearch :: Vector Int -> Int -> Set (Vector Int) -> Set (Vector Int) -breadthFirstSearch v c a = +-- Variable hom is true when solving a homogeneous equation +breadthFirstSearch :: Bool -> Vector Int -> + Set (Vector Int) -> Set (Vector Int) +breadthFirstSearch hom v a = S.fold f S.empty a where f x acc = foldl (flip S.insert) acc [ x // [(k, x!k + 1)] | - k <- indices x, - (prod v x - c) * v!k < 0 ] -- Fortenbacher contribution + k <- indices x, -- When not hom, bound first element + hom || k > 0 || x!k == 0, -- of x to be no more than one + prod v x * v!k < 0 ] -- Fortenbacher contribution -- Inner product prod :: Vector Int -> Vector Int -> Int