diff --git a/benchmark/Main.lhs b/benchmark/Main.lhs
--- a/benchmark/Main.lhs
+++ b/benchmark/Main.lhs
@@ -1,26 +1,26 @@
 We aim to benchmark each implementation of Lloyd's algorithm:
 
+> import Prelude hiding (take, zipWith)
 > import Algorithms.Lloyd.Sequential (Point(..), Cluster(..))
 > import qualified Algorithms.Lloyd.Sequential as Sequential (kmeans)
 > import qualified Algorithms.Lloyd.Strategies as Strategies (kmeans)
 > import Data.Metric (Metric(..), Euclidean(..))
 > import Data.Random.Normal (mkNormals)
-> import Data.Vector (fromList)
+> import Data.Vector (Vector, generate, fromList, take, zipWith)
 > import Control.Monad (forM)
 > import Control.DeepSeq (NFData(..))
 > import Criterion.Main (defaultMain, bench, nf)
 
 We draw 2e3 normally distributed 2D points:
 
-> points :: [Point]
-> points = concat $ [0..2000-1] `forM` \n -> do
->   return . Point . fromList $ [normals!!n, normals!!(n*2)]
->   where normals = take 4000 $ mkNormals 0x29a
+> points :: Vector Point
+> points = generate 2000 $ \n -> Point $ fromList [normals !! n, normals !! (n*2)]
+>   where normals = mkNormals 0x29a
 
 Three of which form the centroids of our initial clusters:
  
-> clusters :: [Cluster]
-> clusters = map (uncurry Cluster) $ zip [0..] (take 3 points)
+> clusters :: Vector Cluster
+> clusters = zipWith Cluster (fromList [0..2]) (take 3 points)
 
 To correctly benchmark the result of a pure function, we need to be able to
 evaluate it to normal form:
diff --git a/kmeans-par.cabal b/kmeans-par.cabal
--- a/kmeans-par.cabal
+++ b/kmeans-par.cabal
@@ -1,5 +1,5 @@
 name:                kmeans-par
-version:             1.3.0
+version:             1.4.0
 synopsis:            Sequential and parallel implementations of Lloyd's algorithm.
 license:             MIT
 license-file:        LICENSE
@@ -12,9 +12,10 @@
 library
   exposed-modules:     Algorithms.Lloyd.Sequential,
                        Algorithms.Lloyd.Strategies
-  other-modules:       Data.Functor.Extras
+  other-modules:       Data.Functor.Extras,
+                       Data.Vector.Split
   other-extensions:    ViewPatterns 
-  build-depends:       base == 4.*, vector, semigroups, parallel, split, metric
+  build-depends:       base == 4.*, vector, semigroups, parallel, metric
   hs-source-dirs:      src
   ghc-options:         -threaded -O2 -with-rtsopts=-N -rtsopts
   default-language:    Haskell2010
diff --git a/src/Algorithms/Lloyd/Sequential.lhs b/src/Algorithms/Lloyd/Sequential.lhs
--- a/src/Algorithms/Lloyd/Sequential.lhs
+++ b/src/Algorithms/Lloyd/Sequential.lhs
@@ -3,19 +3,25 @@
 A sequential implementation of Lloyd's algorithm for k-means clustering,
 adapted from Marlow's _Parallel and Concurrent Programming in Haskell_:
 
-> module Algorithms.Lloyd.Sequential where
+> module Algorithms.Lloyd.Sequential (
+>   Point(..),
+>   Cluster(..), 
+>   kmeans,
+>   PointSum(..),
+>   makeNewClusters,
+>   assign,
+>   assignPS
+> )where
 > 
-> import Prelude hiding (zipWith, map, foldr, replicate)
+> import Prelude hiding (zipWith, map, foldr, replicate, length, zip, head)
 > import Control.Applicative ((<$>))
-> import Control.Monad (guard, forM_)
+> import Control.Monad (guard)
 > import Data.Foldable (Foldable(foldr))
 > import Data.Function (on)
 > import Data.Functor.Extras ((..:),(...:))
-> import Data.List (minimumBy)
 > import Data.Metric (Metric(..))
 > import Data.Semigroup (Semigroup(..))
-> import Data.Vector (Vector(..), toList, create, zipWith, map, empty, replicate, cons)
-> import qualified Data.Vector as V (length)
+> import Data.Vector (Vector(..), toList, fromList, create, zip, zipWith, map, empty, replicate, cons, minimumBy, length, head, forM_)
 > import qualified Data.Vector.Mutable as MV (replicate, read, write)
  
 A data point is represented by an n-dimensional real vector:
@@ -27,8 +33,8 @@
 Point isn't parametrised, so we can't define a functor (or any functor-family)
 instances. Instead, we define a specialised map, `pmap`:
 
-> pmap :: (Double -> Double) -> Point -> Point
-> pmap f = Point . map f . point
+> pmap :: (Vector Double -> Vector Double) -> Point -> Point
+> pmap f = Point . f . point
 
 To define distance between data-points, we admit arbitrary metric spaces
 over real vectors; **However:** it should be noted that convergence is
@@ -80,7 +86,7 @@
 > toCluster :: Int -> PointSum -> Cluster
 > toCluster cid (PointSum count point) = Cluster
 >   { identifier = cid
->   , centroid   = pmap (// count) point
+>   , centroid   = map (//count) `pmap` point
 >   }
 >
 > (//) :: Double -> Int -> Double
@@ -89,34 +95,35 @@
 After each iteration, points are associated with the cluster having the nearest
 centroid:
 
-> closestCluster :: Metric a => (Vector Double -> a) -> [Cluster] -> Point -> Cluster
+> closestCluster :: Metric a => (Vector Double -> a) -> Vector Cluster -> Point -> Cluster
 > closestCluster (useMetric -> d) clusters point = fst . minimumBy (compare `on` snd) $ do
 >   cluster <- clusters
 >   return (cluster, point `d` centroid cluster)
 >
-> assign :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> Vector (Vector Point)
+> assign :: Metric a => (Vector Double -> a) -> Vector Cluster -> Vector Point -> Vector (Vector Point)
 > assign metric clusters points = let nc = length clusters in create $ do
->   vector <- MV.replicate nc empty
+>   vector <- MV.replicate (nc+1) empty
 >   points `forM_` \point -> do
->     let position = identifier $ closestCluster metric clusters point
+>     let cluster  = closestCluster metric clusters point
+>         position = identifier cluster
 >     points' <- MV.read vector position
->     MV.write vector position $! point `cons` points'
+>     MV.write vector position $ point `cons` points'
 >   return vector
 >
-> assignPS :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> Vector PointSum
+> assignPS :: Metric a => (Vector Double -> a) -> Vector Cluster -> Vector Point -> Vector PointSum
 > assignPS metric clusters points = reduce <$> assign metric clusters points
->   where reduce  = foldr (<>) (emptyPointSum length') . map (PointSum 1)
->         length' = V.length . point $ head points
+>   where reduce  = foldr (<>) (emptyPointSum length') . fmap (PointSum 1)
+>         length' = length . point $ head points
 >
-> makeNewClusters :: Vector PointSum -> [Cluster]
+> makeNewClusters :: Vector PointSum -> Vector Cluster
 > makeNewClusters vector = do
->   (pointSum@(PointSum count _), index) <- zip (toList vector) [0..]
+>   (pointSum@(PointSum count _), index) <- zip vector $ fromList [0..length vector]
 >   guard $ count > 0 -- We don't want an empty PointSum: amongst other 
 >                     -- things, this'd lead to division by zero when 
 >                     -- computing the centroid.
 >   return $ toCluster index pointSum
 >
-> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> [Cluster]
+> step :: Metric a => (Vector Double -> a) -> Vector Cluster -> Vector Point -> Vector Cluster
 > step = makeNewClusters ...: assignPS
 >
 
@@ -125,19 +132,19 @@
 convergence. As the algorithm isn't guaranteed to converge, we cut execution if
 convergence hasn't been observed after eighty iterations:
 
-> computeClusters :: Metric a => Int -> (Vector Double -> a) -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters :: Metric a => Int -> (Vector Double -> a) -> Vector Point -> Vector Cluster -> Vector Cluster
 > computeClusters expectDivergent metric = computeClusters' expectDivergent metric 0 
 >
-> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> Vector Point -> Vector Cluster -> Vector Cluster
 > computeClusters' expectDivergent metric iterations points clusters 
 >   | iterations >= expectDivergent = clusters
 >   | clusters' == clusters         = clusters 
 >   | otherwise                     = computeClusters' expectDivergent metric (succ iterations) points clusters'
 >   where clusters' = step metric clusters points
 >
-> kmeans :: Metric a => Int -> (Vector Double -> a) -> [Point] -> [Cluster] -> Vector (Vector Point)
+> kmeans :: Metric a => Int -> (Vector Double -> a) -> Vector Point -> Vector Cluster -> Vector (Vector Point)
 > kmeans expectDivergent metric points initial = assign metric clusters points
->   where clusters = computeClusters 80 metric points initial
+>   where clusters = computeClusters expectDivergent metric points initial
 
 A note on initialisation: typically clusters are randomly assigned to data
 points prior to the first update step.
diff --git a/src/Algorithms/Lloyd/Strategies.lhs b/src/Algorithms/Lloyd/Strategies.lhs
--- a/src/Algorithms/Lloyd/Strategies.lhs
+++ b/src/Algorithms/Lloyd/Strategies.lhs
@@ -3,15 +3,20 @@
 Here we use Evaluation Strategies to parallelise the assignment of
 points to clusters:
 
-> module Algorithms.Lloyd.Strategies where
+> module Algorithms.Lloyd.Strategies (
+>   Point(..),
+>   Cluster(..), 
+>   kmeans
+> ) where
 >
-> import Prelude hiding (zipWith)
+> import Prelude hiding (zipWith, foldr1, map)
 > import Control.Parallel.Strategies (Strategy(..), parTraversable, using, rseq)
+> import Data.Foldable (Foldable(foldr1))
 > import Data.Functor.Extras ((..:))
-> import Data.List.Split (chunksOf)
 > import Data.Metric (Metric(..))
 > import Data.Semigroup (Semigroup(..))
-> import Data.Vector (Vector(..), zipWith)
+> import Data.Vector (Vector(..), zipWith, map)
+> import Data.Vector.Split (chunksOf)
 > import Algorithms.Lloyd.Sequential (Cluster(..), Point(..), PointSum(..), makeNewClusters, assignPS, assign)
 
 We can combine two vectors of some same type $t$ provided we know how to
@@ -23,8 +28,8 @@
 Step is modified to, given a partitioned list of points, perform
 classification in parallel:
 
-> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [[Point]] -> [Cluster]
-> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: map ..: assignPS
+> step :: Metric a => (Vector Double -> a) -> Vector Cluster -> Vector (Vector Point) -> Vector Cluster
+> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: fmap ..: assignPS
 >
 > with :: Strategy a -> a -> a
 > with = flip using
@@ -38,16 +43,16 @@
 are too few items, and those items vary in cost, some of our cores may
 be unused for part of the computation.
 
-> computeClusters :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters :: Metric a => Int -> (Vector Double -> a) -> Int -> Vector Point -> Vector Cluster -> Vector Cluster
 > computeClusters expectDivergent metric = computeClusters' expectDivergent metric 0  ..: chunksOf
 >
-> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> [[Point]] -> [Cluster] -> [Cluster]
+> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> Vector (Vector Point) -> Vector Cluster -> Vector Cluster
 > computeClusters' expectDivergent metric iterations points clusters 
 >   | iterations >= expectDivergent = clusters
 >   | clusters' == clusters         = clusters 
 >   | otherwise                     = computeClusters' expectDivergent metric (succ iterations) points clusters'
 >   where clusters' = step metric clusters points
 >
-> kmeans :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> Vector (Vector Point)
+> kmeans :: Metric a => Int -> (Vector Double -> a) -> Int -> Vector Point -> Vector Cluster -> Vector (Vector Point)
 > kmeans expectDivergent metric iterations points initial = assign metric clusters points
 >   where clusters = computeClusters 80 metric iterations points initial
diff --git a/src/Data/Vector/Split.hs b/src/Data/Vector/Split.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Vector/Split.hs
@@ -0,0 +1,9 @@
+module Data.Vector.Split (chunksOf) where
+
+import Prelude hiding (take, drop, null)
+import Data.Vector (Vector(..), cons, take, drop, empty, null)
+
+chunksOf :: Int -> Vector a -> Vector (Vector a)
+chunksOf n vector 
+  | null vector = empty
+  | otherwise   = take n vector `cons` chunksOf n (drop n vector)