-- | Representation of coherent rings. Traditionally a ring is coherent if every
-- finitely generated ideal is finitely presented. This means that it is 
-- possible to solve homogenous linear equations in them.
module Algebra.Structures.Coherent
  ( Coherent(solve)
  , propCoherent
  , solveMxN, propSolveMxN
  , solveWithIntersection
  , solveGeneralEquation, propSolveGeneralEquation
  , solveGeneral, propSolveGeneral
  ) where

import Test.QuickCheck

import Algebra.Structures.IntegralDomain
import Algebra.Structures.StronglyDiscrete
import Algebra.Matrix
import Algebra.Ideal


-------------------------------------------------------------------------------
-- | Definition of coherent rings.
--
-- We say that R is coherent iff for any M, we can find L such that ML=0 and
--
--   MX=0   \<-\>  \exists Y. X=LY
--
-- that is, iff we can generate the solutions of any linear homogeous system 
-- of equations.
--
-- The main point here is that ML=0, it is not clear how to represent the
-- equivalence in Haskell. This would probably be possible in type theory.
--
class IntegralDomain a => Coherent a where
  solve :: Vector a -> Matrix a

isSolution :: (CommutativeRing a, Eq a) => Matrix a -> Matrix a -> Bool
isSolution m sol = all (==zero) (concat (unMVec (m `mulM` sol)))

propCoherent :: (Coherent a, Eq a) => Vector a -> Bool
propCoherent m = isSolution (vectorToMatrix m) (solve m)

-------------------------------------------------------------------------------
-- | Solves a system of equations.
solveMxN :: (Coherent a, Eq a) => Matrix a -> Matrix a
solveMxN (M (l:ls)) = solveMxN' (solve l) ls
  where
  -- Inductively solve all subsystems. If the computed solution is in fact a 
  -- solution to the next set of equations then don't do anything. 
  -- This solves the problems with having many identical rows in the system,
  -- like [[1,1],[1,1]].
  solveMxN' :: (Coherent a, Eq a) => Matrix a -> [Vector a] -> Matrix a
  solveMxN' m []      = m
  solveMxN' m1 (x:xs) = if isSolution (vectorToMatrix x) m1 
                           then solveMxN' m1 xs 
                           else solveMxN' (m1 `mulM` m2) xs
    where m2 = solve (matrixToVector (mulM (vectorToMatrix x) m1))


-- | Test that the solution of an MxN system is in fact a solution of the system
propSolveMxN :: (Coherent a, Eq a) => Matrix a -> Bool
propSolveMxN m = isSolution m (solveMxN m)


-------------------------------------------------------------------------------
-- | Intersection computable -> Coherence.
-- 
-- Proof that if there is an algorithm to compute a f.g. set of generators for 
-- the intersection of two f.g. ideals then the ring is coherent.
--
-- Takes the vector to solve, \[x1,...,xn\], and a function (int) that computes 
-- the intersection of two ideals. 
--
-- If
--       
--     \[ x_1, ..., x_n \] \`int\` \[ y_1, ..., y_m \] = \[ z_1, ..., z_l \]
--
-- then int should give witnesses us and vs such that:
--
--     z_k = n_k1 * x_1 + ... + u_kn * x_n
--
--         = u_k1 * y_1 + ... + n_km * y_m
--
solveWithIntersection :: (IntegralDomain a, Eq a)
                      => Vector a 
                      -> (Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]]))
                      -> Matrix a
solveWithIntersection (Vec xs) int = transpose $ matrix $ solveInt xs 
  where
  solveInt []     = error "solveInt: Can't solve an empty system"
  solveInt [x]    = [[zero]] -- Base case, could be [x,y] also...
                             -- That wouldn't give the trivial solution... 
--  solveInt [x,y]  | x == zero || y == zero = [[zero,zero]]
--                  | otherwise = 
--    let (Id ts,us,vs) = (Id [x]) `int` (Id [neg y])
--    in [ u ++ v | (u,v) <- zip us vs ]
  solveInt (x:xs)
    | x == zero             = map (zero:) $ solveInt xs 
    | isSameIdeal int as bs = s ++ m'
    | otherwise             = error "solveInt: This does not compute the intersection"
      where
      as            = Id [x]
      bs            = Id (map neg xs)

      -- Compute the intersection of <x1> and <-x2,...,-xn>
      (Id ts,us,vs) = as `int` bs
      s             = [ u ++ v | (u,v) <- zip us vs ]
      
      -- Solve <0,x2,...,xn> recursively
      m             = solveInt xs
      m'            = map (zero:) m


-------------------------------------------------------------------------------
-- | Strongly discrete coherent rings.
--
-- If the ring is strongly discrete and coherent then we can solve arbitrary 
-- equations of the type AX=b. 
--
solveGeneralEquation :: (Coherent a, StronglyDiscrete a) => Vector a -> a -> Maybe (Matrix a)
solveGeneralEquation v@(Vec xs) b = 
  let sol = solve v
  in case b `member` (Id xs) of
    Just as -> Just $ transpose (M (replicate (length (head (unMVec sol))) (Vec as)))
                      `addM` sol
    Nothing -> Nothing


propSolveGeneralEquation :: (Coherent a, StronglyDiscrete a, Eq a) 
                         => Vector a 
                         -> a 
                         -> Bool
propSolveGeneralEquation v b = case solveGeneralEquation v b of
  Just sol -> all (==b) $ concat $ unMVec $ vectorToMatrix v `mulM` sol
  Nothing  -> True


isSolutionB v sol b = all (==b) $ concat $ unMVec $ vectorToMatrix v `mulM` sol


-- | Solves general linear systems of the kind AX = B.
--
-- A is given as a matrix and B is given as a row vector (it should be column 
-- vector).
--
solveGeneral :: (Coherent a, StronglyDiscrete a, Eq a) 
             => Matrix a   -- M
             -> Vector a   -- B
             -> Maybe (Matrix a, Matrix a)  -- (L,X0)
solveGeneral (M (l:ls)) (Vec (a:as)) = 
  case solveGeneral' (solveGeneralEquation l a) ls as [(l,a)] of
    Just x0 -> Just (solveMxN (M (l:ls)), x0)
    Nothing -> Nothing
  where
  -- Compute a new solution inductively and check that the new solution 
  -- satisfies all the previous equations.
  solveGeneral' Nothing _ _ _              = Nothing
  solveGeneral' (Just m) [] [] old         = Just m 
  solveGeneral' (Just m) (l:ls) (a:as) old = 
    if isSolutionB l m a 
       then solveGeneral' (Just m) ls as old
       else case solveGeneralEquation (matrixToVector (vectorToMatrix l `mulM` m)) a of
         Just m' -> let m'' = m `mulM` m'
                    in if all (\(x,y) -> isSolutionB x m'' y) old
                          then solveGeneral' (Just m'') ls as ((l,a):old)
                          else Nothing
         Nothing -> Nothing
  solveGeneral' _ _ _ _ = error "solveGeneral: Bad input"      

-- It would be great to only generate solvable systems...
-- propSolveGeneral :: (Coherent a, StronglyDiscrete a, Eq a) => Matrix a -> Vector a -> Property 
propSolveGeneral m b = length (unM m) == length (unVec b) ==> case solveGeneral m b of
  Just (l,x) -> all (==b) (unM (transpose (m `mulM` x))) &&
                isSolution m l
  Nothing -> True