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.1.5
+Version:             0.1.6
 
 Synopsis:            A library of constructive algebra.
 Description:         
@@ -65,6 +65,8 @@
                        Algebra.Matrix,
                        Algebra.PLM,
                        Algebra.UPoly,
+                       Algebra.EllipticCurve,
+                       Algebra.ZSqrt5,
                        Algebra.Zn,
                        Algebra.Z,
                        Algebra.Q
diff --git a/src/Algebra/EllipticCurve.hs b/src/Algebra/EllipticCurve.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/EllipticCurve.hs
@@ -0,0 +1,112 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+-- | The elliptic curve y^2 = 1 - x^4 in Q[x,y].
+module Algebra.EllipticCurve (EllipticCurve(..)) where
+
+import Test.QuickCheck
+
+import Algebra.Structures.Field hiding ((<*), (*>))
+import Algebra.Structures.EuclideanDomain (quotient, genEuclidAlg)
+import Algebra.Structures.BezoutDomain (toPrincipal)
+import Algebra.Structures.PruferDomain
+import Algebra.Structures.Coherent
+import Algebra.FieldOfRationalFunctions
+import Algebra.Ideal
+import Algebra.UPoly
+
+-- | The elliptic curve y^2=1-x^4 over Q[x,y].
+newtype EllipticCurve = C (Qx,Qx)
+  deriving (Eq,Arbitrary)
+
+instance Show EllipticCurve where
+  show (C (a,b)) | a == zero && b == zero = "0"
+                 | a == zero              = show b ++ "*y"
+                 | b == zero              = show a
+                 | otherwise = case show b of 
+                   ['-','1'] -> show a ++ "-y"
+                   ('-':xs)  -> show a ++ "-" ++ xs ++ "*y"
+                   xs        -> show a ++ "+" ++ xs ++ "*y" 
+
+-- Arithmetical properties
+instance Ring EllipticCurve where
+  (C (a,b)) <+> (C (c,d)) = C (a + c, b + d)
+  (C (a,b)) <*> (C (c,d)) = C (a*c + b*d*(1-x^4), a*d + b*c)
+  neg (C (a,b))           = C (neg a, neg b)
+  zero                    = C (zero,zero)
+  one                     = C (one,zero)
+  
+instance CommutativeRing EllipticCurve where
+
+instance IntegralDomain EllipticCurve where
+
+propIntDomEC :: EllipticCurve -> EllipticCurve-> EllipticCurve -> Property
+propIntDomEC = propIntegralDomain
+
+--------------------------------------------------------------------------------
+-- Useful auxiliary functions:
+
+(*>), (+>) :: Qx -> EllipticCurve -> EllipticCurve
+r *> (C (a,b)) = C (r*a,r*b)
+r +> (C (a,b)) = C (r+a,b)
+
+(<*), (<+) :: EllipticCurve -> Qx -> EllipticCurve
+(C (a,b)) <* r = C (a*r,b*r)
+(C (a,b)) <+ r = C (a+r,b)
+
+infixl 7 *>, <*
+infixl 6 +>, <+
+
+--------------------------------------------------------------------------------
+
+instance PruferDomain EllipticCurve where
+  calcUVW (C (a,b)) (C (c,d)) = (u,v,w)
+    where
+    p = toQX (a * c - b * d * (1 - x^4)) </> toQX (c^2 - d^2 * (1 - x^4))
+    q = toQX (b * c - a * d) </> toQX (c^2 - d^2 * (1 - x^4))
+
+    s :: (QX,QX) 
+    s = (p,q)
+
+    -- a0's^2 + a1's + a2' = 0
+    a0' = (c^2 - d^2 * (1 - x^4))^2
+    a1' = -2 * (a * c - b * d * (1 - x^4)) * (c^2 - d^2 * (1 - x^4))
+    a2' = (a * c - b * d * (1 - x^4))^2 - ((b * c - a * d)^2 * (1-x^4))
+
+    -- Make <a0,a1,a2> = 1
+    g  = genEuclidAlg [a0',a1',a2']
+    a0 = a0' `quotient` g
+    a1 = a1' `quotient` g
+    a2 = a2' `quotient` g
+
+    -- n0 * a0 + n1 * a1 + n2 * a2 = 1
+    (Id [g'],[n0,n1,n2],_) = toPrincipal (Id [a0,a1,a2])
+
+    a0s    = case s of
+      (p,q) -> C (toQx (a0' <*> p), toQx (a0' <*> q))
+      where a0' = toQX a0
+    a0sa1  = a0s <+ a1
+    a0sa1s = C (neg a2,zero)
+
+    alpha = a0s 
+    beta  = a0sa1s
+
+    m0 = n0
+    m1 = -n1
+    m2 = n1
+    m3 = -n2
+
+    u = m0 * a0 +> m2 *> a0sa1
+    v = m0 *> alpha <+> m2 *> beta
+    w = m1 * a0 +> m3 *> a0sa1
+
+instance Coherent EllipticCurve where
+  solve = solvePD
+
+-- Properties:
+propPruferDomEC :: EllipticCurve -> EllipticCurve -> EllipticCurve -> Property
+propPruferDomEC x@(C (a,b)) y@(C (c,d)) z@(C (e,f)) = 
+  a /= zero && b /= zero && c /= zero && d /= zero && e /= zero && f /= zero 
+  ==> propPruferDomain x y z
+
+propIntersectionPEC :: Ideal EllipticCurve -> Ideal EllipticCurve -> Property
+propIntersectionPEC i@(Id is) j@(Id js) = 
+  length is <= 5 && length js <= 5 ==> isSameIdeal intersectionPDWitness i j
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
@@ -2,7 +2,7 @@
 -- principal ideal domains. This means that all finitely generated ideals are
 -- principal.
 --
-{-# LANGUAGE FlexibleInstances, UndecidableInstances #-}
+{-# LANGUAGE FlexibleInstances, UndecidableInstances, OverlappingInstances #-}
 module Algebra.Structures.BezoutDomain
   ( BezoutDomain(..)
   , propBezoutDomain
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
@@ -9,10 +9,11 @@
   , calcUVWT, propCalcUVWT, fromUVWTtoUVW
   , computePLM_PD
   , invertIdeal
+  , intersectionPD, intersectionPDWitness, solvePD
   ) where
 
 import Test.QuickCheck
-import Data.List (nub)
+import Data.List (nub, (\\))
 
 import Algebra.Structures.IntegralDomain
 import Algebra.Structures.Coherent
@@ -120,13 +121,65 @@
       a_njs = [ head (a !! j) | j <- [0..length a - 1]]
   in Id a_njs
 
--- XXX: This is buggy at the moment... Witnesses is not correctly computed!
 -- | Compute the intersection of I and J by:
 --       
 --       (I ∩ J)(I + J) = IJ  => (I ∩ J)(I + J)(I + J)' = IJ(I + J)'
 --
-intersectionP :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
-intersectionP (Id is) (Id js) = case foldr combine ([],[],[]) int of
+intersectionPDWitness :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
+intersectionPDWitness (Id is) (Id js) = (int,wis,wjs)
+  where
+  lj  = length js 
+  li  = length is
+
+  ij  = Id (is ++ js)
+
+  plm = computePLM_PD ij
+
+  as  = take li $ unMVec $ transpose plm
+  as' = drop li $ unMVec $ transpose plm
+
+  int = Id $ concat [ map (j <*>) a | j <- js , a <- as ]
+
+  wis = concat [ [ addZ i li a | a <- as ] | as <- as', i <- [0..li-1] ]
+  wjs = [ addZ i lj a | i <- [0..lj-1], a <- concat as ]
+
+  addZ n l x = replicate n zero ++ x : replicate (l-n-1) zero
+
+{-
+intersectionPD :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
+intersectionPD (Id xs) (Id ys) 
+  | xs' == [] || ys' == [] = zeroIdealWitnesses xs ys 
+  | otherwise              = (Id k, [handleZero xs as], [handleZero ys bs])
+  where
+  -- Compute <x1,...,xn> and <y1,...,ym>
+  xs' = filter (/= zero) xs
+  ys' = filter (/= zero) ys
+
+  -- Compute <z_1...z_k>, k = n+m
+  ij  = Id xs' `addId` Id ys'
+
+  -- Compute <a_11,...,a_k1>
+  inv = fromId $ invertIdeal ij
+
+
+  k  = undefined
+  as = undefined
+  bs = undefined
+
+-- Handle the zeroes specially. If the first element in xs is a zero
+-- then the witness should be zero otherwise use the computed witness. 
+handleZero :: (Ring a, Eq a) => [a] -> [a] -> [a]
+handleZero xs [] 
+  | all (==zero) xs = xs
+  | otherwise       = error "intersectionPD: This should be impossible"
+handleZero (x:xs) (a:as) 
+  | x == zero = zero : handleZero xs (a:as)
+  | otherwise = a    : handleZero xs as
+handleZero [] _  = error "intersectionPD: This should be impossible"
+-}
+{-
+intersectionPDWitness :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> (Ideal a,[[a]],[[a]])
+intersectionPDWitness (Id is) (Id js) = case foldr combine ([],[],[]) int of
   ([],_,_)   -> zeroIdealWitnesses is js
   (xs,ys,zs) -> (Id xs,ys,zs)
   where
@@ -151,20 +204,29 @@
 
   combine (x,y,z) (xs,ys,zs) = (x:xs,y:ys,z:zs)
 
--- intersectionPD :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> Ideal a
-intersectionPD i@(Id is) j@(Id js) = i `mulId` k
-  where
-  plm = unMVec $ computePLM_PD (i `addId` j)
-
-  n = length is - 1 
-  m = n + length js
+  as = filter (/= zero) $ concat 
+     $ drop (length (is \\ js)) 
+     $ unMVec 
+     $ transpose 
+     $ computePLM_PD 
+     $ Id is `addId` Id js
 
-  k = Id [ plm !! i !! j | j <- [n+1..m], i <- [0..m]]
---  k = [ "a" ++ show i  ++ show j | j <- [n+1..m], i <- [0..m]]
+  -- concatMap (replicate (length ys)) as of
+  
+  asdf ys = [ addZ i li a | (i,_) <- zip [0..] is, a <- as ]
 
+  -- case  of 
+  -- case [ concatMap (addZ i (length is)) a | (i,a) <- zip [0..] (replicate (length ys) as) ] of
+--    [[]] -> [ is ]
+--    x    -> x -- map (filter (/= zero)) x
+-}
+intersectionPD :: (PruferDomain a, Eq a) => Ideal a -> Ideal a -> Ideal a
+intersectionPD i j = fst3 (intersectionPDWitness i j)
+  where fst3 (x,_,_) = x
 
+-- | Coherence of Prufer domains.
 solvePD :: (PruferDomain a, Eq a) => Vector a -> Matrix a
-solvePD x = solveWithIntersection x intersectionP
+solvePD x = solveWithIntersection x intersectionPDWitness
 
 -- instance (PruferDomain a, Eq a) => Coherent a where
 --   solve x = solveWithIntersection x intersectIdeals
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
@@ -2,7 +2,7 @@
 module Algebra.Structures.Ring 
   ( Ring(..)
   , propRing
-  , (<->), (<^>), (*>)
+  , (<->), (<^>), (*>), (<*)
   , sumRing, productRing
   ) where
 
@@ -12,6 +12,7 @@
 infixl 8 <^>
 infixl 7 <*>
 infixl 7 *>
+infixl 7 <*
 infixl 6 <+>
 infixl 6 <->
 
@@ -110,9 +111,13 @@
              else x <*> x <^> (y-1)
 
 -- | Multiply from left with an integer; n *> x means x + x + ... + x, n times.
-(*>) :: Ring a => Int -> a -> a
-n *> x = sumRing $ replicate n x
+(*>) :: Ring a => Integer -> a -> a
+0 *> _ = zero
+n *> x | n > 0     = x <+> x <* (n-1)
+       | otherwise = neg (abs n *> x) -- error "<*: Negative input"
 
 -- Multiply from right with an integer.
--- (<*) :: Ring a => a -> Integer -> a
--- x <* n = sumRing $ replicate n x
+(<*) :: Ring a => a -> Integer -> a
+_ <* 0 = zero
+x <* n | n > 0     = x <+> x <* (n-1)
+       | otherwise = neg (x <* abs n) -- error "<*: Negative input"
diff --git a/src/Algebra/Z.hs b/src/Algebra/Z.hs
--- a/src/Algebra/Z.hs
+++ b/src/Algebra/Z.hs
@@ -8,8 +8,9 @@
 
 import Algebra.Structures.IntegralDomain
 import Algebra.Structures.EuclideanDomain
-import Algebra.Structures.BezoutDomain
 import Algebra.Structures.StronglyDiscrete
+import Algebra.Structures.BezoutDomain
+import Algebra.Structures.PruferDomain
 import Algebra.Structures.Coherent
 import Algebra.Ideal
 import Algebra.Matrix
@@ -68,3 +69,10 @@
 -- PLM
 propPLMZ :: Ideal Z -> Bool
 propPLMZ id = propPLM id (computePLM_B id)
+
+-- Prufer domain
+propPruferDomainZ :: Z -> Z -> Z -> Property
+propPruferDomainZ = propPruferDomain
+
+propIsSameIdealPruferDomain :: Ideal Z -> Ideal Z -> Bool
+propIsSameIdealPruferDomain = isSameIdeal intersectionPDWitness
diff --git a/src/Algebra/ZSqrt5.hs b/src/Algebra/ZSqrt5.hs
new file mode 100644
--- /dev/null
+++ b/src/Algebra/ZSqrt5.hs
@@ -0,0 +1,104 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+-- | Proof that Z[sqrt(-5)] is a Prufer domain. This implies that it is 
+-- possible to solve systems of equations over Z[sqrt(-5)].
+module Algebra.ZSqrt5 (ZSqrt5(..)) where
+
+import Test.QuickCheck
+
+import Algebra.Structures.IntegralDomain
+import Algebra.Structures.EuclideanDomain (quotient, genEuclidAlg)
+import Algebra.Structures.BezoutDomain (toPrincipal)
+import Algebra.Structures.PruferDomain
+import Algebra.Structures.Coherent
+import Algebra.Ideal
+import Algebra.Z
+import Algebra.Q
+
+-- | Z[sqrt(-5)] is a pair such that (a,b) = a + b*sqrt(-5)
+newtype ZSqrt5 = ZSqrt5 (Z,Z)
+  deriving (Eq,Ord,Arbitrary)
+
+instance Show ZSqrt5 where
+  show (ZSqrt5 (a,b)) = show a ++ " + " ++ show b ++ " * sqrt(-5)"
+
+-- Arithmetical properties
+instance Ring ZSqrt5 where
+  (ZSqrt5 (a,b)) <+> (ZSqrt5 (c,d)) = ZSqrt5 (a + c, b + d)
+  (ZSqrt5 (a,b)) <*> (ZSqrt5 (c,d)) = ZSqrt5 (a*c - 5*b*d, a*d + b*c)
+  neg (ZSqrt5 (a,b))                = ZSqrt5 (neg a, neg b)
+  zero                              = ZSqrt5 (0,0)
+  one                               = ZSqrt5 (1,0)
+  
+instance CommutativeRing ZSqrt5 where
+
+instance IntegralDomain ZSqrt5 where
+
+propIntDomZSqrt5 :: ZSqrt5 -> ZSqrt5 -> ZSqrt5 -> Property
+propIntDomZSqrt5 = propIntegralDomain
+
+--------------------------------------------------------------------------------
+-- Useful auxiliary functions:
+
+
+(+>) :: Z -> ZSqrt5 -> ZSqrt5
+r +> (ZSqrt5 (a,b)) = ZSqrt5 (r+a,b)
+
+(<+) :: ZSqrt5 -> Z -> ZSqrt5
+(ZSqrt5 (a,b)) <+ r = ZSqrt5 (a+r,b)
+
+infixl 6 +>, <+
+
+--------------------------------------------------------------------------------
+
+instance PruferDomain ZSqrt5 where
+  -- Assume /= 0
+  calcUVW (ZSqrt5 (a,b)) (ZSqrt5 (c,d)) = (u,v,w) 
+    where
+    -- Let s = (a+b*sqrt(-5))/(c+d*sqrt(-5))
+    -- Compute p and q such that: s = p + q*sqrt(-5)
+    p = toQ (a*c+5*b*d) / n
+    q = toQ (b*c-a*d)   / n
+    n = toQ (c^2+5*d^2)
+
+    s :: (Q,Q)
+    s = (p,q)
+
+    -- Rewrite: 
+    -- s = p + q*sqrt(-5)  <==>  s^2 - 2*sp + p^2 + 5*q^2 = 0
+    --                     <==> a0*s^2 + a1*s + a2 = 0
+    --
+    -- Some computations give:
+    a0' = toZ n
+    a1' = -2 * toZ (p*n)
+    a2' = a^2 + 5*b^2 
+    
+    -- Normalize:
+    g  = genEuclidAlg [a0',a1',a2']
+    a0 = a0' `quotient` g
+    a1 = a1' `quotient` g
+    a2 = a2' `quotient` g
+
+    -- Compute m0, m1, m2, m3 such that:
+    -- m0a0 + m1a0s + m2(a0s+a1) + m3(a0s+a1)s = 1
+    a0s    = ZSqrt5 (toZ (a0 *> p), toZ (a0 *> q))
+    a0sa1  = a0s <+ a1
+    a0sa1s = ZSqrt5 (-a2,0)
+
+    (Id [1],[n0,n1,n2],_) = toPrincipal (Id [a0,a1,a2]) 
+
+    m0 = n0
+    m1 = -n1
+    m2 = n1
+    m3 = -n2
+
+    -- Finally we get u, v and w:
+    u = m0 * a0 +> m2 *> a0sa1
+    v = m0 *> a0s <+> m2 *> a0sa1s
+    w = m1 * a0 +> m3 *> a0sa1
+
+propPruferDomZSqrt5 :: ZSqrt5 -> ZSqrt5 -> ZSqrt5 -> Property
+propPruferDomZSqrt5 x@(ZSqrt5 (a,b)) y@(ZSqrt5 (c,d)) z@(ZSqrt5 (e,f)) = 
+  a /= 0 && b /= 0 && c /= 0 && d /= 0 && e /= 0 && f /= 0 ==> propPruferDomain x y z
+
+instance Coherent ZSqrt5 where
+  solve = solvePD