diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,5 @@
+#!/usr/bin/env runhaskell
+
+> import Distribution.Simple
+> main = defaultMain
+
diff --git a/splines.cabal b/splines.cabal
new file mode 100644
--- /dev/null
+++ b/splines.cabal
@@ -0,0 +1,36 @@
+name:                   splines
+version:                0.1
+stability:              provisional
+
+cabal-version:          >= 1.6
+build-type:             Simple
+
+author:                 James Cook <mokus@deepbondi.net>
+maintainer:             James Cook <mokus@deepbondi.net>
+license:                PublicDomain
+
+category:               Graphics, Numerical, Math
+synopsis:               B-Splines, other splines, and NURBS.
+description:            This is a fairly simple implementation of a 
+                        general-purpose spline library, just to get the code
+                        out there.  Its interface is still mildly unstable and
+                        may change (hopefully not drastically) as new needs or
+                        better style ideas come up.  Patches, suggestions
+                        and/or feature requests are welcome.
+
+source-repository head
+    type: darcs
+    location: http://code.haskell.org/~mokus/splines/
+
+Library
+  hs-source-dirs:       src
+  exposed-modules:      Math.Spline
+                        Math.Spline.BezierCurve
+                        Math.Spline.BSpline
+                        Math.Spline.Class
+                        Math.Spline.ISpline
+                        Math.Spline.Knots
+                        Math.Spline.MSpline
+                        Math.NURBS
+  other-modules:        Math.Spline.BSpline.Internal
+  build-depends:        base >= 3 && < 5, containers, vector-space
diff --git a/src/Math/NURBS.hs b/src/Math/NURBS.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/NURBS.hs
@@ -0,0 +1,79 @@
+{-# LANGUAGE StandaloneDeriving, FlexibleContexts, UndecidableInstances, TypeFamilies #-}
+module Math.NURBS
+    ( NURBS
+    , nurbs, toNURBS
+    , evalNURBS, nurbsDomain
+    , nurbsDegree, nurbsKnotVector, nurbsControlPoints
+    , splitNURBS
+    ) where
+
+import Data.VectorSpace
+import Math.Spline.Class (Spline, toBSpline)
+import Math.Spline.BSpline.Internal
+import Math.Spline.BSpline
+import Math.Spline.Knots
+
+newtype NURBS v = NURBS (BSpline (Scalar v, v))
+
+deriving instance (Eq   v, Eq   (Scalar v), Eq   (Scalar (Scalar v))) => Eq   (NURBS v)
+deriving instance (Ord  v, Ord  (Scalar v), Ord  (Scalar (Scalar v))) => Ord  (NURBS v)
+instance (Show v, Show (Scalar v), Show (Scalar (Scalar v))) => Show (NURBS v) where
+    showsPrec p (NURBS spline) = showParen (p>11)
+        ( showString "nurbs "
+        . showsPrec 11 spline
+        )
+
+toNURBS :: (Spline s v, Scalar v ~ Scalar (Scalar v)) => s v -> NURBS v
+toNURBS = NURBS . mapControlPoints (\p -> (1,p)) . toBSpline
+
+nurbs :: (VectorSpace v, Scalar v ~ w,
+          VectorSpace w, Scalar w ~ w)
+       => Knots (Scalar v) -> [(w, v)] -> NURBS v
+nurbs kts cps = NURBS (bSpline kts cps)
+
+-- |Constructs the homogeneous-coordinates B-spline that corresponds to this
+-- NURBS curve
+nurbsAsSpline (NURBS spline) = spline 
+    { controlPoints = map homogenize (controlPoints spline) }
+    where
+        homogenize (w,v) = (w, w *^ v)
+
+-- |Constructs the NURBS curve corresponding to a homogeneous-coordinates B-spline
+splineAsNURBS spline = NURBS spline 
+    { controlPoints = map unHomogenize (controlPoints spline) }
+    where
+        unHomogenize (w,v) = (w, recip w *^ v)
+
+
+evalNURBS
+  :: (VectorSpace v, Scalar v ~ w,
+      VectorSpace w, Scalar w ~ w,
+      Fractional w, Ord w) =>
+     NURBS v -> w -> v
+evalNURBS nurbs = project . evalBSpline (nurbsAsSpline nurbs)
+    where
+        project (w,v) = recip w *^ v
+
+
+-- |Returns the domain of a NURBS - that is, the range of parameter values
+-- over which a spline with this degree and knot vector has a full basis set.
+nurbsDomain :: Scalar v ~ Scalar (Scalar v) => 
+    NURBS v -> Maybe (Scalar v, Scalar v)
+nurbsDomain (NURBS spline) = knotDomain (knotVector spline) (degree spline)
+
+nurbsDegree :: NURBS v -> Int
+nurbsDegree (NURBS spline) = degree spline
+
+nurbsKnotVector :: Scalar v ~ Scalar (Scalar v) => NURBS v -> Knots (Scalar v)
+nurbsKnotVector (NURBS spline) = knotVector spline
+
+nurbsControlPoints :: NURBS v -> [(Scalar v, v)]
+nurbsControlPoints (NURBS spline) = controlPoints spline
+
+splitNURBS :: (VectorSpace v, Scalar v ~ w,
+               VectorSpace w, Scalar w ~ w,
+               Ord w, Fractional w)
+    => NURBS v -> Scalar v -> Maybe (NURBS v, NURBS v)
+splitNURBS nurbs t = do
+    (s0, s1) <- splitBSpline (nurbsAsSpline nurbs) t
+    return (splineAsNURBS s0, splineAsNURBS s1)
diff --git a/src/Math/Spline.hs b/src/Math/Spline.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline.hs
@@ -0,0 +1,17 @@
+module Math.Spline
+    ( Spline(..)
+    
+    , Knots, mkKnots, knots
+    
+    , BezierCurve, bezierCurve
+    , BSpline, bSpline
+    , MSpline, mSpline, toMSpline
+    , ISpline, iSpline, toISpline
+    ) where
+
+import Math.Spline.Class
+import Math.Spline.Knots
+import Math.Spline.BezierCurve
+import Math.Spline.BSpline
+import Math.Spline.MSpline
+import Math.Spline.ISpline
diff --git a/src/Math/Spline/BSpline.hs b/src/Math/Spline/BSpline.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/BSpline.hs
@@ -0,0 +1,91 @@
+{-# LANGUAGE MultiParamTypeClasses, StandaloneDeriving, FlexibleContexts, UndecidableInstances, TypeFamilies, ParallelListComp #-}
+module Math.Spline.BSpline
+    ( BSpline
+    , bSpline
+    , evalBSpline
+    , insertKnot
+    , splitBSpline
+    , differentiateBSpline, integrateBSpline
+    ) where
+
+import Math.Spline.Knots
+import Math.Spline.BSpline.Internal
+
+import Data.Maybe (fromMaybe)
+import Data.VectorSpace
+
+-- |@bSpline kts cps@ creates a B-spline with the given knot vector and control 
+-- points.  The degree is automatically inferred as the difference between the 
+-- number of spans in the knot vector (@numKnots kts - 1@) and the number of 
+-- control points (@length cps@).
+bSpline :: Knots (Scalar a) -> [a] -> BSpline a
+bSpline   _  [] = error "bSpline: no control points"
+bSpline kts cps = fromMaybe (error "bSpline: too few knots") (maybeSpline kts cps)
+
+maybeSpline :: Knots (Scalar a) -> [a] -> Maybe (BSpline a)
+maybeSpline kts cps 
+    | n > m     = Nothing
+    | otherwise = Just (Spline (m - n) kts cps)
+    where
+        n = length cps
+        m = numKnots kts - 1
+
+deriving instance (Eq   (Scalar v), Eq   v) => Eq   (BSpline v)
+deriving instance (Ord  (Scalar v), Ord  v) => Ord  (BSpline v)
+instance (Show (Scalar v), Show v) => Show (BSpline v) where
+    showsPrec p (Spline _ kts cps) = showParen (p>10) 
+        ( showString "bSpline "
+        . showsPrec 11 kts
+        . showChar ' '
+        . showsPrec 11 cps
+        )
+
+differentiateBSpline
+  :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v
+differentiateBSpline spline
+    | numKnots ks  < 2  = spline
+    | numKnots ks == 2  = bSpline ks [zeroV]
+    | otherwise         = bSpline ks' ds'
+    where
+        ks' = mkKnots . init . tail $ ts
+        ds' = zipWith (*^) (tail cs) (zipWith (^-^) (tail ds) ds)
+        
+        ks = knotVector spline; ts = knots ks
+        ds = controlPoints spline
+        
+        p  = degree spline
+        cs = [fromIntegral p / (t1 - t0) | (t0,t1) <- spans p ts]
+
+integrateBSpline
+  :: (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => BSpline v -> BSpline v
+integrateBSpline spline = bSpline (mkKnots ts') (scanl (^+^) zeroV ds')
+    where
+        ds' = zipWith (*^) cs (controlPoints spline)
+        ts = knots (knotVector spline)
+        ts' = head ts : ts ++ [last ts]
+        p = degree spline + 1
+        cs = [(t1 - t0) / fromIntegral p | (t0,t1) <- spans p ts]
+
+spans n xs = zip xs (drop n xs)
+
+splitBSpline
+  :: (VectorSpace v, Ord (Scalar v), Fractional (Scalar v)) =>
+     BSpline v -> Scalar v -> Maybe (BSpline v, BSpline v)
+splitBSpline spline@(Spline p kv ds) t 
+    | inDomain  = Just (Spline p (mkKnots us0) ds0, Spline p (mkKnots us1) ds1)
+    | otherwise = Nothing
+    where
+        inDomain = case knotDomain kv p of
+            Nothing         -> False
+            Just (t0, t1)   -> t >= t0 || t <= t1
+        
+        us = knots kv
+        dss = deBoor spline t
+        
+        us0 = takeWhile (<t) us ++ replicate (p+1) t
+        ds0 = trimTo (drop (p+1) us0) (map head dss)
+        
+        us1 = replicate (p+1) t ++ dropWhile (<=t) us
+        ds1 = reverse (trimTo (drop (p+1) us1) (map last dss))
+
+        trimTo list  xs = zipWith const xs list
diff --git a/src/Math/Spline/BSpline/Internal.hs b/src/Math/Spline/BSpline/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/BSpline/Internal.hs
@@ -0,0 +1,70 @@
+{-# LANGUAGE FlexibleContexts #-}
+module Math.Spline.BSpline.Internal
+    (BSpline(..), mapControlPoints, evalBSpline, insertKnot, deBoor) where
+
+import Math.Spline.Knots
+
+import Data.List (zipWith4)
+import Data.Monoid
+import Data.VectorSpace
+
+data BSpline v = Spline
+    { degree        :: !Int
+    , knotVector    :: Knots (Scalar v)
+    , controlPoints :: [v]
+    }
+
+mapControlPoints f spline = spline
+    { controlPoints = map f (controlPoints spline)
+    , knotVector = knotVector spline
+    }
+
+evalBSpline spline = head . last . deBoor spline
+
+-- |Insert one knot into a 'BSpline'
+insertKnot
+  :: (VectorSpace a, Ord (Scalar a), Fractional (Scalar a)) =>
+     BSpline a -> Scalar a -> BSpline a
+insertKnot spline x = spline
+    { knotVector    = knotVector spline `mappend` knot x
+    , controlPoints = zipWith4 (interp x) us (drop p us) ds (tail ds)
+    }
+    where
+        us = knots (knotVector spline)
+        p  = degree spline
+        ds = extend (controlPoints spline)
+
+
+-- duplicate the endpoints of a list; for example,
+-- extend [1..5] -> [1,1,2,3,4,5,5]
+extend []       = []
+extend (x:xs)   = x : extend' x xs
+    where   extend' x []      = [x,x]
+            extend' x (x':xs) = x:   extend' x' xs
+
+deBoor spline x = go us (controlPoints spline)
+    where
+        us = knots (knotVector spline)
+        
+        -- Upper endpoints of the intervals are the same for
+        -- each row in the table (they just line up differently
+        -- with the lower endpoints):
+        uHi = drop (degree spline + 1) us
+        
+        -- On each pass, the lower endpoints of the 
+        -- interpolation intervals advance and the new 
+        -- coefficients are given by linear interpolation
+        -- on the current intervals:
+        go       _ [] = []
+        go (_:uLo) ds = ds : go uLo ds'
+            where
+                ds' = zipWith4 (interp x) uLo uHi
+                                          ds (tail ds)
+
+interp x x0 x1 y0 y1
+    |  x <  x0  = y0
+    |  x >= x1  = y1
+    | otherwise = lerp y0 y1 a
+    where
+        a = (x - x0) / (x1 - x0)
+
diff --git a/src/Math/Spline/BezierCurve.hs b/src/Math/Spline/BezierCurve.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/BezierCurve.hs
@@ -0,0 +1,57 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}
+module Math.Spline.BezierCurve
+    ( BezierCurve, bezierCurve, splitBezierCurve
+    , evalSpline
+    ) where
+
+import Math.Spline.BSpline
+import Math.Spline.Class
+import Math.Spline.Knots
+
+import Control.Applicative
+import Data.VectorSpace
+
+-- |A BezierCurve curve on @0 <= x <= 1@.
+data BezierCurve v = BezierCurve !Int [v] deriving (Eq, Ord)
+
+-- |Construct a Bezier curve from a list of control points.  The degree
+-- of the curve is one less than the number of control points.
+bezierCurve :: [v] -> BezierCurve v
+bezierCurve cs
+    | null cs   = error "bezierCurve: no control points given"
+    | otherwise = BezierCurve (length cs - 1) cs
+
+instance Show v => Show (BezierCurve v) where
+    showsPrec p (BezierCurve _ cs) = showParen (p>10)
+        ( showString "bezierCurve "
+        . showsPrec 11 cs
+        )
+
+instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BezierCurve v where
+    splineDomain (BezierCurve _  _) = Just (0,1)
+    evalSpline   (BezierCurve _ cs) = head . last . deCasteljau cs
+    splineDegree (BezierCurve p  _) = p
+    knotVector   (BezierCurve p  _) = fromList [(0, p+1), (1, p+1)]
+    toBSpline = bSpline <$> knotVector <*> controlPoints
+
+instance Spline BezierCurve v => ControlPoints BezierCurve v where
+    controlPoints (BezierCurve _ cs) = cs
+
+deCasteljau :: VectorSpace v => [v] -> Scalar v -> [[v]]
+deCasteljau [] t = []
+deCasteljau cs t = cs : deCasteljau (zipWith interp cs (tail cs)) t
+    where
+        interp x0 x1 = lerp x0 x1 t
+
+-- |Split and rescale a Bezier curve.  Given a 'BezierCurve' @b@ and a point 
+-- @t@, @splitBezierCurve b t@ creates 2 curves @(b1, b2)@ such that (up to 
+-- reasonable numerical accuracy expectations):
+-- 
+-- > evalSpline b1  x    == evalSpline b (x * t)
+-- > evalSpline b2 (x-t) == evalSpline b (x * (1-t))
+-- 
+splitBezierCurve :: VectorSpace v => BezierCurve v -> Scalar v -> (BezierCurve v, BezierCurve v)
+splitBezierCurve (BezierCurve n cs) t = 
+    ( BezierCurve n (map head css)
+    , BezierCurve n (reverse (map last css))
+    ) where css = deCasteljau cs t
diff --git a/src/Math/Spline/Class.hs b/src/Math/Spline/Class.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/Class.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, UndecidableInstances #-}
+module Math.Spline.Class where
+
+import Control.Applicative
+import Math.Spline.Knots
+import qualified Math.Spline.BSpline.Internal as BSpline
+
+import Data.VectorSpace
+
+-- |A spline is a piecewise polynomial vector-valued function.  The necessary
+-- and sufficient instance definition is 'toBSpline'.
+class (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline s v where
+    -- |Returns the domain of a spline.  In the case of B-splines, this is
+    -- the domain on which a spline with this degree and knot vector has a 
+    -- full basis set.  In other cases, it should be no larger than 
+    -- @splineDomain . toBSpline@, but may be smaller.  Within this domain,
+    -- 'evalSpline' should agree with @'evalSpline' . 'toBSpline'@ (not 
+    -- necessarily exactly, but up to reasonable expectations of numerical 
+    -- accuracy).
+    splineDomain :: s v -> Maybe (Scalar v, Scalar v)
+    splineDomain = knotDomain <$> knotVector <*> splineDegree
+    
+    evalSpline :: s v -> Scalar v -> v
+    evalSpline = evalSpline . toBSpline
+    
+    splineDegree :: s v -> Int
+    splineDegree = splineDegree . toBSpline
+    
+    knotVector :: s v -> Knots (Scalar v)
+    knotVector = knotVector . toBSpline
+    
+    toBSpline :: s v -> BSpline.BSpline v
+
+class Spline s v => ControlPoints s v where
+    controlPoints :: s v -> [v]
+
+instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline BSpline.BSpline v where
+    evalSpline spline = head . last . BSpline.deBoor spline
+    splineDegree = BSpline.degree
+    knotVector = BSpline.knotVector
+    toBSpline = id
+
+instance Spline BSpline.BSpline v => ControlPoints BSpline.BSpline v where
+    controlPoints = BSpline.controlPoints
diff --git a/src/Math/Spline/ISpline.hs b/src/Math/Spline/ISpline.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/ISpline.hs
@@ -0,0 +1,74 @@
+{-# LANGUAGE
+        MultiParamTypeClasses,
+        FlexibleInstances, FlexibleContexts, UndecidableInstances,
+        ParallelListComp,
+        StandaloneDeriving
+  #-}
+module Math.Spline.ISpline
+    ( ISpline, iSpline, toISpline
+    , evalSpline
+    ) where
+
+import Math.Spline.BSpline
+import Math.Spline.Class
+import Math.Spline.Knots
+
+import Data.VectorSpace
+
+-- |The I-Spline basis functions are the integrals of the M-splines, or
+-- alternatively the integrals of the B-splines normalized to the range
+-- [0,1].  Every I-spline basis function increases monotonically from 0 to 1,
+-- thus it is useful as a basis for monotone functions.  An I-Spline curve
+-- is monotone if and only if every non-zero control point has the same sign.
+data ISpline v = ISpline
+    { iSplineDegree        :: !Int
+    , iSplineKnotVector    :: Knots (Scalar v)
+    , iSplineControlPoints :: [v]
+    }
+
+deriving instance (Eq   (Scalar v), Eq   v) => Eq   (ISpline v)
+deriving instance (Ord  (Scalar v), Ord  v) => Ord  (ISpline v)
+instance (Show (Scalar v), Show v) => Show (ISpline v) where
+    showsPrec p (ISpline _ kts cps) = showParen (p>10) 
+        ( showString "iSpline "
+        . showsPrec 11 kts
+        . showChar ' '
+        . showsPrec 11 cps
+        )
+
+
+-- |@iSpline kts cps@ creates an I-spline with the given knot vector and control 
+-- points.  The degree is automatically inferred as the difference between the 
+-- number of spans in the knot vector (@numKnots kts - 1@) and the number of 
+-- control points (@length cps@).
+iSpline :: Knots (Scalar a) -> [a] -> ISpline a
+iSpline kts cps 
+    | n > m     = error "iSpline: too few knots"
+    | otherwise = ISpline (m - n) kts cps
+    where
+        n = length cps
+        m = numKnots kts - 1
+
+instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline ISpline v where
+    splineDegree = (1 +) . iSplineDegree
+    knotVector spline = mkKnots (head ts : ts ++ [last ts])
+        where ts = knots (iSplineKnotVector spline)
+    toBSpline spline = bSpline (knotVector spline) (scanl (^+^) zeroV cs)
+        where cs = iSplineControlPoints spline
+
+instance Spline ISpline v => ControlPoints ISpline v where
+    controlPoints      (ISpline _  _ cs) = cs
+
+toISpline :: (Spline s v, Eq v) => s v -> ISpline v
+toISpline = fromBSpline . toBSpline
+
+fromBSpline spline
+    | head ds == zeroV 
+    && numKnots ks >= 2 = iSpline (mkKnots (init (tail ts))) (tail ds')
+    | otherwise         = iSpline (mkKnots (init       ts )) ds'
+    where
+        ks = knotVector spline
+        ts = knots ks
+        ds = controlPoints spline
+        
+        ds' = zipWith (^-^) ds (zeroV:ds)
diff --git a/src/Math/Spline/Knots.hs b/src/Math/Spline/Knots.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/Knots.hs
@@ -0,0 +1,108 @@
+module Math.Spline.Knots
+    ( Knots
+    , knot, multipleKnot
+    , mkKnots, fromList
+    
+    , knots, numKnots
+    , toList, distinctKnots, numDistinctKnots
+    
+    , knotMultiplicity, setKnotMultiplicity
+    
+    , knotDomain
+    ) where
+
+import Prelude hiding (sum)
+import Data.Foldable (Foldable(foldMap), sum)
+import qualified Data.Map as M
+import Data.Monoid (Monoid(..))
+import Data.Maybe (fromMaybe)
+
+-- |Knot vectors - multisets of points in a 1-dimensional space.
+data Knots a = Knots !Int (M.Map a Int) deriving (Eq, Ord)
+
+instance Show a => Show (Knots a) where
+    showsPrec p ks@(Knots 1 _) = showParen (p > 10)
+        ( showString "knot "
+        . showsPrec 11 (head $ knots ks)
+        )
+    showsPrec p ks = showParen (p > 10)
+        ( showString "mkKnots "
+        . showsPrec 11 (knots ks)
+        )
+
+instance (Ord a) => Monoid (Knots a) where
+    mempty = Knots 0 M.empty
+    mappend (Knots n1 v1) (Knots n2 v2) =
+        Knots (n1 + n2) (M.filter (/=0) (M.unionWith (+) v1 v2))
+
+instance Foldable Knots where
+    foldMap f = foldMap f . knots
+
+
+-- |Create a knot vector consisting of one knot.
+knot :: Ord a => a -> Knots a
+knot x = multipleKnot x 1
+
+-- |Create a knot vector consisting of one knot with the specified multiplicity.
+multipleKnot :: Ord a => a -> Int -> Knots a
+multipleKnot k n 
+    | n <= 0    = Knots 0 (M.empty)
+    | otherwise = Knots n (M.singleton k n)
+
+-- |Create a knot vector consisting of all the knots in a list.
+mkKnots :: (Ord a) => [a] -> Knots a
+mkKnots ks = fromList (map (\k -> (k,1)) ks)
+
+-- |Create a knot vector consisting of all the knots and corresponding 
+-- multiplicities in a list.
+fromList :: (Ord k) => [(k, Int)] -> Knots k
+fromList ks = Knots (sum kMap) kMap
+    where kMap = M.fromListWith (+) (filter ((>0).snd) ks)
+
+-- |Returns a list of all distinct knots in ascending order along with
+-- their multiplicities.
+toList :: Knots k -> [(k, Int)]
+toList (Knots _ ks) = M.toList ks
+
+-- |Returns the number of knots (not necessarily distinct) in a knot vector.
+numKnots :: Knots t -> Int
+numKnots (Knots n _) = n
+
+-- |Returns the number of distinct knots in a knot vector.
+numDistinctKnots :: Knots t -> Int
+numDistinctKnots (Knots _ ks) = M.size ks
+
+-- |Returns a list of all knots (not necessarily distinct) of a knot vector in ascending order
+knots :: Knots t -> [t]
+knots (Knots _ ks) = concat [replicate n k | (k,n) <- M.toAscList ks]
+
+-- |Returns a list of all distinct knots of a knot vector in ascending order
+distinctKnots :: Knots t -> [t]
+distinctKnots (Knots _ ks) = M.keys ks
+
+-- |Looks up the multiplicity of a knot (which is 0 if the point is not a knot)
+knotMultiplicity :: (Ord k) => k -> Knots k -> Int
+knotMultiplicity k (Knots _ ks) = fromMaybe 0 (M.lookup k ks)
+
+-- |Returns a new knot vector with the given knot set to the specified 
+-- multiplicity and all other knots unchanged.
+setKnotMultiplicity :: Ord k => k -> Int -> Knots k -> Knots k
+setKnotMultiplicity k n (Knots m ks)
+    | n <= 0    = Knots (m     - n') (M.delete k ks)
+    | otherwise = Knots (m + n - n') (M.insert k n ks)
+    where
+        n' = knotMultiplicity k (Knots m ks)
+
+-- |@knotDomain kts p@ return the domain of a B-spline or NURBS with knot
+-- vector @kts@ and degree @p@.  This is the subrange spanned by all
+-- except the first and last @p@ knots.  Outside this domain, the spline
+-- does not have a complete basis set.  De Boor's algorithm assumes that
+-- the basis functions sum to 1, which is only true on this range, and so
+-- this is also precisely the domain on which de Boor's algorithm is valid.
+knotDomain :: Knots a -> Int -> Maybe (a,a)
+knotDomain ks@(Knots n _) p 
+    | n > 2*p   = Just (head (drop p kts), head (drop p (reverse kts)))
+    | otherwise = Nothing
+    where
+        kts = knots ks
+
diff --git a/src/Math/Spline/MSpline.hs b/src/Math/Spline/MSpline.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/Spline/MSpline.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE
+        MultiParamTypeClasses,
+        FlexibleInstances, FlexibleContexts, UndecidableInstances,
+        ParallelListComp,
+        StandaloneDeriving
+  #-}
+module Math.Spline.MSpline
+    ( MSpline, mSpline, toMSpline
+    , evalSpline
+    ) where
+
+import Math.Spline.BSpline
+import Math.Spline.Class
+import Math.Spline.Knots
+
+import Data.VectorSpace
+
+-- |M-Splines are B-splines normalized so that the integral of each basis 
+-- function over the spline domain is 1.
+data MSpline v = MSpline
+    { mSplineDegree        :: !Int
+    , mSplineKnotVector    :: Knots (Scalar v)
+    , mSplineControlPoints :: [v]
+    }
+
+deriving instance (Eq   (Scalar v), Eq   v) => Eq   (MSpline v)
+deriving instance (Ord  (Scalar v), Ord  v) => Ord  (MSpline v)
+instance (Show (Scalar v), Show v) => Show (MSpline v) where
+    showsPrec p (MSpline _ kts cps) = showParen (p>10) 
+        ( showString "mSpline "
+        . showsPrec 11 kts
+        . showChar ' '
+        . showsPrec 11 cps
+        )
+
+
+-- |@mSpline kts cps@ creates a M-spline with the given knot vector and control 
+-- points.  The degree is automatically inferred as the difference between the 
+-- number of spans in the knot vector (@numKnots kts - 1@) and the number of 
+-- control points (@length cps@).
+mSpline :: Knots (Scalar a) -> [a] -> MSpline a
+mSpline kts cps
+    | n > m     = error "mSpline: too few knots"
+    | otherwise = MSpline (m - n) kts cps
+    where
+        n = length cps
+        m = numKnots kts - 1
+
+spans n xs = zip xs (drop n xs)
+
+instance (VectorSpace v, Fractional (Scalar v), Ord (Scalar v)) => Spline MSpline v where
+    splineDegree = mSplineDegree
+    knotVector   = mSplineKnotVector
+    toBSpline (MSpline p ks cs) = bSpline ks cs'
+        where
+            n = p + 1; n' = fromIntegral n
+            cs' = [ (n' / (t1 - t0)) *^ c 
+                  | c <- cs
+                  | (t0, t1) <- spans n (knots ks)
+                  ]
+
+instance Spline MSpline v => ControlPoints MSpline v where
+    controlPoints = mSplineControlPoints
+
+toMSpline :: Spline s v => s v -> MSpline v
+toMSpline = fromBSpline . toBSpline
+
+fromBSpline spline = mSpline ks cs
+    where
+        n = splineDegree spline + 1; n' = fromIntegral n
+        ks = knotVector spline
+        cs =  [ ((t1 - t0) / n') *^ c
+              | c <- controlPoints spline
+              | (t0, t1) <- spans n (knots ks)
+              ]