diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,26 @@
+Copyright (c) Henning Thielemann 2014
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+3. Neither the name of the University nor the names of its contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,3 @@
+#! /usr/bin/env runhaskell
+> import Distribution.Simple
+> main = defaultMain
diff --git a/example/Fit.hs b/example/Fit.hs
new file mode 100644
--- /dev/null
+++ b/example/Fit.hs
@@ -0,0 +1,71 @@
+module Main where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import qualified Numeric.Interpolation.Piecewise as Piecewise
+import qualified Numeric.Interpolation.Type as Type
+
+import qualified Data.Packed.Matrix as Matrix
+import qualified Data.Packed.Vector as Vector
+import qualified Numeric.Container as Container
+import Numeric.Container ((<\>))
+
+import qualified Graphics.Gnuplot.Advanced as GP
+import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
+import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
+
+import System.Random (randomRs, mkStdGen)
+import Control.Monad.HT (void)
+import Data.Monoid ((<>))
+
+
+noisy :: [(Double, Double)]
+noisy =
+   take 100 $
+   zipWith
+      (\x d -> (x, sin x + d))
+      (randomRs (0,7) (mkStdGen 23))
+      (randomRs (-0.2,0.2) (mkStdGen 42))
+
+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
+   in  Type.coefficientsToInterpolator typ xs $
+       Vector.toList $ matrix <\> tys
+
+plotBasisFunctions ::
+   Type.T Double Double ny -> [Double] -> Plot2D.T Double Double
+plotBasisFunctions nodeType xs =
+   let abscissa = Plot2D.linearScale 1000 (minimum xs, maximum xs)
+   in  Plot2D.functions Graph2D.lines abscissa $
+       map (Piecewise.interpolateConstantExt nodeType) $
+       Type.basisFunctions nodeType xs
+
+
+main :: IO ()
+main = do
+   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.cubicLinear exs
+   void $ GP.plotDefault $ plotBasisFunctions Type.cubicParabola exs
+   let linearNodes = fit Type.linear xs noisy
+       hermite1Nodes = fit Type.cubic xs noisy
+       cubicLinearNodes = fit Type.cubicLinear exs noisy
+       cubicParabolaNodes = fit Type.cubicParabola exs noisy
+   void $ GP.plotDefault $
+      Plot2D.list Graph2D.points noisy
+      <>
+      (Plot2D.functions Graph2D.lines (Plot2D.linearScale 1000 (-2,10)) $
+       Piecewise.interpolateConstantExt Type.linear linearNodes :
+       Piecewise.interpolateConstantExt Type.cubic hermite1Nodes :
+       Piecewise.interpolateConstantExt Type.cubicLinear cubicLinearNodes :
+       Piecewise.interpolateConstantExt Type.cubicParabola cubicParabolaNodes :
+       [])
diff --git a/example/Plot.hs b/example/Plot.hs
new file mode 100644
--- /dev/null
+++ b/example/Plot.hs
@@ -0,0 +1,33 @@
+module Main where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import qualified Numeric.Interpolation.Piecewise as Piecewise
+import qualified Numeric.Interpolation.Basis as Basis
+import qualified Numeric.Interpolation.Type as Type
+
+
+import qualified Graphics.Gnuplot.Advanced as GP
+
+import qualified Graphics.Gnuplot.Plot.TwoDimensional as Plot2D
+import qualified Graphics.Gnuplot.Graph.TwoDimensional as Graph2D
+
+import Control.Monad.HT (void)
+
+
+xs :: [Double]
+xs = [0, 1, 3, 4, 6, 7, 9, 10, 11, 13]
+
+main :: IO ()
+main = do
+   let linearNodes = Nodes.fromList $ map (\x -> (x, sin x)) xs
+       hermite1Nodes = Nodes.fromList $ map (\x -> (x, (sin x, cos x))) xs
+       cubicLinearNodes = Basis.coefficientsToCubicLinear xs $ map sin xs
+       cubicParabolaNodes = Basis.coefficientsToCubicParabola xs $ map sin xs
+   void $ GP.plotDefault $
+      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,
+          sin]
diff --git a/interpolation.cabal b/interpolation.cabal
new file mode 100644
--- /dev/null
+++ b/interpolation.cabal
@@ -0,0 +1,112 @@
+Name:             interpolation
+Version:          0.0
+License:          BSD3
+License-File:     LICENSE
+Author:           Henning Thielemann
+Maintainer:       Henning Thielemann <haskell@henning-thielemann.de>
+Homepage:         http://code.haskell.org/~thielema/interpolation/
+Category:         Math
+Synopsis:         piecewise linear and cubic Hermite interpolation
+Description:
+  Represent real functions by linear or cubic segments.
+  The package provides both data structures
+  for efficient lookup of interpolation intervals,
+  and computation of basis functions.
+  .
+  There are two examples that can be built with
+  .
+  > cabal install -fbuildExamples
+  .
+  * @example/Plot.hs@:
+    Interpolate a sinus curve using piecewise linear interpolation
+    and piecewise Hermite cubic interpolation.
+    For the latter one we provide the derivatives of the sinus function
+    at the interpolation nodes.
+  .
+  * @example/Fit.hs@:
+    Demonstrates how to use the basis functions
+    for fitting an interpolation function to a given function
+    using a linear least squares solver like @<\>@ from @hmatrix@.
+    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
+Cabal-Version:    >=1.8
+Build-Type:       Simple
+
+Flag buildExamples
+  description: Build example executables
+  default:     False
+
+Flag hmatrix
+  description: Build examples that depend on hmatrix
+  default:     True
+
+Source-Repository this
+  Tag:         0.0
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/interpolation/
+
+Source-Repository head
+  Type:        darcs
+  Location:    http://code.haskell.org/~thielema/interpolation/
+
+Library
+  Build-Depends:
+    utility-ht >=0.0.1 && <0.1,
+    base >=4 && <5
+
+  GHC-Options:      -Wall
+  Hs-Source-Dirs:   src, private
+  Exposed-Modules:
+    Numeric.Interpolation.NodeList
+    Numeric.Interpolation.Piece
+    Numeric.Interpolation.Piecewise
+    Numeric.Interpolation.Type
+    Numeric.Interpolation.Basis
+    Numeric.Interpolation.Basis.Compact
+    Numeric.Interpolation.Basis.Full
+  Other-Modules:
+    Numeric.Interpolation.Private.Piece
+    Numeric.Interpolation.Private.Basis
+
+Executable interpolation-plot
+  Main-Is:             Plot.hs
+  Hs-Source-Dirs:      example
+  GHC-Options:         -Wall
+  If flag(buildExamples)
+    Build-Depends:
+      interpolation,
+      gnuplot >=0.5.2 && <0.6,
+      utility-ht >=0.0.9 && <0.1,
+      base >=4.5 && <4.8
+  Else
+    Buildable: False
+
+Executable interpolation-fit
+  Main-Is:             Fit.hs
+  Hs-Source-Dirs:      example
+  GHC-Options:         -Wall
+  If flag(buildExamples) && flag(hmatrix)
+    Build-Depends:
+      interpolation,
+      hmatrix >=0.15 && <0.16,
+      random >=1.0 && <1.1,
+      gnuplot >=0.5.2 && <0.6,
+      utility-ht >=0.0.9 && <0.1,
+      base >=4.5 && <4.8
+  Else
+    Buildable: False
+
+Test-Suite interpolation-test
+  Type:                exitcode-stdio-1.0
+  Main-Is:             Test.hs
+  Other-Modules:
+    Test.Piece
+  Hs-Source-Dirs:      test, private
+  GHC-Options:         -Wall
+  Build-Depends:
+    interpolation,
+    QuickCheck >=2.4 && <2.8,
+    utility-ht >=0.0.9 && <0.1,
+    base >=4.5 && <4.8
diff --git a/private/Numeric/Interpolation/Private/Basis.hs b/private/Numeric/Interpolation/Private/Basis.hs
new file mode 100644
--- /dev/null
+++ b/private/Numeric/Interpolation/Private/Basis.hs
@@ -0,0 +1,53 @@
+module Numeric.Interpolation.Private.Basis where
+
+import Numeric.Interpolation.Private.Piece (sqr)
+
+import qualified Data.List.Match as Match
+
+
+hermite1Split :: [a] -> [b] -> [(b, b)]
+hermite1Split xs = uncurry zip . Match.splitAt xs
+
+
+parabolaDerivative ::
+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> a -> (a,a)
+parabolaDerivative (x0,y0) (x1,y1) (x2,y2) x =
+   let l0 = (x-x1)*(x-x2)/((x0-x1)*(x0-x2))
+       l1 = (x-x0)*(x-x2)/((x1-x0)*(x1-x2))
+       l2 = (x-x0)*(x-x1)/((x2-x0)*(x2-x1))
+       dl0 = (2*x-x1-x2)/((x0-x1)*(x0-x2))
+       dl1 = (2*x-x0-x2)/((x1-x0)*(x1-x2))
+       dl2 = (2*x-x0-x1)/((x2-x0)*(x2-x1))
+   in  (y0*l0 + y1*l1 + y2*l2, y0*dl0 + y1*dl1 + y2*dl2)
+
+parabolaBasisDerivativeLeft,
+   parabolaBasisDerivativeCenter,
+   parabolaBasisDerivativeRight ::
+      (Fractional a) => a -> a -> a -> a
+parabolaBasisDerivativeLeft   x0 x1 x2 = (x1-x2)/((x0-x1)*(x0-x2))
+parabolaBasisDerivativeCenter x0 x1 x2 = 1/(x1-x0) + 1/(x1-x2)
+parabolaBasisDerivativeRight  x0 x1 x2 = (x1-x0)/((x2-x0)*(x2-x1))
+
+parabolaDerivativeCenterNode ::
+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> a
+parabolaDerivativeCenterNode (x0,y0) (x1,y1) (x2,y2) =
+   y0 * parabolaBasisDerivativeLeft   x0 x1 x2 +
+   y1 * parabolaBasisDerivativeCenter x0 x1 x2 +
+   y2 * parabolaBasisDerivativeRight  x0 x1 x2
+
+
+parabola2ndDerivativeCenterNode ::
+   (Fractional a) => (a,a) -> (a,a) -> (a,a) -> (a,a) -> a
+parabola2ndDerivativeCenterNode (xl,yl) (x0,y0) (x1,y1) (x2,y2) =
+   let dy0 =
+          yl * (x0-x1)/((xl-x0)*(xl-x1)) +
+          y0 * (1/(x0-xl) + 1/(x0-x1)) +
+          y1 * (x0-xl)/((x1-xl)*(x1-x0))
+       dy1 =
+          y0 * (x1-x2)/((x0-x1)*(x0-x2)) +
+          y1 * (1/(x1-x0) + 1/(x1-x2)) +
+          y2 * (x1-x0)/((x2-x0)*(x2-x1))
+       d = (y1-y0)/(x1-x0)
+       x = x0
+   in  2*(dy0-d) / sqr (x0-x1) * (3*x-2*x1-x0) +
+       2*(dy1-d) / sqr (x1-x0) * (3*x-2*x0-x1)
diff --git a/private/Numeric/Interpolation/Private/Piece.hs b/private/Numeric/Interpolation/Private/Piece.hs
new file mode 100644
--- /dev/null
+++ b/private/Numeric/Interpolation/Private/Piece.hs
@@ -0,0 +1,31 @@
+module Numeric.Interpolation.Private.Piece where
+
+sqr :: (Num a) => a -> a
+sqr x = x*x
+
+
+type T x y ny = (x, ny) -> (x, ny) -> x -> y
+
+linear :: (Fractional a) => T a a a
+linear (x0,y0) (x1,y1) x =
+   (y0*(x1-x) + y1*(x-x0)) / (x1-x0)
+
+hermite1 :: (Fractional a) => T a a (a, a)
+hermite1 (x0,(y0,dy0)) (x1,(y1,dy1)) x =
+   let d = (y1-y0)/(x1-x0)
+   in  linear (x0,y0) (x1,y1) x +
+       (dy0-d) * sqr ((x-x1)/(x0-x1)) * (x-x0) +
+       (dy1-d) * sqr ((x-x0)/(x1-x0)) * (x-x1)
+
+hermite1' :: (Fractional a) => T a a (a, a)
+hermite1' (x0,(y0,dy0)) (x1,(y1,dy1)) x =
+   let d = (y1-y0)/(x1-x0)
+   in  d +
+       (dy0-d) / sqr (x0-x1) * (2*(x-x1) * (x-x0) + sqr (x-x1)) +
+       (dy1-d) / sqr (x1-x0) * (2*(x-x0) * (x-x1) + sqr (x-x0))
+
+hermite1'' :: (Fractional a) => T a a (a, a)
+hermite1'' (x0,(y0,dy0)) (x1,(y1,dy1)) x =
+   let d = (y1-y0)/(x1-x0)
+   in  2*(dy0-d) / sqr (x0-x1) * (3*x-2*x1-x0) +
+       2*(dy1-d) / sqr (x1-x0) * (3*x-2*x0-x1)
diff --git a/src/Numeric/Interpolation/Basis.hs b/src/Numeric/Interpolation/Basis.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Basis.hs
@@ -0,0 +1,62 @@
+module Numeric.Interpolation.Basis (
+   Compact.linear,
+   Compact.hermite1,
+   Compact.cubicLinear,
+   Compact.cubicParabola,
+   coefficientsToLinear,
+   coefficientsToHermite1,
+   coefficientsToCubicLinear,
+   coefficientsToCubicParabola,
+   ) where
+
+import qualified Numeric.Interpolation.Basis.Compact as Compact
+import qualified Numeric.Interpolation.NodeList as Nodes
+import Numeric.Interpolation.Private.Basis
+          (parabolaDerivativeCenterNode, hermite1Split)
+
+import qualified Data.List as List
+
+
+{- |
+@coefficientsToLinear nodes coefficients@
+creates an interpolation function for @nodes@,
+where the @coefficients@ correspond to the basis functions
+constructed with @Basis.linear nodes@.
+-}
+coefficientsToLinear :: [a] -> [b] -> Nodes.T a b
+coefficientsToLinear xs = Nodes.fromList . zip xs
+
+{- |
+Cf. 'coefficientsToLinear'
+-}
+coefficientsToHermite1 :: [a] -> [b] -> Nodes.T a (b, b)
+coefficientsToHermite1 xs =
+   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'
+-}
+coefficientsToCubicLinear :: (Fractional a) => [a] -> [a] -> Nodes.T a (a, a)
+coefficientsToCubicLinear xs =
+   Nodes.fromList .
+   mapAdjacent3 (\(xl,yl) (xn,yn) (xr,yr) -> (xn, (yn, (yr-yl)/(xr-xl)))) .
+   zip xs
+
+{- |
+Cf. 'coefficientsToLinear'
+-}
+coefficientsToCubicParabola :: (Fractional a) => [a] -> [a] -> Nodes.T a (a, a)
+coefficientsToCubicParabola xs =
+   Nodes.fromList .
+   mapAdjacent3
+      (\pl pn@(xn,yn) pr ->
+         (xn, (yn, parabolaDerivativeCenterNode pl pn pr))) .
+   zip xs
diff --git a/src/Numeric/Interpolation/Basis/Compact.hs b/src/Numeric/Interpolation/Basis/Compact.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Basis/Compact.hs
@@ -0,0 +1,109 @@
+module Numeric.Interpolation.Basis.Compact (
+   linear, hermite1, cubicLinear, cubicParabola,
+   ) where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import Numeric.Interpolation.Private.Basis (
+   parabolaBasisDerivativeRight,
+   parabolaBasisDerivativeCenter,
+   parabolaBasisDerivativeLeft,
+   )
+
+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
+      (\l n r ->
+          Nodes.Node (n,ny)
+             (maybe Nodes.Interval (flip Nodes.singleton nz) l)
+             (maybe Nodes.Interval (flip Nodes.singleton nz) r))
+
+
+linear :: (Num b) => [a] -> [Nodes.T a b]
+linear = generic 0 1
+
+hermite1 :: (Num b) => [a] -> [Nodes.T a (b, b)]
+hermite1 xs =
+   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) =>
+   (a -> a -> a -> b) ->
+   (a -> a -> a -> b) ->
+   (a -> a -> a -> b) ->
+   [a] -> [Nodes.T a (b, b)]
+cubicAutoGeneric dl dn dr =
+   mapAdjacentMaybe5
+      (\ml2 ml1 n mr1 mr2 ->
+         let node x y y' = (x, (y,y'))
+         in  Nodes.fromList $ catMaybes $
+             liftM (\l2 -> node l2 0 0) ml2 :
+             liftM2 (\l2 l1 -> node l1 0 (dl l2 l1 n)) ml2 ml1 :
+             liftM2 (\l1 r1 -> node n  1 (dn l1 n r1)) ml1 mr1 :
+             liftM2 (\r1 r2 -> node r1 0 (dr n r1 r2)) mr1 mr2 :
+             liftM (\r2 -> node r2 0 0) mr2 :
+             [])
+
+
+{- |
+Cubic interpolation
+where the derivative at a node is set to the slope of the two adjacent nodes.
+-}
+cubicLinear :: (Fractional a) => [a] -> [Nodes.T a (a, a)]
+cubicLinear =
+   cubicAutoGeneric
+      (\ll _l n -> recip $ n-ll)
+      (\_l _n _r -> 0)
+      (\n _r rr -> recip $ n-rr)
+
+
+{- |
+Cubic interpolation
+where the derivative at a node is set to the slope of the parabola
+through the current and the two adjacent nodes.
+-}
+cubicParabola :: (Fractional a) => [a] -> [Nodes.T a (a, a)]
+cubicParabola =
+   cubicAutoGeneric
+      parabolaBasisDerivativeRight
+      parabolaBasisDerivativeCenter
+      parabolaBasisDerivativeLeft
+
+
+{- |
+Experimental interpolation
+which is mean of 'cubicLinear' and 'cubicParabola'.
+The result looks reasonable, too.
+-}
+_cubicMean :: (Fractional a) => [a] -> [Nodes.T a (a, a)]
+_cubicMean =
+   cubicAutoGeneric
+      (\ll l n -> (parabolaBasisDerivativeRight ll l n + recip (n-ll))/2)
+      (\l n r -> parabolaBasisDerivativeCenter l n r / 2)
+      (\n r rr -> (parabolaBasisDerivativeLeft n r rr + recip (n-rr))/2)
diff --git a/src/Numeric/Interpolation/Basis/Full.hs b/src/Numeric/Interpolation/Basis/Full.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Basis/Full.hs
@@ -0,0 +1,20 @@
+module Numeric.Interpolation.Basis.Full (linear, hermite1) where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+
+import qualified Data.List.Match as Match
+
+
+generic :: ny -> ny -> [x] -> [Nodes.T x ny]
+generic nz ny xs =
+   map (Nodes.fromList . zip xs) $
+   Match.take xs $ iterate (nz:) $ ny : repeat nz
+
+linear :: (Num b) => [a] -> [Nodes.T a b]
+linear = generic 0 1
+
+hermite1 :: (Num b) => [a] -> [Nodes.T a (b, b)]
+hermite1 xs =
+   generic (0,0) (1,0) xs
+   ++
+   generic (0,0) (0,1) xs
diff --git a/src/Numeric/Interpolation/NodeList.hs b/src/Numeric/Interpolation/NodeList.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/NodeList.hs
@@ -0,0 +1,48 @@
+module Numeric.Interpolation.NodeList (
+   T(Interval, Node),
+   fromList,
+   toList,
+   singleton,
+   lookup,
+   ) where
+
+import Data.Tuple.HT (mapFst)
+
+import Prelude hiding (lookup)
+
+
+data T x y = Interval | Node (x, y) (T x y) (T x y)
+   deriving (Eq, Ord, Show)
+
+{- |
+list must be sorted with respect to first element
+-}
+fromList :: [(x,y)] -> T x y
+fromList =
+   let merge n0 xys0 =
+          case xys0 of
+             (xy0,n1):(xy1,n2):xys ->
+                (Node xy0 n0 n1,
+                 uncurry (:) $ mapFst ((,) xy1) $ merge n2 xys)
+             (xy0,n1):[] -> (Node xy0 n0 n1, [])
+             [] -> (n0, [])
+       rep (n,xyns) = if null xyns then n else rep $ merge n xyns
+   in  rep . merge Interval . map (flip (,) Interval)
+
+singleton :: x -> y -> T x y
+singleton x y = Node (x,y) Interval Interval
+
+toList :: T x y -> [(x,y)]
+toList =
+   let go Interval = []
+       go (Node p l r) = go l ++ p : go r
+   in  go
+
+lookup :: Ord x => T x y -> x -> (Maybe (x,y), Maybe (x,y))
+lookup nodes0 x0 =
+   let go lb rb Interval = (lb, rb)
+       go lb rb (Node n@(x,_y) ln rn) =
+          if x0>=x
+            then go (Just n) rb rn
+            else go lb (Just n) ln
+   in  go Nothing Nothing nodes0
diff --git a/src/Numeric/Interpolation/Piece.hs b/src/Numeric/Interpolation/Piece.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Piece.hs
@@ -0,0 +1,22 @@
+module Numeric.Interpolation.Piece (
+   Piece.T,
+   Piece.linear,
+   hermite1,
+   ) where
+
+import qualified Numeric.Interpolation.Private.Piece as Piece
+
+
+{- |
+Hermite interpolation with one derivative per node.
+That is, the interpolating polynomial is cubic.
+-}
+hermite1 :: (Fractional a) => Piece.T a a (a, a)
+hermite1 (x0,(y0,dy0)) (x1,(y1,dy1)) x =
+   let d = (y1-y0) / dx10
+       dx0 = x-x0
+       dx1 = x1-x
+       dx10 = x1-x0
+   in  (y0*dx1 + y1*dx0 +
+        ((dy0-d) * dx1 - (dy1-d) * dx0) * dx0 * dx1 / dx10)
+          / dx10
diff --git a/src/Numeric/Interpolation/Piecewise.hs b/src/Numeric/Interpolation/Piecewise.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Piecewise.hs
@@ -0,0 +1,30 @@
+module Numeric.Interpolation.Piecewise (
+   interpolate,
+   interpolateConstantExt,
+   ) where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import qualified Numeric.Interpolation.Type as Type
+
+
+{- |
+It is a checked error to interpolate outside of the range of nodes.
+-}
+interpolate :: (Ord x) => Type.T x y ny -> Nodes.T x ny -> x -> y
+interpolate typ ns x =
+   case Nodes.lookup ns x of
+      (Just p0, Just p1) -> Type.interpolatePiece typ p0 p1 x
+      _ -> error "interpolate: argument outside range"
+
+{- |
+Outside the range of nodes the interpolation function
+takes the value of the respective border.
+-}
+interpolateConstantExt ::
+   (Ord x) => Type.T x y ny -> Nodes.T x ny -> x -> y
+interpolateConstantExt typ ns x =
+   case Nodes.lookup ns x of
+      (Just p0, Just p1) -> Type.interpolatePiece typ p0 p1 x
+      (Just p, Nothing) -> Type.valueFromNode typ $ snd p
+      (Nothing, Just p) -> Type.valueFromNode typ $ snd p
+      (Nothing, Nothing) -> error "interpolateConstantExt: empty node list"
diff --git a/src/Numeric/Interpolation/Type.hs b/src/Numeric/Interpolation/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/Interpolation/Type.hs
@@ -0,0 +1,81 @@
+module Numeric.Interpolation.Type (
+   T(..),
+   linear,
+   cubic,
+   cubicLinear,
+   cubicParabola,
+   ) where
+
+import qualified Numeric.Interpolation.NodeList as Nodes
+import qualified Numeric.Interpolation.Piece as Piece
+import qualified Numeric.Interpolation.Basis as Basis
+import Numeric.Interpolation.Private.Basis (hermite1Split)
+
+
+data T x y ny =
+   Cons {
+      ssvFromNodes :: [x] -> [y] -> String,
+      interpolatePiece :: Piece.T x y ny,
+      basisFunctions :: [x] -> [Nodes.T x ny],
+      coefficientsToInterpolator :: [x] -> [y] -> Nodes.T x ny,
+      valueFromNode :: ny -> y
+   }
+
+linear :: T Double Double Double
+linear =
+   Cons {
+      ssvFromNodes =
+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
+      interpolatePiece = Piece.linear,
+      basisFunctions = Basis.linear,
+      coefficientsToInterpolator = Basis.coefficientsToLinear,
+      valueFromNode = id
+   }
+
+cubic :: T Double Double (Double, Double)
+cubic =
+   Cons {
+      ssvFromNodes =
+         \xs ys ->
+            unlines .
+            zipWith (\x (y,dy) -> show x ++ " " ++ show y ++ " " ++ show dy) xs $
+            hermite1Split xs ys,
+      interpolatePiece = Piece.hermite1,
+      basisFunctions = Basis.hermite1,
+      coefficientsToInterpolator = Basis.coefficientsToHermite1,
+      valueFromNode = fst
+   }
+
+cubicLinear :: T Double Double (Double, Double)
+cubicLinear =
+   Cons {
+      ssvFromNodes =
+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
+      interpolatePiece = Piece.hermite1,
+      basisFunctions = Basis.cubicLinear,
+      coefficientsToInterpolator = Basis.coefficientsToCubicLinear,
+      valueFromNode = fst
+   }
+
+cubicParabola :: T Double Double (Double, Double)
+cubicParabola =
+   Cons {
+      ssvFromNodes =
+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
+      interpolatePiece = Piece.hermite1,
+      basisFunctions = Basis.cubicParabola,
+      coefficientsToInterpolator = Basis.coefficientsToCubicParabola,
+      valueFromNode = fst
+   }
+
+
+_cubicMean :: T Double Double (Double, Double)
+_cubicMean =
+   Cons {
+      ssvFromNodes =
+         \xs ys -> unlines $ zipWith (\x y -> show x ++ " " ++ show y) xs ys,
+      interpolatePiece = Piece.hermite1,
+      basisFunctions = Basis.cubicParabola, -- Basis.cubicMean,
+      coefficientsToInterpolator = Basis.coefficientsToCubicParabola, -- not correct
+      valueFromNode = fst
+   }
diff --git a/test/Test.hs b/test/Test.hs
new file mode 100644
--- /dev/null
+++ b/test/Test.hs
@@ -0,0 +1,12 @@
+module Main where
+
+import qualified Test.Piece as Piece
+
+
+run :: String -> [(String, IO ())] -> IO ()
+run prefix =
+   mapM_ (\(msg,act) -> putStr (prefix ++ '.' : msg ++ ": ") >> act)
+
+main :: IO ()
+main =
+   run "Piece" Piece.tests
diff --git a/test/Test/Piece.hs b/test/Test/Piece.hs
new file mode 100644
--- /dev/null
+++ b/test/Test/Piece.hs
@@ -0,0 +1,60 @@
+module Test.Piece where
+
+import qualified Numeric.Interpolation.Piece as Piece
+import qualified Numeric.Interpolation.Private.Piece as PiecePriv
+
+import Test.QuickCheck (Property, quickCheck, (==>), )
+
+
+type Point = (Rational, Rational)
+
+linearCommutative ::
+   Point -> Point -> Rational -> Property
+linearCommutative p1@(x1,_) p2@(x2,_) x =
+   x1/=x2
+   ==>
+   Piece.linear p1 p2 x
+   ==
+   Piece.linear p2 p1 x
+
+
+type PointSlope = (Rational, (Rational, Rational))
+
+hermite1Commutative ::
+   PointSlope -> PointSlope -> Rational -> Property
+hermite1Commutative p1@(x1,_) p2@(x2,_) x =
+   x1/=x2
+   ==>
+   Piece.hermite1 p1 p2 x
+   ==
+   Piece.hermite1 p2 p1 x
+
+
+linearHermite1 ::
+   Point -> Point -> Rational -> Property
+linearHermite1 p1@(x1,y1) p2@(x2,y2) x =
+   x1/=x2
+   ==>
+   Piece.linear p1 p2 x
+   ==
+   let slope = (y2-y1)/(x2-x1)
+   in  Piece.hermite1 (x1, (y1,slope)) (x2, (y2, slope)) x
+
+
+hermite1Alternative ::
+   PointSlope -> PointSlope -> Rational -> Property
+hermite1Alternative p1@(x1,_) p2@(x2,_) x =
+   x1/=x2
+   ==>
+   Piece.hermite1 p1 p2 x
+   ==
+   PiecePriv.hermite1 p1 p2 x
+
+
+tests :: [(String, IO ())]
+tests =
+   ("linearCommutative", quickCheck linearCommutative) :
+   ("hermite1Commutative", quickCheck hermite1Commutative) :
+   ("linearHermite1", quickCheck linearHermite1) :
+   ("hermite1Alternative", quickCheck hermite1Alternative) :
+   []