diff --git a/README.hs b/README.hs
deleted file mode 100644
--- a/README.hs
+++ /dev/null
@@ -1,80 +0,0 @@
--------------------------------------------------------------------------------
--- | Constructive Algebra Library 
--- 
--- Anders Mortberg    <mortberg@student.chalmers.se>
--- Bassel Mannaa      <mannaa@student.chalmers.se>
---
--- Abstract:
--- This is a library written as part of our master theses. It focuses mainly
--- on the theory of commutative rings from a constructive point of view. 
---
--------------------------------------------------------------------------------
-
-module README where
-
-
---------------------------------------------------------------------------------
--- Structures
-
--- Rings with basic operations. 
-import Algebra.Structures.Ring
-
--- Commutative rings.
-import Algebra.Structures.CommutativeRing
-
--- Integral domains.
-import Algebra.Structures.IntegralDomain
-
--- Fields.
-import Algebra.Structures.Field
-
--- Strongly discrete rings - Rings with decidable ideal membership.
-import Algebra.Structures.StronglyDiscrete
-
--- EuclideanDomains - Integral domains with decidable division and and Euclidean
--- function. Contains lots of functions that are possible at the level of 
--- Euclidean domain like the Euclidean algorithm and extended Euclidean 
--- algorithm.
-import Algebra.Structures.EuclideanDomain
-
--- Bezout domains - Non-Noetherian analogues of principal ideal domains. All 
--- finitely generated ideals are principal.
-import Algebra.Structures.BezoutDomain
-
--- GCD domains - Non-Noetherian analogues of unique factorization domains. 
--- All pairs of nonzero elements have a greatest common divisor.
-import Algebra.Structures.GCDDomain
-
--- Field of fractions of a GCD domain.
-import Algebra.Structures.FieldOfFractions
-
--- Coherent rings. That is rings in which it is possible to solve homogenous
--- linear equations. 
-import Algebra.Structures.Coherent
-
-
--------------------------------------------------------------------------------
--- Special constructions.
-
--- Finitely generated ideals over commutative rings. 
-import Algebra.Ideal
-
--- Simple matrix library
-import Algebra.Matrix
-
--- Principle localization matrices
-import Algebra.PLM
-
-
--------------------------------------------------------------------------------
--- Instances.
-
--- The integers.
-import Algebra.Z
-
--- The rational numbers as the field of fractions of Z. 
-import Algebra.Q
-
-
--------------------------------------------------------------------------------
--- The end.
diff --git a/constructive-algebra.cabal b/constructive-algebra.cabal
--- a/constructive-algebra.cabal
+++ b/constructive-algebra.cabal
@@ -7,7 +7,7 @@
 -- The package version. See the Haskell package versioning policy
 -- (http://www.haskell.org/haskellwiki/Package_versioning_policy) for
 -- standards guiding when and how versions should be incremented.
-Version:             0.2.0
+Version:             0.3.0
 
 Synopsis:            A library of constructive algebra.
 Description:         
@@ -19,7 +19,8 @@
         generated. For example, instead of principal ideal domains one gets 
         Bezout domains which are integral domains in which all finitely 
         generated ideals are principal (and not necessarily that all ideals are
-        principal).
+        principal). This give a good framework for implementing many 
+        interesting algorithms.
 
 License:             BSD3
 License-file:        LICENSE
@@ -40,7 +41,7 @@
 
 -- Extra files to be distributed with the package, such as examples or
 -- a README.
-Extra-source-files:  README.hs, examples/Z_Examples.hs
+Extra-source-files:  examples/Z_Examples.hs
 
 -- Constraint on the version of Cabal needed to build this package.
 Cabal-version:       >=1.2
@@ -48,34 +49,36 @@
 
 Library
   -- Modules exported by the library.
-  Exposed-modules:     Algebra.Structures.Group, 
-                       Algebra.Structures.Ring,
+  Exposed-modules:     Algebra.Structures.BezoutDomain,
+                       Algebra.Structures.Coherent,
                        Algebra.Structures.CommutativeRing,
-                       Algebra.Structures.IntegralDomain, 
+                       Algebra.Structures.EuclideanDomain,
+                       Algebra.Structures.ExplicitUnits,
                        Algebra.Structures.Field,
+                       Algebra.Structures.FieldOfFractions,
+                       Algebra.Structures.GCDDomain, 
+                       Algebra.Structures.Group, 
+                       Algebra.Structures.IntegralDomain, 
                        Algebra.Structures.Module,
-                       Algebra.Structures.BezoutDomain,
                        Algebra.Structures.PruferDomain,
-                       Algebra.Structures.EuclideanDomain,
+                       Algebra.Structures.Ring,
                        Algebra.Structures.StronglyDiscrete,
-                       Algebra.Structures.FieldOfFractions,
-                       Algebra.Structures.GCDDomain, 
-                       Algebra.Structures.Coherent,
                        Algebra.TypeChar.Char,
+                       Algebra.EllipticCurve,
                        Algebra.FieldOfRationalFunctions,
                        Algebra.Ideal,
                        Algebra.Matrix,
                        Algebra.PLM,
+                       Algebra.Q,
+--                       Algebra.SmithNormalForm,  <- Should be added to separate package
                        Algebra.UPoly,
-                       Algebra.EllipticCurve,
+                       Algebra.Z,                  
                        Algebra.ZSqrt5,
-                       Algebra.Zn,
-                       Algebra.Z,
-                       Algebra.Q
+                       Algebra.Zn
                        
 
   -- Packages needed in order to build this package.
-  Build-depends:       base >= 3 && <= 4, QuickCheck >= 2, type-level >= 0.2
+  Build-depends:       base >= 3 && <= 4.3.1.0, QuickCheck >= 2, type-level >= 0.2
   
   -- Modules not exported by this package.
   -- Other-modules:       
diff --git a/src/Algebra/Matrix.hs b/src/Algebra/Matrix.hs
--- a/src/Algebra/Matrix.hs
+++ b/src/Algebra/Matrix.hs
@@ -4,16 +4,23 @@
   ( Vector(Vec)
   , unVec, lengthVec
   , Matrix(M), matrix
-  , matrixToVector, vectorToMatrix, unMVec, unM 
-  , identity, mulM, addM, transpose, isSquareMatrix, dimension
+  , matrixToVector, vectorToMatrix, unMVec, unM, (!!!)
+  , identity, propLeftIdentity, propRightIdentity
+  , mulM, addM, transpose, isSquareMatrix, dimension
+  , scale, swap, pivot
+  , addRow, subRow, addCol, subCol
+  , findPivot, forwardElim, gaussElim, gaussElimCorrect
   ) where
 
 import qualified Data.List as L
+import Data.Function (on)
 import Control.Monad (liftM)
+import Control.Arrow hiding ((<+>))
 import Test.QuickCheck
 
-import Algebra.Structures.IntegralDomain
+import Algebra.Structures.Field
 
+import Debug.Trace
 
 -------------------------------------------------------------------------------
 -- | Row vectors
@@ -23,15 +30,19 @@
 instance Show r => Show (Vector r) where
   show (Vec vs) = show vs
 
-instance (Ring r, Arbitrary r, Eq r) => Arbitrary (Vector r) where
+-- Generate vector of length 1-10
+instance Arbitrary r => Arbitrary (Vector r) where
   arbitrary = do n <- choose (1,10) :: Gen Int
                  liftM Vec $ gen n
     where
     gen 0 = return []
     gen n = do x <- arbitrary
                xs <- gen (n-1)
-               if x == zero then return (one:xs) else return (x:xs)
+               return (x:xs)
 
+instance Functor Vector where 
+  fmap f = Vec . map f . unVec
+
 {-
 instance Ring r => Ring (Vector r) where
   (Vec xs) <+> (Vec ys) | length xs == length ys = Vec (zipWith (<+>) xs ys)
@@ -62,7 +73,8 @@
     [] -> "[]" 
     xs -> init xs ++ "\n"
 
-instance (Eq r, Arbitrary r, Ring r) => Arbitrary (Matrix r) where
+-- Generate matrices with at most 10 rows
+instance Arbitrary r => Arbitrary (Matrix r) where
   arbitrary = do n <- choose (1,10) :: Gen Int
                  m <- choose (1,10) :: Gen Int
                  xs <- sequence [ liftM Vec (gen n) | _ <- [1..m]]
@@ -71,8 +83,10 @@
     gen 0 = return []
     gen n = do x <- arbitrary
                xs <- gen (n-1)
-               if x == zero then return (one:xs) else return (x:xs)
+               return (x:xs)
 
+instance Functor Matrix where
+  fmap f = M . map (fmap f) . unM
 
 -- | Construct a mxn matrix.
 matrix :: [[r]] -> Matrix r
@@ -96,11 +110,18 @@
 matrixToVector m | fst (dimension m) == 1 = head (unM m)
                  | otherwise              = error "matrixToVector: Bad dimension"
 
+(!!!) :: Matrix a -> (Int,Int) -> a
+m !!! (r,c) | r >= 0 && r < rows && c >= 0 && c < cols = unMVec m !! r !! c
+            | otherwise = error "!!!: Out of bounds"
+  where
+  (rows,cols) = dimension m 
+
 -- | Compute the dimension of a matrix.
 dimension :: Matrix r -> (Int, Int)
 dimension (M xs) | null xs   = (0,0)
                  | otherwise = (length xs, length (unVec (head xs)))
 
+
 isSquareMatrix :: Matrix r -> Bool
 isSquareMatrix (M xs) = all (== length xs) (map lengthVec xs)
 
@@ -110,7 +131,7 @@
 
 -- | Matrix addition.
 addM :: Ring r => Matrix r -> Matrix r -> Matrix r
-addM (M xs) (M ys) 
+addM (M xs) (M ys)
   | dimension (M xs) == dimension (M ys) = m
   | otherwise = error "Bad dimensions in matrix addition"
   where
@@ -118,33 +139,33 @@
 
 -- | Matrix multiplication.
 mulM :: Ring r => Matrix r -> Matrix r -> Matrix r
-mulM (M xs) (M ys) 
+mulM (M xs) (M ys)
   | snd (dimension (M xs)) == fst (dimension (M ys)) = m
   | otherwise = error "Bad dimensions in matrix multiplication"
     where
-    m = matrix [ [ foldr1 (<+>) (zipWith (<*>) x y) 
-                 | y <- L.transpose (map unVec ys) ]
+    m = matrix [ [ mulVec x y | y <- L.transpose (map unVec ys) ]
                | x <- map unVec xs ]
 
-
+mulVec xs ys | length xs == length ys = foldr (<+>) zero $ zipWith (<*>) xs ys
+             | otherwise = error "mulVec: Bad dimension"
 
 {-
 -- In order to do this the size of the matrix need to be encoded in the type
--- There is also a problem with the fact that it is not possible to add or 
+-- There is also a problem with the fact that it is not possible to add or
 -- multiply matrices with bad dimensions, so the generation of matrices has to be better...
 instance Ring r => Ring (Matrix r) where
   (<+>) = add
   (<*>) = mul
   neg (Vec xs d) = Vec [ map neg x | x <- xs ] d
-  zero  = undefined 
+  zero  = undefined
 -}
 
 -- | Construct a nxn identity matrix.
 identity :: IntegralDomain r => Int -> Matrix r
 identity n = matrix (xs 0)
   where
-  xs x | x == n    = [] 
-       | otherwise = (replicate x zero ++ [one] ++ 
+  xs x | x == n    = []
+       | otherwise = (replicate x zero ++ [one] ++
                       replicate (n-x-1) zero) : xs (x+1)
 
 -- Specification of identity.
@@ -155,3 +176,163 @@
 propRightIdentity :: (IntegralDomain r, Eq r) => Matrix r -> Bool
 propRightIdentity a = a == a `mulM` identity m
   where m = snd (dimension a)
+
+
+-------------------------------------------------------------------------------
+-- Operations on matrices.
+
+-- | Scale a row in a matrix.
+scale :: CommutativeRing a => Matrix a -> Int -> a -> Matrix a
+scale m r s
+  | 0 <= r && r < rows = matrix $ take r m' ++ map (s <*>) (m' !! r) : drop (r+1) m'
+  | otherwise = error "scale: Index out of bounds"
+  where
+  (rows,_) = dimension m
+  m'       = unMVec m
+
+-- Scaling does not affect dimension
+propScaleDimension :: (Arbitrary r, CommutativeRing r) => Matrix r -> Int -> r -> Bool
+propScaleDimension m r s = d == dimension (scale m (mod r rows) s)
+  where d@(rows,_) = dimension m
+
+-- | Swap two rows of a matrix.
+swap :: Matrix a -> Int -> Int -> Matrix a
+swap m i j
+  | 0 <= i && i <= r && 0 <= j && j <= r = matrix $ swap' m' i j
+  | otherwise = error "swap: Index out of bounds"
+  where
+  (r,_) = dimension m
+  m'    = unMVec m
+
+  swap' xs 0 0     = xs
+  swap' (x:xs) 0 j = (x:xs) !! j : take (j-1) xs ++ x : drop j xs
+  swap' xs i 0     = swap' xs 0 i
+  swap' (x:xs) i j = x : swap' xs (i-1) (j-1)
+
+-- Swapping does not affect dimension
+propSwapDimension :: Matrix () -> Int -> Int -> Bool
+propSwapDimension m i j = d == dimension (swap m (mod i r) (mod j r))
+  where d@(r,_) = dimension m
+
+-- Swap is itselfs identity.
+propSwapIdentity :: Matrix () -> Int -> Int -> Bool
+propSwapIdentity m i j = m == swap (swap m i' j') i' j'
+  where
+  d@(r,_) = dimension m
+  i'      = mod i r
+  j'      = mod j r
+
+
+-- Add the row-vector to the specified row of the matrix.
+addRow :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a
+addRow m row@(Vec xs) x
+  | 0 <= x && x < r = matrix $ take x m' ++
+                               zipWith (<+>) (m' !! x) xs :
+                               drop (x+1) m'
+  | c /= length xs  = error "addRow: Bad length of row"
+  | otherwise       = error "addRow: Bad row number"
+    where
+    (r,c) = dimension m
+    m'    = unMVec m
+
+propAddRowDimension :: (CommutativeRing a, Arbitrary a)
+                    => Matrix a -> Vector a -> Int -> Property
+propAddRowDimension m row@(Vec xs) r =
+  length xs == c ==> d == dimension (addRow m row (mod r r'))
+  where d@(r',c) = dimension m
+
+addCol :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a
+addCol m c x = transpose $ addRow (transpose m) c x
+
+subRow, subCol :: CommutativeRing a => Matrix a -> Vector a -> Int -> Matrix a
+subRow m (Vec xs) x = addRow m (Vec (map neg xs)) x
+subCol m (Vec xs) x = addCol m (Vec (map neg xs)) x
+
+-- Multiply the pivot row and add it to the target row.
+pivot :: CommutativeRing a => Matrix a -> a -> Int -> Int -> Matrix a
+pivot m s p t = addRow m (fmap (s <*>) (unM m !! p)) t
+
+-- Find first non-zero number below the pivot and return its value and row number
+-- given that it exists
+findPivot :: (CommutativeRing a, Eq a) => Matrix a -> (Int,Int) -> Maybe (a,Int)
+findPivot m (r,c) = safeHead $ filter ((/= zero) . fst) $ drop (r+1) $ zip (head $ drop c $ unMVec $ transpose m) [0..]
+  where
+  m' = unMVec m
+
+  safeHead []     = Nothing
+  safeHead (x:xs) = Just x
+
+fE :: (Field a, Eq a) => Matrix a -> Matrix a
+fE (M [])         = M []
+fE (M (Vec []:_)) = M []
+fE m     = case L.findIndices (/= zero) (map head xs) of
+  (i:is) -> case fE (cancelOut m [ (i,map head xs !! i) | i <- is ] (i,map head xs !! i)) of
+    ys -> matrix (xs !! i : map (zero :) (unMVec ys))
+  []     -> case fE (matrix (map tail xs)) of
+    ys -> matrix (map (zero:) (unMVec ys))
+  where
+  cancelOut :: (Field a, Eq a) => Matrix a -> [(Int,a)] -> (Int,a) -> Matrix a
+  cancelOut m [] (i,_)    = let xs = unMVec m in matrix $ map tail (L.delete (xs !! i) xs)
+  cancelOut m ((t,x):xs) (i,p) = cancelOut (pivot m (neg (x </> p)) i t) xs (i,p)
+
+  xs = unMVec m
+
+
+-- | Compute row echelon form of a system Ax=b.
+forwardElim :: (Field a, Eq a) => (Matrix a,Vector a) -> (Matrix a,Vector a)
+forwardElim (m,v) = fE m' (0,0)
+  where
+  -- fE takes the matrix to eliminate and the current row and column
+  fE :: (Field a, Eq a) => Matrix a -> (Int,Int) -> (Matrix a,Vector a)
+  fE (M []) _  = error "forwardElim: Empty input matrix"
+  fE m rc@(r,c)
+      -- The algorithm is done when it reaches the last column or row.
+    | c == mc || r == mr =
+      -- Decompose the matrix into A and b again
+      (matrix *** Vec) $ unzip $ map (init &&& last) $ unMVec m
+
+    | m !!! rc == zero   = case findPivot m rc of
+      -- If the pivot element is zero swap the pivot row with the first row
+      -- with a nonzero element in the pivot column.
+      Just (_,r') -> fE (swap m r r') rc
+      -- If all elements in the pivot column is zero the move right.
+      Nothing     -> fE m (r,c+1)
+
+    | m !!! rc /= one    =
+      -- Make the pivot element 1.
+      fE (scale m r (inv (m !!! rc))) rc
+
+    | otherwise          = case findPivot m rc of
+      -- Make the first nonzero element in the pivot row 0.
+      Just (v,r') -> fE (pivot m (neg v) r r') (r,c)
+      -- If all elements in the pivot column is zero then move down and right.
+      Nothing     -> fE m (r+1,c+1)
+
+  (mr,mc) = dimension m
+
+  -- Combine A and b to a matrix where the last column is b
+  m' = matrix $ [ r ++ [x] | (r,x) <- zip (unMVec m) (unVec v) ]
+
+
+-- | Perform "jordan"-step in Gauss-Jordan elimination. That is make every
+-- element above the diagonal zero. In other words compute the reduced
+-- echelon form of a matrix given that the input is in row echelon form.
+jordan :: (Field a, Eq a) => (Matrix a, Vector a) -> (Matrix a, Vector a)
+jordan (m, Vec ys) = case L.unzip (jordan' (zip (unMVec m) ys) (r-1)) of
+  (a,b) -> (matrix a, Vec b)
+  where
+  (r,_) = dimension m
+
+  jordan' [] _ = []
+  jordan' xs c =
+    jordan' [ (take c x ++ zero : drop (c+1) x, v <-> x !! c <*> snd (last xs))
+            | (x,v) <- init xs ] (c-1) ++ [last xs]
+
+
+-- | Gauss-Jordan elimination: Given A and B solve Ax=B.
+gaussElim :: (Field a, Eq a, Show a) => (Matrix a, Vector a) -> (Matrix a, Vector a)
+gaussElim = jordan . forwardElim
+
+gaussElimCorrect :: (Field a, Eq a, Arbitrary a, Show a) => (Matrix a, Vector a) -> Property
+gaussElimCorrect m@(a,b) = fst (dimension a) == lengthVec b && isSquareMatrix a ==>
+  matrixToVector (transpose (a `mulM` transpose (M [snd (gaussElim m)]))) == b
diff --git a/src/Algebra/Q.hs b/src/Algebra/Q.hs
--- a/src/Algebra/Q.hs
+++ b/src/Algebra/Q.hs
@@ -1,16 +1,13 @@
 {-# LANGUAGE TypeSynonymInstances #-}
 -- | Representation of rational numbers as the field of fractions of Z.
-module Algebra.Q 
-  ( Q
-  , toQ, toZ
-  ) where
+module Algebra.Q ( Q, toQ, toZ ) where
 
+import Data.Ratio (numerator, denominator)
 import Test.QuickCheck
--- import qualified Math.Algebra.Field.Base as A (Q(..)) 
--- import Data.Ratio 
 
 import Algebra.Structures.Field
-import Algebra.Structures.FieldOfFractions
+import Algebra.Structures.FieldOfFractions hiding (numerator, denominator)
+import Algebra.Structures.ExplicitUnits
 import Algebra.Z
 
 -------------------------------------------------------------------------------
@@ -26,12 +23,18 @@
   fromInteger      = toQ
 
 instance Fractional Q where
-  (/) = (</>)
-  fromRational = undefined
---   fromRational (a :% b) = reduce $ F (a,b)
+  (/)            = (</>)
+  fromRational x = 
+    reduce $ F (fromIntegral (numerator x), fromIntegral (denominator x))
 
 toQ :: Z -> Q
 toQ = toFieldOfFractions
 
 toZ :: Q -> Z
-toZ = fromFieldOfFractions   
+toZ = fromFieldOfFractions
+
+propFieldQ :: Q -> Q -> Q -> Property
+propFieldQ = propField
+
+instance ExplicitUnits Q where
+  unit a = if a == 0 then Nothing else Just (inv a)
diff --git a/src/Algebra/Structures/BezoutDomain.hs b/src/Algebra/Structures/BezoutDomain.hs
--- a/src/Algebra/Structures/BezoutDomain.hs
+++ b/src/Algebra/Structures/BezoutDomain.hs
@@ -1,18 +1,16 @@
 {-# LANGUAGE FlexibleInstances, UndecidableInstances #-}
--- | Representation of Bezout domains. That is non-Noetherian analogues of 
+-- | Representation of Bezout domains. That is non-Noetherian analogues of
 -- principal ideal domains. This means that all finitely generated ideals are
 -- principal.
 --
 module Algebra.Structures.BezoutDomain
-  ( BezoutDomain(..)
-  , propToPrincipal, propIsSameIdeal, propBezoutDomain
-  , dividesB, gcdB
-  , intersectionB, intersectionBWitness
-  , solveB
-  , crt
+  ( BezoutDomain(..), propBezoutDomain
+  , toPrincipal, propToPrincipal, propIsSameIdeal
+  , gcdB, intersectionB, intersectionBWitness
+  , solveB, crt
   ) where
 
-import Test.QuickCheck 
+import Test.QuickCheck
 
 import Algebra.Structures.IntegralDomain
 import Algebra.Structures.Coherent
@@ -23,21 +21,49 @@
 
 
 -------------------------------------------------------------------------------
--- | Bezout domains
--- 
--- Compute a principal ideal from another ideal. Also give witness that the
--- principal ideal is equal to the first ideal.
---
--- toPrincipal \<a_1,...,a_n> = (\<a>,u_i,v_i)
---   where
---
---   sum (u_i * a_i) = a
---
---   a_i = v_i * a
---
+{- | Bezout domains
+
+Has a Bezout function which given a and b give g, a1, b1, x and y such that:
+
+ - g = gcd(a,b)
+
+ - a = g * a1 and b = g * b1
+
+ - g = a * x + b * y
+
+-}
 class IntegralDomain a => BezoutDomain a where
-  toPrincipal :: Ideal a -> (Ideal a,[a],[a])
+  bezout :: a -> a -> (a,a,a,a,a)
 
+propBezoutDomain :: (BezoutDomain a, Eq a) => a -> a -> Property
+propBezoutDomain a b = 
+  let (g,a1,b1,x,y) = bezout a b
+  in if a == g <*> a1 && b == g <*> b1 && a <*> x <+> b <*> y == g
+        then propIntegralDomain a b b
+        else whenFail (print "propBezoutDomain") False
+
+
+{- | Compute a principal ideal from another ideal. Also give witness that the
+principal ideal is equal to the first ideal.
+
+toPrincipal \<a_1,...,a_n> = (\<a>,u_i,v_i)
+  where
+
+  sum (u_i * a_i) = a
+
+  a_i = v_i * a
+-}
+toPrincipal :: BezoutDomain a => Ideal a -> (Ideal a,[a],[a])
+toPrincipal (Id [])    = error "toPrincipal: Empty input"
+toPrincipal (Id [a])   = (Id [a],[one],[one])
+toPrincipal (Id [a,b]) = 
+  let (g,a1,b1,x,y) = bezout a b
+  in (Id [g],[x,y],[a1,b1])
+toPrincipal (Id (a:xs)) = 
+  let (Id [g],us,vs) = toPrincipal (Id xs)
+      (g',a1,b1,x,y) = bezout a g
+  in (Id [g'],x : map (y <*>) us,a1 : map (b1 <*>) vs)
+
 -- | Test that the generated ideal is principal.
 propToPrincipal :: (BezoutDomain a, Eq a) => Ideal a -> Bool
 propToPrincipal = isPrincipal . (\(a,_,_) -> a) . toPrincipal
@@ -45,24 +71,12 @@
 -- | Test that the generated ideal generate the same elements as the given.
 propIsSameIdeal :: (BezoutDomain a, Eq a) => Ideal a -> Bool
 propIsSameIdeal (Id as) =
-  let (Id [a], us, vs) = toPrincipal (Id as) 
-  in a == foldr1 (<+>) (zipWith (<*>) as us) 
+  let (Id [a], us, vs) = toPrincipal (Id as)
+  in a == foldr1 (<+>) (zipWith (<*>) as us)
   && and [ ai == a <*> vi | (ai,vi) <- zip as vs ]
   && length us == l_as && length vs == l_as
   where l_as = length as
 
-propBezoutDomain :: (BezoutDomain a, Eq a) => Ideal a -> a -> a -> a -> Property
-propBezoutDomain id@(Id xs) a b c = zero `notElem` xs ==> 
-  if propToPrincipal id
-     then if propIsSameIdeal id
-             then propIntegralDomain a b c 
-             else whenFail (print "propIsSameIdeal") False
-     else whenFail (print "propToPrincipal") False
-
-dividesB :: (BezoutDomain a, Eq a) => a -> a -> Bool
-dividesB a b = a == x || a == neg x
-    where (Id [x],_,_) = toPrincipal (Id [a,b])
-
 -- TODO: Add error cases...
 gcdB :: BezoutDomain a => a -> a -> a
 gcdB a b = g
@@ -72,24 +86,28 @@
 -- Euclidean domain -> Bezout domain
 
 instance (EuclideanDomain a, Eq a) => BezoutDomain a where
-  toPrincipal (Id [x]) = (Id [x], [one], [one])
-  toPrincipal (Id xs)  = (Id [a], as, [ quotient ai a | ai <- xs ])
-    where
-    a  = genEuclidAlg xs
-    as = genExtendedEuclidAlg xs
+  bezout a b
+    | b == zero        = (a,one,zero,one,zero)
+    | a == zero        = (b,zero,one,zero,one)
+    | norm a <= norm b = let (q,r)           = quotientRemainder b a
+                             (g,a1,r1,u',v') = bezout a r
+                         in (g,a1,r1<+>(q<*>a1),u'<->(q<*>v'),v')
+    | otherwise        = let (q,r)           = quotientRemainder a b
+                             (g,b1,r1,u',v') = bezout b r
+                         in (g,r1<+>(q<*>b1),b1,v',u'<->(v'<*>q))
 
 
 -------------------------------------------------------------------------------
 -- | Intersection of ideals with witness.
--- 
--- If one of the ideals is the zero ideal then the intersection is the zero 
+--
+-- If one of the ideals is the zero ideal then the intersection is the zero
 -- ideal.
--- 
-intersectionBWitness :: (BezoutDomain a, Eq a) 
-              => Ideal a 
-              -> Ideal a 
+
+intersectionBWitness :: (BezoutDomain a, Eq a)
+              => Ideal a
+              -> Ideal a
               -> (Ideal a, [[a]], [[a]])
-intersectionBWitness (Id xs) (Id ys) 
+intersectionBWitness (Id xs) (Id ys)
   | xs' == [] = zeroIdealWitnesses xs ys
   | ys' == [] = zeroIdealWitnesses xs ys
   | otherwise = (Id [l], [handleZero xs as], [handleZero ys bs])
@@ -97,7 +115,7 @@
   xs'            = filter (/= zero) xs
   ys'            = filter (/= zero) ys
 
-  (Id [a],us1,vs1) = toPrincipal (Id xs') 
+  (Id [a],us1,vs1) = toPrincipal (Id xs')
   (Id [b],us2,vs2) = toPrincipal (Id ys')
 
   (Id [g],[u1,u2],[v1,v2]) = toPrincipal (Id [a,b])
@@ -105,14 +123,14 @@
   l  = g <*> v1 <*> v2
   as = map (v2 <*>) us1
   bs = map (v1 <*>) us2
-  
+
 -- Handle the zeroes specially. If the first element in xs is a zero
--- then the witness should be zero otherwise use the computed witness. 
+-- then the witness should be zero otherwise use the computed witness.
 handleZero :: (Ring a, Eq a) => [a] -> [a] -> [a]
-handleZero xs [] 
+handleZero xs []
   | all (==zero) xs = xs
   | otherwise       = error "intersectionB: This should be impossible"
-handleZero (x:xs) (a:as) 
+handleZero (x:xs) (a:as)
   | x == zero = zero : handleZero xs (a:as)
   | otherwise = a    : handleZero xs as
 handleZero [] _  = error "intersectionB: This should be impossible"
@@ -134,13 +152,13 @@
 
 -------------------------------------------------------------------------------
 -- | Strongly discreteness for Bezout domains
--- 
+--
 -- Given x, compute as such that x = sum (a_i * x_i)
 --
 instance (BezoutDomain a, Eq a) => StronglyDiscrete a where
   member x (Id xs) | x == zero = Just (replicate (length xs) zero)
-                   | otherwise = if a == g 
-                                    then Just witness 
+                   | otherwise = if a == g
+                                    then Just witness
                                     else Nothing
     where
     -- (<g>, as, bs)   = <x1,...,xn>
@@ -148,30 +166,27 @@
     -- x_i             = b_i * g
     (Id [g], as, bs) = toPrincipal (Id (filter (/= zero) xs))
     (Id [a], _,[q1,q2]) = toPrincipal (Id [x,g])
-    
+
     -- x = qg = q (sum (ai * xi)) = sum (q * ai * xi)
     witness = handleZero xs (map (q1 <*>) as)
 
 
 -------------------------------------------------------------------------------
 -- | Chinese remainder theorem
--- 
+--
 -- Given a_1,...,a_n and m_1,...,m_n such that gcd(m_i,m_j) = 1.
 -- Let m = m_1*...*m_n compute a such that:
--- 
--- (1) a = a_i (mod m_i) 
--- 
--- (2) If b is such that
---      
---        b = a_i (mod m_i)
--- 
---     then a = b (mod m)
 --
+-- (1) a = a_i (mod m_i)
+--
+-- (2) If b is such that b = a_i (mod m_i) then a = b (mod m)
+--
 -- The function return (a,m).
+
 crt :: (BezoutDomain a, Eq a) => [a] -> [a] -> (a,a)
 crt as ms
   | length as /= length ms = error "crt: Input lists need to have same length"
-  | not (and [ gcdB m1 m2 == one | m1 <- ms, m2 <- ms, m1 /= m2 ]) = 
+  | not (and [ gcdB m1 m2 == one | m1 <- ms, m2 <- ms, m1 /= m2 ]) =
       error "crt: All ms need to be relatively prime"
   | otherwise = crt' as ms
   where
diff --git a/src/Algebra/Structures/Coherent.hs b/src/Algebra/Structures/Coherent.hs
--- a/src/Algebra/Structures/Coherent.hs
+++ b/src/Algebra/Structures/Coherent.hs
@@ -1,9 +1,9 @@
 -- | Representation of coherent rings. Traditionally a ring is coherent if every
--- finitely generated ideal is finitely presented. This means that it is 
+-- 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
+  , propCoherent, isSolution
   , solveMxN, propSolveMxN
   , solveWithIntersection
   , solveGeneralEquation, propSolveGeneralEquation
@@ -25,7 +25,7 @@
 --
 --   MX=0   \<-\>  \exists Y. X=LY
 --
--- that is, iff we can generate the solutions of any linear homogeous system 
+-- 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
@@ -34,6 +34,7 @@
 class IntegralDomain a => Coherent a where
   solve :: Vector a -> Matrix a
 
+-- | Test that the second matrix is a solution to the first.
 isSolution :: (CommutativeRing a, Eq a) => Matrix a -> Matrix a -> Bool
 isSolution m sol = all (==zero) (concat (unMVec (m `mulM` sol)))
 
@@ -45,14 +46,14 @@
 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. 
+  -- 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 
+  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))
 
@@ -64,38 +65,34 @@
 
 -------------------------------------------------------------------------------
 -- | Intersection computable -> Coherence.
--- 
--- Proof that if there is an algorithm to compute a f.g. set of generators for 
+--
+-- 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. 
+-- 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 \]
+-- 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
+--     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 
+                      => Vector a
                       -> (Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]]))
                       -> Matrix a
-solveWithIntersection (Vec xs) int = transpose $ matrix $ solveInt xs 
+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 ]
+                             -- 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 
+    | x == zero             = (one : replicate (length xs) zero) : (map (zero:) $ solveInt xs)
     | isSameIdeal int as bs = s ++ m'
     | otherwise             = error "solveInt: This does not compute the intersection"
       where
@@ -105,7 +102,7 @@
       -- 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
@@ -114,11 +111,11 @@
 -------------------------------------------------------------------------------
 -- | Strongly discrete coherent rings.
 --
--- If the ring is strongly discrete and coherent then we can solve arbitrary 
--- equations of the type AX=b. 
+-- 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 = 
+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)))
@@ -126,9 +123,9 @@
     Nothing -> Nothing
 
 
-propSolveGeneralEquation :: (Coherent a, StronglyDiscrete a, Eq a) 
-                         => Vector a 
-                         -> a 
+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
@@ -140,24 +137,24 @@
 
 -- | 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 
+-- 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) 
+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)) = 
+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 
+  -- 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 
+  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'
@@ -165,10 +162,10 @@
                           then solveGeneral' (Just m'') ls as ((l,a):old)
                           else Nothing
          Nothing -> Nothing
-  solveGeneral' _ _ _ _ = error "solveGeneral: Bad input"      
+  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 :: (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
diff --git a/src/Algebra/Structures/CommutativeRing.hs b/src/Algebra/Structures/CommutativeRing.hs
--- a/src/Algebra/Structures/CommutativeRing.hs
+++ b/src/Algebra/Structures/CommutativeRing.hs
@@ -1,3 +1,5 @@
+-- | Structure for commutative rings. 
+--
 module Algebra.Structures.CommutativeRing
   ( module Algebra.Structures.Ring
   , CommutativeRing(..)
diff --git a/src/Algebra/Structures/EuclideanDomain.hs b/src/Algebra/Structures/EuclideanDomain.hs
--- a/src/Algebra/Structures/EuclideanDomain.hs
+++ b/src/Algebra/Structures/EuclideanDomain.hs
@@ -1,10 +1,10 @@
--- | Representation of Euclidean domains. That is integral domains with an 
+-- | Representation of Euclidean domains. That is integral domains with an
 -- Euclidean functions and decidable division.
 --
-module Algebra.Structures.EuclideanDomain 
+module Algebra.Structures.EuclideanDomain
   ( EuclideanDomain(..)
-  , propD, propQuotRem, propEuclideanDomain
-  , modulo, quotient, divides 
+  , propNorm, propQuotRem, propEuclideanDomain
+  , modulo, quotient, divides
   , euclidAlg, genEuclidAlg
   , lcmE, genLcmE
   , extendedEuclidAlg, genExtendedEuclidAlg
@@ -20,25 +20,27 @@
 -------------------------------------------------------------------------------
 -- | Euclidean domains
 --
--- Given a and b compute (q,r) such that a = bq + r and r = 0 || d r < d b. 
--- Where d is the Euclidean function.
+-- Given a and b compute (q,r) such that a = bq + r and r = 0 || norm r < norm b.
+-- Where norm is the Euclidean function.
 
 class IntegralDomain a => EuclideanDomain a where
-  d :: a -> Integer
+  norm :: a -> Integer
   quotientRemainder :: a -> a -> (a,a)
 
 -- | Check both that |a| <= |ab| and |a| >= 0 for all a,b.
-propD :: (EuclideanDomain a, Eq a) => a -> a -> Bool
-propD a b = 
-  a == zero || b == zero || (d a <= d (a <*> b) && d a >= 0 && d b >= 0)
+propNorm :: (EuclideanDomain a, Eq a) => a -> a -> Bool
+propNorm a b =
+  a == zero || b == zero || 
+  (norm a <= norm (a <*> b) && norm a >= 0 && norm b >= 0)
 
 propQuotRem :: (EuclideanDomain a, Eq a) => a -> a -> Bool
-propQuotRem a b = b == zero || (a == b <*> q <+> r && (r == zero || d r < d b))
-  where (q,r) = quotientRemainder a b 
+propQuotRem a b = 
+  b == zero || (a == b <*> q <+> r && (r == zero || norm r < norm b))
+    where (q,r) = quotientRemainder a b
 
 propEuclideanDomain :: (EuclideanDomain a, Eq a) => a -> a -> a -> Property
 propEuclideanDomain a b c =
-  if propD a b 
+  if propNorm a b
      then if propQuotRem a b
              then propIntegralDomain a b c
              else whenFail (print "propQuotRem") False
@@ -59,7 +61,7 @@
 
 -- | The Euclidean algorithm for calculating the GCD of a and b.
 euclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> a
-euclidAlg a b | a == zero && b == zero = error "GCD of 0 and 0 is undefined"
+euclidAlg a b | a == zero && b == zero = zero
               | b == zero = a
               | otherwise = euclidAlg b (a `modulo` b)
 
@@ -75,18 +77,20 @@
 genLcmE :: (EuclideanDomain a, Eq a) => [a] -> a
 genLcmE xs = quotient (foldr1 (<*>) xs) (genEuclidAlg xs)
 
--- | The extended Euclidean algorithm. 
--- 
+-- | The extended Euclidean algorithm.
+--
 -- Computes x and y in ax + by = gcd(a,b).
--- 
+--
 extendedEuclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> (a,a)
-extendedEuclidAlg a b | modulo a b == zero = (zero,one)
+extendedEuclidAlg a b | a == zero = (zero,one)
+                      | b == zero = (one,zero)
+--                      | modulo a b == zero = (one,zero)
                       | otherwise          = (y, x <-> y <*> (a `quotient` b))
   where (x,y) = extendedEuclidAlg b (a `modulo` b)
 
 -- Specification of extended Euclidean algorithm.
 propExtendedEuclidAlg :: (EuclideanDomain a, Eq a) => a -> a -> Property
-propExtendedEuclidAlg a b = a /= zero && b /= zero ==> 
+propExtendedEuclidAlg a b = a /= zero && b /= zero ==>
   let (x,y) = extendedEuclidAlg a b in a <*> x <+> b <*> y == euclidAlg a b
 
 -- | Generalized extended Euclidean algorithm.
@@ -99,7 +103,7 @@
   let (x,y) = extendedEuclidAlg (genEuclidAlg (init xs)) (last xs)
   in map (x<*>) (genExtendedEuclidAlg (init xs)) ++ [y]
 
--- Specification of generalized extended Euclidean algorithm. 
+-- Specification of generalized extended Euclidean algorithm.
 propGenExtEuclidAlg :: (EuclideanDomain a, Eq a) => [a] -> Property
-propGenExtEuclidAlg xs = all (/= zero) xs && length xs >= 2 ==> 
+propGenExtEuclidAlg xs = all (/= zero) xs && length xs >= 2 ==>
   foldr (<+>) zero (zipWith (<*>) (genExtendedEuclidAlg xs) xs) == genEuclidAlg xs
diff --git a/src/Algebra/Structures/ExplicitUnits.hs b/src/Algebra/Structures/ExplicitUnits.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/Structures/ExplicitUnits.hs
@@ -0,0 +1,34 @@
+-- | Structure of rings with explicit units. 
+module Algebra.Structures.ExplicitUnits
+  ( ExplicitUnits(..)
+  , propUnit, isUnit, (%|), (~=)
+  ) where 
+
+import Algebra.Structures.IntegralDomain
+import Algebra.Structures.GCDDomain
+
+infix 5 %|
+infix 4 ~=
+
+-- | A ring has explicit units if there is a function that can test if an
+-- element is invertible and if this is the case give the inverse. 
+class IntegralDomain a => ExplicitUnits a where
+  unit :: a -> Maybe a
+
+propUnit :: (ExplicitUnits a, Eq a) => a -> Bool
+propUnit a = case unit a of
+  Just a' -> a <*> a' == one
+  Nothing -> True
+
+-- | An element is a unit if it is invertible. 
+isUnit :: ExplicitUnits a => a -> Bool
+isUnit = maybe False (const True) . unit
+
+-- | Decidable units is sufficient to decide divisibility in GCD domains. 
+(%|) :: (ExplicitUnits a, GCDDomain a) => a -> a -> Bool
+a %| b = let (g,x,y) = gcd' a b
+         in isUnit x
+
+-- | Test for associatedness, i.e. a ~ b iff a | b /\\ b | a.
+(~=) :: (ExplicitUnits a, GCDDomain a) => a -> a -> Bool
+a ~= b = a %| b && b %| a 
diff --git a/src/Algebra/Structures/Field.hs b/src/Algebra/Structures/Field.hs
--- a/src/Algebra/Structures/Field.hs
+++ b/src/Algebra/Structures/Field.hs
@@ -1,3 +1,4 @@
+-- | Structure for fields.
 module Algebra.Structures.Field
   ( module Algebra.Structures.IntegralDomain
   , Field(inv)
@@ -9,7 +10,6 @@
 
 import Algebra.Structures.Ring
 import Algebra.Structures.IntegralDomain
-
 
 infixl 7 </>
 
diff --git a/src/Algebra/Structures/FieldOfFractions.hs b/src/Algebra/Structures/FieldOfFractions.hs
--- a/src/Algebra/Structures/FieldOfFractions.hs
+++ b/src/Algebra/Structures/FieldOfFractions.hs
@@ -2,6 +2,7 @@
 -- domain is that we only want to work over reduced quotients.
 module Algebra.Structures.FieldOfFractions
   ( FieldOfFractions(..)
+  , numerator, denominator
   , toFieldOfFractions, fromFieldOfFractions
   , reduce, propReduce
   ) where
@@ -17,7 +18,11 @@
 
 newtype GCDDomain a => FieldOfFractions a = F (a,a)
 
+numerator, denominator :: GCDDomain a => FieldOfFractions a -> a
+numerator (F (x,_))   = x
+denominator (F (_,x)) = x
 
+
 --------------------------------------------------------------------------------
 -- Instances
 
@@ -33,13 +38,13 @@
     b <- arbitrary
     if b == zero 
        then return $ F (a,one)
-       else return $ F (a,b)
+       else return $ reduce $ F (a,b)
 
 instance (GCDDomain a, Eq a) => Eq (FieldOfFractions a) where
-  f == g = a <*> d == b <*> c
-    where
-    F (a,b) = reduce f
-    F (c,d) = reduce g
+  (F (a,b)) == (F (c,d)) = a <*> d == b <*> c
+--    where
+--    F (a,b) = reduce f
+--    F (c,d) = reduce g
 
 instance (GCDDomain a, Eq a) => Ring (FieldOfFractions a) where
   (F (a,b)) <+> (F (c,d)) = reduce (F (a <*> d <+> c <*> b,b <*> d))
@@ -52,7 +57,7 @@
 instance (GCDDomain a, Eq a) => IntegralDomain (FieldOfFractions a)
 
 instance (GCDDomain a, Eq a) => Field (FieldOfFractions a) where
-  inv (F (a,b)) | b /= zero && a /= zero = reduce $ F (b,a)
+  inv (F (a,b)) | a /= zero && b /= zero = reduce $ F (b,a)
                 | otherwise = error "FieldOfFraction: Division by zero"
 
 
@@ -66,17 +71,15 @@
 -- | Extract a value from the field of fractions. This is only possible if the
 -- divisor is one.
 fromFieldOfFractions :: (GCDDomain a, Eq a) => FieldOfFractions a -> a
-fromFieldOfFractions (F (a,b)) 
-  | b == one  = a
-  | otherwise = error "FieldOfFractions: Can't extract value"
+fromFieldOfFractions x | b' == one = a'
+                       | otherwise = error "fromFieldOfFractions: Division by zero"
+  where F (a',b') = reduce x 
 
 -- | Reduce an element.
 reduce :: (GCDDomain a, Eq a) => FieldOfFractions a -> FieldOfFractions a
-reduce (F (a,b)) | b == zero = error "FieldOfFractions: Division by zero"
+reduce (F (a,b)) | b == zero = error "reduce: Division by zero"
                  | a == zero = F (zero,one)
-                 | otherwise = if g == one
-                                  then F (a,b)
-                                  else F (x,y)
+                 | otherwise = F (x,y)
   where
   (g,x,y) = gcd' a b
 
diff --git a/src/Algebra/Structures/GCDDomain.hs b/src/Algebra/Structures/GCDDomain.hs
--- a/src/Algebra/Structures/GCDDomain.hs
+++ b/src/Algebra/Structures/GCDDomain.hs
@@ -8,6 +8,7 @@
 module Algebra.Structures.GCDDomain 
   ( GCDDomain(gcd')
   , propGCD, propGCDDomain
+  , ggcd
   ) where
 
 import Test.QuickCheck
@@ -16,6 +17,7 @@
 import Algebra.Structures.BezoutDomain
 import Algebra.Ideal
 
+-- infix 4 ~~
 
 -------------------------------------------------------------------------------
 -- | GCD domains
@@ -41,10 +43,27 @@
                          then propIntegralDomain a b c
                          else whenFail (print "propGCD") False
 
--- This can be used to compute gcd of a list of non-zero elements
--- genGCD :: ?
--- genGCD = ?
+fst3 :: (a,b,c) -> a
+fst3 (x,_,_) = x
 
+-- Generalized greatest common divisor, computes the gcd of a list of elements.
+ggcd :: GCDDomain a => [a] -> a
+ggcd []     = error "ggcd: Can't compute ggcd of the empty list"
+ggcd (x:xs) = foldr (\a b -> fst3 (gcd' a b)) x xs
+
 instance BezoutDomain a => GCDDomain a where
   gcd' a b = (g,x,y)
-    where (Id [g],_,[x,y]) = toPrincipal (Id [a,b])
+   where (Id [g],_,[x,y]) = toPrincipal (Id [a,b])
+
+{-
+class IntegralDomain a => DecidableUnits a where
+  unit :: a -> Maybe a -- Just x = the inverse
+
+-- Divisibility is decidable if A is a gcd domain with decidable units
+divides :: (GCDDomain a, DecidableUnits a) => a -> a -> Bool
+divides a b = unit g
+  where (g,_,_) = gcd' a b
+
+(~~) :: (GCDDomain a, DecidableUnits a) => a -> a -> Bool
+(~~) = divides
+-}
diff --git a/src/Algebra/Structures/Module.hs b/src/Algebra/Structures/Module.hs
--- a/src/Algebra/Structures/Module.hs
+++ b/src/Algebra/Structures/Module.hs
@@ -1,24 +1,26 @@
 {-# LANGUAGE MultiParamTypeClasses, TypeSynonymInstances, FlexibleInstances #-}
+-- | R-modules. 
 module Algebra.Structures.Module 
-  ( Module(..), (<*)
+  ( Module((*>)), (<*)
   , propScalarMul, propScalarAdd, propScalarAssoc, propModule
   ) where
 
 import Algebra.Structures.Group as G
-import Algebra.Structures.CommutativeRing as R -- hiding ((<*),(*>))
+import Algebra.Structures.CommutativeRing as R
 import Algebra.Z
 import Algebra.Zn
 
 import Test.QuickCheck
 
 infixl 7 *>
-infixl 7 <*
+infixr 7 <*
 
 -- Consider only the commutative case, it would be possible to implement left
 -- and right modules instead.
 
--- A module over a commutative ring r.
+-- | Module over a commutative ring r.
 class (CommutativeRing r, AbelianGroup m) => Module r m where
+  -- | Scalar multiplication.
   (*>) :: r -> m -> m
 
 propScalarMul :: (Module r m, Eq m) => r -> m -> m -> Bool
@@ -38,11 +40,11 @@
     (_,False,_)      -> whenFail (print "propScalarAdd") False
     (_,_,False)      -> whenFail (print "propScalarAssoc") False
 
--- Since the ring is commutative we can turn this around.
+-- | Since the ring is commutative we can turn the scalar multiplication around.
 (<*) :: Module r m => m -> r -> m
 (<*) = flip (*>)
 
--- Z-module
+-- | Z-module structure.
 instance AbelianGroup m => Module Z m where
   n *> x | n > 0  = sumGroup (replicate (fromInteger n) x)
          | n == 0 = G.zero
diff --git a/src/Algebra/Structures/PruferDomain.hs b/src/Algebra/Structures/PruferDomain.hs
--- a/src/Algebra/Structures/PruferDomain.hs
+++ b/src/Algebra/Structures/PruferDomain.hs
@@ -25,7 +25,12 @@
 
 
 -------------------------------------------------------------------------------
--- | Prufer domain
+-- | Given a and b it computes u, v and we such that:
+-- 
+--  (1) au = bv
+--
+--  (2) b(1-u) = aw
+--
 class IntegralDomain a => PruferDomain a where
   --         a    b     u v w
   calcUVW :: a -> a -> (a,a,a)
@@ -59,7 +64,7 @@
 
 
 --------------------------------------------------------------------------------
--- | Bezout domain -> Prüfer domain
+-- | Bezout domain -> Prufer domain
 --
 {-
 Prufer: forall a b exists u v w t.  u+t = 1 &  ua = vb & wa = tb
@@ -179,7 +184,7 @@
 
 -- | Compute the intersection of I and J by:
 --       
---       (I ∩ J)(I + J) = IJ  => (I ∩ J)(I + J)(I + J)' = IJ(I + J)'
+--       (I \\cap J)(I + J) = IJ  => (I \\cap J)(I + J)(I + J)' = IJ(I + J)'
 --
 intersectionPDWitness :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
 intersectionPDWitness (Id is) (Id js) = (int,wis,wjs)
diff --git a/src/Algebra/Structures/Ring.hs b/src/Algebra/Structures/Ring.hs
--- a/src/Algebra/Structures/Ring.hs
+++ b/src/Algebra/Structures/Ring.hs
@@ -6,6 +6,7 @@
   , propRing
   , (<->), (<^>) -- , (*>), (<*)
   , sumRing, productRing
+--  , (~~)
   ) where
 
 import Test.QuickCheck
@@ -17,7 +18,7 @@
 -- infixl 7 <*
 infixl 6 <+>
 infixl 6 <->
-
+--infix  4 ~~
 
 -------------------------------------------------------------------------------
 -- | Definition of rings.
@@ -111,6 +112,10 @@
 x <^> y = if y < 0 
              then error "<^>: Input should be positive"
              else x <*> x <^> (y-1)
+
+-- | Check if a == b or -a == b or a == -b or -a == -b
+-- (~~) :: (Ring a, Eq a) => a -> a -> Bool
+-- x ~~ y = x == y || neg x == y || x == neg y || neg x == neg y
 
 {-
 -- | Multiply from left with an integer; n *> x means x + x + ... + x, n times.
diff --git a/src/Algebra/UPoly.hs b/src/Algebra/UPoly.hs
--- a/src/Algebra/UPoly.hs
+++ b/src/Algebra/UPoly.hs
@@ -1,11 +1,14 @@
 {-# LANGUAGE ScopedTypeVariables, FlexibleContexts #-}
 -- | Univariate polynomials parametrised by the variable name.
-module Algebra.UPoly 
+module Algebra.UPoly
   ( UPoly(..)
   , deg
   , Qx, x
   , toUPoly, monomial
   , lt, deriv
+  , cont, isPrimitive
+  , toPrimitive, propToPrimitive, gaussLemma
+  , gcdUPolyWitness
   , sqfr, sqfrDec
   ) where
 
@@ -13,17 +16,22 @@
 import Test.QuickCheck
 import Control.Monad (liftM)
 
-import Algebra.TypeChar.Char hiding (Q)
+import Algebra.TypeChar.Char hiding (Z,Q)
 import Algebra.Structures.Field
-import Algebra.Structures.BezoutDomain
+import Algebra.Structures.FieldOfFractions
 import Algebra.Structures.EuclideanDomain
+import Algebra.Structures.ExplicitUnits
+import Algebra.Structures.BezoutDomain
+import Algebra.Structures.GCDDomain
 import Algebra.Structures.PruferDomain
-import Algebra.Structures.StronglyDiscrete
+-- import Algebra.Structures.StronglyDiscrete
+
 import Algebra.Ideal
+import Algebra.Z
 import Algebra.Q
 
--- | Polynomials over a commutative ring, indexed by a phantom type x that 
--- denote the name of the variable that the polynomial is over. For example 
+-- | Polynomials over a commutative ring, indexed by a phantom type x that
+-- denote the name of the variable that the polynomial is over. For example
 -- UPoly Q X_ is Q[x] and UPoly Q T_ is Q[t].
 newtype CommutativeRing r => UPoly r x = UP [r]
   deriving (Eq,Ord)
@@ -44,7 +52,11 @@
 toUPoly :: (CommutativeRing r, Eq r) => [r] -> UPoly r x
 toUPoly = UP . reverse . dropWhile (==zero) . reverse
 
--- | Take an element of the ring and the degree of the desired monomial, for 
+isZeroPoly :: CommutativeRing r => UPoly r x -> Bool
+isZeroPoly (UP []) = True
+isZeroPoly _       = False
+
+-- | Take an element of the ring and the degree of the desired monomial, for
 -- example: monomial 3 7 = 3x^7
 monomial :: CommutativeRing r => r -> Integer -> UPoly r x
 monomial a i = UP $ replicate (fromInteger i) zero ++ [a]
@@ -67,7 +79,7 @@
 
 instance (CommutativeRing r, Eq r, Show r, Show x) => Show (UPoly r x) where
   show (UP []) = "0"
-  show (UP ps) = init $ fixSign $ concat 
+  show (UP ps) = init $ fixSign $ concat
                   [ show' (show (undefined :: x)) p n
                   | (p,n) <- zip ps [0..]
                   , p /= zero ]
@@ -75,10 +87,10 @@
     show' :: (CommutativeRing r, Show r) => String -> r -> Integer -> String
     show' x p 0 = show p ++ "+"
     show' x p 1 = if p == one then x ++ "+" else show p ++ x ++ "+"
-    show' x p n = if p == one  
-                     then x ++ "^" ++ show n ++ "+" 
+    show' x p n = if p == one
+                     then x ++ "^" ++ show n ++ "+"
                      else show p ++ x ++ "^" ++ show n ++ "+"
-    
+
     fixSign []  = []
     fixSign [x] = [x]
     fixSign ('+':'-':xs) = '-' : fixSign xs
@@ -89,9 +101,9 @@
 
 -- Addition of polynomials.
 addUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x
-addUP (UP ps) (UP qs) | length ps >= length qs = add' ps qs 
+addUP (UP ps) (UP qs) | length ps >= length qs = add' ps qs
                       | otherwise              = add' qs ps
-  where add' a b = toUPoly $ zipWith (<+>) a b ++ drop (length b) a 
+  where add' a b = toUPoly $ zipWith (<+>) a b ++ drop (length b) a
 
 -- Multiplication of polynomials.
 mulUP :: (CommutativeRing r, Eq r) => UPoly r x -> UPoly r x -> UPoly r x
@@ -99,7 +111,7 @@
   where
   m ps qs r | r > length ps + length qs - 2 = []
             | otherwise = c r 0 (length ps-1) (length qs-1) : m ps qs (r+1)
-  
+
   c (-1) _ _ _ = zero
   c r k m n | r > m || k > n = c (r-1) (k+1) m n
             | otherwise      = ps !! r <*> qs !! k <+> c (r-1) (k+1) m n
@@ -115,7 +127,7 @@
   (+)    = (<+>)
   (-)    = (<->)
   (*)    = (<*>)
-  abs    = fromInteger . d
+  abs    = fromInteger . norm
   signum = undefined -- Is it possible to define this?
   fromInteger x = UP [fromInteger x]
 
@@ -124,31 +136,151 @@
 
 -- Polynomial rings are Euclidean.
 instance (Field k, Eq k) => EuclideanDomain (UPoly k x) where
-  d (UP ps)             = fromIntegral (length ps) - 1
+  norm (UP ps)             = fromIntegral (length ps) - 1
   quotientRemainder f g = qr zero f
     where
     -- This is the division algorithm in k[x]. Page 39 in Cox.
-    qr q r | d g <= d r = qr (q <+> monomial (lt r </> lt g) (d r - d g))
-                            (r <-> monomial (lt r </> lt g) (d r - d g) <*> g)
+    qr q r | norm g <= norm r = 
+              qr (q <+> monomial (lt r </> lt g) (norm r - norm g))
+                 (r <-> monomial (lt r </> lt g) (norm r - norm g) <*> g)
            | otherwise = (q,r)
 
 instance (Field k, Eq k) => PruferDomain (UPoly k x) where
   calcUVW = calcUVW_B
 
+instance (ExplicitUnits a, Eq a) => ExplicitUnits (UPoly a x) where
+  unit (UP [a]) = case unit a of
+    Just a' -> Just (UP [a'])
+    Nothing -> Nothing
+  unit _        = Nothing
+
 -- Now that we know that the polynomial ring k[x] is a Bezout domain it is
 -- possible to implement membership in an ideal of k[x]. f is a member of the
 -- ideal <f1,...,fn> if the rest is zero when dividing f with the principal
 -- ideal <h>.
 -- instance (Field k, Eq k, Show x) => StronglyDiscrete (UPoly k x) where
---  member p ps = modulo p h == zero 
+--  member p ps = modulo p h == zero
 --    where Id [h] = (\(a,_,_) -> a) $ toPrincipal ps
 
+
+-------------------------------------------------------------------------------
+-- Proof that if A is a GCD domain then A[x] is a GCD domain following
+-- section 4.4 in A course in constructive algebra.
+
+-- Some test polynomials:
+
+test1 :: UPoly Z X_
+test1 = toUPoly [1,2,3]
+
+test2 :: UPoly Z X_
+test2 = toUPoly [2,4,6,8,10]
+
+test3 :: UPoly Q X_
+test3 = toUPoly [inv 2, inv 3, inv 4]
+
+
+-- | Compute the content of a polynomial, i.e. the gcd of the coefficients.
+cont :: (GCDDomain a, Eq a) => UPoly a x -> a
+cont (UP xs) = case filter (/= zero) xs of
+  []  -> error "cont: Can't compute the content of the zero polynomial"
+  xs' -> ggcd xs'
+
+-- *Algebra.UPoly> cont test1
+-- 1
+-- *Algebra.UPoly> cont test2
+-- 2
+
+
+-- | If all coefficients are relatively prime then the polynomial is primitive.
+isPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> Bool
+isPrimitive = isUnit . cont
+
+-- *Algebra.UPoly> isPrimitive test1
+-- True
+-- *Algebra.UPoly> isPrimitive test2
+-- False
+
+
+-- | Lemma 4.2: Given a polynomial p in K[x] where K=Quot(A) we can find c in K
+-- and q primitive in A[x] such that p = cq.
+toPrimitive :: (GCDDomain a, Eq a)
+            => UPoly (FieldOfFractions a) x
+            -> (FieldOfFractions a, UPoly a x)
+toPrimitive p@(UP xs) = (c,q)
+  where
+  c0' = toFieldOfFractions $ productRing $ map denominator xs
+  c0  = inv c0'
+  g0  = map (fromFieldOfFractions . (c0' <*>)) xs
+  cg0 = toFieldOfFractions $ cont $ UP g0
+  c   = c0 <*> cg0
+  q   = toUPoly $ map (\x -> fromFieldOfFractions (toFieldOfFractions x </> cg0)) g0
+
+-- *Algebra.UPoly> toPrimitive test3
+-- (1/12,6+4x+3x^2)
+
+-- Specification of toPrimitive.
+propToPrimitive :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly (FieldOfFractions a) x -> Property
+propToPrimitive p =
+  not (isZeroPoly p) ==>
+    p == toUPoly (map (\x -> c <*> toFieldOfFractions x) q) && isPrimitive (UP q)
+  where (c,UP q) = toPrimitive p
+
+-- *Algebra.UPoly> quickCheck (propToPrimitive :: UPoly Q X_ -> Property)
+-- +++ OK, passed 100 tests.
+
+
+-- | Gauss lemma says that if p and q are polynomials over a GCD domain then
+-- cont(pq) = cont(p) * cont(q).
+gaussLemma :: (ExplicitUnits a, GCDDomain a, Eq a) => UPoly a x -> UPoly a x -> Property
+gaussLemma p q =
+  not (isZeroPoly p) && not (isZeroPoly q) ==>
+    cont (p <*> q) ~= cont p <*> cont q
+
+-- *Algebra.UPoly> quickCheck (gaussLemma :: UPoly Z X_ -> UPoly Z X_ -> Property)
+-- +++ OK, passed 100 tests.
+
+liftUPoly :: (GCDDomain a, Eq a) => UPoly a x -> UPoly (FieldOfFractions a) x
+liftUPoly (UP xs) = toUPoly $ map toFieldOfFractions xs
+
+-- | Proof that if A is a GCD domain then A[x] also is a GCD domain. This also
+-- computes witnesses that the computed GCD divides the given polynomials.
+gcdUPolyWitness :: (GCDDomain a, Eq a)
+                => UPoly a x
+                -> UPoly a x
+                -> (UPoly a x, UPoly a x, UPoly a x)
+gcdUPolyWitness p q = (constUPoly d <*> h, constUPoly x <*> a, constUPoly y <*> b)
+  where
+  (h',a',b') = gcd' (liftUPoly p) (liftUPoly q)
+
+  (_,h) = toPrimitive h'
+  (_,a) = toPrimitive a'
+  (_,b) = toPrimitive b'
+
+  (d,x,y) = gcd' (cont p) (cont q)
+
+test4, test5 :: UPoly Z X_
+test4 = toUPoly [6,7,1]
+test5 = toUPoly [-6,-5,1]
+
+-- *Algebra.UPoly> gcdUPolyWitness test4 test5
+-- (1+x,6+x,-6+x)
+-- *Algebra.UPoly> gcdUPolyWitness test4 test4
+-- (6+7x+x^2,1,1)
+-- *Algebra.UPoly> gcdUPolyWitness test1 test2
+-- (1,1+2x+3x^2,2+4x+6x^2+8x^3+10x^4)
+
+-- This does not work:
+-- instance (GCDDomain a, Eq a) => GCDDomain (UPoly a x) where
+--   gcd' = gcdUPolyWitness
+
+
+-------------------------------------------------------------------------------
 -- Square free decomposition.
 -- Teo Mora; Solving Polynomial Equations Systems I. pg 69-70
 -- Works only for Char 0
 --TODO: Add check for char
 --square free associate of f
--- | Square free decomposition of a polynomial. 
+-- | Square free decomposition of a polynomial.
 sqfr :: (Num k, Field k) => UPoly k x -> UPoly k x
 sqfr f = f `quotient` euclidAlg f f'
   where f' = deriv f
@@ -156,24 +288,24 @@
 -- | Distinct power factorization, aka square free decomposition
 sqfrDec :: (Num k, Field k) => UPoly k x -> [UPoly k x]
 sqfrDec f = help p q
-  where 
+  where
   p = euclidAlg f (deriv f)
-  q = f `quotient` p 
-  
-  help p q | d q < 1    = []
+  q = f `quotient` p
+
+  help p q | norm q < 1    = []
            | otherwise  = t : help (p `quotient` s) s
-    where 
+    where
     s = euclidAlg p q
     t = q `quotient` s
 
 -- | Pseudo-division of polynomials.
--- 
+--
 -- Given s(x) and p(x) compute c, q(x) and r(x) such that:
---   
+--
 --   cs(x) = p(x)q(x)+r(x), deg r < deg p.
-pseudoDivide :: (CommutativeRing a, Eq a) 
+pseudoDivide :: (CommutativeRing a, Eq a)
              => UPoly a x -> UPoly a x -> (a, UPoly a x, UPoly a x)
-pseudoDivide s p 
+pseudoDivide s p
   | m < n     = (one,zero,s)
   | otherwise = pD (a' <*> s' <-> b' <*> xmn <*> p') 1 (b' <*> xmn) s2
   where
@@ -189,7 +321,7 @@
   p'  = p <-> monomial a n
   s2  = a' <*> s' <-> b' <*> xmn <-> p'
 
-  pD s k out1 out2 
+  pD s k out1 out2
     | deg s < n = (a <^> k,out1,out2)
     | otherwise = pD s3 (k+1) (b2xm2n <+> a' <*> out1) s3
     where
@@ -200,6 +332,6 @@
     s3 = (a' <*> s) <-> (b2xm2n <*> p)
 
 
-propPD :: Qx -> Qx -> Property
+propPD :: (CommutativeRing a, Eq a) => UPoly a x -> UPoly a x -> Property
 propPD s p = deg s > 1 && deg p > 1 ==> constUPoly c <*> s == p <*> q <+> r
-  where (c,q,r) = pseudoDivide s p 
+  where (c,q,r) = pseudoDivide s p
diff --git a/src/Algebra/Z.hs b/src/Algebra/Z.hs
--- a/src/Algebra/Z.hs
+++ b/src/Algebra/Z.hs
@@ -1,15 +1,17 @@
 {-# LANGUAGE TypeSynonymInstances #-}
 module Algebra.Z 
   ( Z
-  , Ring(..)
+  , IntegralDomain(..)
   ) where
 
 import Test.QuickCheck
 
 import Algebra.Structures.IntegralDomain
 import Algebra.Structures.EuclideanDomain
+import Algebra.Structures.ExplicitUnits
 import Algebra.Structures.StronglyDiscrete
 import Algebra.Structures.BezoutDomain
+import Algebra.Structures.GCDDomain
 import Algebra.Structures.PruferDomain
 import Algebra.Structures.Coherent
 import Algebra.Ideal
@@ -34,17 +36,28 @@
 propIntegralDomainZ :: Z -> Z -> Z -> Property
 propIntegralDomainZ = propIntegralDomain
 
+instance ExplicitUnits Z where
+  unit 1    = Just 1
+  unit (-1) = Just (-1)
+  unit _    = Nothing
+
 instance EuclideanDomain Z where
-  d = abs
+  norm = abs
   quotientRemainder = quotRem
 
 propEuclideanDomainZ :: Z -> Z -> Z -> Property
 propEuclideanDomainZ = propEuclideanDomain
 
 -- Euclidean domain -> Bezout domain
-propBezoutDomainZ :: Ideal Z -> Z -> Z -> Z -> Property
+propBezoutDomainZ :: Z -> Z -> Property
 propBezoutDomainZ = propBezoutDomain
 
+propToPrincipalZ :: Ideal Z -> Bool
+propToPrincipalZ = propToPrincipal
+
+propIsSameIdealZ :: Ideal Z -> Bool
+propIsSameIdealZ = propIsSameIdeal
+
 -- Bezout domain -> Strongly discrete
 propStronglyDiscreteZ :: Z -> Ideal Z -> Bool
 propStronglyDiscreteZ = propStronglyDiscrete
@@ -65,6 +78,11 @@
 -- Not working perfectly...
 propSolveGeneralZ :: Matrix Z -> Vector Z -> Property
 propSolveGeneralZ = propSolveGeneral
+
+
+-- GCD Domain
+propGCDDomainZ :: Z -> Z -> Z -> Property
+propGCDDomainZ = propGCDDomain
 
 -- PLM
 propPLMZ :: Ideal Z -> Bool
diff --git a/src/Algebra/Zn.hs b/src/Algebra/Zn.hs
--- a/src/Algebra/Zn.hs
+++ b/src/Algebra/Zn.hs
@@ -2,7 +2,7 @@
              FunctionalDependencies, FlexibleContexts,    UndecidableInstances,
              FlexibleInstances #-}
 -- | Integers modulo n parametrised by the n. This also has type-level primality
--- testing used for instantiating integral domain and field type classes. The 
+-- testing used for instantiating integral domain and field type classes. The
 -- primality testing is very slow, but it seem to be working fine for relatively
 -- small numbers.
 module Algebra.Zn (Zn(..), Z3) where
@@ -28,12 +28,12 @@
   signum (Zn x) = Zn $ signum x
   negate (Zn x) = Zn $ negate x `mod` toNum (undefined :: n)
   fromInteger x = Zn $ fromInteger $ x `mod` toNum (undefined :: n)
-  
+
 instance Nat n => Arbitrary (Zn n) where
   arbitrary = liftM Zn (choose (0,toNum (undefined :: n) - 1))
- 
+
 instance Nat n => Ring (Zn n) where
-  (<+>) = (+) 
+  (<+>) = (+)
   zero  = Zn 0
   one   = Zn 1
   neg   = negate
@@ -61,13 +61,13 @@
 
 type Z17 = Zn D17
 
-test2 :: Z17
-test2 = inv 13
+-- test2 :: Z17
+-- test2 = inv 13
 
 -- Test that all elements of Z17 get correct inverses
-test3 :: Prelude.Bool
-test3 = all (==1) [ inv x * x | x <- xs ]
-  where xs :: [Z17] = map fromInteger [1..16]
+-- test3 :: Prelude.Bool
+-- test3 = all (==1) [ inv x * x | x <- xs ]
+--   where xs :: [Z17] = map fromInteger [1..16]
 
 -----------------------------------------------------------------------
 -- Lots of crazy type-level stuff:
@@ -78,7 +78,7 @@
 class Sqrt' x y r sqrt | x y r -> sqrt
 instance Sub y D2 y' => Sqrt' x y LT y'
 instance Pred y y'   => Sqrt' x y EQ y'
-instance (ExpBase y D2 square, Succ y y', Trich x square r, 
+instance (ExpBase y D2 square, Succ y y', Trich x square r,
           Sqrt' x y' r sqrt) => Sqrt' x y GT sqrt
 
 sqrt :: Sqrt x sqrt => x -> sqrt
@@ -89,7 +89,7 @@
 
 class Data.TypeLevel.Bool b => Prime' x y r b | x y r -> b
 instance Prime' x D1 EQ True
-instance (Pred y z, Trich z D1 r1, Mod x y rest, IsZero rest b1, 
+instance (Pred y z, Trich z D1 r1, Mod x y rest, IsZero rest b1,
           Not b1 b', Prime' x z r1 b2, And b' b2 b3) => Prime' x y GT b3
 
 prime :: Prime x b => x -> b