diff --git a/Geom2D.hs b/Geom2D.hs
--- a/Geom2D.hs
+++ b/Geom2D.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE BangPatterns #-}
+
 -- | Basic 2 dimensional geometry functions.
 module Geom2D where
 
@@ -6,8 +8,8 @@
 infixr 5 $*
 
 data Point = Point {
-  pointX :: Double,
-  pointY :: Double}
+  pointX :: {-# UNPACK #-} !Double,
+  pointY :: {-# UNPACK #-} !Double}
 
 instance Show Point where
   show (Point x y) =
@@ -15,12 +17,12 @@
 
 -- | A transformation (x, y) -> (ax + by + c, dx + ey + d)
 data Transform = Transform {
-  xformA :: Double,
-  xformB :: Double,
-  xformC :: Double,
-  xformD :: Double,
-  xformE :: Double,
-  xformF :: Double }
+  xformA :: {-# UNPACK #-} !Double,
+  xformB :: {-# UNPACK #-} !Double,
+  xformC :: {-# UNPACK #-} !Double,
+  xformD :: {-# UNPACK #-} !Double,
+  xformE :: {-# UNPACK #-} !Double,
+  xformF :: {-# UNPACK #-} !Double }
                deriving Show
 
 data Line = Line Point Point
@@ -49,7 +51,7 @@
 inverse :: Transform -> Maybe Transform
 inverse (Transform a b c d e f) = case a*e - b*d of
   0 -> Nothing
-  det -> Just $ Transform (a/det) (d/det) (-(a*c + d*f)/det) (b/det) (e/det)
+  det -> Just $! Transform (a/det) (d/det) (-(a*c + d*f)/det) (b/det) (e/det)
          (-(b*c + e*f)/det)
 
 -- | Return the parameters (a, b, c) for the normalised equation
@@ -71,15 +73,18 @@
 -- | The lenght of the vector.
 vectorMag :: Point -> Double
 vectorMag (Point x y) = sqrt(x*x + y*y)
+{-# INLINE vectorMag #-}
 
 -- | The angle of the vector, in the range @(-'pi', 'pi']@.
 vectorAngle :: Point -> Double
 vectorAngle (Point 0.0 0.0) = 0.0
 vectorAngle (Point x y) = atan2 y x
+{-# INLINE vectorAngle #-}
 
 -- | The unitvector with the given angle.
 dirVector :: Double -> Point
 dirVector angle = Point (cos angle) (sin angle)
+{-# INLINE dirVector #-}
 
 -- | The unit vector with the same direction.
 normVector :: Point -> Point
@@ -89,38 +94,47 @@
 -- | Scale vector by constant.
 (*^) :: Double -> Point -> Point
 s *^ (Point x y) = Point (s*x) (s*y)
+{-# INLINE (*^) #-}
 
 -- | Scale vector by reciprocal of constant.
 (^/) :: Point -> Double -> Point
 (Point x y) ^/ s = Point (x/s) (y/s)
+{-# INLINE (^/) #-}
 
 -- | Scale vector by constant, with the arguments swapped.
 (^*) :: Point -> Double -> Point
 p ^* s = s *^ p
+{-# INLINE (^*) #-}
 
 -- | Add two vectors.
 (^+^) :: Point -> Point -> Point
 (Point x1 y1) ^+^ (Point x2 y2) = Point (x1+x2) (y1+y2)
+{-# INLINE (^+^) #-}
 
 -- | Subtract two vectors.
 (^-^) :: Point -> Point -> Point
 (Point x1 y1) ^-^ (Point x2 y2) = Point (x1-x2) (y1-y2)
+{-# INLINE (^-^) #-}
 
 -- | Dot product of two vectors.
 (^.^) :: Point -> Point -> Double
 (Point x1 y1) ^.^ (Point x2 y2) = x1*x2 + y1*y2
+{-# INLINE (^.^) #-}
 
 -- | Cross product of two vectors.
 vectorCross :: Point -> Point -> Double
 vectorCross (Point x1 y1) (Point x2 y2) = x1*y2 - y1*x2
+{-# INLINE vectorCross #-}
 
 -- | Distance between two vectors.
 vectorDistance :: Point -> Point -> Double
 vectorDistance p q = vectorMag (p^-^q)
+{-# INLINE vectorDistance #-}
 
 -- | Interpolate between two vectors.
 interpolateVector :: Point -> Point -> Double -> Point
 interpolateVector a b t = t*^b ^+^ (1-t)*^a
+{-# INLINE interpolateVector #-}
 
 -- | Create a transform that rotates by the angle of the given vector
 -- with the x-axis
diff --git a/Geom2D/CubicBezier.hs b/Geom2D/CubicBezier.hs
--- a/Geom2D/CubicBezier.hs
+++ b/Geom2D/CubicBezier.hs
@@ -1,18 +1,22 @@
--- | Export all the cubic bezier functions, but not the numeric helper functions
+-- | Export all the cubic bezier functions.
 
 module Geom2D.CubicBezier
        (module Geom2D.CubicBezier.Basic,
         module Geom2D.CubicBezier.Approximate,
         module Geom2D.CubicBezier.Outline,
         module Geom2D.CubicBezier.Curvature,
-        module Geom2D.CubicBezier.Intersection
+        module Geom2D.CubicBezier.Intersection,
+        module Geom2D.CubicBezier.MetaPath,
+        module Geom2D
        ) where
 
+import Geom2D
 import Geom2D.CubicBezier.Basic
 import Geom2D.CubicBezier.Approximate
 import Geom2D.CubicBezier.Outline
 import Geom2D.CubicBezier.Curvature
 import Geom2D.CubicBezier.Intersection
+import Geom2D.CubicBezier.MetaPath
        
        
         
diff --git a/Geom2D/CubicBezier/Approximate.hs b/Geom2D/CubicBezier/Approximate.hs
--- a/Geom2D/CubicBezier/Approximate.hs
+++ b/Geom2D/CubicBezier/Approximate.hs
@@ -1,5 +1,6 @@
+{-# LANGUAGE BangPatterns #-}
 module Geom2D.CubicBezier.Approximate (
-  approximateCurve, approximateCurveWithParams)
+  approximatePath, approximatePathMax, approximateCurve, approximateCurveWithParams)
        where
 import Geom2D
 import Geom2D.CubicBezier.Numeric
@@ -7,7 +8,92 @@
 import Data.Function
 import Data.List
 import Data.Maybe
+import qualified Data.Map as M
 
+interpolate :: Double -> Double -> Double -> Double
+interpolate a b x = (1-x)*a + x*b
+
+-- | Approximate a function with piecewise cubic bezier splines using
+-- a least-squares fit, within the given tolerance.  Each subcurve is
+-- approximated by using a finite number of samples.  It is recommended
+-- to avoid changes in direction by subdividing the original function
+-- at points of inflection.
+
+approximatePath :: (Double -> (Point, Point)) -- ^ The function to approximate and it's derivative
+                -> Double                     -- ^ The number of discrete samples taken to approximate each subcurve
+                -> Double                     -- ^ The tolerance
+                -> Double                     -- ^ The lower parameter of the function      
+                -> Double                     -- ^ The upper parameter of the function
+                -> [CubicBezier]
+approximatePath f n tol tmin tmax
+  | err <= tol = [cb_out]
+  | otherwise = approximatePath f n tol tmin terr ++
+                approximatePath f n tol terr tmax
+  where
+    (cb_out, terr', err) = approximateCurveWithParams curveCb
+                           points ts tol
+    terr = interpolate tmin tmax terr'
+    ts        = [i/(n+1) | i <- [1..n]]
+    points    = map (fst . f . interpolate tmin tmax) ts
+    (t0, t0') = f tmin
+    (t1, t1') = f tmax
+    curveCb = CubicBezier t0 (t0^+^t0') (t1^-^t1') t1
+
+
+-- | Like approximatePath, but limit the number of subcurves.
+approximatePathMax :: Int                        -- ^ The maximum number of subcurves
+                   -> (Double -> (Point, Point)) -- ^ The function to approximate and it's derivative
+                   -> Double                     -- ^ The number of discrete samples taken to approximate each subcurve
+                   -> Double                     -- ^ The tolerance
+                   -> Double                     -- ^ The lower parameter of the function      
+                   -> Double                     -- ^ The upper parameter of the function
+                   -> [CubicBezier]
+approximatePathMax m f n tol tmin tmax =
+  approxMax f tol m ts segments
+  where segments              = M.singleton err (FunctionSegment tmin tmax t_err outline)
+        (p0, p0') = f tmin
+        (p1, p1') = f tmax
+        ts = [i/(n+1) | i <- [1..n]]
+        points = map (fst . f . interpolate tmin tmax) ts
+        curveCb = CubicBezier p0 (p0^+^p0') (p1^-^p1') p1
+        (outline, t_err', err) = approximateCurveWithParams curveCb
+                                 points ts tol
+        t_err = interpolate tmin tmax t_err'
+
+data FunctionSegment = FunctionSegment {
+  fs_t_min :: {-# UNPACK #-} !Double,  -- the least t param of the segment in the original curve
+  _fs_t_max :: {-# UNPACK #-} !Double,  -- the max t param of the segment in the original curve
+  _fs_t_err :: {-# UNPACK #-} !Double,  -- the param where the error is maximal
+  fs_curve :: CubicBezier -- the curve segment
+  }
+
+-- Keep a map from maxError to FunctionSegment for each subsegment to keep
+-- track of the segment with the maximum error.  This ensures a n
+-- log(n) execution time, rather than n^2 when a list is used.
+approxMax :: (Double -> (Point, Point)) -> Double -> Int
+          -> [Double] -> M.Map Double FunctionSegment -> [CubicBezier]
+approxMax f tol n ts segments
+  | n < 1 = error "Minimum number of segments is one."
+  | (n == 1) || (err < tol) =
+    map fs_curve $ sortBy (compare `on` fs_t_min) $ map snd $ M.toList segments
+  | otherwise = approxMax f tol (n-1) ts $
+                M.insert err_l (FunctionSegment t_min t_err t_err_l curve_l) $
+                M.insert err_r (FunctionSegment t_err t_max t_err_r curve_r)
+                newSegments
+  where
+    ((err, FunctionSegment t_min t_max t_err _), newSegments) = M.deleteFindMax segments
+    (fmin, fmin') = f t_min
+    (fmid, fmid') = f t_err
+    (fmax, fmax') = f t_max
+    fcurve_l = CubicBezier fmin (fmin^+^fmin') (fmid^-^fmid') fmid
+    fcurve_r = CubicBezier fmid (fmid^+^fmid') (fmax^-^fmax') fmax
+    pointsl = map (fst . f . interpolate t_min t_err) ts
+    pointsr = map (fst . f . interpolate t_err t_max) ts
+    t_err_l = interpolate t_min t_err t_err_l'
+    t_err_r = interpolate t_err t_max t_err_r'
+    (curve_l, t_err_l', err_l)  = approximateCurveWithParams fcurve_l pointsl ts tol
+    (curve_r, t_err_r', err_r)  = approximateCurveWithParams fcurve_r pointsr ts tol
+
 -- | @approximateCurve b pts eps@ finds the least squares fit of a bezier
 -- curve to the points @pts@.  The resulting bezier has the same first
 -- and last control point as the curve @b@, and have tangents colinear with @b@.
@@ -31,9 +117,16 @@
       (t, maxError) = maximumBy (compare `on` snd) (zip ts distances)
   in (c, t, maxError)
 
-add6 (a, b, c, d, e, f) (a', b', c', d', e', f') =
-  (a+a', b+b', c+c', d+d', e+e', f+f')
+data LSParams = LSParams {-# UNPACK #-} !Double
+                {-# UNPACK #-} !Double
+                {-# UNPACK #-} !Double
+                {-# UNPACK #-} !Double
+                {-# UNPACK #-} !Double
+                {-# UNPACK #-} !Double
 
+addParams :: LSParams -> LSParams -> LSParams
+addParams (LSParams a b c d e f) (LSParams a' b' c' d' e' f') =
+  LSParams (a+a') (b+b') (c+c') (d+d') (e+e') (f+f')
 
 -- find the least squares between the points p_i and B(t_i) for
 -- bezier curve B, where pts contains the points p_i and ts
@@ -46,8 +139,8 @@
 -- with two unknown values (alpha1 and alpha2), which can be
 -- solved easily
 leastSquares :: CubicBezier -> [Point] -> [Double] -> Maybe CubicBezier
-leastSquares (CubicBezier (Point p1x p1y) (Point p2x p2y) (Point p3x p3y) (Point p4x p4y)) pts ts = let
-  calcParams t (Point px py)  = let
+leastSquares (CubicBezier (Point !p1x !p1y) (Point !p2x !p2y) (Point !p3x !p3y) (Point !p4x !p4y)) pts ts = let
+  calcParams t (Point px py) = let
     t2 = t * t; t3 = t2 * t
     ax = 3 * (p2x - p1x) * (t3 - 2 * t2 + t)
     ay = 3 * (p2y - p1y) * (t3 - 2 * t2 + t)
@@ -55,13 +148,14 @@
     by = 3 * (p3y - p4y) * (t2 - t3)
     cx = (p4x - p1x) * (3 * t2 - 2 * t3) + p1x - px
     cy = (p4y - p1y) * (3 * t2 - 2 * t3) + p1y - py
-    in (ax * ax + ay * ay,
-        ax * bx + ay * by,
-        ax * cx + ay * cy,
-        bx * ax + by * ay,
-        bx * bx + by * by,
-        bx * cx + by * cy)
-  (a, b, c, d, e, f) = foldl1' add6 $ zipWith calcParams ts pts
+    in LSParams
+       (ax * ax + ay * ay)
+       (ax * bx + ay * by)
+       (ax * cx + ay * cy)
+       (bx * ax + by * ay)
+       (bx * bx + by * by)
+       (bx * cx + by * cy)
+  LSParams !a !b !c !d !e !f = foldl1' addParams (zipWith calcParams ts pts)
   in do (alpha1, alpha2) <- solveLinear2x2 a b c d e f
         let cp1 = Point (alpha1 * (p2x - p1x) + p1x) (alpha1 * (p2y - p1y) + p1y)
             cp2 = Point (alpha2 * (p3x - p4x) + p4x) (alpha2 * (p3y - p4y) + p4y)
@@ -76,8 +170,8 @@
   newCurve <- leastSquares curve pts ts
   let deltaTs = zipWith (calcDeltaT newCurve) pts ts
       ts' = map (max 0 . min 1) $ zipWith (-) ts deltaTs
-  newCurve <- leastSquares curve pts ts'
-  let deltaTs' = zipWith (calcDeltaT newCurve) pts ts'
+  newerCurve <- leastSquares curve pts ts'
+  let deltaTs' = zipWith (calcDeltaT newerCurve) pts ts'
       newTs = interpolateTs ts ts' deltaTs deltaTs'
       thisDeltaT = maximum $ map abs $ zipWith (-) newTs ts
   if maxiter < 1 ||
@@ -122,7 +216,8 @@
 -- the reduction of t is one iteration of Newton Raphson:  f'(t)/f''(t)
 -- using more iterations doesn't appear to give an improvement
 -- See Curve Fitting with Piecewise Parametric Cubics by Stone & Plass
-calcDeltaT curve (Point ptx pty) t = let
+calcDeltaT :: CubicBezier -> Point -> Double -> Double
+calcDeltaT curve (Point !ptx !pty) t = let
   [Point bezx bezy, Point dbezx dbezy, Point ddbezx ddbezy, _] = evalBezierDerivs curve t
   in ((bezx - ptx) * dbezx + (bezy - pty) * dbezy) /
      (dbezx * dbezx + dbezy * dbezy + (bezx - ptx) * ddbezx + (bezy - pty) * ddbezy)
diff --git a/Geom2D/CubicBezier/Basic.hs b/Geom2D/CubicBezier/Basic.hs
--- a/Geom2D/CubicBezier/Basic.hs
+++ b/Geom2D/CubicBezier/Basic.hs
@@ -19,7 +19,10 @@
   bezierC3 :: Point} deriving Show
 
 data PathJoin = JoinLine | JoinCurve Point Point
-data Path = Path Point [(PathJoin, Point)]
+              deriving Show
+data Path = OpenPath [(Point, PathJoin)] Point
+          | ClosedPath [(Point, PathJoin)]
+          deriving Show
 
 instance AffineTransform CubicBezier where
   transform t (CubicBezier c0 c1 c2 c3) =
@@ -39,11 +42,11 @@
 -- can use the maximum of the convex hull of the derivative, and double it to
 -- have some margin for larger values.
 bezierParamTolerance :: CubicBezier -> Double -> Double
-bezierParamTolerance (CubicBezier p1 p2 p3 p4) eps = eps / maxDist
+bezierParamTolerance (CubicBezier !p1 !p2 !p3 !p4) eps = eps / maxDist
   where 
-    maxDist = 6 * maximum [vectorDistance p1 p2,
-                           vectorDistance p2 p3,
-                           vectorDistance p3 p4]
+    maxDist = 6 * (max (vectorDistance p1 p2) $
+                   max (vectorDistance p2 p3)
+                   (vectorDistance p3 p4))
 
 -- | Reorient to the curve B(1-t).
 reorient :: CubicBezier -> CubicBezier
@@ -57,23 +60,28 @@
 
 -- | Calculate a value on the curve.
 evalBezier :: CubicBezier -> Double -> Point
-evalBezier b t = Point (bernsteinEval x t) (bernsteinEval y t)
-  where (x, y) = bezierToBernstein b
+evalBezier b = fst . evalBezierDeriv b 
 
 -- | Calculate a value and the first derivative on the curve.
 evalBezierDeriv :: CubicBezier -> Double -> (Point, Point)
-evalBezierDeriv b =
-  let (px, py) = bezierToBernstein b
-      px' = bernsteinDeriv px
-      py' = bernsteinDeriv py
-  in \t -> (Point (bernsteinEval px t) (bernsteinEval py t),
-            Point (bernsteinEval px' t) (bernsteinEval py' t))
-
+evalBezierDeriv (CubicBezier !p0 !p1 !p2 !p3) t = (bt, bt')
+  where
+    b0' = 3*^(p1^-^p0)
+    b0'' = 2*^(3*^(p2^-^p1) ^-^ b0')
+    b0''' = 6*^(p3^-^ 2*^p2 ^+^ p1) ^-^ b0''
+    bt' = b0'^+^(b0''^+^ t*^b0'''^/2)^*t
+    bt = p0 ^+^ t*^(b0' ^+^ t*^(b0''^/2 ^+^ t*^(b0'''^/6)))
+    
 -- | Calculate a value and all derivatives on the curve.
 evalBezierDerivs :: CubicBezier -> Double -> [Point]
-evalBezierDerivs b t = zipWith Point (bernsteinEvalDerivs px t)
-                       (bernsteinEvalDerivs py t)
-  where (px, py) = bezierToBernstein b
+evalBezierDerivs (CubicBezier !p0 !p1 !p2 !p3) t = [bt, bt', bt'', b0''']
+  where
+    b0' = 3*^(p1^-^p0)
+    b0'' = 2*^(3*^(p2^-^p1) ^-^ b0')
+    b0''' = 6*^(p3^-^ 2*^p2 ^+^ p1) ^-^ b0''
+    bt'' = b0''^+^ b0'''^*t
+    bt' = b0'^+^(b0''^+^ t*^b0'''^/2)^*t
+    bt = p0 ^+^ t*^(b0' ^+^ t*^(b0''^/2 ^+^ t*^(b0'''^/6)))
 
 -- | @findBezierTangent p b@ finds the parameters where
 -- the tangent of the bezier curve @b@ has the same direction as vector p.
@@ -140,10 +148,11 @@
   where distDeriv t' = vectorMag $ snd $ evalD t'
         evalD = evalBezierDeriv b 
 
+outline :: CubicBezier -> Double
 outline (CubicBezier c0 c1 c2 c3) =
-  sum [vectorDistance c0 c1,
-       vectorDistance c1 c2,
-       vectorDistance c2 c3]
+  vectorDistance c0 c1 +
+  vectorDistance c1 c2 +
+  vectorDistance c2 c3
 
 arcLengthEstimate :: CubicBezier -> Double -> (Double, (Double, Double))
 arcLengthEstimate b eps = (arclen, (estimate, ol))
@@ -153,17 +162,20 @@
     ol = outline b
     (arcL, (estL, olL)) = arcLengthEstimate bl eps
     (arcR, (estR, olR)) = arcLengthEstimate br eps
-    arclen | (abs(estL + estR - estimate) < eps) = estL + estR
+    arclen | abs(estL + estR - estimate) < eps = estL + estR
            | otherwise = arcL + arcR
 
 -- | arcLengthParam c len tol finds the parameter where the curve c has the arclength len,
 -- within tolerance tol.
+arcLengthParam :: CubicBezier -> Double -> Double -> Double
 arcLengthParam b len eps =
   arcLengthP b len ol (len/ol) 1 eps
   where ol = outline b
 
 -- Use the Newton rootfinding method.  Start with large tolerance
 -- values, and decrease tolerance as we go closer to the root.
+arcLengthP :: CubicBezier -> Double -> Double ->
+              Double -> Double -> Double -> Double
 arcLengthP !b !len !tot !t !dt !eps
   | abs diff < eps = t - newDt
   | otherwise = arcLengthP b len tot (t - newDt) newDt eps
@@ -202,7 +214,7 @@
 
 -- | Return True if all the control points are colinear within tolerance.
 colinear :: CubicBezier -> Double -> Bool
-colinear (CubicBezier a b c d) eps =
+colinear (CubicBezier !a !b !c !d) eps =
   abs (ld b) < eps && abs (ld c) < eps
   where ld = lineDistance (Line a d)
 
diff --git a/Geom2D/CubicBezier/Curvature.hs b/Geom2D/CubicBezier/Curvature.hs
--- a/Geom2D/CubicBezier/Curvature.hs
+++ b/Geom2D/CubicBezier/Curvature.hs
@@ -6,7 +6,8 @@
 import Geom2D.CubicBezier.Intersection
 import Math.BernsteinPoly
 
--- | Curvature of the Bezier curve.
+-- | Curvature of the Bezier curve.  A negative curvature means the curve
+-- curves to the right.
 curvature :: CubicBezier -> Double -> Double
 curvature b t
   | t == 0 = curvature' b
@@ -14,7 +15,8 @@
   | t < 0.5 = curvature' $ snd $ splitBezier b t
   | otherwise = negate $ curvature' $ reorient $ fst $ splitBezier b t
 
-curvature' (CubicBezier c0 c1 c2 c3) = 2/3 * b/a^3
+curvature' :: CubicBezier -> Double
+curvature' (CubicBezier c0 c1 c2 _c3) = 2/3 * b/a^(3::Int)
   where 
     a = vectorDistance c1 c0
     b = (c1^-^c0) `vectorCross` (c2^-^c1)
@@ -25,29 +27,29 @@
 radiusOfCurvature b t = 1 / curvature b t
 
 extrema :: CubicBezier -> BernsteinPoly
-extrema (CubicBezier p0 p1 p2 p3) =
-  let bez = [p0, p1, p2, p3]
-      x' = bernsteinDeriv $ listToBernstein $ map pointX bez
-      y' = bernsteinDeriv $ listToBernstein $ map pointY bez
+extrema bez = 
+  let (x, y) = bezierToBernstein bez
+      x' = bernsteinDeriv x
+      y' = bernsteinDeriv y
       x'' = bernsteinDeriv x'
       y'' = bernsteinDeriv y'
       x''' = bernsteinDeriv x''
       y''' = bernsteinDeriv y''
-  in -- (y'^2 + x'^2) * (x'*y''' - y'*x''') -
-     -- 3 * (x'*y'' - y'*x'') * (y'*y'' + x'*x'')
-   (y'~*y' ~+ x'~*x') ~* (x'~*y''' ~- y'~*x''') ~-
-   3 *~ (x'~*y'' ~- y'~*x'') ~* (y'~*y'' ~+ x'~*x'')
+  in (y'~*y' ~+ x'~*x') ~* (x'~*y''' ~- y'~*x''') ~-
+     3 *~ (x'~*y'' ~- y'~*x'') ~* (y'~*y'' ~+ x'~*x'')
 
 -- | Find extrema of the curvature, but not inflection points.
 curvatureExtrema :: CubicBezier -> Double -> [Double]
-curvatureExtrema b eps = bezierFindRoot (extrema b) 0 1 $
-                         bezierParamTolerance b eps
+curvatureExtrema b eps
+  | colinear b eps = []
+  | otherwise = bezierFindRoot (extrema b) 0 1 $
+                bezierParamTolerance b eps
 
 radiusSquareEq :: CubicBezier -> Double -> BernsteinPoly
-radiusSquareEq (CubicBezier p0 p1 p2 p3) d =
-  let bez = [p0, p1, p2, p3]
-      x' = bernsteinDeriv $ listToBernstein $ map pointX bez
-      y' = bernsteinDeriv $ listToBernstein $ map pointY bez
+radiusSquareEq bez d =
+  let (x, y) = bezierToBernstein bez
+      x' = bernsteinDeriv x
+      y' = bernsteinDeriv y
       x'' = bernsteinDeriv x'
       y'' = bernsteinDeriv y'
       a =  x'~*x' ~+  y'~*y'
@@ -55,9 +57,12 @@
   in (a~*a~*a) ~- (d*d) *~ b~*b
 
 -- | Find points on the curve that have a certain radius of curvature.
+-- Values to the left and to the right of the curve are returned.
 findRadius :: CubicBezier  -- ^ the curve
            -> Double       -- ^ distance
            -> Double       -- ^ tolerance
            -> [Double]     -- ^ result
-findRadius b d eps = bezierFindRoot (radiusSquareEq b d) 0 1 $
-                     bezierParamTolerance b eps
+findRadius b d eps
+  | colinear b eps = []  -- either empty or a huge list of t's
+  | otherwise = bezierFindRoot (radiusSquareEq b d) 0 1 $
+                bezierParamTolerance b eps
diff --git a/Geom2D/CubicBezier/Intersection.hs b/Geom2D/CubicBezier/Intersection.hs
--- a/Geom2D/CubicBezier/Intersection.hs
+++ b/Geom2D/CubicBezier/Intersection.hs
@@ -1,7 +1,7 @@
 {-# LANGUAGE BangPatterns #-}
 -- | Intersection routines using Bezier Clipping.  Provides also functions for finding the roots of onedimensional bezier curves.  This can be used as a general polynomial root solver by converting from the power basis to the bernstein basis.
 module Geom2D.CubicBezier.Intersection
-       (bezierIntersection, bezierLineIntersections, bezierFindRoot)
+       (bezierIntersection, bezierLineIntersections, bezierFindRoot, closest)
        where
 import Geom2D
 import Geom2D.CubicBezier.Basic
@@ -11,19 +11,21 @@
 
 -- find the convex hull by comparing the angles of the vectors with
 -- the cross product and backtracking if necessary.
-findOuter' upper !dir !p1 l@(p2:rest)
+findOuter' :: Bool -> Point -> Point -> [Point] -> Either [Point] [Point]
+findOuter' !upper !dir !p1 l@(p2:rest)
   -- backtrack if the direction is outward
   | if upper
     then dir `vectorCross` (p2^-^p1) > 0 -- left turn
-    else dir `vectorCross` (p2^-^p1) < 0 = Left l
+    else dir `vectorCross` (p2^-^p1) < 0 = Left $! l
   -- succeed
   | otherwise = case findOuter' upper (p2^-^p1) p2 rest of
     Left m -> findOuter' upper dir p1 m
     Right m -> Right (p1:m)
 
-findOuter' _ _ p1 p = Right (p1:p)
+findOuter' _ _ p1 p = Right $! (p1:p)
 
 -- find the outermost point.  It doesn't look at the x values.
+findOuter :: Bool -> [Point] -> [Point]
 findOuter upper (p1:p2:rest) =
   case findOuter' upper (p2^-^p1) p2 rest of
     Right l -> p1:l
@@ -40,38 +42,40 @@
       findOuter False points)
 
 -- test if the chords cross the fat line
--- use continuation passing style
+-- return the continuation if above the line
 testBelow :: Double -> [Point] -> Maybe Double -> Maybe Double
-testBelow dmin [] _ = Nothing
-testBelow dmin [_] _ = Nothing
-testBelow dmin (p:q:rest) cont
+testBelow _    [] _ = Nothing
+testBelow _    [_] _ = Nothing
+testBelow !dmin (p:q:rest) cont
   | pointY p >= dmin = cont
   | pointY p > pointY q = Nothing
   | pointY q < dmin = testBelow dmin (q:rest) cont
-  | otherwise = Just $ intersectPt dmin p q
+  | otherwise = Just $! intersectPt dmin p q
 
 testBetween :: Double -> Point -> Maybe Double -> Maybe Double
-testBetween dmax (Point x y) cont
+testBetween !dmax (Point !x !y) cont
   | y <= dmax = Just x
   | otherwise = cont
 
 -- test if the chords cross the line y=dmax somewhere
 testAbove :: Double -> [Point] -> Maybe Double
-testAbove dmax [] = Nothing
-testAbove dmax [_] = Nothing
+testAbove _    [] = Nothing
+testAbove _    [_] = Nothing
 testAbove dmax (p:q:rest)
   | pointY p < pointY q = Nothing
   | pointY q > dmax = testAbove dmax (q:rest)
-  | otherwise = Just $ intersectPt dmax p q
+  | otherwise = Just $! intersectPt dmax p q
 
 -- find the x value where the line through the two points
 -- intersect the line y=d
+intersectPt :: Double -> Point -> Point -> Double
 intersectPt d (Point x1 y1) (Point x2 y2) =
   x1 + (d  - y1) * (x2 - x1) / (y2 - y1)
 
 -- make a hull and test over which interval the
 -- curve is garuanteed to lie inside the fat line
-chopHull dmin dmax ds = do
+chopHull :: Double -> Double -> [Double] -> Maybe (Double, Double)
+chopHull !dmin !dmax ds = do
   let (upper, lower) = makeHull ds
   left_t <- testBelow dmin upper $
             testBetween dmax (head upper) $
@@ -81,8 +85,11 @@
              testAbove dmax (reverse lower)
   Just (left_t, right_t)
 
+bezierClip :: CubicBezier -> CubicBezier -> Double -> Double
+           -> Double -> Double -> Double -> Double -> Bool
+           -> [(Double, Double)]
 bezierClip p@(CubicBezier !p0 !p1 !p2 !p3) q@(CubicBezier !q0 !q1 !q2 !q3)
-  tmin tmax umin umax prevClip eps reverse
+  tmin tmax umin umax prevClip eps revCurves
 
   -- no intersection
   | isNothing chop_interval = []
@@ -94,17 +101,17 @@
     then let
       (pl, pr) = splitBezier newP 0.5
       half_t = new_tmin + (new_tmax - new_tmin) / 2
-      in bezierClip q pl umin umax new_tmin half_t newClip eps (not reverse) ++
-         bezierClip q pr umin umax half_t new_tmax newClip eps (not reverse)
+      in bezierClip q pl umin umax new_tmin half_t newClip eps (not revCurves) ++
+         bezierClip q pr umin umax half_t new_tmax newClip eps (not revCurves)
     else let
       (ql, qr) = splitBezier q 0.5
       half_t = umin + (umax - umin) / 2
-      in bezierClip ql newP umin half_t new_tmin new_tmax newClip eps (not reverse) ++
-         bezierClip qr newP half_t umax new_tmin new_tmax newClip eps (not reverse)
+      in bezierClip ql newP umin half_t new_tmin new_tmax newClip eps (not revCurves) ++
+         bezierClip qr newP half_t umax new_tmin new_tmax newClip eps (not revCurves)
 
   -- within tolerance      
   | max (umax - umin) (new_tmax - new_tmin) < eps =
-    if reverse
+    if revCurves
     then [ (umin + (umax-umin)/2,
             new_tmin + (new_tmax-new_tmin)/2) ]
     else [ (new_tmin + (new_tmax-new_tmin)/2,
@@ -112,7 +119,7 @@
 
   -- iterate with the curves reversed.
   | otherwise =
-      bezierClip q newP umin umax new_tmin new_tmax newClip eps (not reverse)
+      bezierClip q newP umin umax new_tmin new_tmax newClip eps (not revCurves)
 
   where
     d = lineDistance (Line q0 q3)
@@ -188,3 +195,14 @@
   bezierParamTolerance b eps
   where (CubicBezier p0 p1 p2 p3) = 
           fromJust (inverse $ translate p $* rotateVec (q ^-^ p)) $* b
+
+-- | Find the closest value(s) on the bezier to the given point, within tolerance.
+closest :: CubicBezier -> Point -> Double -> [Double]
+closest cb (Point px py) eps = bezierFindRoot poly 0 1 eps
+  where
+    (bx, by) = bezierToBernstein cb
+    bx' = bernsteinDeriv bx
+    by' = bernsteinDeriv by
+    poly = (bx ~- listToBernstein [px, px, px, px]) ~* bx' ~+
+           (by ~- listToBernstein [py, py, py, py]) ~* by'
+
diff --git a/Geom2D/CubicBezier/MetaPath.hs b/Geom2D/CubicBezier/MetaPath.hs
new file mode 100644
--- /dev/null
+++ b/Geom2D/CubicBezier/MetaPath.hs
@@ -0,0 +1,430 @@
+{-# LANGUAGE BangPatterns #-}
+-- | This module implements an extension to paths as used in
+-- D.E.Knuth's /Metafont/.  Metafont gives a more intuitive method to
+-- specify paths than bezier curves.  I'll give a brief overview of
+-- the metafont curves.  For a more in depth explanation look at
+-- /The MetafontBook/.
+-- 
+-- Each spline has a tension parameter, which is a relative measure of
+-- the length of the curve.  You can specify the tension for the left
+-- side and the right side of the spline separately.  By default
+-- metafont gives a tension of 1, which gives a good looking curve.
+-- Tensions shouldn't be less than 3/4, but this implementation
+-- doesn't check for it.  If you want to avoid points of inflection on
+-- the spline, you can use @TensionAtLeast@ instead of @Tension@,
+-- which will adjust the length of the control points so they fall
+-- into the /bounding triangle/, if such a triangle exist.
+--
+-- You can either give directions for each node, or let metafont find
+-- them.  Metafont will solve a set of equations to find the
+-- directions.  You can also let metafont find directions at corner
+-- points by setting the /curl/, which is how much the point /curls/
+-- at that point.  At endpoints a curl of 1 is implied when it is not
+-- given.
+--
+-- Metafont will then find the control points from the path for you.
+-- You can also specify the control points explicitly.
+--
+-- Here is an example path from the metafont program text:
+-- 
+-- @
+-- z0..z1..tension atleast 1..{curl 2}z2..z3{-1,-2}..tension 3 and 4..z4..controls z45 and z54..z5
+-- @
+-- 
+-- This path is equivalent to:
+--
+-- @
+-- z0{curl 1}..tension atleast 1 and atleast 1..{curl 2}z2{curl 2}..tension 1 and 1..
+-- {-1,-2}z3{-1,-2}..tension 3 and 4..z4..controls z45 and z54..z5
+-- @
+--
+-- This path can be used with the following datatype:
+-- 
+-- @
+-- OpenMetaPath [ (z0, MetaJoin Open (Tension 1) (Tension 1) Open)
+--              , (z1, MetaJoin Open (TensionAtLeast 1) (TensionAtLeast 1) (Curl 2))
+--              , (z2, MetaJoin Open (Tension 1) (Tension 1) Open)
+--              , (z3, MetaJoin (Direction (Point (-1) (-2))) (Tension 3) (Tension 4) Open)
+--              , (z4, Controls z45 z54)
+--              ] z5
+-- @
+--
+-- Cyclic paths are similar, but use the @CyclicMetaPath@ contructor.
+-- There is no ending point, since the ending point will be the same
+-- as the first point.
+
+module Geom2D.CubicBezier.MetaPath
+       (unmeta, MetaPath (..), MetaJoin (..), MetaNodeType (..), Tension (..))
+where
+import Geom2D
+import Geom2D.CubicBezier.Basic
+import Data.List
+import Text.Printf
+
+data MetaPath = OpenMetaPath [(Point, MetaJoin)] Point
+              | CyclicMetaPath [(Point, MetaJoin)]
+
+data MetaJoin = MetaJoin { metaTypeL :: MetaNodeType
+                         , tensionL :: Tension
+                         , tensionR :: Tension
+                         , metaTypeR :: MetaNodeType
+                         }
+              | Controls Point Point
+              deriving Show
+
+data MetaNodeType = Open
+                  | Curl {curlgamma :: Double}
+                  | Direction {nodedir :: Point}
+                  deriving Show
+
+data Tension = Tension {tensionValue :: Double}
+             | TensionAtLeast {tensionValue :: Double}
+             deriving (Eq, Show)
+
+instance Show MetaPath where
+  show (CyclicMetaPath nodes) =
+    showPath nodes ++ "cycle"
+  show (OpenMetaPath nodes lastpoint) =
+    showPath nodes ++ showPoint lastpoint
+
+showPath :: [(Point, MetaJoin)] -> [Char]
+showPath = concatMap showNodes
+  where
+    showNodes (p, Controls u v) =
+      showPoint p ++ "..controls " ++ showPoint u ++ "and " ++ showPoint v ++ ".."
+    showNodes (p, MetaJoin m1 t1 t2 m2) =
+      showPoint p ++ typename m1 ++ ".." ++ tensions ++ typename m2
+      where
+        tensions
+          | t1 == t2 && t1 == Tension 1 = ""
+          | t1 == t2 = printf "tension %s.." (showTension t1)
+          | otherwise = printf "tension %s and %s.."
+                        (showTension t1) (showTension t2)
+    showTension (TensionAtLeast t) = printf "atleast %.3f" t :: String
+    showTension (Tension t) = printf "%.3f" t :: String
+    typename Open = ""
+    typename (Curl g) = printf "{curl %.3f}" g :: String
+    typename (Direction dir) = printf "{%s}" (showPoint dir) :: String
+    
+showPoint :: Point -> String
+showPoint (Point x y) = printf "(%.3f, %.3f)" x y
+
+-- | Create a normal path from a metapath.
+unmeta :: MetaPath -> Path
+unmeta (OpenMetaPath nodes endpoint) =
+  unmetaOpen (sanitizeOpen nodes) endpoint
+
+unmeta (CyclicMetaPath nodes) =
+  case span (bothOpen . snd) nodes of
+    (l, []) -> unmetaCyclic l
+    (l, (m:n)) ->
+      if leftOpen $ snd m
+      then unmetaAsOpen (l++[m]) n
+      else unmetaAsOpen l (m:n)
+
+unmetaOpen :: [(Point, MetaJoin)] -> Point -> Path
+unmetaOpen nodes endpoint =
+  let subsegs = openSubSegments nodes endpoint
+      path = joinSegments $ map unmetaSubSegment subsegs
+  in OpenPath path endpoint
+
+-- decompose into a list of subsegments that need to be solved.
+openSubSegments :: [(Point, MetaJoin)] -> Point -> [MetaPath]
+openSubSegments l p = openSubSegments' (tails l) p
+
+openSubSegments' :: [[(Point, MetaJoin)]] -> Point -> [MetaPath]
+openSubSegments' [[]] _ = []
+openSubSegments' [] _   = []
+openSubSegments' l lastPoint = case break breakPoint l of
+  (m, n:o) ->
+    let point = case o of
+          (((p,_):_):_) -> p
+          _ -> lastPoint
+    in OpenMetaPath (map head (m ++ [n])) point :
+       openSubSegments' o lastPoint
+  _ -> error "openSubSegments': unexpected end of segments"
+
+-- join subsegments into one segment
+joinSegments :: [Path] -> [(Point, PathJoin)]
+joinSegments = concatMap nodes
+  where nodes (OpenPath n _) = n
+        nodes (ClosedPath n) = n
+
+-- solve a cyclic metapath where all angles depend on the each other.
+unmetaCyclic :: [(Point, MetaJoin)] -> Path
+unmetaCyclic nodes =
+  let points = map fst nodes
+      chords = zipWith (^-^) points (last points : points)
+      tensionsA = (map (tensionL . snd) nodes)
+      tensionsB = (map (tensionR . snd) nodes)
+      turnAngles = zipWith turnAngle chords (tail $ cycle chords)
+      thetas = solveCyclicTriD $
+               eqsCycle tensionsA
+               points
+               tensionsB
+               turnAngles
+      phis = zipWith (\x y -> -(x+y)) turnAngles (tail thetas ++ [head thetas])
+  in ClosedPath $ zip points $
+     zipWith6 unmetaJoin points (tail points ++ [head points])
+     thetas phis tensionsA tensionsB
+
+-- solve a cyclic metapath as an open path if possible.
+-- rotate to the defined node, and rotate back after
+-- solving the path.
+unmetaAsOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)] -> Path
+unmetaAsOpen l m = ClosedPath (l'++m') 
+  where n = length m
+        OpenPath o _ = unmetaOpen (sanitizeCycle (m++l)) (fst $ head m)
+        (m',l') = splitAt n o
+
+-- solve a subsegment
+unmetaSubSegment :: MetaPath -> Path
+
+-- the simple case where the control points are already given.
+unmetaSubSegment (OpenMetaPath [(p, Controls u v)] q) =
+  OpenPath [(p, JoinCurve u v)] q
+
+-- otherwise solve the angles, and find the control points
+unmetaSubSegment (OpenMetaPath nodes lastpoint) =
+  let points = map fst nodes ++ [lastpoint]
+      joins = map snd nodes
+      chords = zipWith (^-^) (tail points) points
+      tensionsA = map tensionL  joins
+      tensionsB = map tensionR joins
+      turnAngles = zipWith turnAngle chords (tail chords) ++ [0]
+      thetas = solveTriDiagonal $
+               eqsOpen points joins chords turnAngles
+               (map tensionValue tensionsA)
+               (map tensionValue tensionsB)
+      phis = zipWith (\x y -> -x-y) turnAngles (tail thetas)
+      pathjoins = zipWith6 unmetaJoin points (tail points) thetas phis tensionsA tensionsB
+  in OpenPath (zip points pathjoins) lastpoint
+
+unmetaSubSegment _ = error "unmetaSubSegment: subsegment should not be cyclic"
+
+bothOpen :: MetaJoin -> Bool
+bothOpen (MetaJoin Open _ _ Open) = True
+bothOpen _ = False
+
+leftOpen :: MetaJoin -> Bool
+leftOpen (MetaJoin Open _ _ _) = True
+leftOpen _ = False
+
+replaceLast :: [a] -> a -> [a]
+replaceLast [] _ = []
+replaceLast [_] n = [n]
+replaceLast (l:ls) n = l:replaceLast ls n
+
+sanitizeCycle :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
+sanitizeCycle l = replaceLast ls l'
+  where
+    (l':ls) = sanitizeRest (last l: l)
+
+-- replace open nodetypes with more defined nodetypes if possible
+sanitizeOpen :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
+sanitizeOpen [] = []
+
+-- starting open => curl
+sanitizeOpen ((p, MetaJoin Open t1 t2 m):rest) =
+  sanitizeRest ((p, MetaJoin (Curl 1) t1 t2 m):rest)
+sanitizeOpen l = sanitizeRest l
+   
+sanitizeRest :: [(Point, MetaJoin)] -> [(Point, MetaJoin)]
+sanitizeRest [] = []
+
+-- ending open => curl
+sanitizeRest [(p, MetaJoin m t1 t2 Open)] =
+  [(p, MetaJoin m t1 t2 (Curl 1))]
+
+sanitizeRest (node1@(p, MetaJoin m1 tl tr m2): node2@(q, MetaJoin n1 sl sr n2): rest) =
+  case (m2, n1) of
+    (Curl g, Open) -> -- curl, open => curl, curl
+      node1 : sanitizeRest ((q, MetaJoin (Curl g) sl sr n2):rest)
+    (Open, Curl g) -> -- open, curl => curl, curl
+      (p, MetaJoin m1 tl tr (Curl g)) : sanitizeRest (node2:rest)
+    (Direction dir, Open) ->   -- given, open => given, given
+      node1 : sanitizeRest ((q, (MetaJoin (Direction dir) sl sr n2)) : rest)
+    (Open, Direction dir) ->   -- open, given => given, given
+      (p, MetaJoin m1 tl tr (Direction dir)) : sanitizeRest (node2:rest)
+    _ -> node1 : sanitizeRest (node2:rest)
+
+sanitizeRest ((p, m): (q, n): rest) =
+  case (m, n) of
+    (Controls _u v, MetaJoin Open t1 t2 mt2) ->  -- explicit, open => explicit, given
+      (p, m) : sanitizeRest ((q, MetaJoin (Direction (q^-^v)) t1 t2 mt2): rest)
+    (MetaJoin mt1 tl tr Open, Controls u _v) ->  -- open, explicit => given, explicit
+      (p, MetaJoin mt1 tl tr (Direction (u^-^p))) : sanitizeRest ((q, n): rest)
+    _ -> (p, m) : sanitizeRest ((q, n) : rest)
+
+sanitizeRest (n:l) = n:sanitizeRest l
+
+-- break the subsegment if the angle to the left or the right is defined or a curl.
+breakPoint :: [(Point, MetaJoin)] -> Bool
+breakPoint ((_,  MetaJoin _ _ _ Open):(_, MetaJoin Open _ _ _):_) = False
+breakPoint _ = True
+
+-- solve the tridiagonal system for t[i]:
+-- a[n] t[i-1] + b[i] t[i] + c[b] t[i+1] = d[i]
+-- where a[0] = c[n] = 0
+-- by first rewriting it into
+-- the system t[i] + u[i] t[i+1] = v[i]
+-- where u[n] = 0
+-- then solving for t[n]
+-- see metafont the program: ¶ 283
+solveTriDiagonal :: [(Double, Double, Double, Double)] -> [Double]
+solveTriDiagonal [] = error "solveTriDiagonal: not enough equations"
+solveTriDiagonal ((_, b0, c0, d0): rows) = solutions
+  where
+    ((_, vn): twovars) =
+      reverse $ scanl nextrow (c0/b0, d0/b0) rows
+    nextrow (u, v) (ai, bi, ci, di) =
+      (ci/(bi - u*ai), (di - v*ai)/(bi - u*ai))
+    solutions = reverse $ scanl nextsol vn twovars
+    nextsol ti (u, v) = v - u*ti
+
+-- test = ((80.0,58.0,51.0),[(-432.0,78.0,102.0,503.0),(71.0,-82.0,20.0,2130.0),(52.39,-10.43,4.0,56.0),(34.0,38.0,0.0,257.0)])
+
+-- solve the cyclic tridiagonal system.
+-- see metafont the program: ¶ 286
+solveCyclicTriD :: [(Double, Double, Double, Double)] -> [Double]
+solveCyclicTriD rows = solutions
+  where
+    (!un, !vn, !wn): threevars =
+      reverse $ tail $ scanl nextrow (0, 0, 1) rows
+    nextrow (!u, !v, !w) (!ai, !bi, !ci, !di) =
+      (ci/(bi - ai*u), (di - ai*v)/(bi - ai*u), -ai*w/(bi - ai*u))
+    (totvn, totwn) = foldl (\(v', w') (u, v, w) ->
+                             (v - u*v', w - u*w'))
+                     (0, 1) threevars
+    t0 = (vn - un*totvn) / (1 - (wn - un*totwn))
+    solutions = scanl nextsol t0
+                ((un, vn, wn) : reverse (tail threevars))
+    nextsol t (!u, !v, !w) = (v + w*t0 - t)/u
+
+turnAngle :: Point -> Point -> Double
+turnAngle (Point x y) q = vectorAngle $ rotateVec p $* q
+  where p = Point x (-y)
+
+zipPrev :: [a] -> [(a, a)]
+zipPrev [] = []
+zipPrev l = zip (last l : l) l
+
+-- find the equations for a cycle containing only open points
+eqsCycle :: [Tension] -> [Point] -> [Tension]
+         -> [Double] -> [(Double, Double, Double, Double)]
+eqsCycle tensionsA points tensionsB turnAngles = 
+  zipWith4 eqTension
+  (zipPrev (map tensionValue tensionsA))
+  (zipPrev dists)
+  (zipPrev turnAngles)
+  (zipPrev (map tensionValue tensionsB))
+  where 
+    dists = zipWith vectorDistance points (tail $ cycle points)
+
+-- find the equations for an path with open points.
+-- The first and last node should be a curl or a given angle
+
+eqsOpen :: [Point] -> [MetaJoin] -> [Point] -> [Double]
+        -> [Double] -> [Double] -> [(Double, Double, Double, Double)]
+eqsOpen _ [join] [delta] _ _ _ =
+  case join of
+    MetaJoin (Curl _) _ _ (Curl _) ->
+      [(0, 1, 0, 0), (0, 1, 0, 0)]
+    MetaJoin (Curl g) t1 t2 (Direction dir) ->
+      [eqCurl0 g (tensionValue t1) (tensionValue t2) 0,
+       (0, 1, 0, turnAngle delta dir)]
+    MetaJoin (Direction dir) t1 t2 (Curl g) ->
+      [(0, 1, 0, turnAngle delta dir),
+       eqCurlN g (tensionValue t1) (tensionValue t2)]
+    MetaJoin (Direction dir) _ _ (Direction dir2) ->
+      [(0, 1, 0, turnAngle delta dir),
+       (0, 1, 0, turnAngle delta dir2)]
+    _ -> error "eqsOpen: illegal nodetype in subsegment"
+
+eqsOpen points joins chords turnAngles tensionsA tensionsB =
+  eq0 : restEquations joins tensionsA dists turnAngles tensionsB
+  where
+    dists = zipWith vectorDistance points (tail points)      
+    eq0 = case head joins of
+      (MetaJoin (Curl g) _ _ _) -> eqCurl0 g (head tensionsA) (head tensionsB) (head turnAngles)
+      (MetaJoin (Direction dir) _ _ _) -> (0, 1, 0, turnAngle (head chords) dir)
+      _ -> error "eqsOpen: illegal subsegment first nodetype"
+
+    restEquations [lastnode] (tensionA:_) _ _ (tensionB:_) =
+      case lastnode of
+        MetaJoin _ _ _ (Curl g) -> [eqCurlN g tensionA tensionB]
+        MetaJoin _ _ _ (Direction dir) -> [(0, 1, 0, turnAngle (last chords) dir)]
+        _  -> error "eqsOpen: illegal subsegment last nodetype"
+
+    restEquations (_:othernodes) (tensionA:restTA) (d:restD) (turn:restTurn) (tensionB:restTB) =
+      eqTension (tensionA, head restTA) (d, head restD) (turn, head restTurn) (tensionB, head restTB) :
+      restEquations othernodes restTA restD restTurn restTB
+
+    restEquations _ _ _ _ _ = error "eqsOpen: illegal rest equations"
+
+-- the equation for an open node
+eqTension :: (Double, Double) -> (Double, Double)
+          -> (Double, Double) -> (Double, Double)
+          -> (Double, Double, Double, Double)
+eqTension (tensionA', tensionA) (dist', dist) (psi', psi) (tensionB', tensionB) =
+  (a, b+c, d, -b*psi' - d*psi)
+  where
+    a = (tensionB' * tensionB' / (tensionA' * dist'))
+    b = (3 - 1/tensionA') * tensionB' * tensionB' / dist'
+    c = (3 - 1/tensionB) * tensionA * tensionA / dist
+    d = tensionA * tensionA / (tensionB * dist)
+
+-- the equation for a starting curl
+eqCurl0 :: Double -> Double -> Double -> Double -> (Double, Double, Double, Double)
+eqCurl0 gamma tensionA tensionB psi = (0, c, d, r)
+  where
+    c = chi/tensionA + 3 - 1/tensionB
+    d = (3 - 1/tensionA)*chi + 1/tensionB
+    chi = gamma*tensionB*tensionB / (tensionA*tensionA)
+    r = -d*psi
+
+-- the equation for an ending curl
+eqCurlN :: Double -> Double -> Double -> (Double, Double, Double, Double)
+eqCurlN gamma tensionA tensionB = (a, b, 0, 0)
+  where
+    a = (3 - 1/tensionB)*chi + 1/tensionA
+    b = chi/tensionB + 3 - 1/tensionA
+    chi = gamma*tensionA*tensionA / (tensionB*tensionB)
+
+-- magic formula for getting the control points by John Hobby
+unmetaJoin :: Point -> Point -> Double -> Double -> Tension -> Tension -> PathJoin
+unmetaJoin !z0 !z1 !theta !phi !alpha !beta
+  | abs phi < 1e-4 && abs theta < 1e-4 = JoinLine
+  | otherwise = JoinCurve u v
+  where Point dx dy = z1^-^z0
+        bounded = (sf <= 0 && st <= 0 && sf <= 0) ||
+                  (sf >= 0 && st >= 0 && sf >= 0)
+        rr' = velocity st sf ct cf alpha
+        ss' = velocity sf st cf ct beta
+        stf = st*cf + sf*ct -- sin (theta + phi)
+        st = sin theta
+        sf = sin phi
+        ct = cos theta
+        cf = cos phi
+        rr = case alpha of
+          TensionAtLeast _ | bounded ->
+            min rr' (sf/stf)
+          _ -> rr'
+        ss = case beta of
+          TensionAtLeast _ | bounded ->
+            min ss' (st/stf)
+          _ -> ss'
+        u = z0 ^+^ rr *^ Point (dx*ct - dy*st) (dy*ct + dx*st)  -- z0 + rr * (rotate theta chord)
+        v = z1 ^-^ ss *^ Point (dx*cf + dy*sf) (dy*cf - dx*sf)  -- z1 - ss * (rotate (-phi) chord)
+
+constant1, constant2, sqrt2 :: Double
+constant1 = 3/2*(sqrt 5 - 1)
+constant2 = 3/2*(3 - sqrt 5)
+sqrt2 = sqrt 2
+
+-- another magic formula by John Hobby.
+velocity :: Double -> Double -> Double
+         -> Double -> Tension -> Double
+velocity st sf ct cf t =
+  (2 + sqrt2 * (st - sf/16)*(sf - st/16)*(ct - cf)) /
+  ((3 + constant1*ct + constant2*cf) * tensionValue t)
diff --git a/Geom2D/CubicBezier/Outline.hs b/Geom2D/CubicBezier/Outline.hs
--- a/Geom2D/CubicBezier/Outline.hs
+++ b/Geom2D/CubicBezier/Outline.hs
@@ -7,79 +7,24 @@
 import Geom2D.CubicBezier.Basic
 import Geom2D.CubicBezier.Approximate
 import Geom2D.CubicBezier.Curvature
-import qualified Data.Map as M
-import Data.Function
-import Data.List
 
 offsetPoint :: Double -> Point -> Point -> Point
 offsetPoint dist start tangent =
   start ^+^ (rotate90L $* dist *^ normVector tangent)
 
-bezierOffsetPoint :: CubicBezier -> Double -> Double -> Point
-bezierOffsetPoint cb dist t =
-  uncurry (offsetPoint dist) $
-  evalBezierDeriv cb t
+bezierOffsetPoint :: CubicBezier -> Double -> Double -> (Point, Point)
+bezierOffsetPoint cb dist t = (offsetPoint dist p p', p')
+  where (p, p') = evalBezierDeriv cb t
 
 -- Approximate the bezier curve offset by dist.  A positive value
 -- means to the left, a negative to the right.
-approximateOffset :: CubicBezier -> Double -> Double -> (CubicBezier, Double, Double)
-approximateOffset cb@(CubicBezier p1 p2 p3 p4) dist tol =
-  approximateCurveWithParams offsetCb points ts tol
-  where tan1     = p2 ^-^ p1
-        tan2     = p4 ^-^ p3
-        offsetCb = CubicBezier
-                   (offsetPoint dist p1 tan1)
-                   (offsetPoint dist p2 tan1)
-                   (offsetPoint dist p3 tan2)
-                   (offsetPoint dist p4 tan2)
-        points   = map (bezierOffsetPoint cb dist) ts
-        ts = [i/16 | i <- [1..15]]
-
--- subdivide the original curve and approximate the offset until
--- the maximum error is below tolerance
 offsetSegment :: Double -> Double -> CubicBezier -> [CubicBezier]
-offsetSegment dist tol cb
-  | err <= tol = [cb_out]
-  | otherwise     = offsetSegment dist tol cb_l ++
-                    offsetSegment dist tol cb_r 
-  where
-    (cb_out, t, err) = approximateOffset cb dist tol
-    (cb_l, cb_r) = splitBezier cb t
-
-data OutlineSegment = OutlineSegment {
-  os_t_min :: Double,  -- the least t param of the segment in the original curve
-  os_t_err :: Double,  -- the param where the error is maximal
-  os_curve :: CubicBezier, -- the segment on the original curve
-  os_outline :: CubicBezier } -- the outline of the segment
-
--- Keep a map from maxError to OutlineSegment for each subsegment to keep
--- track of the segment with the maximum error.  This ensures a n
--- log(n) execution time, rather than n^2 when a list is used.
-offsetMax :: Double -> Double -> Int ->
-             M.Map Double OutlineSegment ->
-             [CubicBezier]
-offsetMax dist tol n segments
-  | n <= 1 = error "minimum segments to offset is 1"
-  | (n == 1) || (err < tol) = map os_outline $
-                              sortBy (compare `on` os_t_min) $
-                              M.elems segments
+offsetSegment dist tol cb =
+  approximatePath (bezierOffsetPoint cb dist) 15 tol 0 1
 
-    -- split the maximum curve in two and add the two segments to the map
-  | otherwise = offsetMax dist tol (n-1) $
-                M.insert err_l (OutlineSegment t_min t_err_l cb_l outline_l) $
-                M.insert err_r (OutlineSegment t_err t_err_r cb_r outline_r) $
-                newSegments
-  where
-    ((err, OutlineSegment t_min t_err curve _), newSegments) = M.deleteFindMax segments
-    (cb_l, cb_r) = splitBezier curve t_err
-    (outline_l, t_err_l, err_l)  = approximateOffset cb_l dist tol
-    (outline_r, t_err_r, err_r)  = approximateOffset cb_r dist tol
-    
 offsetSegmentMax :: Int -> Double -> Double -> CubicBezier -> [CubicBezier]
-offsetSegmentMax n dist tol cb =
-  offsetMax dist tol n segments
-  where segments              = M.singleton err (OutlineSegment 0 t_err cb outline)
-        (outline, t_err, err) = approximateOffset cb dist tol
+offsetSegmentMax m dist tol cb =
+  approximatePathMax m (bezierOffsetPoint cb dist) 15 tol 0 1
 
 -- | Calculate an offset path from the bezier curve to within
 -- tolerance.  If the distance is positive offset to the left,
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,339 +1,28 @@
-GNU GENERAL PUBLIC LICENSE
-                       Version 2, June 1991
-
- Copyright (C) 1989, 1991 Free Software Foundation, Inc.,
- 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
- Everyone is permitted to copy and distribute verbatim copies
- of this license document, but changing it is not allowed.
-
-                            Preamble
-
-  The licenses for most software are designed to take away your
-freedom to share and change it.  By contrast, the GNU General Public
-License is intended to guarantee your freedom to share and change free
-software--to make sure the software is free for all its users.  This
-General Public License applies to most of the Free Software
-Foundation's software and to any other program whose authors commit to
-using it.  (Some other Free Software Foundation software is covered by
-the GNU Lesser General Public License instead.)  You can apply it to
-your programs, too.
-
-  When we speak of free software, we are referring to freedom, not
-price.  Our General Public Licenses are designed to make sure that you
-have the freedom to distribute copies of free software (and charge for
-this service if you wish), that you receive source code or can get it
-if you want it, that you can change the software or use pieces of it
-in new free programs; and that you know you can do these things.
-
-  To protect your rights, we need to make restrictions that forbid
-anyone to deny you these rights or to ask you to surrender the rights.
-These restrictions translate to certain responsibilities for you if you
-distribute copies of the software, or if you modify it.
-
-  For example, if you distribute copies of such a program, whether
-gratis or for a fee, you must give the recipients all the rights that
-you have.  You must make sure that they, too, receive or can get the
-source code.  And you must show them these terms so they know their
-rights.
-
-  We protect your rights with two steps: (1) copyright the software, and
-(2) offer you this license which gives you legal permission to copy,
-distribute and/or modify the software.
-
-  Also, for each author's protection and ours, we want to make certain
-that everyone understands that there is no warranty for this free
-software.  If the software is modified by someone else and passed on, we
-want its recipients to know that what they have is not the original, so
-that any problems introduced by others will not reflect on the original
-authors' reputations.
-
-  Finally, any free program is threatened constantly by software
-patents.  We wish to avoid the danger that redistributors of a free
-program will individually obtain patent licenses, in effect making the
-program proprietary.  To prevent this, we have made it clear that any
-patent must be licensed for everyone's free use or not licensed at all.
-
-  The precise terms and conditions for copying, distribution and
-modification follow.
-
-                    GNU GENERAL PUBLIC LICENSE
-   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
-
-  0. This License applies to any program or other work which contains
-a notice placed by the copyright holder saying it may be distributed
-under the terms of this General Public License.  The "Program", below,
-refers to any such program or work, and a "work based on the Program"
-means either the Program or any derivative work under copyright law:
-that is to say, a work containing the Program or a portion of it,
-either verbatim or with modifications and/or translated into another
-language.  (Hereinafter, translation is included without limitation in
-the term "modification".)  Each licensee is addressed as "you".
-
-Activities other than copying, distribution and modification are not
-covered by this License; they are outside its scope.  The act of
-running the Program is not restricted, and the output from the Program
-is covered only if its contents constitute a work based on the
-Program (independent of having been made by running the Program).
-Whether that is true depends on what the Program does.
-
-  1. You may copy and distribute verbatim copies of the Program's
-source code as you receive it, in any medium, provided that you
-conspicuously and appropriately publish on each copy an appropriate
-copyright notice and disclaimer of warranty; keep intact all the
-notices that refer to this License and to the absence of any warranty;
-and give any other recipients of the Program a copy of this License
-along with the Program.
-
-You may charge a fee for the physical act of transferring a copy, and
-you may at your option offer warranty protection in exchange for a fee.
-
-  2. You may modify your copy or copies of the Program or any portion
-of it, thus forming a work based on the Program, and copy and
-distribute such modifications or work under the terms of Section 1
-above, provided that you also meet all of these conditions:
-
-    a) You must cause the modified files to carry prominent notices
-    stating that you changed the files and the date of any change.
-
-    b) You must cause any work that you distribute or publish, that in
-    whole or in part contains or is derived from the Program or any
-    part thereof, to be licensed as a whole at no charge to all third
-    parties under the terms of this License.
-
-    c) If the modified program normally reads commands interactively
-    when run, you must cause it, when started running for such
-    interactive use in the most ordinary way, to print or display an
-    announcement including an appropriate copyright notice and a
-    notice that there is no warranty (or else, saying that you provide
-    a warranty) and that users may redistribute the program under
-    these conditions, and telling the user how to view a copy of this
-    License.  (Exception: if the Program itself is interactive but
-    does not normally print such an announcement, your work based on
-    the Program is not required to print an announcement.)
-
-These requirements apply to the modified work as a whole.  If
-identifiable sections of that work are not derived from the Program,
-and can be reasonably considered independent and separate works in
-themselves, then this License, and its terms, do not apply to those
-sections when you distribute them as separate works.  But when you
-distribute the same sections as part of a whole which is a work based
-on the Program, the distribution of the whole must be on the terms of
-this License, whose permissions for other licensees extend to the
-entire whole, and thus to each and every part regardless of who wrote it.
-
-Thus, it is not the intent of this section to claim rights or contest
-your rights to work written entirely by you; rather, the intent is to
-exercise the right to control the distribution of derivative or
-collective works based on the Program.
-
-In addition, mere aggregation of another work not based on the Program
-with the Program (or with a work based on the Program) on a volume of
-a storage or distribution medium does not bring the other work under
-the scope of this License.
-
-  3. You may copy and distribute the Program (or a work based on it,
-under Section 2) in object code or executable form under the terms of
-Sections 1 and 2 above provided that you also do one of the following:
-
-    a) Accompany it with the complete corresponding machine-readable
-    source code, which must be distributed under the terms of Sections
-    1 and 2 above on a medium customarily used for software interchange; or,
-
-    b) Accompany it with a written offer, valid for at least three
-    years, to give any third party, for a charge no more than your
-    cost of physically performing source distribution, a complete
-    machine-readable copy of the corresponding source code, to be
-    distributed under the terms of Sections 1 and 2 above on a medium
-    customarily used for software interchange; or,
-
-    c) Accompany it with the information you received as to the offer
-    to distribute corresponding source code.  (This alternative is
-    allowed only for noncommercial distribution and only if you
-    received the program in object code or executable form with such
-    an offer, in accord with Subsection b above.)
-
-The source code for a work means the preferred form of the work for
-making modifications to it.  For an executable work, complete source
-code means all the source code for all modules it contains, plus any
-associated interface definition files, plus the scripts used to
-control compilation and installation of the executable.  However, as a
-special exception, the source code distributed need not include
-anything that is normally distributed (in either source or binary
-form) with the major components (compiler, kernel, and so on) of the
-operating system on which the executable runs, unless that component
-itself accompanies the executable.
-
-If distribution of executable or object code is made by offering
-access to copy from a designated place, then offering equivalent
-access to copy the source code from the same place counts as
-distribution of the source code, even though third parties are not
-compelled to copy the source along with the object code.
-
-  4. You may not copy, modify, sublicense, or distribute the Program
-except as expressly provided under this License.  Any attempt
-otherwise to copy, modify, sublicense or distribute the Program is
-void, and will automatically terminate your rights under this License.
-However, parties who have received copies, or rights, from you under
-this License will not have their licenses terminated so long as such
-parties remain in full compliance.
-
-  5. You are not required to accept this License, since you have not
-signed it.  However, nothing else grants you permission to modify or
-distribute the Program or its derivative works.  These actions are
-prohibited by law if you do not accept this License.  Therefore, by
-modifying or distributing the Program (or any work based on the
-Program), you indicate your acceptance of this License to do so, and
-all its terms and conditions for copying, distributing or modifying
-the Program or works based on it.
-
-  6. Each time you redistribute the Program (or any work based on the
-Program), the recipient automatically receives a license from the
-original licensor to copy, distribute or modify the Program subject to
-these terms and conditions.  You may not impose any further
-restrictions on the recipients' exercise of the rights granted herein.
-You are not responsible for enforcing compliance by third parties to
-this License.
-
-  7. If, as a consequence of a court judgment or allegation of patent
-infringement or for any other reason (not limited to patent issues),
-conditions are imposed on you (whether by court order, agreement or
-otherwise) that contradict the conditions of this License, they do not
-excuse you from the conditions of this License.  If you cannot
-distribute so as to satisfy simultaneously your obligations under this
-License and any other pertinent obligations, then as a consequence you
-may not distribute the Program at all.  For example, if a patent
-license would not permit royalty-free redistribution of the Program by
-all those who receive copies directly or indirectly through you, then
-the only way you could satisfy both it and this License would be to
-refrain entirely from distribution of the Program.
-
-If any portion of this section is held invalid or unenforceable under
-any particular circumstance, the balance of the section is intended to
-apply and the section as a whole is intended to apply in other
-circumstances.
-
-It is not the purpose of this section to induce you to infringe any
-patents or other property right claims or to contest validity of any
-such claims; this section has the sole purpose of protecting the
-integrity of the free software distribution system, which is
-implemented by public license practices.  Many people have made
-generous contributions to the wide range of software distributed
-through that system in reliance on consistent application of that
-system; it is up to the author/donor to decide if he or she is willing
-to distribute software through any other system and a licensee cannot
-impose that choice.
-
-This section is intended to make thoroughly clear what is believed to
-be a consequence of the rest of this License.
-
-  8. If the distribution and/or use of the Program is restricted in
-certain countries either by patents or by copyrighted interfaces, the
-original copyright holder who places the Program under this License
-may add an explicit geographical distribution limitation excluding
-those countries, so that distribution is permitted only in or among
-countries not thus excluded.  In such case, this License incorporates
-the limitation as if written in the body of this License.
-
-  9. The Free Software Foundation may publish revised and/or new versions
-of the General Public License from time to time.  Such new versions will
-be similar in spirit to the present version, but may differ in detail to
-address new problems or concerns.
-
-Each version is given a distinguishing version number.  If the Program
-specifies a version number of this License which applies to it and "any
-later version", you have the option of following the terms and conditions
-either of that version or of any later version published by the Free
-Software Foundation.  If the Program does not specify a version number of
-this License, you may choose any version ever published by the Free Software
-Foundation.
-
-  10. If you wish to incorporate parts of the Program into other free
-programs whose distribution conditions are different, write to the author
-to ask for permission.  For software which is copyrighted by the Free
-Software Foundation, write to the Free Software Foundation; we sometimes
-make exceptions for this.  Our decision will be guided by the two goals
-of preserving the free status of all derivatives of our free software and
-of promoting the sharing and reuse of software generally.
-
-                            NO WARRANTY
-
-  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
-FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
-OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
-PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
-OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
-MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
-TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
-PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
-REPAIR OR CORRECTION.
-
-  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
-WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
-REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
-INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
-OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
-TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
-YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
-PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
-POSSIBILITY OF SUCH DAMAGES.
-
-                     END OF TERMS AND CONDITIONS
-
-            How to Apply These Terms to Your New Programs
-
-  If you develop a new program, and you want it to be of the greatest
-possible use to the public, the best way to achieve this is to make it
-free software which everyone can redistribute and change under these terms.
-
-  To do so, attach the following notices to the program.  It is safest
-to attach them to the start of each source file to most effectively
-convey the exclusion of warranty; and each file should have at least
-the "copyright" line and a pointer to where the full notice is found.
-
-    Haskell library for manipulating cubic bezier curves
-    Copyright (C) 2013  Kristof Bastiaensen
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-Also add information on how to contact you by electronic and paper mail.
-
-If the program is interactive, make it output a short notice like this
-when it starts in an interactive mode:
-
-    Gnomovision version 69, Copyright (C) year name of author
-    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
-    This is free software, and you are welcome to redistribute it
-    under certain conditions; type `show c' for details.
-
-The hypothetical commands `show w' and `show c' should show the appropriate
-parts of the General Public License.  Of course, the commands you use may
-be called something other than `show w' and `show c'; they could even be
-mouse-clicks or menu items--whatever suits your program.
-
-You should also get your employer (if you work as a programmer) or your
-school, if any, to sign a "copyright disclaimer" for the program, if
-necessary.  Here is a sample; alter the names:
+Copyright (c) 2013, Kristof Bastiaensen
+All rights reserved.
 
-  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
-  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
 
-  {signature of Ty Coon}, 1 April 1989
-  Ty Coon, President of Vice
+    Redistributions of source code must retain the above copyright
+    notice, this list of conditions and the following disclaimer.
+    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.  Neither the name of the organisation nor the
+    names of its contributors may be used to endorse or promote
+    products derived from this software without specific prior written
+    permission.
 
-This General Public License does not permit incorporating your program into
-proprietary programs.  If your program is a subroutine library, you may
-consider it more useful to permit linking proprietary applications with the
-library.  If this is what you want to do, use the GNU Lesser General
-Public License instead of this License.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS 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 COPYRIGHT
+HOLDER 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/Math/BernsteinPoly.hs b/Math/BernsteinPoly.hs
--- a/Math/BernsteinPoly.hs
+++ b/Math/BernsteinPoly.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE BangPatterns #-}
 module Math.BernsteinPoly
        (BernsteinPoly(..), bernsteinSubsegment, listToBernstein, zeroPoly, (~*), (*~), (~+),
         (~-), degreeElevate, bernsteinSplit, bernsteinEval,
@@ -7,8 +8,8 @@
 import Data.List
 
 data BernsteinPoly = BernsteinPoly {
-  bernsteinDegree :: Int,
-  bernsteinCoeffs :: [Double] }
+  bernsteinDegree :: !Int,
+  bernsteinCoeffs :: ![Double] }
                    deriving Show
 
 infixl 7 ~*, *~
@@ -49,32 +50,39 @@
         a' = zipWith (*) a (binCoeff la)
         b' = zipWith (*) b (binCoeff lb)
 
-degreeElevate' :: BernsteinPoly -> Int -> BernsteinPoly
-degreeElevate' b 0 = b
-degreeElevate' (BernsteinPoly lp p) times =
-  degreeElevate' (BernsteinPoly (lp+1) (head p:inner p 1)) (times-1)
-  where
-    inner [a] _ = [a]
-    inner (a:b:rest) i =
-      (i*a/fromIntegral lp + b*(1 - i/fromIntegral lp))
-      : inner (b:rest) (i+1)
 
 -- find the binomial coefficients of degree n.
 binCoeff :: Int -> [Double]
 binCoeff n = map fromIntegral $
              scanl (\x m -> x * (n-m+1) `quot` m) 1 [1..n]
 
--- | Degree elevate a bernstein polynomail.
+-- | Degree elevate a bernstein polynomail a number of times.
 degreeElevate :: BernsteinPoly -> Int -> BernsteinPoly
-degreeElevate l times = degreeElevate' l times
+degreeElevate b 0 = b
+degreeElevate (BernsteinPoly lp p) times =
+  degreeElevate (BernsteinPoly (lp+1) (head p:inner p 1)) (times-1)
+  where
+    n = fromIntegral lp
+    inner []  _ = error "empty bernstein coefficients"
+    inner [a] _ = [a]
+    inner (a:b:rest) i =
+      (i*a/(n+1) + b*(1 - i/(n+1)))
+      : inner (b:rest) (i+1)
 
+
 -- | Evaluate the bernstein polynomial.
 bernsteinEval :: BernsteinPoly -> Double -> Double
-bernsteinEval (BernsteinPoly lp p) t = foldl' addcoeff 0 $
-                                       zip3 ts (binCoeff lp) p
-  where ts = iterate (*t) 1
-        u = 1-t
-        addcoeff a (s, d, b) = (a*u + b*s*d)
+bernsteinEval (BernsteinPoly lp [b]) _ = b
+bernsteinEval (BernsteinPoly lp (b':bs)) t = go t n (b'*u) 2 bs
+  where u = 1-t
+        n = fromIntegral lp
+        go !tn !bc !tmp _  [b] = tmp + tn*bc*b
+        go !tn !bc !tmp !i (b:rest) =
+          go (tn*t)         -- tn
+          (bc*(n-i+1)/i)    -- bc
+          ((tmp + tn*bc*b)*u) -- tmp
+          (i+1)             -- i
+          rest
 
 -- | Evaluate the bernstein polynomial and its derivatives.
 bernsteinEvalDerivs :: BernsteinPoly -> Double -> [Double]
@@ -88,8 +96,7 @@
 bernsteinDeriv (BernsteinPoly 0 _) = zeroPoly
 bernsteinDeriv (BernsteinPoly lp p) =
   BernsteinPoly (lp-1) $
-  map (* fromIntegral lp) $
-  zipWith (-) (tail p) p
+  map (* fromIntegral lp) $ zipWith (-) (tail p) p
 
 -- | Split a bernstein polynomial
 bernsteinSplit :: BernsteinPoly -> Double -> (BernsteinPoly, BernsteinPoly)
diff --git a/cubicbezier.cabal b/cubicbezier.cabal
--- a/cubicbezier.cabal
+++ b/cubicbezier.cabal
@@ -1,10 +1,10 @@
 Name:		cubicbezier
-Version: 	0.1.0
+Version: 	0.2.0
 Synopsis:	Efficient manipulating of 2D cubic bezier curves.
 Category: 	Graphics, Geometry, Typography
 Copyright: 	Kristof Bastiaensen (2013)
 Stability:	Unstable
-License:	GPL-2
+License:	BSD3
 License-file:	LICENSE
 Author:		Kristof Bastiaensen
 Maintainer:	Kristof Bastiaensen
@@ -26,7 +26,8 @@
   location:	https://github.com/kuribas/cubicbezier
 
 Library
-  Build-depends: base >= 3 && < 5, containers > 0.4, integration >= 0.1.1
+  Ghc-options: -Wall
+  Build-depends: base >= 3 && < 5, containers > 0.4, integration >= 0.1.1, deepseq >= 1.3.0
   Exposed-Modules:
     Geom2D
     Geom2D.CubicBezier
@@ -35,6 +36,7 @@
     Geom2D.CubicBezier.Outline
     Geom2D.CubicBezier.Curvature
     Geom2D.CubicBezier.Intersection
+    Geom2D.CubicBezier.MetaPath
     Math.BernsteinPoly
   Other-Modules:
     Geom2D.CubicBezier.Numeric