diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+The MIT License (MIT)
+
+Copyright (c) vi
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/benchmark/Main.lhs b/benchmark/Main.lhs
new file mode 100644
--- /dev/null
+++ b/benchmark/Main.lhs
@@ -0,0 +1,39 @@
+We aim to benchmark each implementation of Lloyd's algorithm:
+
+> import Algorithms.Lloyd.Sequential (Point(..), Cluster(..))
+> import qualified Algorithms.Lloyd.Sequential as Sequential (kmeans)
+> import qualified Algorithms.Lloyd.Strategies as Strategies (kmeans)
+> import Data.Random.Normal (mkNormals)
+> import Data.Vector (fromList)
+> 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
+
+Three of which form the centroids of our initial clusters:
+ 
+> clusters :: [Cluster]
+> clusters = map (uncurry Cluster) $ zip [0..] (take 3 points)
+
+To correctly benchmark the result of a pure function, we need to be able to
+evaluate it to normal form:
+
+> instance NFData Cluster where
+>   rnf (Cluster i c) = rnf i `seq` rnf c
+>
+> instance NFData Point where
+>   rnf (Point v) = rnf v
+ 
+Together the subject of our benchmarks:
+ 
+> main :: IO ()
+> main = defaultMain
+>   [ bench "Sequential" $ nf (Sequential.kmeans points) clusters
+>   , bench "Strategies" $ nf (Strategies.kmeans 64 points) clusters
+>   ] 
diff --git a/kmeans-par.cabal b/kmeans-par.cabal
new file mode 100644
--- /dev/null
+++ b/kmeans-par.cabal
@@ -0,0 +1,26 @@
+name:                kmeans-par
+version:             1.0.2
+synopsis:            Sequential and parallel implementations of Lloyd's algorithm.
+license:             MIT
+license-file:        LICENSE
+author:              vi
+maintainer:          me@vikramverma.com
+category:            Algorithm
+build-type:          Simple
+cabal-version:       >=1.10
+
+library
+  exposed-modules:     Algorithms.Lloyd.Sequential,
+                       Algorithms.Lloyd.Strategies
+  build-depends:       base == 2.*, vector, semigroups, parallel, split
+  hs-source-dirs:      src
+  ghc-options:         -threaded
+  default-language:    Haskell2010
+
+benchmark kmeans-benchmark
+  type:             exitcode-stdio-1.0
+  main-is:          Main.lhs
+  hs-source-dirs:   benchmark
+  build-depends:    base == 2.*, random, criterion, normaldistribution, kmeans-par, deepseq, vector
+  ghc-options:      -threaded
+  default-language: Haskell2010
diff --git a/src/Algorithms/Lloyd/Sequential.lhs b/src/Algorithms/Lloyd/Sequential.lhs
new file mode 100644
--- /dev/null
+++ b/src/Algorithms/Lloyd/Sequential.lhs
@@ -0,0 +1,148 @@
+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
+> 
+> import Prelude hiding (zipWith, map, foldr1, replicate)
+> import Control.Monad (guard, forM_)
+> import Data.Function (on)
+> import Data.List (minimumBy)
+> import Data.Semigroup (Semigroup(..))
+> import Data.Vector (Vector(..), toList, create, zipWith, map, foldr1, empty, replicate)
+> import qualified Data.Vector.Mutable as MV (replicate, read, write)
+ 
+A data point is represented by an n-dimensional real vector:
+
+> data Point = Point (Vector Double) deriving (Eq, Show)
+
+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 vector) = Point $ map f vector
+ 
+The `Metric` typeclass, as defined here, is intended to contain types that
+admit metrics (i.e. are metric spaces.) Instances can be defined in terms of 
+`distance` or the infix `<->`:
+
+> class Metric a where
+>   distance :: a -> a -> Double
+>   distance = (<->)
+>   (<->) :: a -> a -> Double
+>   (<->) = distance
+
+We define distance between data-points as the sum of squares of differences of
+corresponding coordinates. The actual (Euclidean) distance is given by the
+square root of this value, but as we're only interested in comparing distances,
+we omit that function for efficiency:
+
+> instance Metric Point where
+>   Point v0 <-> Point v1 = foldr1 (+) . map (**2) $ zipWith (-) v0 v1
+ 
+`Point`s may be combined using simple vector addition:
+
+> instance Semigroup Point where
+>   Point v0 <> Point v1 = Point $ zipWith (+) v0 v1
+
+Clusters are represented by an identifier and a centroid:
+
+> data Cluster = Cluster
+>   { identifier :: Int
+>   , centroid   :: Point
+>   } deriving (Show)
+
+Clusters are treated as equivalent if they share a centroid:
+
+> instance Eq Cluster where
+>   (Cluster _ c) == (Cluster _ d) = c == d
+
+`PointSum` is a point sum; this is an intermediate type used to contain the sum
+of points in a set when constructing new clusters. It consists of a count of
+included points, as well as the vector sum of the constituent points:
+
+> data PointSum = PointSum Int Point deriving (Show)
+
+Two `PointSum`s can be combined by summing corresponding values:
+
+> instance Semigroup PointSum where
+>   PointSum c0 p0 <> PointSum c1 p1 = PointSum (c0+c1) (p0<>p1)
+
+Without dependent types, we can't define a valid Monoid instance for PointSum,
+as the identity varies with the dimensions of the point vector. But we can
+construct one at runtime:
+
+> emptyPointSum :: Int -> PointSum
+> emptyPointSum length = PointSum 0 . Point $ replicate length 0
+
+A point sum is constructed by incrementally adding points likeso:
+
+> addPoint :: PointSum -> Point -> PointSum
+> addPoint sum point = sum <> pure point
+>   where pure = PointSum 1 
+
+A point sum can be turned into a cluster by computing its centroid, the average
+of all points associated with the cluster:
+
+> toCluster :: Int -> PointSum -> Cluster
+> toCluster cid (PointSum count point) = Cluster
+>   { identifier = cid
+>   , centroid   = pmap (//count) point
+>   }
+>
+> (//) :: Double -> Int -> Double
+> x // y = x / fromIntegral y
+         
+After each iteration, points are associated with the cluster having the nearest
+centroid:
+
+> closestCluster :: [Cluster] -> Point -> Cluster
+> closestCluster clusters point = fst . minimumBy (compare `on` snd) $ do
+>   cluster <- clusters
+>   return (cluster, point <-> centroid cluster)
+>
+> assign :: [Cluster] -> [Point] -> Vector PointSum
+> assign clusters points = let nc = length clusters in create $ do
+>   vector <- MV.replicate nc $ emptyPointSum nc
+>   points `forM_` \point -> do
+>     let cluster  = closestCluster clusters point 
+>         position = identifier cluster
+>     sum <- MV.read vector position
+>     MV.write vector position $! addPoint sum point
+>   return vector
+>
+> makeNewClusters :: Vector PointSum -> [Cluster]
+> makeNewClusters vector = do
+>   (pointSum@(PointSum count _), index) <- zip (toList vector) [0..]
+>   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 :: [Cluster] -> [Point] -> [Cluster]
+> step = makeNewClusters ..: assign
+>
+> (..:) :: (Functor f, Functor g) => (a -> b) -> f (g a) -> f (g b)
+> (..:) = fmap . fmap
+> 
+> infixr 8 ..:
+
+The algorithm consists of iteratively finding the centroid of each existing
+cluster, then reallocating points according to which centroid is closest, until
+convergence. As the algorithm isn't guaranteed to converge, we cut execution if
+convergence hasn't been observed after eighty iterations:
+
+> kmeans :: [Point] -> [Cluster] -> [Cluster]
+> kmeans = kmeans' 0
+>
+> kmeans' :: Int -> [Point] -> [Cluster] -> [Cluster]
+> kmeans' iterations points clusters 
+>   | iterations >= expectDivergent = clusters
+>   | clusters' == clusters         = clusters 
+>   | otherwise                     = kmeans' (succ iterations) points clusters'
+>   where clusters' = step clusters points
+>
+> expectDivergent :: Int
+> expectDivergent = 80
+
+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
new file mode 100644
--- /dev/null
+++ b/src/Algorithms/Lloyd/Strategies.lhs
@@ -0,0 +1,47 @@
+A parallel implementation of Lloyd's algorithm for k-means clustering,
+adapted from Marlow's _Parallel and Concurrent Programming in Haskell_.
+Here we use Evaluation Strategies to parallelise the assignment of
+points to clusters:
+
+> module Algorithms.Lloyd.Strategies where
+>
+> import Prelude hiding (zipWith)
+> import Control.Parallel.Strategies (Strategy(..), parTraversable, using, rseq)
+> import Data.List.Split (chunksOf)
+> import Data.Semigroup (Semigroup(..))
+> import Data.Vector (Vector(..), zipWith)
+> import Algorithms.Lloyd.Sequential (Cluster(..), Point(..), PointSum(..), makeNewClusters, assign, (..:), expectDivergent)
+
+We can combine two vectors of some same type $t$ provided we know how to
+combine two $t$s:
+
+> instance Semigroup a => Semigroup (Vector a) where
+>   (<>) = zipWith (<>)
+
+Step is modified to, given a partitioned list of points, perform
+classification in parallel:
+
+> step :: [Cluster] -> [[Point]] -> [Cluster]
+> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: map . assign 
+>
+> with :: Strategy a -> a -> a
+> with = flip using
+
+This version of k-means takes an additional arguments -- the number of
+partitions the set of points'll be divided into. This needn't equal the
+number of processors: if there are more spark than cores, the runtime
+can be trusted to schedule unallocated sparks so soon as a core becomes
+available. That said: if there are too many work items, the overhead of
+recombination may exceed the speed-up provided by parallellism; if there
+are too few items, and those items vary in cost, some of our cores may
+be unused for part of the computation.
+
+> kmeans :: Int -> [Point] -> [Cluster] -> [Cluster]
+> kmeans = kmeans' 0 ..: chunksOf
+>
+> kmeans' :: Int -> [[Point]] -> [Cluster] -> [Cluster]
+> kmeans' iterations points clusters 
+>   | iterations >= expectDivergent = clusters
+>   | clusters' == clusters         = clusters 
+>   | otherwise                     = kmeans' (succ iterations) points clusters'
+>   where clusters' = step clusters points