diff --git a/ChangeLog b/ChangeLog
new file mode 100644
--- /dev/null
+++ b/ChangeLog
@@ -0,0 +1,9 @@
+0.1:
+
+* Hermite1 interpolation: different order of coefficients
+
+  Interleave node values and derivatives
+  in order to get a narrow banded Gramian matrix
+  from a sampled interpolation basis.
+
+  Rename from 'cubic' to 'hermite1'.
diff --git a/example/Fit.hs b/example/Fit.hs
--- a/example/Fit.hs
+++ b/example/Fit.hs
@@ -4,8 +4,13 @@
 import qualified Numeric.Interpolation.Piecewise as Piecewise
 import qualified Numeric.Interpolation.Type as Type
 
+import qualified Data.Packed.ST as PackST
 import qualified Data.Packed.Matrix as Matrix
 import qualified Data.Packed.Vector as Vector
+import Data.Packed.Matrix (Matrix)
+import Data.Packed.Vector (Vector)
+
+import qualified Numeric.LinearAlgebra.Banded as Banded
 import qualified Numeric.Container as Container
 import Numeric.Container ((<\>))
 
@@ -15,6 +20,9 @@
 
 import System.Random (randomRs, mkStdGen)
 import Control.Monad.HT (void)
+import Control.Monad (when, zipWithM_, forM_)
+
+import qualified Data.Foldable as Fold
 import Data.Monoid ((<>))
 
 
@@ -26,19 +34,93 @@
       (randomRs (0,7) (mkStdGen 23))
       (randomRs (-0.2,0.2) (mkStdGen 42))
 
+basisMatrixFull ::
+   Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double
+basisMatrixFull typ xs txs0 =
+   let txs = Vector.fromList txs0
+   in  Matrix.fromColumns $
+       map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $
+       Type.basisFunctions typ xs
+
+basisMatrixSparse ::
+   Type.T Double Double ny -> [Double] -> [Double] -> Matrix Double
+basisMatrixSparse typ xs txs = PackST.runSTMatrix $ do
+   mat <- PackST.newMatrix 0 (length txs) (length $ Type.basisFunctions typ xs)
+   zipWithM_
+      (\k -> mapM_ (uncurry (PackST.writeMatrix mat k))) [0..] $
+      map (Type.sampleBasisFunctions typ xs) txs
+   return mat
+
 fit ::
    Type.T Double Double ny ->
    [Double] -> [(Double, Double)] -> Nodes.T Double ny
 fit typ xs target =
-   let txs = Vector.fromList $ map fst target
-       tys = Vector.fromList $ map snd target
-       matrix =
-          Matrix.fromColumns $
-          map (flip Container.cmap txs . Piecewise.interpolateConstantExt typ) $
-          Type.basisFunctions typ xs
+   let (txs, tys) = unzip target
+       matrix = basisMatrixSparse typ xs txs
    in  Type.coefficientsToInterpolator typ xs $
-       Vector.toList $ matrix <\> tys
+       Vector.toList $ matrix <\> Vector.fromList tys
 
+matrixDiff ::
+   Type.T Double Double ny ->
+   [Double] -> [(Double, Double)] -> Double
+matrixDiff typ xs target =
+   let txs = map fst target
+   in  Container.maxElement $ Container.cmap abs $
+       Container.sub
+          (basisMatrixFull typ xs txs)
+          (basisMatrixSparse typ xs txs)
+
+
+mulSparseMatrixVector ::
+   Int -> [[(Int, Double)]] -> [Double] -> Vector Double
+mulSparseMatrixVector size samples tys = PackST.runSTVector $ do
+   vec <- PackST.newVector 0 size
+   forM_ (zip samples tys) $ \(row,ty) ->
+      forM_ row $ \(k,y) ->
+         PackST.modifyVector vec k (+y*ty)
+   return vec
+
+bandedGramian ::
+   Int -> Int -> [[(Int, Double)]] -> Banded.SymmetricMatrix Double
+bandedGramian size width samples =
+      Banded.SymmetricMatrix $ PackST.runSTMatrix $ do
+   mat <- PackST.newMatrix 0 size width
+   forM_ samples $ \row ->
+      forM_ row $ \(k,yk) ->
+      forM_ row $ \(j,yj) ->
+         when (k<=j) $ PackST.modifyMatrix mat k (j-k) (+yk*yj)
+   return mat
+
+fitBanded ::
+   Type.T Double Double ny ->
+   [Double] -> [(Double, Double)] -> Nodes.T Double ny
+fitBanded typ xs target =
+   let size = length $ Type.basisFunctions typ xs
+       (txs, tys) = unzip target
+       samples = map (Type.sampleBasisFunctions typ xs) txs
+       matrix =
+          Banded.choleskyDecompose $
+          bandedGramian size (Type.basisOverlap typ) samples
+   in  Type.coefficientsToInterpolator typ xs $ Vector.toList $
+       Banded.choleskySolve matrix $ mulSparseMatrixVector size samples tys
+
+bandedDiff ::
+   (ny -> ny -> Double) ->
+   Type.T Double Double ny ->
+   [Double] -> [(Double, Double)] -> Double
+bandedDiff absDiff typ xs target =
+   maximum $
+   zipWith absDiff
+      (Fold.toList $ fit typ xs target)
+      (Fold.toList $ fitBanded typ xs target)
+
+absDiffSingle :: Double -> Double -> Double
+absDiffSingle x y = abs (x-y)
+
+absDiffPair :: (Double,Double) -> (Double,Double) -> Double
+absDiffPair (x,dx) (y,dy) = max (abs (x-y)) (abs (dx-dy))
+
+
 plotBasisFunctions ::
    Type.T Double Double ny -> [Double] -> Plot2D.T Double Double
 plotBasisFunctions nodeType xs =
@@ -53,11 +135,11 @@
    let xs = [0, 1, 3, 4, 6, 7]
        exs = (-1) : xs ++ [8]
    void $ GP.plotDefault $ plotBasisFunctions Type.linear xs
-   void $ GP.plotDefault $ plotBasisFunctions Type.cubic xs
+   void $ GP.plotDefault $ plotBasisFunctions Type.hermite1 xs
    void $ GP.plotDefault $ plotBasisFunctions Type.cubicLinear exs
    void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs
    let linearNodes = fit Type.linear xs noisy
-       hermite1Nodes = fit Type.cubic xs noisy
+       hermite1Nodes = fit Type.hermite1 xs noisy
        cubicLinearNodes = fit Type.cubicLinear exs noisy
        cubicParabolaNodes = fit Type.cubicParabola exs noisy
    void $ GP.plotDefault $
@@ -65,7 +147,19 @@
       <>
       (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $
        Piecewise.interpolateConstantExt Type.linear linearNodes :
-       Piecewise.interpolateConstantExt Type.cubic hermite1Nodes :
+       Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes :
        Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :
        Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :
        [])
+
+   putStrLn "differences between matrices should be almost zero:"
+   putStrLn $ "linear: " ++ show (matrixDiff Type.linear xs noisy)
+   putStrLn $ "hermite1: " ++ show (matrixDiff Type.hermite1 xs noisy)
+   putStrLn $ "cubicLinear: " ++ show (matrixDiff Type.cubicLinear exs noisy)
+   putStrLn $ "cubicParabola: " ++ show (matrixDiff Type.cubicParabola exs noisy)
+
+   putStrLn "differences between samples should be almost zero:"
+   putStrLn $ "linear: " ++ show (bandedDiff absDiffSingle Type.linear xs noisy)
+   putStrLn $ "hermite1: " ++ show (bandedDiff absDiffPair Type.hermite1 xs noisy)
+   putStrLn $ "cubicLinear: " ++ show (bandedDiff absDiffPair Type.cubicLinear exs noisy)
+   putStrLn $ "cubicParabola: " ++ show (bandedDiff absDiffPair Type.cubicParabola exs noisy)
diff --git a/example/Plot.hs b/example/Plot.hs
--- a/example/Plot.hs
+++ b/example/Plot.hs
@@ -27,7 +27,7 @@
       Plot2D.functions Graph2D.lines
          (Plot2D.linearScale 1000 (-2,15))
          [Piecewise.interpolateConstantExt Type.linear linearNodes,
-          Piecewise.interpolateConstantExt Type.cubic hermite1Nodes,
-          Piecewise.interpolateConstantExt Type.cubic cubicLinearNodes,
-          Piecewise.interpolateConstantExt Type.cubic cubicParabolaNodes,
+          Piecewise.interpolateConstantExt Type.hermite1 hermite1Nodes,
+          Piecewise.interpolateConstantExt Type.hermite1 cubicLinearNodes,
+          Piecewise.interpolateConstantExt Type.hermite1 cubicParabolaNodes,
           sin]
diff --git a/interpolation.cabal b/interpolation.cabal
--- a/interpolation.cabal
+++ b/interpolation.cabal
@@ -1,5 +1,5 @@
 Name:             interpolation
-Version:          0.0
+Version:          0.1
 License:          BSD3
 License-File:     LICENSE
 Author:           Henning Thielemann
@@ -8,7 +8,7 @@
 Category:         Math
 Synopsis:         piecewise linear and cubic Hermite interpolation
 Description:
-  Represent real functions by linear or cubic segments.
+  Represent real functions by linear or cubic polynomial segments.
   The package provides both data structures
   for efficient lookup of interpolation intervals,
   and computation of basis functions.
@@ -30,9 +30,13 @@
     We use a distorted sinus as target.
   .
   The package needs only Haskell 98.
-Tested-With:      GHC==7.4.2, GHC==7.6.3, GHC==7.8.2
+  Most of the package dependencies are only needed for the examples
+  and are only installed if you enable to build them.
+Tested-With:      GHC==7.4.2, GHC==7.6.3, GHC==7.8.4, GHC==7.10.1
 Cabal-Version:    >=1.8
 Build-Type:       Simple
+Extra-Source-Files:
+  ChangeLog
 
 Flag buildExamples
   description: Build example executables
@@ -43,7 +47,7 @@
   default:     True
 
 Source-Repository this
-  Tag:         0.0
+  Tag:         0.1
   Type:        darcs
   Location:    http://code.haskell.org/~thielema/interpolation/
 
@@ -66,9 +70,11 @@
     Numeric.Interpolation.Basis
     Numeric.Interpolation.Basis.Compact
     Numeric.Interpolation.Basis.Full
+    Numeric.Interpolation.Sample
   Other-Modules:
     Numeric.Interpolation.Private.Piece
     Numeric.Interpolation.Private.Basis
+    Numeric.Interpolation.Private.List
 
 Executable interpolation-plot
   Main-Is:             Plot.hs
@@ -90,6 +96,7 @@
   If flag(buildExamples) && flag(hmatrix)
     Build-Depends:
       interpolation,
+      hmatrix-banded >=0.0 && <0.1,
       hmatrix >=0.15 && <0.16,
       random >=1.0 && <1.1,
       gnuplot >=0.5.2 && <0.6,
@@ -103,10 +110,14 @@
   Main-Is:             Test.hs
   Other-Modules:
     Test.Piece
+    Test.Sample
+    Test.Overlap
   Hs-Source-Dirs:      test, private
   GHC-Options:         -Wall
   Build-Depends:
     interpolation,
     QuickCheck >=2.4 && <2.8,
     utility-ht >=0.0.9 && <0.1,
+    array >=0.4 && <0.6,
+    containers >=0.4 && <0.6,
     base >=4.5 && <4.8
diff --git a/private/Numeric/Interpolation/Private/Basis.hs b/private/Numeric/Interpolation/Private/Basis.hs
--- a/private/Numeric/Interpolation/Private/Basis.hs
+++ b/private/Numeric/Interpolation/Private/Basis.hs
@@ -5,8 +5,16 @@
 import qualified Data.List.Match as Match
 
 
+_hermite1Split :: [a] -> [b] -> [(b, b)]
+_hermite1Split xs = uncurry zip . Match.splitAt xs
+
 hermite1Split :: [a] -> [b] -> [(b, b)]
-hermite1Split xs = uncurry zip . Match.splitAt xs
+hermite1Split _ = pairs
+
+pairs :: [a] -> [(a,a)]
+pairs (x0:x1:xs) = (x0,x1) : pairs xs
+pairs [] = []
+pairs _ = error "pairs: odd number of elements"
 
 
 parabolaDerivative ::
diff --git a/private/Numeric/Interpolation/Private/List.hs b/private/Numeric/Interpolation/Private/List.hs
new file mode 100644
--- /dev/null
+++ b/private/Numeric/Interpolation/Private/List.hs
@@ -0,0 +1,27 @@
+module Numeric.Interpolation.Private.List where
+
+import qualified Data.List as List
+
+
+mapAdjacent3 :: (a -> a -> a -> b) -> [a] -> [b]
+mapAdjacent3 f xs0 =
+   let xs1 = drop 1 xs0
+       xs2 = drop 1 xs1
+   in  List.zipWith3 f xs0 xs1 xs2
+
+mapAdjacentMaybe3 :: (Maybe a -> a -> Maybe a -> b) -> [a] -> [b]
+mapAdjacentMaybe3 f xs =
+   let jxs = map Just xs
+   in  zipWith3 f (Nothing : jxs) xs (drop 1 jxs ++ [Nothing])
+
+mapAdjacentMaybe5 ::
+   (Maybe a -> Maybe a -> a -> Maybe a -> Maybe a -> b) ->
+   [a] -> [b]
+mapAdjacentMaybe5 f xs =
+   let jxs = map Just xs
+       lxs1 = Nothing : jxs
+       lxs2 = Nothing : lxs1
+       rxs1 = drop 1 $ jxs ++ repeat Nothing
+       rxs2 = drop 1 $ rxs1
+   in  List.zipWith5 f lxs2 lxs1 xs rxs1 rxs2
+
diff --git a/src/Numeric/Interpolation/Basis.hs b/src/Numeric/Interpolation/Basis.hs
--- a/src/Numeric/Interpolation/Basis.hs
+++ b/src/Numeric/Interpolation/Basis.hs
@@ -13,8 +13,7 @@
 import qualified Numeric.Interpolation.NodeList as Nodes
 import Numeric.Interpolation.Private.Basis
           (parabolaDerivativeCenterNode, hermite1Split)
-
-import qualified Data.List as List
+import Numeric.Interpolation.Private.List (mapAdjacent3, )
 
 
 {- |
@@ -34,12 +33,6 @@
    Nodes.fromList . zip xs . hermite1Split xs
 
 
-
-mapAdjacent3 :: (a -> a -> a -> b) -> [a] -> [b]
-mapAdjacent3 f xs0 =
-   let xs1 = drop 1 xs0
-       xs2 = drop 1 xs1
-   in  List.zipWith3 f xs0 xs1 xs2
 
 {- |
 Cf. 'coefficientsToLinear'
diff --git a/src/Numeric/Interpolation/Basis/Compact.hs b/src/Numeric/Interpolation/Basis/Compact.hs
--- a/src/Numeric/Interpolation/Basis/Compact.hs
+++ b/src/Numeric/Interpolation/Basis/Compact.hs
@@ -8,18 +8,16 @@
    parabolaBasisDerivativeCenter,
    parabolaBasisDerivativeLeft,
    )
+import Numeric.Interpolation.Private.List (
+   mapAdjacentMaybe3,
+   mapAdjacentMaybe5,
+   )
 
 import Control.Monad (liftM, liftM2)
 
-import qualified Data.List as List
 import Data.Maybe (catMaybes)
 
 
-mapAdjacentMaybe3 :: (Maybe a -> a -> Maybe a -> b) -> [a] -> [b]
-mapAdjacentMaybe3 f xs =
-   let jxs = map Just xs
-   in  zipWith3 f (Nothing : jxs) xs (drop 1 jxs ++ [Nothing])
-
 generic :: ny -> ny -> [x] -> [Nodes.T x ny]
 generic nz ny =
    mapAdjacentMaybe3
@@ -34,23 +32,13 @@
 
 hermite1 :: (Num b) => [a] -> [Nodes.T a (b, b)]
 hermite1 xs =
-   generic (0,0) (1,0) xs
-   ++
-   generic (0,0) (0,1) xs
-
+   concat $
+   zipWith (\f df -> [f,df])
+      (generic (0,0) (1,0) xs)
+      (generic (0,0) (0,1) xs)
 
 
 
-mapAdjacentMaybe5 ::
-   (Maybe a -> Maybe a -> a -> Maybe a -> Maybe a -> b) ->
-   [a] -> [b]
-mapAdjacentMaybe5 f xs =
-   let jxs = map Just xs
-       lxs1 = Nothing : jxs
-       lxs2 = Nothing : lxs1
-       rxs1 = drop 1 $ jxs ++ repeat Nothing
-       rxs2 = drop 1 $ rxs1
-   in  List.zipWith5 f lxs2 lxs1 xs rxs1 rxs2
 
 cubicAutoGeneric ::
    (Num b) =>
diff --git a/src/Numeric/Interpolation/NodeList.hs b/src/Numeric/Interpolation/NodeList.hs
--- a/src/Numeric/Interpolation/NodeList.hs
+++ b/src/Numeric/Interpolation/NodeList.hs
@@ -8,11 +8,37 @@
 
 import Data.Tuple.HT (mapFst)
 
+import Data.Traversable (Traversable, traverse)
+import Data.Foldable (Foldable, foldMap)
+import Data.Monoid (mempty, (<>))
+
+import Control.Applicative (liftA3, pure)
+
 import Prelude hiding (lookup)
 
 
 data T x y = Interval | Node (x, y) (T x y) (T x y)
    deriving (Eq, Ord, Show)
+
+instance Functor (T x) where
+   fmap f =
+      let go Interval = Interval
+          go (Node (x,y) l r) = Node (x, f y) (go l) (go r)
+      in  go
+
+instance Foldable (T x) where
+   foldMap f =
+      let go Interval = mempty
+          go (Node (_x,y) l r) = go l <> f y <> go r
+      in  go
+
+instance Traversable (T x) where
+   traverse f =
+      let go Interval = pure Interval
+          go (Node (x,y) l0 r0) =
+             liftA3 (\l m r -> Node (x,m) l r) (go l0) (f y) (go r0)
+      in  go
+
 
 {- |
 list must be sorted with respect to first element
diff --git a/src/Numeric/Interpolation/Sample.hs b/src/Numeric/Interpolation/Sample.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Sample.hs
@@ -0,0 +1,113 @@
+module Numeric.Interpolation.Sample (
+   T,
+   linear,
+   hermite1,
+   cubicLinear,
+   cubicParabola,
+   ) where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import qualified Numeric.Interpolation.Piece as Piece
+import Numeric.Interpolation.Private.List (mapAdjacentMaybe3, )
+import Numeric.Interpolation.Private.Basis (
+   parabolaBasisDerivativeRight,
+   parabolaBasisDerivativeCenter,
+   parabolaBasisDerivativeLeft,
+   )
+
+
+type T x y = [x] -> x -> [(Int, y)]
+
+linear :: (Fractional a, Ord a) => T a a
+linear nodeXs =
+   let nodes = Nodes.fromList $ zip nodeXs [0..]
+   in  \x ->
+          case Nodes.lookup nodes x of
+             (Just (l,nl), Just (r,nr)) ->
+                [(nl, Piece.linear (l,1) (r,0) x),
+                 (nr, Piece.linear (l,0) (r,1) x)]
+             (Just (_l,nl), Nothing) -> [(nl, 1)]
+             (Nothing, Just (_r,nr)) -> [(nr, 1)]
+             (Nothing, Nothing) -> []
+
+hermite1 :: (Fractional a, Ord a) => T a a
+hermite1 nodeXs =
+   let nodes = Nodes.fromList $ zip nodeXs [0..]
+   in  \x ->
+          case Nodes.lookup nodes x of
+             (Just (l,nl), Just (r,nr)) ->
+                [(2*nl+0, Piece.hermite1 (l,(1,0)) (r,(0,0)) x),
+                 (2*nl+1, Piece.hermite1 (l,(0,1)) (r,(0,0)) x),
+                 (2*nr+0, Piece.hermite1 (l,(0,0)) (r,(1,0)) x),
+                 (2*nr+1, Piece.hermite1 (l,(0,0)) (r,(0,1)) x)]
+             (Just (_l,nl), Nothing) -> [(2*nl, 1)]
+             (Nothing, Just (_r,nr)) -> [(2*nr, 1)]
+             (Nothing, Nothing) -> []
+
+cubicLinear :: (Fractional a, Ord a) => T a a
+cubicLinear nodeXs =
+   let nodes =
+          Nodes.fromList $ zip nodeXs $ zip [0..] $
+          mapAdjacentMaybe3 (\l _ r -> (l,r)) nodeXs
+   in  \x ->
+          case Nodes.lookup nodes x of
+             (Nothing, Nothing) -> []
+             (Just (_l,(nl,_)), Nothing) -> [(nl-1, 1)]
+             (Nothing, Just (_r,(nr,_))) -> [(nr+1, 1)]
+             (Just (l,(nl,(mll,_))), Just (r,(nr,(_,mrr)))) ->
+                let interL ll =
+                       (nl-1, Piece.hermite1 (l,(0,recip(ll-r))) (r,(0,0)) x)
+                    interR rr =
+                       (nr+1, Piece.hermite1 (l,(0,0)) (r,(0,recip(rr-l))) x)
+                in  case (mll,mrr) of
+                       (Just ll, Just rr) ->
+                          interL ll :
+                          (nl, Piece.hermite1 (l,(1,0)) (r,(0,recip(l-rr))) x) :
+                          (nr, Piece.hermite1 (l,(0,recip(r-ll))) (r,(1,0)) x) :
+                          interR rr :
+                          []
+                       (Just ll, Nothing) -> interL ll : [(nl, 1)]
+                       (Nothing, Just rr) -> interR rr : [(nr, 1)]
+                       (Nothing, Nothing) -> []
+
+cubicParabola :: (Fractional a, Ord a) => T a a
+cubicParabola nodeXs =
+   let nodes =
+          Nodes.fromList $ zip nodeXs $ zip [0..] $
+          mapAdjacentMaybe3 (\l _ r -> (l,r)) nodeXs
+   in  \x ->
+          case Nodes.lookup nodes x of
+             (Nothing, Nothing) -> []
+             (Just (_l,(nl,_)), Nothing) -> [(nl-1, 1)]
+             (Nothing, Just (_r,(nr,_))) -> [(nr+1, 1)]
+             (Just (l,(nl,(mll,_))), Just (r,(nr,(_,mrr)))) ->
+                let interL ll =
+                       (nl-1,
+                        Piece.hermite1
+                           (l,(0, parabolaBasisDerivativeLeft ll l r))
+                           (r,(0, 0))
+                           x)
+                    interR rr =
+                       (nr+1,
+                        Piece.hermite1
+                           (l,(0, 0))
+                           (r,(0, parabolaBasisDerivativeRight l r rr))
+                           x)
+                in  case (mll,mrr) of
+                       (Just ll, Just rr) ->
+                          interL ll :
+                          (nl,
+                           Piece.hermite1
+                              (l, (1, parabolaBasisDerivativeCenter ll l r))
+                              (r, (0, parabolaBasisDerivativeLeft l r rr))
+                              x) :
+                          (nr,
+                           Piece.hermite1
+                              (l, (0, parabolaBasisDerivativeRight ll l r))
+                              (r, (1, parabolaBasisDerivativeCenter l r rr))
+                              x) :
+                          interR rr :
+                          []
+                       (Just ll, Nothing) -> interL ll : [(nl, 1)]
+                       (Nothing, Just rr) -> interR rr : [(nr, 1)]
+                       (Nothing, Nothing) -> []
diff --git a/src/Numeric/Interpolation/Type.hs b/src/Numeric/Interpolation/Type.hs
--- a/src/Numeric/Interpolation/Type.hs
+++ b/src/Numeric/Interpolation/Type.hs
@@ -1,7 +1,7 @@
 module Numeric.Interpolation.Type (
    T(..),
    linear,
-   cubic,
+   hermite1,
    cubicLinear,
    cubicParabola,
    ) where
@@ -9,6 +9,7 @@
 import qualified Numeric.Interpolation.NodeList as Nodes
 import qualified Numeric.Interpolation.Piece as Piece
 import qualified Numeric.Interpolation.Basis as Basis
+import qualified Numeric.Interpolation.Sample as Sample
 import Numeric.Interpolation.Private.Basis (hermite1Split)
 
 
@@ -16,24 +17,31 @@
    Cons {
       ssvFromNodes :: [x] -> [y] -> String,
       interpolatePiece :: Piece.T x y ny,
+      basisOverlap :: Int
+         {- ^
+         maximum difference of indices of basis functions that overlap plus one
+         -},
       basisFunctions :: [x] -> [Nodes.T x ny],
+      sampleBasisFunctions :: [x] -> x -> [(Int, y)],
       coefficientsToInterpolator :: [x] -> [y] -> Nodes.T x ny,
       valueFromNode :: ny -> y
    }
 
-linear :: T Double Double Double
+linear :: (Fractional a, Ord a, Show a) => T a a a
 linear =
    Cons {
       ssvFromNodes =
          \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
       interpolatePiece = Piece.linear,
+      basisOverlap = 2,
       basisFunctions = Basis.linear,
+      sampleBasisFunctions = Sample.linear,
       coefficientsToInterpolator = Basis.coefficientsToLinear,
       valueFromNode = id
    }
 
-cubic :: T Double Double (Double, Double)
-cubic =
+hermite1 :: (Fractional a, Ord a, Show a) => T a a (a, a)
+hermite1 =
    Cons {
       ssvFromNodes =
          \xs ys ->
@@ -41,41 +49,49 @@
             zipWith (\x (y,dy) -> show x ++ " " ++ show y ++ " " ++ show dy) xs $
             hermite1Split xs ys,
       interpolatePiece = Piece.hermite1,
+      basisOverlap = 4,
       basisFunctions = Basis.hermite1,
+      sampleBasisFunctions = Sample.hermite1,
       coefficientsToInterpolator = Basis.coefficientsToHermite1,
       valueFromNode = fst
    }
 
-cubicLinear :: T Double Double (Double, Double)
+cubicLinear :: (Fractional a, Ord a, Show a) => T a a (a, a)
 cubicLinear =
    Cons {
       ssvFromNodes =
          \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
       interpolatePiece = Piece.hermite1,
+      basisOverlap = 4,
       basisFunctions = Basis.cubicLinear,
+      sampleBasisFunctions = Sample.cubicLinear,
       coefficientsToInterpolator = Basis.coefficientsToCubicLinear,
       valueFromNode = fst
    }
 
-cubicParabola :: T Double Double (Double, Double)
+cubicParabola :: (Fractional a, Ord a, Show a) => T a a (a, a)
 cubicParabola =
    Cons {
       ssvFromNodes =
          \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
       interpolatePiece = Piece.hermite1,
+      basisOverlap = 4,
       basisFunctions = Basis.cubicParabola,
+      sampleBasisFunctions = Sample.cubicParabola,
       coefficientsToInterpolator = Basis.coefficientsToCubicParabola,
       valueFromNode = fst
    }
 
 
-_cubicMean :: T Double Double (Double, Double)
+_cubicMean :: (Fractional a, Ord a, Show a) => T a a (a, a)
 _cubicMean =
    Cons {
       ssvFromNodes =
          \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
       interpolatePiece = Piece.hermite1,
+      basisOverlap = 4,
       basisFunctions = Basis.cubicParabola, -- Basis.cubicMean,
+      sampleBasisFunctions = Sample.cubicParabola, -- Sample.cubicMean,
       coefficientsToInterpolator = Basis.coefficientsToCubicParabola, -- not correct
       valueFromNode = fst
    }
diff --git a/test/Test.hs b/test/Test.hs
--- a/test/Test.hs
+++ b/test/Test.hs
@@ -1,6 +1,8 @@
 module Main where
 
 import qualified Test.Piece as Piece
+import qualified Test.Sample as Sample
+import qualified Test.Overlap as Overlap
 
 
 run :: String -> [(String, IO ())] -> IO ()
@@ -8,5 +10,7 @@
    mapM_ (\(msg,act) -> putStr (prefix ++ '.' : msg ++ ": ") >> act)
 
 main :: IO ()
-main =
+main = do
    run "Piece" Piece.tests
+   run "Sample" Sample.tests
+   run "Overlap" Overlap.tests
diff --git a/test/Test/Overlap.hs b/test/Test/Overlap.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Overlap.hs
@@ -0,0 +1,25 @@
+module Test.Overlap where
+
+import qualified Numeric.Interpolation.Type as Type
+
+import Test.QuickCheck (quickCheck, )
+
+
+
+test :: Type.T Double y ny -> IO ()
+test typ =
+   quickCheck $ \xs xi ->
+      let samples = map fst $ Type.sampleBasisFunctions typ xs xi
+          {- not total:
+          maximum samples - minimum samples < Type.basisOverlap typ
+          -}
+      in  all (< minimum samples + Type.basisOverlap typ) samples
+
+
+tests :: [(String, IO ())]
+tests =
+   ("linear", test Type.linear) :
+   ("hermite1", test Type.hermite1) :
+   ("cubicLinear", test Type.cubicLinear) :
+   ("cubicParabola", test Type.cubicParabola) :
+   []
diff --git a/test/Test/Sample.hs b/test/Test/Sample.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Sample.hs
@@ -0,0 +1,62 @@
+module Test.Sample where
+
+import qualified Numeric.Interpolation.Type as Type
+import qualified Numeric.Interpolation.Piecewise as Piecewise
+
+import qualified Data.Set as Set
+import Data.Array (accumArray, listArray, )
+import Data.List.HT (lengthAtLeast, )
+
+import Test.QuickCheck (Property, quickCheck, (==>), )
+
+
+
+withSortedRatios ::
+   ([Rational] -> Rational -> a) ->
+   ([Integer] -> Integer -> a)
+withSortedRatios f nodeXs x =
+   f (map fromInteger $ Set.toAscList $ Set.fromList nodeXs) (fromInteger x)
+
+checkEq ::
+   (Ord x, Eq y, Num y) =>
+   Type.T x y ny -> [x] -> x -> Bool
+checkEq typ nodeXs x =
+   let ys =
+          map
+             (flip (Piecewise.interpolateConstantExt typ) x)
+             (Type.basisFunctions typ nodeXs)
+       bounds = (0, length ys - 1)
+   in  listArray bounds ys
+       ==
+       accumArray (flip const) 0 bounds
+          (Type.sampleBasisFunctions typ nodeXs x)
+
+
+linear :: [Integer] -> Integer -> Bool
+linear = withSortedRatios $ checkEq Type.linear
+
+hermite1 :: [Integer] -> Integer -> Bool
+hermite1 = withSortedRatios $ checkEq Type.hermite1
+
+
+derivativeFree ::
+   Type.T Rational Rational ny ->
+   [Integer] -> Integer -> Property
+derivativeFree typ =
+   withSortedRatios $ \nodeXs x ->
+      lengthAtLeast 4 nodeXs ==> checkEq typ nodeXs x
+
+cubicLinear :: [Integer] -> Integer -> Property
+cubicLinear = derivativeFree Type.cubicLinear
+
+cubicParabola :: [Integer] -> Integer -> Property
+cubicParabola = derivativeFree Type.cubicParabola
+
+
+tests :: [(String, IO ())]
+tests =
+   ("linear", quickCheck linear) :
+   ("hermite1", quickCheck hermite1) :
+   ("cubicLinear", quickCheck cubicLinear) :
+   ("cubicParabola", quickCheck cubicParabola) :
+   []