diff --git a/Graphics/Typography.lhs b/Graphics/Typography.lhs
deleted file mode 100644
--- a/Graphics/Typography.lhs
+++ /dev/null
@@ -1,89 +0,0 @@
-\begin{code}
-{-# OPTIONS -XFlexibleInstances -XNamedFieldPuns #-}
--- | This module contains basic tools for geometric types and functions.
-module Graphics.Typography (Matrix2(..),
-                            inverse,rotation,
-                            Geometric(..),
-                            leftMost,rightMost,topMost,bottomMost)
-       where
-
-import Algebra.Polynomials.Numerical
-
--- | The type of the transformation matrices used in all geometrical applications.
-data Matrix2 a=
-  -- | The application of @Matrix2 a b c d@ to vector @(x,y)@ should be
-  -- @(ax+by,cx+dy)@.
-  Matrix2 a a a a deriving (Show, Read, Eq)
-
--- | Inverses an inversible matrix. If it is not inversible,
--- The behaviour is undefined.
-inverse::(Fractional a, Num a)=>Matrix2 a->Matrix2 a
-inverse (Matrix2 a b c d)=
-  let det=a*d-c*b in
-  Matrix2 (d/det) (-b/det) (-c/det) (a/det)
-
-     
-
-instance Num a=>Num (Matrix2 a) where
-  (+) (Matrix2 a b c d) (Matrix2 e f g h)=
-    Matrix2 (a+e) (b+f) (c+g) (d+h)
-  (*) (Matrix2 a b c d) (Matrix2 e f g h)=
-    Matrix2 (a*e+b*g) (a*f+b*h) (c*e+d*g) (c*f+d*h)
-  fromInteger a=Matrix2 (fromInteger a) 0 0 (fromInteger a)
-  abs=undefined
-  signum=undefined
-
-instance Intervalize Matrix2 where
-  intervalize (Matrix2 a b c d)=
-    Matrix2 (interval a) (interval b) (interval c) (interval d)
-
-  intersects (Matrix2 a b c d) (Matrix2 a' b' c' d')=
-    (intersectsd a a') &&
-    (intersectsd b b') &&
-    (intersectsd c c') &&
-    (intersectsd d d')
-    
--- | A class for applying geometric applications to objects
-class Geometric g where
-  translate::Double->Double->g->g
-  apply::Matrix2 Double->g->g
-
--- | The matrix of a rotation
-rotation::Floating a=>a->Matrix2 a
-rotation theta=
-  let ct=cos theta
-      st=sin theta
-  in
-   Matrix2 ct (-st) st ct
-
-instance Geometric g=>Geometric [g] where
-  
-  translate x y cur=map (translate x y) cur
-  apply m cur=map (apply m) cur
-  
-
--- | @'leftMost' a b@ is the leftmost point between @a@ and @b@.
-leftMost::(Double,Double)->(Double,Double)->(Double,Double)
-leftMost u@(a,_) v@(b,_)
-  | a<b = u
-  | otherwise = v
--- | @'rightMost' a b@ is the rightmost point between @a@ and @b@.
-rightMost::(Double,Double)->(Double,Double)->(Double,Double)
-rightMost u@(a,_) v@(b,_)
-  | a<b = v
-  | otherwise = u
--- | @'bottomMost' a b@ is the lower point between @a@ and @b@.
-bottomMost::(Double,Double)->(Double,Double)->(Double,Double)
-bottomMost u@(_,a) v@(_,b)
-  | a<b = u
-  | otherwise = v
--- | @'topMost' a b@ is the upper point between @a@ and @b@.
-topMost::(Double,Double)->(Double,Double)->(Double,Double)
-topMost u@(_,a) v@(_,b)
-  | a<b = v
-  | otherwise = u
-
-
-
-
-\end{code}
diff --git a/Graphics/Typography/Approximation.lhs b/Graphics/Typography/Approximation.lhs
deleted file mode 100644
--- a/Graphics/Typography/Approximation.lhs
+++ /dev/null
@@ -1,170 +0,0 @@
-\begin{code}
-{-# OPTIONS -XRecordWildCards -XNamedFieldPuns #-}
--- | This module contains the function to approximate a list of curves with
--- degree 3 Bezier curves, using a least squares method.
-
-module Graphics.Typography.Approximation(approximate) where
-
-import qualified Data.Vector.Unboxed as UV
-import Graphics.Typography.Bezier
-import Graphics.Typography
-import Algebra.Polynomials.Bernstein
-
-import Algebra.Polynomials.Numerical
--- import Debug.Trace
-rnd::Interval->Double
-rnd (Interval a b)=(a+b)/2
-
-
--- | Approximates a list of 'Curves' with a list of degree 3 Bernstein curves.
-approximate::[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
-approximate []=[]
-approximate l0@(h0:_)= -- traceShow "starting" $
-  let approx::Double->Double->[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
-      approx _ _ []=[]
-      approx x0 y0 (cc@(Circle {..}):s)= -- traceShow "circle" $
-        let theta=abs $ t1-t0 in
-        if theta <= pi/2 then
-          let x0_=cos $ theta/2
-              y0_=sin $ theta/2
-              x1_=(4-x0_)/3
-              y1_=(1-x0_)*(3-x0_)/(3*y0_)
-        
-              c0=cos $! theta/2+t0
-              s0=sin $! theta/2+t0
-        
-        
-              px0=c0*x0_ - s0*y0_
-              py0=s0*x0_ + c0*y0_
-        
-              px1=c0*x1_ - s0*y1_
-              py1=s0*x1_ + c0*y1_
-        
-              px2=c0*x1_ + s0*y1_
-              py2=s0*x1_ - c0*y1_
-        
-              -- px3=c0*x0_ + s0*y0_
-              -- py3=s0*x0_ - c0*y0_
-        
-              x1=cx0+(a*px0+b*py0)
-              y1=cy0+(c*px0+d*py0)
-              
-              (Matrix2 a b c d)=matrix
-              x=Bernsteinp 4 $ UV.fromList
-                [ x0, -- cx0+(a*px3+b*py3),
-                  cx0+(a*px2+b*py2),
-                  cx0+(a*px1+b*py1),
-                  x1]
-              
-              y=Bernsteinp 4 $ UV.fromList
-                [ y0, -- cy0+(c*px3+d*py3),
-                  cy0+(c*px2+d*py2),
-                  cy0+(c*px1+d*py1),
-                  y1 ]
-          in
-           (x,y):(approx x1 y1 s)
-     
-        else
-          let t1'=(t1+t0)/2 in
-          approx x0 y0 $ (cc { t1=t1' }):(cc { t0=t1' }):s
-{-
-      approx x0 y0 (h@(Bezier{}):s)=
-        -- incorrect !
-        (restriction (cx h) (t0 h) (t1 h),
-         restriction (cy h) (t0 h) (t1 h)):
-        (approx (UV.last $ coefs $ cx h)
-         (UV.last $ coefs $ cy h) s)
--}
-      -- Ce qui suit est une methode de moindres carres
-      approx x0 y0 (off_:s)= -- traceShow ("offset") $
-        -- On commence par chercher les points ou la derivee de la norme
-        -- de la tangente est maximale. C'est la qu'on va couper s'il y
-        -- a un probleme.
-        let bx=restriction (cx off_) (t0 off_) (t1 off_)
-            by=restriction (cy off_) (t0 off_) (t1 off_)
-            off=off_ { cx=bx,cy=by,t0=0,t1=1 }
-            ibx=elevate (bounds by-bounds bx) $ intervalize bx
-            iby=elevate (bounds bx-bounds by) $ intervalize by
-            points=
-              let np=10 in
-              map (\x->(x/np,x/np)) [0..np]
-            -- Ensuite, moindres carres standard, comme dans Hoschek 88.
-      
-            vx0=ibx?1-ibx?0
-            vy0=iby?1-iby?0
-            vx1=ibx?(bounds ibx-2)-ibx?(bounds ibx-1)
-            vy1=iby?(bounds iby-2)-iby?(bounds iby-1)
-    
-            (wx0,wy0)=evalCurve off 0
-            (wx1,wy1)=evalCurve off 1
-
-            wx=Bernsteinp 4 $ UV.fromList [wx0,wx0,wx1,wx1] :: Bernsteinp Int Interval
-            wy=Bernsteinp 4 $ UV.fromList [wy0,wy0,wy1,wy1] :: Bernsteinp Int Interval
-
-            bern1=Bernsteinp 4 $ UV.fromList [0,1,0,0] :: Bernsteinp Int Interval
-            bern2=Bernsteinp 4 $ UV.fromList [0,0,1,0] :: Bernsteinp Int Interval
-
-            sumAll a b c d x1 y1 ((h1,h2):ss)=
-        
-              let h=Interval h1 h2
-                  (xi,yi)=evalCurve off h
-            
-                  b1=eval bern1 h
-                  b2=eval bern2 h
-            
-                  a'=a + (vx0*vx0+vy0*vy0)*b1*b1
-                  b'=b + (vx0*vx1 + vy0*vy1)*b1*b2
-                  c'=c + (vx0*vx1 + vy0*vy1)*b1*b2
-                  d'=d + (vx1*vx1+vy1*vy1)*b2*b2
-            
-                  dx=xi-(eval wx h)
-                  dy=yi-(eval wy h)
-            
-                  x1'=x1 + (vx0*dx + vy0*dy)*b1
-                  y1'=y1 + (vx1*dx + vy1*dy)*b2
-              in
-               sumAll a' b' c' d' x1' y1' ss
-         
-            sumAll a b c d x1 y1 []=(a,b,c,d,x1,y1)
-                                    
-            (ra,rb,rc,rd,rx1,ry1)=sumAll 0 0 0 0 0 0 points
-            
-            (Matrix2 a_ b_ c_ d_)=inverse $ Matrix2 ra rb rc rd
-            lambda1=a_*rx1+b_*ry1
-            lambda2=c_*rx1+d_*ry1
-      
-            -- On a la courbe optimale. Il faut chercher ou on va couper, maintenant
-            xap=Bernsteinp 4 $ UV.fromList [wx0,
-                                            wx0+lambda1*vx0,
-                                            wx1+lambda2*vx1,
-                                            wx1]
-            yap=Bernsteinp 4 $ UV.fromList [wy0,
-                                            wy0+lambda1*vy0,
-                                            wy1+lambda2*vy1,
-                                            wy1]
-            (err,arg)=foldl (\m (h1,h2)->
-                              let (xi,yi)=evalCurve off (Interval h1 h2)
-                                  xj=eval xap (Interval h1 h2)
-                                  yj=eval yap (Interval h1 h2)
-                              in
-                               max m (iup $ abs $ (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj), (h1+h2)/2))
-                      (0,0) points
-        in
-         if err<=0.01 then
-           (desintervalize xap,desintervalize yap):(approx (rnd wx1) (rnd wy1) s)
-         else
-           approx x0 y0 $
-            (off { cx=restriction (cx off) 0 arg,
-                   cy=restriction (cy off) 0 arg }):
-            (off { cx=restriction (cx off) arg 1,
-                   cy=restriction (cy off) arg 1 }):s
-           
-      (x0h,y0h)=evalCurve h0 $ Interval (t0 h0) (t0 h0)
-      
-  in
-   approx (rnd x0h) (rnd y0h) l0
-           
-desintervalize::(Bernsteinp a Interval)->(Bernsteinp a Double)
-desintervalize b=b { coefs=UV.map rnd $ coefs b}
-  
-\end{code}
diff --git a/Graphics/Typography/Bezier.lhs b/Graphics/Typography/Bezier.lhs
deleted file mode 100644
--- a/Graphics/Typography/Bezier.lhs
+++ /dev/null
@@ -1,749 +0,0 @@
-\documentclass{article}
-%include lhs2TeX.fmt
-\begin{document}
-\begin{code}
-{-# OPTIONS -XUnboxedTuples -XBangPatterns -XNamedFieldPuns -XRecordWildCards -XMagicHash -cpp #-}
--- | This module contains the basic functions for manipulating Bezier curves. It is heavily
--- based on the book by N. M. Patrikalakis and T. Maekawa, Shape Interrogation for Computer
--- Aided Design and Manufacturing.
-
-module Graphics.Typography.Bezier (
-  Curve(..),line,bezier3,
-  offset,
-  inter,
-  evalCurve,distance,
-  left,bottom,right,top)  where
-
-import Algebra.Polynomials.Bernstein
-import Algebra.Polynomials.Numerical
-import Graphics.Typography
-import qualified Data.Vector as V
-import qualified Data.Vector.Unboxed as UV
-import Data.List (partition,sort)
-
--- | The type for representing all types of curves.
-data Curve=
-  Bezier { cx::Bernsteinp Int Double,
-           cy::Bernsteinp Int Double, 
-           t0::Double,
-           t1::Double }
-  
-  | Offset { cx::Bernsteinp Int Double,
-             cy::Bernsteinp Int Double,
-             t0::Double,
-             t1::Double,
-             matrix::Matrix2 Double
-           }
-  | Circle { cx0::Double,
-             cy0::Double,
-             t0::Double,
-             t1::Double,
-             matrix::Matrix2 Double
-           }
-  deriving (Show)
-
-
--- | The basic constructor for lines : a line is a degree 1 Bezier curve
-line::Double->Double->Double->Double->Curve
-line px py px' py'=Bezier { cx=Bernsteinp 2 $ UV.fromList [px,px'], 
-                            cy=Bernsteinp 2 $ UV.fromList [py,py'],
-                            t0=0,t1=1 }
-
--- | A shortcut to define degree 3 Bezier curves from points. If the control
--- points are @a,b,c,d@, the function should be called with
--- @'bezier3' xa ya xb yb xc yc xd yd@.
-bezier3::Double->Double->Double->Double->Double->Double->Double->Double->Curve
-bezier3 px0 py0 px1 py1 px2 py2 px3 py3=
-  Bezier { cx=Bernsteinp 4 $ UV.fromList [px0,px1,px2,px3],
-           cy=Bernsteinp 4 $ UV.fromList [py0,py1,py2,py3],
-           t0=0,t1=1 }
-
-instance Geometric Curve where
-  translate x y cur@(Circle{cx0,cy0})=        
-    cur { cx0=cx0+x,cy0=cy0+y }
-  translate x y cur=
-    cur { cx=(cx cur) { coefs=UV.map (+x) $ coefs $ cx cur},
-          cy=(cy cur) { coefs=UV.map (+y) $ coefs $ cy cur} }
-
-  apply m0@(Matrix2 a b c d) cir@(Circle{cx0,cy0,matrix})=
-    cir { cx0=a*cx0+b*cy0, cy0=c*cx0+d*cy0, matrix=m0*matrix }
-  apply (Matrix2 a b c d) cur=
-    cur { cx=(scale a $ cx cur)+(scale b $ cy cur),
-          cy=(scale c $ cx cur)+(scale d $ cy cur) }
-
-
-
--- | Gives the point corresponding to the given value of the parameter
-evalCurve::Curve->Interval->(Interval,Interval)
-evalCurve (Offset{..}) t=
-  let ix=intervalize cx
-      iy=intervalize cy
-      xt=eval ix t
-      yt=eval iy t
-      m@(Matrix2 a b c d)=intervalize matrix
-      (Matrix2 a_ b_ c_ d_)=inverse m
-      xt0'=eval (derivate ix) t
-      yt0'=eval (derivate iy) t
-      xt'=a_*xt0' + b_*yt0'
-      yt'=c_*xt0' + d_*yt0'
-      dd=sqrt $ xt'*xt' + yt'*yt'
-  in
-   (xt+(a*yt'-b*xt')/dd, yt+(c*yt'-d*xt')/dd)
-
-evalCurve (Circle{..}) alpha=
-  let xx=cos alpha
-      yy=sin alpha
-      (Matrix2 a b c d)=intervalize matrix
-  in
-   (interval cx0+a*xx+b*yy, interval cy0+c*xx+d*yy)
-   
-evalCurve (Bezier{..}) t=
-  let ix=intervalize cx
-      iy=intervalize cy
-      xx=eval ix t
-      yy=eval iy t
-  in
-   (xx,yy)
-   
-data Topo=Dehors | SurLaLigne | Dedans deriving Eq
-
--- | @'inter' c0 c1@ is a list of all possible points of intersection
--- between curves @c0@ and @c1@ : if @(u,v,w,x)@ is returned by 'inter',
--- then curve @c0@ may intersect with @c1@ between parameter values @u@
--- and @v@, which corresponds to parameter values between @w@ and @x@ for
--- @c1@. The implementation guarantees that all actual solutions are found,
--- but possibly false solutions may also be returned.
-
-inter::Curve->Curve->[((Double,Double,Double,Double))]
-inter op@(Offset { cx=bxp_,cy=byp_,matrix=mp,t0=t0a,t1=t1a })
-  (Offset { cx=bxq_,cy=byq_,matrix=mq,t0=t0b,t1=t1b })=
-  
-  -- Attention : verifier si c'est la meme generatrice
-  let thrx=1e-5
-      solutions _ []=[]
-      solutions thr boxes@(_:_)=
-        let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2,eq3])) boxes 
-        
-            (correct,toRefine)=partition (\(u,v,_,_,_,_,_,_)->
-                              let (xu,yu)=evalCurve op (Interval u u)
-                                  (xv,yv)=evalCurve op (Interval v v)
-                              in
-                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
-        in
-         correct++(solutions (thr/2) toRefine)
-  in
-   map (\(u,v,w,x,_,_,_,_)->(u,v,w,x)) $ solutions 1e-2 $
-   [(t0a,t1a,t0b,t1b,0,1,0,1)::
-       (Double,Double,Double,Double,Double,Double,Double,Double)]
-
-     where
-       
-       
-       imp@(Matrix2 ap bp cp dp)=intervalize mp
-       imq@(Matrix2 aq bq cq dq)=intervalize mq
-       (Matrix2 ap_ bp_ cp_ dp_)=inverse imp
-       (Matrix2 aq_ bq_ cq_ dq_)=inverse imq
-
-       bxp=intervalize bxp_
-       byp=intervalize byp_
-       bxq=intervalize bxq_
-       byq=intervalize byq_
-    
-       bxp4=promote 1 bxp
-       byp4=promote 1 byp
-       bxq4=promote 2 bxq
-       byq4=promote 2 byq
-    
-
-       bxp'=derivate bxp
-       byp'=derivate byp
-       bxq'=derivate bxq
-       byq'=derivate byq
-
-       bXp'=promote 1 $ (scale ap_ bxp')+(scale bp_ byp')
-       bYp'=promote 1 $ (scale cp_ bxp')+(scale dp_ byp')
-
-       bXq'=promote 2 $ (scale aq_ bxq')+(scale bq_ byq')
-       bYq'=promote 2 $ (scale cq_ bxq')+(scale dq_ byq')
-
-       bomp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp') :: Bernsteinp (Int,Int,Int,Int) Interval
-       bomq@(Bernsteinp _ omegaq)=(bXq'*bXq')+(bYq'*bYq') :: Bernsteinp (Int,Int,Int,Int) Interval
-       
-       au=
-         let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
-         max 0 $ fpred au_
-       bu=
-         let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
-         fsucc bu_
-       av=
-         let av_=sqrt $ UV.minimum $ UV.map ilow omegaq in
-         max 0 $ fpred av_
-       bv=
-         let bv_=sqrt $ UV.maximum $ UV.map iup omegaq in
-         fsucc bv_
-    
-       alphau=Bernsteinp (1,1,2,1) $ UV.fromList [Interval au au,Interval bu bu]
-       alphav=Bernsteinp (1,1,1,2) $ UV.fromList [Interval av av,Interval bv bv]
-
-       eq0=
-         ((bxp4*alphau*alphav) + (scale ap $ bYp'*alphav) - (scale bp $ bXp'*alphav)
-          -(bxq4*alphau*alphav) - (scale aq $ bYq'*alphau) + (scale bq $ bXq'*alphau))
-       eq1=
-         ((byp4*alphau*alphav) + (scale cp $ bYp'*alphav) - (scale dp $ bXp'*alphav)
-          -(byq4*alphau*alphav) - (scale cq $ bYq'*alphau) + (scale dq $ bXq'*alphau))    
-       eq2=bomp-(alphau*alphau)
-       eq3=bomq-(alphav*alphav)
-      
-    
-
-
-inter b@(Circle{}) a@(Offset{})=
-  map (\(i,j,k,l)->(k,l,i,j)) $ inter a b
-
-inter o@(Offset { cx=bxp, cy=byp, matrix=mp })
-  cir@(Circle{cx0,cy0,matrix=mq})=
-
-  let ix=intervalize bxp
-      iy=intervalize byp
-      m@(Matrix2 a b c d)=intervalize mp
-      (Matrix2 a_ b_ c_ d_)=inverse m
-      x'=derivate ix
-      y'=derivate iy
-      xx'=(scale a_ x')+(scale b_ y')
-      yy'=(scale c_ x')+(scale d_ y')
-      omega@(Bernsteinp _ omegap)=xx'*xx'+yy'*yy'
-      au=
-        let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
-        max 0 $ fpred au_
-      bu=
-        let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
-        fsucc bu_
-      
-      alphau=Bernsteinp (1,2) $ UV.fromList [Interval au au,Interval bu bu]
-      
-      lambda=(promote 1 omega) - alphau*alphau
-      -- Avant multiplication par M_C^-1
-      xx0=(promote 1 $ ix-(intervalize $ constant cx0))*alphau
-          +(promote 1 $ scale a yy'-scale b xx')
-      yy0=(promote 1 $ iy-(intervalize $ constant cy0))*alphau
-          +(promote 1 $ scale c yy'-scale d xx') :: Bernsteinp (Int,Int) Interval
-    
-      (Matrix2 ac_ bc_ cc_ dc_)=inverse $ intervalize mq
-      xx1=(scale ac_ xx0)+(scale bc_ yy0)
-      yy1=(scale cc_ xx0)+(scale dc_ yy0)
-         
-      eqc=xx1*xx1+yy1*yy1-alphau*alphau
-      
-      thrx=1e-5
-      
-      solutions _ []=[]
-      solutions thr boxes@(_:_)=
-        let sol0=concatMap (solve thr (V.fromList [eqc,lambda])) boxes 
-        
-            (correct,toRefine)=partition (\(u,v,_,_)->
-                              let (xu,yu)=evalCurve o (Interval u u)
-                                  (xv,yv)=evalCurve o (Interval v v)
-                              in
-                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
-        in
-         correct++(solutions (thr/2) toRefine)
-         
-
-      -- Removing false positives by computing the distance to the center of
-      -- the circle (this is quite fast).
-      
-      removeFalse cl0 (h@(_,v,_,_):h'@(u',_,_,_):s)=
-        let u''=(v+u')/2
-            (xu,yu)=evalCurve o (Interval u'' u'')
-            Interval dl du=distance xu yu cir
-            cl1
-              | du<1 = Dedans
-              | dl>1 = Dehors
-              | otherwise = SurLaLigne
-        in
-         if cl0/=cl1 then h:(removeFalse cl1 (h':s)) else
-           removeFalse cl1 (h':s)
-      removeFalse _ l=l
-      
-      initCl=
-        let (x0,y0)=evalCurve o (Interval (t0 o) (t0 o)) 
-            Interval dl du=distance x0 y0 cir
-        in
-        if dl>1 then Dehors else if du<1 then Dedans else SurLaLigne
-  in
-   foldl (\l (u,v,_,_)->
-              let (Interval xl xu,Interval yl yu)=evalCurve o (Interval u v) in
-              case angle (Interval xl xu) (Interval yl yu) cir of
-                Just (Interval a0 a1)->
-                  (u,v,a0,a1):l
-                Nothing->l
-            ) [] $ removeFalse initCl $ sort $ solutions 1e-2 [(t0 o,t1 o,0::Double,1::Double)]
-   
-  
-inter a@(Circle{cx0=x0a,cy0=y0a,matrix=ma})
-  b@(Circle{cx0=x0b,cy0=y0b,matrix=mb})=
-  
-  if (intervalize ma)`intersects`(intervalize mb) && x0a==x0b && y0a==y0b then
-    let up ix@(Interval _ x_) tt0 tt1
-          | x_<tt0 =
-            up (ix+(2*interval pi)) tt0 tt1
-          | otherwise = down ix tt0 tt1
-        down ix@(Interval x_ x__) tt0 tt1
-          | x_>tt1 =
-            down (ix-(2*interval pi)) tt0 tt1
-          | x__<tt0 =
-              Nothing
-          | otherwise =
-            Just ix
-            
-        alpha=up (interval $ t0 a) (t0 b) (t1 b)
-        beta=up (interval $ t0 b) (t0 b) (t1 b)
-    in
-     
-     case (alpha,beta) of
-       (Just aa,Just ab)->
-         case (up aa (t0 a) (t1 a),
-               up ab (t0 a) (t1 a)) of
-           
-           (Just ba,Just bb)
-             | ilow aa<=iup ab -> [(ilow ba, iup bb,
-                                    ilow aa, iup ab)]
-             | otherwise->
-               case (up (interval $ t0 b) (t0 a) (t1 a),
-                     up (interval $ t1 b) (t0 a) (t1 a)) of
-                 (Just b0,Just b1)->
-                   [(ilow b0,iup bb,
-                     t0 b, iup ab),
-                    (ilow ba,iup b1,
-                     ilow aa, t1 b)]
-                 _->[]
-           _->[]
-       _->[]
-
-  else
-    let thr=1e-5
-        solutions=solve thr (V.fromList [eq0,eq1]) (fpred u0,fsucc v0,
-                                                    fpred w0,fsucc x0)
-    in
-     foldl (\l (u,v,w,x)->
-             let alpha=angle (Interval u v) (Interval w x) a
-                 beta=angle (Interval u v) (Interval w x) b
-             in
-              case alpha of
-                Just (Interval a0l a0u)->
-                  case beta of
-                    Just (Interval b0l b0u)->(a0l,a0u,b0l,b0u):l
-                    _->l
-                _->l
-           ) [] solutions 
-  where
-    
-    ima@(Matrix2 am bm cm dm)=intervalize ma
-    
-    maxa=max (iup $ abs am+abs bm) (iup $ abs cm+abs dm)
-    (u0,v0,w0,x0)=(x0a-maxa,x0a+maxa,y0a-maxa,y0a+maxa)
-    
-    -- x-x0
-    xxa0=intervalize $ Bernsteinp (2,1) $ UV.fromList [-x0a,1-x0a] :: Bernsteinp (Int,Int) Interval
-    yya0=intervalize $ Bernsteinp (1,2) $ UV.fromList [-y0a,1-y0a] :: Bernsteinp (Int,Int) Interval
-    (Matrix2 aa_ ba_ ca_ da_)=inverse ima
-    xxa=(scale aa_ xxa0)+(scale ba_ yya0)::Bernsteinp (Int,Int) Interval
-    yya=(scale ca_ xxa0)+(scale da_ yya0)
-    
-    xxb0=intervalize $ Bernsteinp (2,1) $ UV.fromList [-x0b,1-x0b]
-    yyb0=intervalize $ Bernsteinp (1,2) $ UV.fromList [-y0b,1-y0b]
-    (Matrix2 ab_ bb_ cb_ db_)=inverse $ intervalize mb
-    xxb=(scale ab_ xxb0)+(scale bb_ yyb0)
-    yyb=(scale cb_ xxb0)+(scale db_ yyb0)
-    
-    c1=Bernsteinp (1,1) $ UV.singleton 1
-    
-    eq0=xxa*xxa+yya*yya-c1
-    eq1=xxb*xxb+yyb*yyb-c1
-
-inter op@(Bezier{cx=bxa,cy=bya,t0=t0a,t1=t1a}) (Bezier{cx=xb,cy=yb,t0=t0b,t1=t1b})=
-  
-  let p0=(promote 1 $ intervalize bxa)-(promote 2 $ intervalize xb) :: Bernsteinp (Int,Int) Interval
-      p1=(promote 1 $ intervalize bya)-(promote 2 $ intervalize yb) :: Bernsteinp (Int,Int) Interval
-      thrx=1e-2
-      solutions _ []=[]
-      solutions thr boxes@(_:_)=
-        let sol0=concatMap (solve thr (V.fromList [p0,p1])) boxes 
-        
-            (correct,toRefine)=partition (\(u,v,_,_)->
-                              let (xu,yu)=evalCurve op (Interval u u)
-                                  (xv,yv)=evalCurve op (Interval v v)
-                              in
-                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
-        in
-         correct++(solutions (thr/2) toRefine)
-  in
-   solutions 1e-2 [(t0a,t1a,t0b,t1b)]
-
-
-inter cir@(Circle{}) bez@(Bezier{})=map (\(u,v,w,x)->(w,x,u,v)) $ inter bez cir
-
-inter bez@(Bezier{}) cir@(Circle{})=
-  let xx=(intervalize $ cx bez)-(intervalize $ constant $ cx0 cir)
-      yy=(intervalize $ cy bez)-(intervalize $ constant $ cy0 cir)
-      (Matrix2 a b c d)=inverse $ intervalize $ matrix cir
-      xx0=scale a xx+scale b yy
-      yy0=scale c xx+scale d yy
-      
-      thrx=1e-5
-      
-      solutions _ []=[]
-      solutions thr boxes@(_:_)=
-        let sol0=concatMap (solve thr (V.singleton (xx0*xx0+yy0*yy0-(constant 1)))) boxes
-        
-            (correct,toRefine)=partition (\(u,v)->
-                              let (xu,yu)=evalCurve bez (Interval u u)
-                                  (xv,yv)=evalCurve bez (Interval v v)
-                              in
-                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
-        in
-         correct++(solutions (thr/2) toRefine)
-  in
-   foldl (\l (u,v)->
-              let (Interval xl xu,Interval yl yu)=evalCurve bez (Interval u v) in
-              case angle (Interval xl xu) (Interval yl yu) cir of
-                Just (Interval a0 a1)->
-                  (u,v,a0,a1):l
-                Nothing->l
-            ) [] $!
-   solutions (1e-2) [(t0 bez,t1 bez)]
-
-inter bez@(Bezier{}) off@(Offset{})=map (\(u,v,w,x)->(w,x,u,v)) $ inter off bez
-
-
-inter off@(Offset{}) bez@(Bezier{})=
-  
-  let thr=1e-2
-      thrx=1e-5
-      solutions _ []=[]
-      solutions thr boxes@(_:_)=
-        let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2])) boxes
-        
-            (correct,toRefine)=partition (\(u,v,_,_,_,_)->
-                              let (xu,yu)=evalCurve off (Interval u u)
-                                  (xv,yv)=evalCurve off (Interval v v)
-                              in
-                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
-        in
-         correct++(solutions (thr/2) toRefine)
-  in
-   map (\(a,b,c,d,_,_)->(a,b,c,d)) $ solutions 1e-2 $
-   [(0,1,0,1,0,1)::(Double,Double,Double,Double,Double,Double)]
-    where
-      
-      bxp=intervalize $ cx off
-      byp=intervalize $ cy off
-      bxq=intervalize $ cx bez
-      byq=intervalize $ cy bez
-      
-      bxp'=derivate bxp
-      byp'=derivate byp
-      bxp3=promote 1 bxp
-      byp3=promote 1 byp
-      bxq3=promote 2 bxq
-      byq3=promote 2 byq
-      
-      mp@(Matrix2 ap bp cp dp)=intervalize $ matrix $ off
-      (Matrix2 ap_ bp_ cp_ dp_)=inverse mp
-      bXp'=promote 1 $ (scale ap_ bxp')+(scale bp_ byp')
-      bYp'=promote 1 $ (scale cp_ bxp')+(scale dp_ byp')
-
-      omp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp')
-      au=
-        let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
-        max 0 $ fpred au_
-      bu=
-        let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
-        fsucc bu_
-    
-      alphau=Bernsteinp (1,1,2) $ UV.fromList [Interval au au,Interval bu bu]
-      eq0=bxp3*alphau + (scale ap bYp') - (scale bp bXp') - bxq3
-      eq1=byp3*alphau + (scale cp bYp') - (scale dp bXp') - byq3
-      eq2=alphau*alphau-omp
-      
-
-
-angle::Interval->Interval->Curve->Maybe Interval
-angle x y (Circle { cx0,cy0,matrix,t0,t1 })=
-  let vx=x-interval cx0
-      vy=y-interval cy0
-      
-      Matrix2 a b c d=inverse $ intervalize matrix
-      -- L'arithmetique d'intervalles fait un peu n'importe quoi
-      -- quand le vecteur est trop long. On le raccourcit.
-      alpha=
-        let co@(Interval col cou)=a*vx+b*vy
-            Interval sil siu=c*vx+d*vy
-            co2=
-              let (col2,cou2)=if col*col<cou*cou then (col*col,cou*cou) else
-                                (cou*cou,col*col)
-              in
-               Interval (fpred col2) (fsucc cou2)
-            si2=
-              let (sil2,siu2)=if sil*sil<siu*siu then (sil*sil,siu*siu) else
-                                (siu*siu,sil*sil)
-              in
-               Interval (fpred sil2) (fsucc siu2)
-            coco=co/(sqrt (co2+si2))
-            ac@(Interval acl acu)=acos $ Interval (max (-1) $ ilow coco) (min 1 $ iup coco)
-        in
-         if siu<0 then negate ac else
-           if sil>=0 then ac else
-             Interval (negate $ min (abs acl) (abs acu))
-             (max (abs acl) (abs acu))
-      up ix
-        | iup ix<t0 =
-          up $ ix+(2*interval pi)
-        | otherwise =
-            down ix
-      down ix
-        | ilow ix>t1 =
-          down $ ix-(2*interval pi)
-        | iup ix<t0 =
-          Nothing
-        | otherwise =
-          Just ix
-  in
-   up alpha
-
-
-angle _ _ _=error "angle"
-
--- | Pseudo-distance from a point to a curve. Is the result is
--- smaller than 1, the point is inside the curve. If it is greater
--- than 1, the point is outside. Else we don't know (as usual with
--- interval arithmetic).
-
-distance::Interval->Interval->Curve->Interval
-distance x0 y0 (Bezier{..})=
-  distance x0 y0 (Offset{cx,cy,t0,t1,matrix=Matrix2 1 0 0 1})
-  
-distance x0 y0 (Offset{..})=
-  let (Matrix2 a b c d)=inverse $ intervalize matrix
-      vx_=intervalize cx-(constant x0)
-      vy_=intervalize cy-(constant y0)
-      vx=scale a vx_+scale b vy_
-      vy=scale c vx_+scale d vy_
-      
-      dist=vx*vx+vy*vy
-  in
-   foldl (\di (u,v)->let di'=eval dist (Interval u v) in
-           if iup di<iup di' then di else di') (Interval (1/0) (1/0)) $
-   (t0,t0):(t1,t1):(solve 1e-5 (V.singleton (derivate dist)) (t0,t1))
-  
-  
-distance x1 y1 (Circle{..})=
-  let (Matrix2 a b c d)=inverse $ intervalize matrix
-      vx_=x1-Interval cx0 cx0
-      vy_=y1-Interval cy0 cy0
-      vx=a*vx_+b*vy_
-      vy=c*vx_+d*vy_
-  in
-   vx*vx+vy*vy
-   
--- | Offsets a given Bezier curve with the given pen matrix. The original
--- pen is a circle of radius one, the matrix, if inversible, is applied to it.
-
-offset::Matrix2 Double->Curve->[Curve]
-offset m (Bezier{cx=x@(Bernsteinp nx bx),cy=y@(Bernsteinp ny by)})=
-  if nx <=1 && ny <=1 then
-    [Circle { cx0=UV.head bx,cy0=UV.head by,t0=ilow 0,t1=iup $ 2*pi,matrix=m }]
-  else
-    [ c0,c1,c2,c3 ]
-    
-  where
-    im=intervalize m
-    (Matrix2 a_ b_ c_ d_)=inverse im
-    
-    ibx=intervalize x
-    iby=intervalize y
-    
-    lastCoef (Bernsteinp n c)
-      | n>=1 = UV.last c
-      | otherwise = 0
-    firstCoef (Bernsteinp n c)
-      | n>=1 = UV.head c
-      | otherwise = 0
-    
-    -- Premiere courbe offset
-    c0=Offset { cx=x, cy=y, t0=0,t1=1,matrix=m }
-    
-    -- Demi-cercle 1
-    
-    ibx'=derivate ibx
-    iby'=derivate iby
-    
-    -- Calcul du vecteur tangent au bout du premier
-
-    alpha0=
-      let xx0=lastCoef ibx'
-          yy0=lastCoef iby'
-
-          xx0'=a_*xx0+b_*yy0
-          yy0'=c_*xx0+d_*yy0
-          norm0=sqrt $ xx0'*xx0'+yy0'*yy0'
-          
-          xx'=xx0'/norm0
-          yy'=yy0'/norm0
-      in
-       if ilow xx'>=0 then
-         -(acos yy')
-       else
-         if iup xx'<=0 then
-           acos yy'
-         else
-           let Interval u v=acos yy' in
-           Interval (negate $ max (abs u) (abs v))
-           (max (abs u) (abs v))
-          
-                    
-    alpha0'=alpha0+interval pi
-    c1=Circle { cx0=lastCoef x,
-                cy0=lastCoef y,
-                t0=ilow alpha0,
-                t1=iup alpha0',
-                matrix=m }
-    
-    
-    -- Deuxieme courbe offset
-    c2=Offset { cx=reorient x,
-                cy=reorient y,
-                t0=0,t1=1,
-                matrix=m }
-       
-    -- Deuxieme demi-cercle
-    alpha1=
-      let xx0=firstCoef ibx'
-          yy0=firstCoef iby'
-
-          xx0'=a_*xx0+b_*yy0
-          yy0'=c_*xx0+d_*yy0
-          norm0=sqrt $ xx0'*xx0'+yy0'*yy0'
-          
-          xx'=xx0'/norm0
-          yy'=yy0'/norm0
-      in
-       if ilow xx'>=0 then
-         -(acos yy')
-       else
-         if iup xx'<=0 then
-           acos yy'
-         else
-           let Interval u v=acos yy' in
-           Interval (negate $ max (abs u) (abs v))
-           (max (abs u) (abs v))
-           
-    alpha1'=alpha1-pi
-    c3=Circle { cx0=firstCoef x, 
-                cy0=firstCoef y,
-                t0=ilow alpha1',
-                t1=iup alpha1,
-                matrix=m }
-    
-
-offset _ _=error "offset : undefined yet for other than Bezier"
-
-
-rnd::Interval->Double
-rnd (Interval a b)=(a+b)/2
-
-derivRoots::Double->Curve->([(Double,Double)],[(Double,Double)])
-derivRoots thr (Bezier{..})=
-  (solve thr (V.singleton $ derivate $ intervalize cx) (t0,t1),
-   solve thr (V.singleton $ derivate $ intervalize cy) (t0,t1))
-derivRoots thr (Offset{..})=
-  let ix=intervalize cx
-      iy=intervalize cy
-      x'=derivate ix
-      y'=derivate iy
-      m@(Matrix2 a b c d)=intervalize matrix
-      (Matrix2 a_ b_ c_ d_)=inverse m
-      
-      xx'=(scale a_ x')+(scale b_ y')
-      yy'=(scale c_ x')+(scale d_ y')
-      xx''=derivate xx'
-      yy''=derivate yy'
-      
-      omega=xx'*xx'+yy'*yy'
-      au=
-        let au_=sqrt $ UV.minimum $ UV.map ilow $ coefs omega in
-        max 0 $ fpred au_
-      bu=
-        let bu_=sqrt $ UV.maximum $ UV.map iup $ coefs omega in
-        fsucc bu_
-      alphau=Bernsteinp (1,2) $ UV.fromList [Interval au au,Interval bu bu]
-      
-      eqx0=(promote 1 yy'')*(alphau^(2::Int))-(promote 1 $ yy'*yy''+xx'*xx'')
-      eqy0=(promote 1 $ yy'*yy''+xx'*xx'')-(promote 1 xx'')*(alphau^(2::Int))
-      
-      eqx=(promote 1 x')*(alphau^(3::Int))+(scale a eqx0)+(scale b eqy0)
-          :: Bernsteinp (Int,Int) Interval
-      eqy=(promote 1 y')*(alphau^(3::Int))+(scale c eqx0)+(scale d eqy0)
-      
-      eq=(promote 1 omega)-alphau^(2::Int) :: Bernsteinp (Int,Int) Interval
-  in
-   (map (\(u,v,_,_)->(u,v)) $ solve thr (V.fromList [eqx,eq])
-    ((t0,t1,0,1)::(Double,Double,Double,Double)),
-    map (\(u,v,_,_)->(u,v)) $ solve thr (V.fromList [eqy,eq])
-    ((t0,t1,0,1)::(Double,Double,Double,Double)))
-derivRoots _ (Circle{..})=
-  let (Matrix2 a b c d)=intervalize matrix
-      aa=sqrt $ a*a+b*b
-      cc=sqrt $ c*c+d*d
-      sx
-        | ilow a>=0 = acos $ b/aa
-        | otherwise = negate $ acos $ b/aa
-      sy
-        | ilow c>=0 = acos $ d/cc
-        | otherwise = negate $ acos $ d/cc
-      Interval ux vx=(-sx)-pi/2
-      Interval uy vy=(-sy)-pi/2
-  in
-   ([(ux,vx)],[(uy,vy)])
-
--- | The leftmost point on a curve
-left::Curve->(Double,Double)
-left cur=
-  let (x,y)=
-        foldl (\m@(xx,_) (s,t)->
-                let m'@(xx',_)=evalCurve cur $ interval $ (s+t)/2 in
-                if ilow xx<ilow xx' then m else m') (1/0,1/0) $
-        (t0 cur,t0 cur):(t1 cur,t1 cur):(fst $ derivRoots 1e-5 cur)
-  in
-   (rnd x,rnd y)
--- | The bottommost point on a curve
-bottom::Curve->(Double,Double)
-bottom cur=
-  let (x,y)=
-        foldl (\m@(_,yy) (s,t)->
-                let m'@(_,yy')=evalCurve cur $ interval $ (s+t)/2 in
-                if ilow yy<ilow yy' then m else m') (1/0,1/0) $
-        (t0 cur,t0 cur):(t1 cur,t1 cur):(snd $ derivRoots 1e-5 cur)
-  in
-   (rnd x,rnd y)
--- | The rightmost point on a curve
-right::Curve->(Double,Double)
-right cur=
-  let (x,y)=
-        foldl (\m@(xx,_) (s,t)->
-                let m'@(xx',_)=evalCurve cur $ interval $ (s+t)/2 in
-                if iup xx>iup xx' then m else m') (-1/0,-1/0) $
-        (t0 cur,t0 cur):(t1 cur,t1 cur):(fst $ derivRoots 1e-5 cur)
-  in
-   (rnd x,rnd y)
--- | The topmost point on a curve
-top::Curve->(Double,Double)
-top cur=
-  let (x,y)=
-        foldl (\m@(_,yy) (s,t)->
-                let m'@(_,yy')=evalCurve cur $ interval $ (s+t)/2 in
-                if iup yy>iup yy' then m else m') (-1/0,-1/0) $
-        (t0 cur,t0 cur):(t1 cur,t1 cur):(snd $ derivRoots 1e-5 cur)
-  in
-   (rnd x,rnd y)
-\end{code}
diff --git a/Graphics/Typography/Geometry.lhs b/Graphics/Typography/Geometry.lhs
new file mode 100644
--- /dev/null
+++ b/Graphics/Typography/Geometry.lhs
@@ -0,0 +1,89 @@
+\begin{code}
+{-# OPTIONS -XFlexibleInstances -XNamedFieldPuns #-}
+-- | This module contains basic tools for geometric types and functions.
+module Graphics.Typography.Geometry (Matrix2(..),
+                            inverse,rotation,
+                            Geometric(..),
+                            leftMost,rightMost,topMost,bottomMost)
+       where
+
+import Algebra.Polynomials.Numerical
+
+-- | The type of the transformation matrices used in all geometrical applications.
+data Matrix2 a=
+  -- | The application of @Matrix2 a b c d@ to vector @(x,y)@ should be
+  -- @(ax+by,cx+dy)@.
+  Matrix2 a a a a deriving (Show, Read, Eq)
+
+-- | Inverses an inversible matrix. If it is not inversible,
+-- The behaviour is undefined.
+inverse::(Fractional a, Num a)=>Matrix2 a->Matrix2 a
+inverse (Matrix2 a b c d)=
+  let det=a*d-c*b in
+  Matrix2 (d/det) (-b/det) (-c/det) (a/det)
+
+     
+
+instance Num a=>Num (Matrix2 a) where
+  (+) (Matrix2 a b c d) (Matrix2 e f g h)=
+    Matrix2 (a+e) (b+f) (c+g) (d+h)
+  (*) (Matrix2 a b c d) (Matrix2 e f g h)=
+    Matrix2 (a*e+b*g) (a*f+b*h) (c*e+d*g) (c*f+d*h)
+  fromInteger a=Matrix2 (fromInteger a) 0 0 (fromInteger a)
+  abs=undefined
+  signum=undefined
+
+instance Intervalize Matrix2 where
+  intervalize (Matrix2 a b c d)=
+    Matrix2 (interval a) (interval b) (interval c) (interval d)
+
+  intersects (Matrix2 a b c d) (Matrix2 a' b' c' d')=
+    (intersectsd a a') &&
+    (intersectsd b b') &&
+    (intersectsd c c') &&
+    (intersectsd d d')
+    
+-- | A class for applying geometric applications to objects
+class Geometric g where
+  translate::Double->Double->g->g
+  apply::Matrix2 Double->g->g
+
+-- | The matrix of a rotation
+rotation::Floating a=>a->Matrix2 a
+rotation theta=
+  let ct=cos theta
+      st=sin theta
+  in
+   Matrix2 ct (-st) st ct
+
+instance Geometric g=>Geometric [g] where
+  
+  translate x y cur=map (translate x y) cur
+  apply m cur=map (apply m) cur
+  
+
+-- | @'leftMost' a b@ is the leftmost point between @a@ and @b@.
+leftMost::(Double,Double)->(Double,Double)->(Double,Double)
+leftMost u@(a,_) v@(b,_)
+  | a<b = u
+  | otherwise = v
+-- | @'rightMost' a b@ is the rightmost point between @a@ and @b@.
+rightMost::(Double,Double)->(Double,Double)->(Double,Double)
+rightMost u@(a,_) v@(b,_)
+  | a<b = v
+  | otherwise = u
+-- | @'bottomMost' a b@ is the lower point between @a@ and @b@.
+bottomMost::(Double,Double)->(Double,Double)->(Double,Double)
+bottomMost u@(_,a) v@(_,b)
+  | a<b = u
+  | otherwise = v
+-- | @'topMost' a b@ is the upper point between @a@ and @b@.
+topMost::(Double,Double)->(Double,Double)->(Double,Double)
+topMost u@(_,a) v@(_,b)
+  | a<b = v
+  | otherwise = u
+
+
+
+
+\end{code}
diff --git a/Graphics/Typography/Geometry/Approximation.lhs b/Graphics/Typography/Geometry/Approximation.lhs
new file mode 100644
--- /dev/null
+++ b/Graphics/Typography/Geometry/Approximation.lhs
@@ -0,0 +1,170 @@
+\begin{code}
+{-# OPTIONS -XRecordWildCards -XNamedFieldPuns #-}
+-- | This module contains the function to approximate a list of curves with
+-- degree 3 Bezier curves, using a least squares method.
+
+module Graphics.Typography.Geometry.Approximation(approximate) where
+
+import qualified Data.Vector.Unboxed as UV
+import Graphics.Typography.Geometry.Bezier
+import Graphics.Typography.Geometry
+import Algebra.Polynomials.Bernstein
+
+import Algebra.Polynomials.Numerical
+-- import Debug.Trace
+rnd::Interval->Double
+rnd (Interval a b)=(a+b)/2
+
+
+-- | Approximates a list of 'Curves' with a list of degree 3 Bernstein curves.
+approximate::[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
+approximate []=[]
+approximate l0@(h0:_)= -- traceShow "starting" $
+  let approx::Double->Double->[Curve]->[(Bernsteinp Int Double,Bernsteinp Int Double)]
+      approx _ _ []=[]
+      approx x0 y0 (cc@(Circle {..}):s)= -- traceShow "circle" $
+        let theta=abs $ t1-t0 in
+        if theta <= pi/2 then
+          let x0_=cos $ theta/2
+              y0_=sin $ theta/2
+              x1_=(4-x0_)/3
+              y1_=(1-x0_)*(3-x0_)/(3*y0_)
+        
+              c0=cos $! theta/2+t0
+              s0=sin $! theta/2+t0
+        
+        
+              px0=c0*x0_ - s0*y0_
+              py0=s0*x0_ + c0*y0_
+        
+              px1=c0*x1_ - s0*y1_
+              py1=s0*x1_ + c0*y1_
+        
+              px2=c0*x1_ + s0*y1_
+              py2=s0*x1_ - c0*y1_
+        
+              -- px3=c0*x0_ + s0*y0_
+              -- py3=s0*x0_ - c0*y0_
+        
+              x1=cx0+(a*px0+b*py0)
+              y1=cy0+(c*px0+d*py0)
+              
+              (Matrix2 a b c d)=matrix
+              x=Bernsteinp 4 $ UV.fromList
+                [ x0, -- cx0+(a*px3+b*py3),
+                  cx0+(a*px2+b*py2),
+                  cx0+(a*px1+b*py1),
+                  x1]
+              
+              y=Bernsteinp 4 $ UV.fromList
+                [ y0, -- cy0+(c*px3+d*py3),
+                  cy0+(c*px2+d*py2),
+                  cy0+(c*px1+d*py1),
+                  y1 ]
+          in
+           (x,y):(approx x1 y1 s)
+     
+        else
+          let t1'=(t1+t0)/2 in
+          approx x0 y0 $ (cc { t1=t1' }):(cc { t0=t1' }):s
+{-
+      approx x0 y0 (h@(Bezier{}):s)=
+        -- incorrect !
+        (restriction (cx h) (t0 h) (t1 h),
+         restriction (cy h) (t0 h) (t1 h)):
+        (approx (UV.last $ coefs $ cx h)
+         (UV.last $ coefs $ cy h) s)
+-}
+      -- Ce qui suit est une methode de moindres carres
+      approx x0 y0 (off_:s)= -- traceShow ("offset") $
+        -- On commence par chercher les points ou la derivee de la norme
+        -- de la tangente est maximale. C'est la qu'on va couper s'il y
+        -- a un probleme.
+        let bx=restriction (cx off_) (t0 off_) (t1 off_)
+            by=restriction (cy off_) (t0 off_) (t1 off_)
+            off=off_ { cx=bx,cy=by,t0=0,t1=1 }
+            ibx=elevate (bounds by-bounds bx) $ intervalize bx
+            iby=elevate (bounds bx-bounds by) $ intervalize by
+            points=
+              let np=10 in
+              map (\x->(x/np,x/np)) [0..np]
+            -- Ensuite, moindres carres standard, comme dans Hoschek 88.
+      
+            vx0=ibx?1-ibx?0
+            vy0=iby?1-iby?0
+            vx1=ibx?(bounds ibx-2)-ibx?(bounds ibx-1)
+            vy1=iby?(bounds iby-2)-iby?(bounds iby-1)
+    
+            (wx0,wy0)=evalCurve off 0
+            (wx1,wy1)=evalCurve off 1
+
+            wx=Bernsteinp 4 $ UV.fromList [wx0,wx0,wx1,wx1] :: Bernsteinp Int Interval
+            wy=Bernsteinp 4 $ UV.fromList [wy0,wy0,wy1,wy1] :: Bernsteinp Int Interval
+
+            bern1=Bernsteinp 4 $ UV.fromList [0,1,0,0] :: Bernsteinp Int Interval
+            bern2=Bernsteinp 4 $ UV.fromList [0,0,1,0] :: Bernsteinp Int Interval
+
+            sumAll a b c d x1 y1 ((h1,h2):ss)=
+        
+              let h=Interval h1 h2
+                  (xi,yi)=evalCurve off h
+            
+                  b1=eval bern1 h
+                  b2=eval bern2 h
+            
+                  a'=a + (vx0*vx0+vy0*vy0)*b1*b1
+                  b'=b + (vx0*vx1 + vy0*vy1)*b1*b2
+                  c'=c + (vx0*vx1 + vy0*vy1)*b1*b2
+                  d'=d + (vx1*vx1+vy1*vy1)*b2*b2
+            
+                  dx=xi-(eval wx h)
+                  dy=yi-(eval wy h)
+            
+                  x1'=x1 + (vx0*dx + vy0*dy)*b1
+                  y1'=y1 + (vx1*dx + vy1*dy)*b2
+              in
+               sumAll a' b' c' d' x1' y1' ss
+         
+            sumAll a b c d x1 y1 []=(a,b,c,d,x1,y1)
+                                    
+            (ra,rb,rc,rd,rx1,ry1)=sumAll 0 0 0 0 0 0 points
+            
+            (Matrix2 a_ b_ c_ d_)=inverse $ Matrix2 ra rb rc rd
+            lambda1=a_*rx1+b_*ry1
+            lambda2=c_*rx1+d_*ry1
+      
+            -- On a la courbe optimale. Il faut chercher ou on va couper, maintenant
+            xap=Bernsteinp 4 $ UV.fromList [wx0,
+                                            wx0+lambda1*vx0,
+                                            wx1+lambda2*vx1,
+                                            wx1]
+            yap=Bernsteinp 4 $ UV.fromList [wy0,
+                                            wy0+lambda1*vy0,
+                                            wy1+lambda2*vy1,
+                                            wy1]
+            (err,arg)=foldl (\m (h1,h2)->
+                              let (xi,yi)=evalCurve off (Interval h1 h2)
+                                  xj=eval xap (Interval h1 h2)
+                                  yj=eval yap (Interval h1 h2)
+                              in
+                               max m (iup $ abs $ (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj), (h1+h2)/2))
+                      (0,0) points
+        in
+         if err<=0.01 then
+           (desintervalize xap,desintervalize yap):(approx (rnd wx1) (rnd wy1) s)
+         else
+           approx x0 y0 $
+            (off { cx=restriction (cx off) 0 arg,
+                   cy=restriction (cy off) 0 arg }):
+            (off { cx=restriction (cx off) arg 1,
+                   cy=restriction (cy off) arg 1 }):s
+           
+      (x0h,y0h)=evalCurve h0 $ Interval (t0 h0) (t0 h0)
+      
+  in
+   approx (rnd x0h) (rnd y0h) l0
+           
+desintervalize::(Bernsteinp a Interval)->(Bernsteinp a Double)
+desintervalize b=b { coefs=UV.map rnd $ coefs b}
+  
+\end{code}
diff --git a/Graphics/Typography/Geometry/Bezier.lhs b/Graphics/Typography/Geometry/Bezier.lhs
new file mode 100644
--- /dev/null
+++ b/Graphics/Typography/Geometry/Bezier.lhs
@@ -0,0 +1,749 @@
+\documentclass{article}
+%include lhs2TeX.fmt
+\begin{document}
+\begin{code}
+{-# OPTIONS -XUnboxedTuples -XBangPatterns -XNamedFieldPuns -XRecordWildCards -XMagicHash -cpp #-}
+-- | This module contains the basic functions for manipulating Bezier curves. It is heavily
+-- based on the book by N. M. Patrikalakis and T. Maekawa, Shape Interrogation for Computer
+-- Aided Design and Manufacturing.
+
+module Graphics.Typography.Geometry.Bezier (
+  Curve(..),line,bezier3,
+  offset,
+  inter,
+  evalCurve,distance,
+  left,bottom,right,top)  where
+
+import Algebra.Polynomials.Bernstein
+import Algebra.Polynomials.Numerical
+import Graphics.Typography.Geometry
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as UV
+import Data.List (partition,sort)
+
+-- | The type for representing all types of curves.
+data Curve=
+  Bezier { cx::Bernsteinp Int Double,
+           cy::Bernsteinp Int Double, 
+           t0::Double,
+           t1::Double }
+  
+  | Offset { cx::Bernsteinp Int Double,
+             cy::Bernsteinp Int Double,
+             t0::Double,
+             t1::Double,
+             matrix::Matrix2 Double
+           }
+  | Circle { cx0::Double,
+             cy0::Double,
+             t0::Double,
+             t1::Double,
+             matrix::Matrix2 Double
+           }
+  deriving (Show)
+
+
+-- | The basic constructor for lines : a line is a degree 1 Bezier curve
+line::Double->Double->Double->Double->Curve
+line px py px' py'=Bezier { cx=Bernsteinp 2 $ UV.fromList [px,px'], 
+                            cy=Bernsteinp 2 $ UV.fromList [py,py'],
+                            t0=0,t1=1 }
+
+-- | A shortcut to define degree 3 Bezier curves from points. If the control
+-- points are @a,b,c,d@, the function should be called with
+-- @'bezier3' xa ya xb yb xc yc xd yd@.
+bezier3::Double->Double->Double->Double->Double->Double->Double->Double->Curve
+bezier3 px0 py0 px1 py1 px2 py2 px3 py3=
+  Bezier { cx=Bernsteinp 4 $ UV.fromList [px0,px1,px2,px3],
+           cy=Bernsteinp 4 $ UV.fromList [py0,py1,py2,py3],
+           t0=0,t1=1 }
+
+instance Geometric Curve where
+  translate x y cur@(Circle{cx0,cy0})=        
+    cur { cx0=cx0+x,cy0=cy0+y }
+  translate x y cur=
+    cur { cx=(cx cur) { coefs=UV.map (+x) $ coefs $ cx cur},
+          cy=(cy cur) { coefs=UV.map (+y) $ coefs $ cy cur} }
+
+  apply m0@(Matrix2 a b c d) cir@(Circle{cx0,cy0,matrix})=
+    cir { cx0=a*cx0+b*cy0, cy0=c*cx0+d*cy0, matrix=m0*matrix }
+  apply (Matrix2 a b c d) cur=
+    cur { cx=(scale a $ cx cur)+(scale b $ cy cur),
+          cy=(scale c $ cx cur)+(scale d $ cy cur) }
+
+
+
+-- | Gives the point corresponding to the given value of the parameter
+evalCurve::Curve->Interval->(Interval,Interval)
+evalCurve (Offset{..}) t=
+  let ix=intervalize cx
+      iy=intervalize cy
+      xt=eval ix t
+      yt=eval iy t
+      m@(Matrix2 a b c d)=intervalize matrix
+      (Matrix2 a_ b_ c_ d_)=inverse m
+      xt0'=eval (derivate ix) t
+      yt0'=eval (derivate iy) t
+      xt'=a_*xt0' + b_*yt0'
+      yt'=c_*xt0' + d_*yt0'
+      dd=sqrt $ xt'*xt' + yt'*yt'
+  in
+   (xt+(a*yt'-b*xt')/dd, yt+(c*yt'-d*xt')/dd)
+
+evalCurve (Circle{..}) alpha=
+  let xx=cos alpha
+      yy=sin alpha
+      (Matrix2 a b c d)=intervalize matrix
+  in
+   (interval cx0+a*xx+b*yy, interval cy0+c*xx+d*yy)
+   
+evalCurve (Bezier{..}) t=
+  let ix=intervalize cx
+      iy=intervalize cy
+      xx=eval ix t
+      yy=eval iy t
+  in
+   (xx,yy)
+   
+data Topo=Dehors | SurLaLigne | Dedans deriving Eq
+
+-- | @'inter' c0 c1@ is a list of all possible points of intersection
+-- between curves @c0@ and @c1@ : if @(u,v,w,x)@ is returned by 'inter',
+-- then curve @c0@ may intersect with @c1@ between parameter values @u@
+-- and @v@, which corresponds to parameter values between @w@ and @x@ for
+-- @c1@. The implementation guarantees that all actual solutions are found,
+-- but possibly false solutions may also be returned.
+
+inter::Curve->Curve->[((Double,Double,Double,Double))]
+inter op@(Offset { cx=bxp_,cy=byp_,matrix=mp,t0=t0a,t1=t1a })
+  (Offset { cx=bxq_,cy=byq_,matrix=mq,t0=t0b,t1=t1b })=
+  
+  -- Attention : verifier si c'est la meme generatrice
+  let thrx=1e-5
+      solutions _ []=[]
+      solutions thr boxes@(_:_)=
+        let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2,eq3])) boxes 
+        
+            (correct,toRefine)=partition (\(u,v,_,_,_,_,_,_)->
+                              let (xu,yu)=evalCurve op (Interval u u)
+                                  (xv,yv)=evalCurve op (Interval v v)
+                              in
+                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
+        in
+         correct++(solutions (thr/2) toRefine)
+  in
+   map (\(u,v,w,x,_,_,_,_)->(u,v,w,x)) $ solutions 1e-2 $
+   [(t0a,t1a,t0b,t1b,0,1,0,1)::
+       (Double,Double,Double,Double,Double,Double,Double,Double)]
+
+     where
+       
+       
+       imp@(Matrix2 ap bp cp dp)=intervalize mp
+       imq@(Matrix2 aq bq cq dq)=intervalize mq
+       (Matrix2 ap_ bp_ cp_ dp_)=inverse imp
+       (Matrix2 aq_ bq_ cq_ dq_)=inverse imq
+
+       bxp=intervalize bxp_
+       byp=intervalize byp_
+       bxq=intervalize bxq_
+       byq=intervalize byq_
+    
+       bxp4=promote 1 bxp
+       byp4=promote 1 byp
+       bxq4=promote 2 bxq
+       byq4=promote 2 byq
+    
+
+       bxp'=derivate bxp
+       byp'=derivate byp
+       bxq'=derivate bxq
+       byq'=derivate byq
+
+       bXp'=promote 1 $ (scale ap_ bxp')+(scale bp_ byp')
+       bYp'=promote 1 $ (scale cp_ bxp')+(scale dp_ byp')
+
+       bXq'=promote 2 $ (scale aq_ bxq')+(scale bq_ byq')
+       bYq'=promote 2 $ (scale cq_ bxq')+(scale dq_ byq')
+
+       bomp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp') :: Bernsteinp (Int,Int,Int,Int) Interval
+       bomq@(Bernsteinp _ omegaq)=(bXq'*bXq')+(bYq'*bYq') :: Bernsteinp (Int,Int,Int,Int) Interval
+       
+       au=
+         let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
+         max 0 $ fpred au_
+       bu=
+         let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
+         fsucc bu_
+       av=
+         let av_=sqrt $ UV.minimum $ UV.map ilow omegaq in
+         max 0 $ fpred av_
+       bv=
+         let bv_=sqrt $ UV.maximum $ UV.map iup omegaq in
+         fsucc bv_
+    
+       alphau=Bernsteinp (1,1,2,1) $ UV.fromList [Interval au au,Interval bu bu]
+       alphav=Bernsteinp (1,1,1,2) $ UV.fromList [Interval av av,Interval bv bv]
+
+       eq0=
+         ((bxp4*alphau*alphav) + (scale ap $ bYp'*alphav) - (scale bp $ bXp'*alphav)
+          -(bxq4*alphau*alphav) - (scale aq $ bYq'*alphau) + (scale bq $ bXq'*alphau))
+       eq1=
+         ((byp4*alphau*alphav) + (scale cp $ bYp'*alphav) - (scale dp $ bXp'*alphav)
+          -(byq4*alphau*alphav) - (scale cq $ bYq'*alphau) + (scale dq $ bXq'*alphau))    
+       eq2=bomp-(alphau*alphau)
+       eq3=bomq-(alphav*alphav)
+      
+    
+
+
+inter b@(Circle{}) a@(Offset{})=
+  map (\(i,j,k,l)->(k,l,i,j)) $ inter a b
+
+inter o@(Offset { cx=bxp, cy=byp, matrix=mp })
+  cir@(Circle{cx0,cy0,matrix=mq})=
+
+  let ix=intervalize bxp
+      iy=intervalize byp
+      m@(Matrix2 a b c d)=intervalize mp
+      (Matrix2 a_ b_ c_ d_)=inverse m
+      x'=derivate ix
+      y'=derivate iy
+      xx'=(scale a_ x')+(scale b_ y')
+      yy'=(scale c_ x')+(scale d_ y')
+      omega@(Bernsteinp _ omegap)=xx'*xx'+yy'*yy'
+      au=
+        let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
+        max 0 $ fpred au_
+      bu=
+        let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
+        fsucc bu_
+      
+      alphau=Bernsteinp (1,2) $ UV.fromList [Interval au au,Interval bu bu]
+      
+      lambda=(promote 1 omega) - alphau*alphau
+      -- Avant multiplication par M_C^-1
+      xx0=(promote 1 $ ix-(intervalize $ constant cx0))*alphau
+          +(promote 1 $ scale a yy'-scale b xx')
+      yy0=(promote 1 $ iy-(intervalize $ constant cy0))*alphau
+          +(promote 1 $ scale c yy'-scale d xx') :: Bernsteinp (Int,Int) Interval
+    
+      (Matrix2 ac_ bc_ cc_ dc_)=inverse $ intervalize mq
+      xx1=(scale ac_ xx0)+(scale bc_ yy0)
+      yy1=(scale cc_ xx0)+(scale dc_ yy0)
+         
+      eqc=xx1*xx1+yy1*yy1-alphau*alphau
+      
+      thrx=1e-5
+      
+      solutions _ []=[]
+      solutions thr boxes@(_:_)=
+        let sol0=concatMap (solve thr (V.fromList [eqc,lambda])) boxes 
+        
+            (correct,toRefine)=partition (\(u,v,_,_)->
+                              let (xu,yu)=evalCurve o (Interval u u)
+                                  (xv,yv)=evalCurve o (Interval v v)
+                              in
+                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
+        in
+         correct++(solutions (thr/2) toRefine)
+         
+
+      -- Removing false positives by computing the distance to the center of
+      -- the circle (this is quite fast).
+      
+      removeFalse cl0 (h@(_,v,_,_):h'@(u',_,_,_):s)=
+        let u''=(v+u')/2
+            (xu,yu)=evalCurve o (Interval u'' u'')
+            Interval dl du=distance xu yu cir
+            cl1
+              | du<1 = Dedans
+              | dl>1 = Dehors
+              | otherwise = SurLaLigne
+        in
+         if cl0/=cl1 then h:(removeFalse cl1 (h':s)) else
+           removeFalse cl1 (h':s)
+      removeFalse _ l=l
+      
+      initCl=
+        let (x0,y0)=evalCurve o (Interval (t0 o) (t0 o)) 
+            Interval dl du=distance x0 y0 cir
+        in
+        if dl>1 then Dehors else if du<1 then Dedans else SurLaLigne
+  in
+   foldl (\l (u,v,_,_)->
+              let (Interval xl xu,Interval yl yu)=evalCurve o (Interval u v) in
+              case angle (Interval xl xu) (Interval yl yu) cir of
+                Just (Interval a0 a1)->
+                  (u,v,a0,a1):l
+                Nothing->l
+            ) [] $ removeFalse initCl $ sort $ solutions 1e-2 [(t0 o,t1 o,0::Double,1::Double)]
+   
+  
+inter a@(Circle{cx0=x0a,cy0=y0a,matrix=ma})
+  b@(Circle{cx0=x0b,cy0=y0b,matrix=mb})=
+  
+  if (intervalize ma)`intersects`(intervalize mb) && x0a==x0b && y0a==y0b then
+    let up ix@(Interval _ x_) tt0 tt1
+          | x_<tt0 =
+            up (ix+(2*interval pi)) tt0 tt1
+          | otherwise = down ix tt0 tt1
+        down ix@(Interval x_ x__) tt0 tt1
+          | x_>tt1 =
+            down (ix-(2*interval pi)) tt0 tt1
+          | x__<tt0 =
+              Nothing
+          | otherwise =
+            Just ix
+            
+        alpha=up (interval $ t0 a) (t0 b) (t1 b)
+        beta=up (interval $ t0 b) (t0 b) (t1 b)
+    in
+     
+     case (alpha,beta) of
+       (Just aa,Just ab)->
+         case (up aa (t0 a) (t1 a),
+               up ab (t0 a) (t1 a)) of
+           
+           (Just ba,Just bb)
+             | ilow aa<=iup ab -> [(ilow ba, iup bb,
+                                    ilow aa, iup ab)]
+             | otherwise->
+               case (up (interval $ t0 b) (t0 a) (t1 a),
+                     up (interval $ t1 b) (t0 a) (t1 a)) of
+                 (Just b0,Just b1)->
+                   [(ilow b0,iup bb,
+                     t0 b, iup ab),
+                    (ilow ba,iup b1,
+                     ilow aa, t1 b)]
+                 _->[]
+           _->[]
+       _->[]
+
+  else
+    let thr=1e-5
+        solutions=solve thr (V.fromList [eq0,eq1]) (fpred u0,fsucc v0,
+                                                    fpred w0,fsucc x0)
+    in
+     foldl (\l (u,v,w,x)->
+             let alpha=angle (Interval u v) (Interval w x) a
+                 beta=angle (Interval u v) (Interval w x) b
+             in
+              case alpha of
+                Just (Interval a0l a0u)->
+                  case beta of
+                    Just (Interval b0l b0u)->(a0l,a0u,b0l,b0u):l
+                    _->l
+                _->l
+           ) [] solutions 
+  where
+    
+    ima@(Matrix2 am bm cm dm)=intervalize ma
+    
+    maxa=max (iup $ abs am+abs bm) (iup $ abs cm+abs dm)
+    (u0,v0,w0,x0)=(x0a-maxa,x0a+maxa,y0a-maxa,y0a+maxa)
+    
+    -- x-x0
+    xxa0=intervalize $ Bernsteinp (2,1) $ UV.fromList [-x0a,1-x0a] :: Bernsteinp (Int,Int) Interval
+    yya0=intervalize $ Bernsteinp (1,2) $ UV.fromList [-y0a,1-y0a] :: Bernsteinp (Int,Int) Interval
+    (Matrix2 aa_ ba_ ca_ da_)=inverse ima
+    xxa=(scale aa_ xxa0)+(scale ba_ yya0)::Bernsteinp (Int,Int) Interval
+    yya=(scale ca_ xxa0)+(scale da_ yya0)
+    
+    xxb0=intervalize $ Bernsteinp (2,1) $ UV.fromList [-x0b,1-x0b]
+    yyb0=intervalize $ Bernsteinp (1,2) $ UV.fromList [-y0b,1-y0b]
+    (Matrix2 ab_ bb_ cb_ db_)=inverse $ intervalize mb
+    xxb=(scale ab_ xxb0)+(scale bb_ yyb0)
+    yyb=(scale cb_ xxb0)+(scale db_ yyb0)
+    
+    c1=Bernsteinp (1,1) $ UV.singleton 1
+    
+    eq0=xxa*xxa+yya*yya-c1
+    eq1=xxb*xxb+yyb*yyb-c1
+
+inter op@(Bezier{cx=bxa,cy=bya,t0=t0a,t1=t1a}) (Bezier{cx=xb,cy=yb,t0=t0b,t1=t1b})=
+  
+  let p0=(promote 1 $ intervalize bxa)-(promote 2 $ intervalize xb) :: Bernsteinp (Int,Int) Interval
+      p1=(promote 1 $ intervalize bya)-(promote 2 $ intervalize yb) :: Bernsteinp (Int,Int) Interval
+      thrx=1e-2
+      solutions _ []=[]
+      solutions thr boxes@(_:_)=
+        let sol0=concatMap (solve thr (V.fromList [p0,p1])) boxes 
+        
+            (correct,toRefine)=partition (\(u,v,_,_)->
+                              let (xu,yu)=evalCurve op (Interval u u)
+                                  (xv,yv)=evalCurve op (Interval v v)
+                              in
+                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
+        in
+         correct++(solutions (thr/2) toRefine)
+  in
+   solutions 1e-2 [(t0a,t1a,t0b,t1b)]
+
+
+inter cir@(Circle{}) bez@(Bezier{})=map (\(u,v,w,x)->(w,x,u,v)) $ inter bez cir
+
+inter bez@(Bezier{}) cir@(Circle{})=
+  let xx=(intervalize $ cx bez)-(intervalize $ constant $ cx0 cir)
+      yy=(intervalize $ cy bez)-(intervalize $ constant $ cy0 cir)
+      (Matrix2 a b c d)=inverse $ intervalize $ matrix cir
+      xx0=scale a xx+scale b yy
+      yy0=scale c xx+scale d yy
+      
+      thrx=1e-5
+      
+      solutions _ []=[]
+      solutions thr boxes@(_:_)=
+        let sol0=concatMap (solve thr (V.singleton (xx0*xx0+yy0*yy0-(constant 1)))) boxes
+        
+            (correct,toRefine)=partition (\(u,v)->
+                              let (xu,yu)=evalCurve bez (Interval u u)
+                                  (xv,yv)=evalCurve bez (Interval v v)
+                              in
+                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
+        in
+         correct++(solutions (thr/2) toRefine)
+  in
+   foldl (\l (u,v)->
+              let (Interval xl xu,Interval yl yu)=evalCurve bez (Interval u v) in
+              case angle (Interval xl xu) (Interval yl yu) cir of
+                Just (Interval a0 a1)->
+                  (u,v,a0,a1):l
+                Nothing->l
+            ) [] $!
+   solutions (1e-2) [(t0 bez,t1 bez)]
+
+inter bez@(Bezier{}) off@(Offset{})=map (\(u,v,w,x)->(w,x,u,v)) $ inter off bez
+
+
+inter off@(Offset{}) bez@(Bezier{})=
+  
+  let thr=1e-2
+      thrx=1e-5
+      solutions _ []=[]
+      solutions thr boxes@(_:_)=
+        let sol0=concatMap (solve thr (V.fromList [eq0,eq1,eq2])) boxes
+        
+            (correct,toRefine)=partition (\(u,v,_,_,_,_)->
+                              let (xu,yu)=evalCurve off (Interval u u)
+                                  (xv,yv)=evalCurve off (Interval v v)
+                              in
+                               (iup $ (xu-xv)^(2::Int)+(yu-yv)^(2::Int))<=thrx) sol0
+        in
+         correct++(solutions (thr/2) toRefine)
+  in
+   map (\(a,b,c,d,_,_)->(a,b,c,d)) $ solutions 1e-2 $
+   [(0,1,0,1,0,1)::(Double,Double,Double,Double,Double,Double)]
+    where
+      
+      bxp=intervalize $ cx off
+      byp=intervalize $ cy off
+      bxq=intervalize $ cx bez
+      byq=intervalize $ cy bez
+      
+      bxp'=derivate bxp
+      byp'=derivate byp
+      bxp3=promote 1 bxp
+      byp3=promote 1 byp
+      bxq3=promote 2 bxq
+      byq3=promote 2 byq
+      
+      mp@(Matrix2 ap bp cp dp)=intervalize $ matrix $ off
+      (Matrix2 ap_ bp_ cp_ dp_)=inverse mp
+      bXp'=promote 1 $ (scale ap_ bxp')+(scale bp_ byp')
+      bYp'=promote 1 $ (scale cp_ bxp')+(scale dp_ byp')
+
+      omp@(Bernsteinp _ omegap)=(bXp'*bXp')+(bYp'*bYp')
+      au=
+        let au_=sqrt $ UV.minimum $ UV.map ilow omegap in
+        max 0 $ fpred au_
+      bu=
+        let bu_=sqrt $ UV.maximum $ UV.map iup omegap in
+        fsucc bu_
+    
+      alphau=Bernsteinp (1,1,2) $ UV.fromList [Interval au au,Interval bu bu]
+      eq0=bxp3*alphau + (scale ap bYp') - (scale bp bXp') - bxq3
+      eq1=byp3*alphau + (scale cp bYp') - (scale dp bXp') - byq3
+      eq2=alphau*alphau-omp
+      
+
+
+angle::Interval->Interval->Curve->Maybe Interval
+angle x y (Circle { cx0,cy0,matrix,t0,t1 })=
+  let vx=x-interval cx0
+      vy=y-interval cy0
+      
+      Matrix2 a b c d=inverse $ intervalize matrix
+      -- L'arithmetique d'intervalles fait un peu n'importe quoi
+      -- quand le vecteur est trop long. On le raccourcit.
+      alpha=
+        let co@(Interval col cou)=a*vx+b*vy
+            Interval sil siu=c*vx+d*vy
+            co2=
+              let (col2,cou2)=if col*col<cou*cou then (col*col,cou*cou) else
+                                (cou*cou,col*col)
+              in
+               Interval (fpred col2) (fsucc cou2)
+            si2=
+              let (sil2,siu2)=if sil*sil<siu*siu then (sil*sil,siu*siu) else
+                                (siu*siu,sil*sil)
+              in
+               Interval (fpred sil2) (fsucc siu2)
+            coco=co/(sqrt (co2+si2))
+            ac@(Interval acl acu)=acos $ Interval (max (-1) $ ilow coco) (min 1 $ iup coco)
+        in
+         if siu<0 then negate ac else
+           if sil>=0 then ac else
+             Interval (negate $ min (abs acl) (abs acu))
+             (max (abs acl) (abs acu))
+      up ix
+        | iup ix<t0 =
+          up $ ix+(2*interval pi)
+        | otherwise =
+            down ix
+      down ix
+        | ilow ix>t1 =
+          down $ ix-(2*interval pi)
+        | iup ix<t0 =
+          Nothing
+        | otherwise =
+          Just ix
+  in
+   up alpha
+
+
+angle _ _ _=error "angle"
+
+-- | Pseudo-distance from a point to a curve. Is the result is
+-- smaller than 1, the point is inside the curve. If it is greater
+-- than 1, the point is outside. Else we don't know (as usual with
+-- interval arithmetic).
+
+distance::Interval->Interval->Curve->Interval
+distance x0 y0 (Bezier{..})=
+  distance x0 y0 (Offset{cx,cy,t0,t1,matrix=Matrix2 1 0 0 1})
+  
+distance x0 y0 (Offset{..})=
+  let (Matrix2 a b c d)=inverse $ intervalize matrix
+      vx_=intervalize cx-(constant x0)
+      vy_=intervalize cy-(constant y0)
+      vx=scale a vx_+scale b vy_
+      vy=scale c vx_+scale d vy_
+      
+      dist=vx*vx+vy*vy
+  in
+   foldl (\di (u,v)->let di'=eval dist (Interval u v) in
+           if iup di<iup di' then di else di') (Interval (1/0) (1/0)) $
+   (t0,t0):(t1,t1):(solve 1e-5 (V.singleton (derivate dist)) (t0,t1))
+  
+  
+distance x1 y1 (Circle{..})=
+  let (Matrix2 a b c d)=inverse $ intervalize matrix
+      vx_=x1-Interval cx0 cx0
+      vy_=y1-Interval cy0 cy0
+      vx=a*vx_+b*vy_
+      vy=c*vx_+d*vy_
+  in
+   vx*vx+vy*vy
+   
+-- | Offsets a given Bezier curve with the given pen matrix. The original
+-- pen is a circle of radius one, the matrix, if inversible, is applied to it.
+
+offset::Matrix2 Double->Curve->[Curve]
+offset m (Bezier{cx=x@(Bernsteinp nx bx),cy=y@(Bernsteinp ny by)})=
+  if nx <=1 && ny <=1 then
+    [Circle { cx0=UV.head bx,cy0=UV.head by,t0=ilow 0,t1=iup $ 2*pi,matrix=m }]
+  else
+    [ c0,c1,c2,c3 ]
+    
+  where
+    im=intervalize m
+    (Matrix2 a_ b_ c_ d_)=inverse im
+    
+    ibx=intervalize x
+    iby=intervalize y
+    
+    lastCoef (Bernsteinp n c)
+      | n>=1 = UV.last c
+      | otherwise = 0
+    firstCoef (Bernsteinp n c)
+      | n>=1 = UV.head c
+      | otherwise = 0
+    
+    -- Premiere courbe offset
+    c0=Offset { cx=x, cy=y, t0=0,t1=1,matrix=m }
+    
+    -- Demi-cercle 1
+    
+    ibx'=derivate ibx
+    iby'=derivate iby
+    
+    -- Calcul du vecteur tangent au bout du premier
+
+    alpha0=
+      let xx0=lastCoef ibx'
+          yy0=lastCoef iby'
+
+          xx0'=a_*xx0+b_*yy0
+          yy0'=c_*xx0+d_*yy0
+          norm0=sqrt $ xx0'*xx0'+yy0'*yy0'
+          
+          xx'=xx0'/norm0
+          yy'=yy0'/norm0
+      in
+       if ilow xx'>=0 then
+         -(acos yy')
+       else
+         if iup xx'<=0 then
+           acos yy'
+         else
+           let Interval u v=acos yy' in
+           Interval (negate $ max (abs u) (abs v))
+           (max (abs u) (abs v))
+          
+                    
+    alpha0'=alpha0+interval pi
+    c1=Circle { cx0=lastCoef x,
+                cy0=lastCoef y,
+                t0=ilow alpha0,
+                t1=iup alpha0',
+                matrix=m }
+    
+    
+    -- Deuxieme courbe offset
+    c2=Offset { cx=reorient x,
+                cy=reorient y,
+                t0=0,t1=1,
+                matrix=m }
+       
+    -- Deuxieme demi-cercle
+    alpha1=
+      let xx0=firstCoef ibx'
+          yy0=firstCoef iby'
+
+          xx0'=a_*xx0+b_*yy0
+          yy0'=c_*xx0+d_*yy0
+          norm0=sqrt $ xx0'*xx0'+yy0'*yy0'
+          
+          xx'=xx0'/norm0
+          yy'=yy0'/norm0
+      in
+       if ilow xx'>=0 then
+         -(acos yy')
+       else
+         if iup xx'<=0 then
+           acos yy'
+         else
+           let Interval u v=acos yy' in
+           Interval (negate $ max (abs u) (abs v))
+           (max (abs u) (abs v))
+           
+    alpha1'=alpha1-pi
+    c3=Circle { cx0=firstCoef x, 
+                cy0=firstCoef y,
+                t0=ilow alpha1',
+                t1=iup alpha1,
+                matrix=m }
+    
+
+offset _ _=error "offset : undefined yet for other than Bezier"
+
+
+rnd::Interval->Double
+rnd (Interval a b)=(a+b)/2
+
+derivRoots::Double->Curve->([(Double,Double)],[(Double,Double)])
+derivRoots thr (Bezier{..})=
+  (solve thr (V.singleton $ derivate $ intervalize cx) (t0,t1),
+   solve thr (V.singleton $ derivate $ intervalize cy) (t0,t1))
+derivRoots thr (Offset{..})=
+  let ix=intervalize cx
+      iy=intervalize cy
+      x'=derivate ix
+      y'=derivate iy
+      m@(Matrix2 a b c d)=intervalize matrix
+      (Matrix2 a_ b_ c_ d_)=inverse m
+      
+      xx'=(scale a_ x')+(scale b_ y')
+      yy'=(scale c_ x')+(scale d_ y')
+      xx''=derivate xx'
+      yy''=derivate yy'
+      
+      omega=xx'*xx'+yy'*yy'
+      au=
+        let au_=sqrt $ UV.minimum $ UV.map ilow $ coefs omega in
+        max 0 $ fpred au_
+      bu=
+        let bu_=sqrt $ UV.maximum $ UV.map iup $ coefs omega in
+        fsucc bu_
+      alphau=Bernsteinp (1,2) $ UV.fromList [Interval au au,Interval bu bu]
+      
+      eqx0=(promote 1 yy'')*(alphau^(2::Int))-(promote 1 $ yy'*yy''+xx'*xx'')
+      eqy0=(promote 1 $ yy'*yy''+xx'*xx'')-(promote 1 xx'')*(alphau^(2::Int))
+      
+      eqx=(promote 1 x')*(alphau^(3::Int))+(scale a eqx0)+(scale b eqy0)
+          :: Bernsteinp (Int,Int) Interval
+      eqy=(promote 1 y')*(alphau^(3::Int))+(scale c eqx0)+(scale d eqy0)
+      
+      eq=(promote 1 omega)-alphau^(2::Int) :: Bernsteinp (Int,Int) Interval
+  in
+   (map (\(u,v,_,_)->(u,v)) $ solve thr (V.fromList [eqx,eq])
+    ((t0,t1,0,1)::(Double,Double,Double,Double)),
+    map (\(u,v,_,_)->(u,v)) $ solve thr (V.fromList [eqy,eq])
+    ((t0,t1,0,1)::(Double,Double,Double,Double)))
+derivRoots _ (Circle{..})=
+  let (Matrix2 a b c d)=intervalize matrix
+      aa=sqrt $ a*a+b*b
+      cc=sqrt $ c*c+d*d
+      sx
+        | ilow a>=0 = acos $ b/aa
+        | otherwise = negate $ acos $ b/aa
+      sy
+        | ilow c>=0 = acos $ d/cc
+        | otherwise = negate $ acos $ d/cc
+      Interval ux vx=(-sx)-pi/2
+      Interval uy vy=(-sy)-pi/2
+  in
+   ([(ux,vx)],[(uy,vy)])
+
+-- | The leftmost point on a curve
+left::Curve->(Double,Double)
+left cur=
+  let (x,y)=
+        foldl (\m@(xx,_) (s,t)->
+                let m'@(xx',_)=evalCurve cur $ interval $ (s+t)/2 in
+                if ilow xx<ilow xx' then m else m') (1/0,1/0) $
+        (t0 cur,t0 cur):(t1 cur,t1 cur):(fst $ derivRoots 1e-5 cur)
+  in
+   (rnd x,rnd y)
+-- | The bottommost point on a curve
+bottom::Curve->(Double,Double)
+bottom cur=
+  let (x,y)=
+        foldl (\m@(_,yy) (s,t)->
+                let m'@(_,yy')=evalCurve cur $ interval $ (s+t)/2 in
+                if ilow yy<ilow yy' then m else m') (1/0,1/0) $
+        (t0 cur,t0 cur):(t1 cur,t1 cur):(snd $ derivRoots 1e-5 cur)
+  in
+   (rnd x,rnd y)
+-- | The rightmost point on a curve
+right::Curve->(Double,Double)
+right cur=
+  let (x,y)=
+        foldl (\m@(xx,_) (s,t)->
+                let m'@(xx',_)=evalCurve cur $ interval $ (s+t)/2 in
+                if iup xx>iup xx' then m else m') (-1/0,-1/0) $
+        (t0 cur,t0 cur):(t1 cur,t1 cur):(fst $ derivRoots 1e-5 cur)
+  in
+   (rnd x,rnd y)
+-- | The topmost point on a curve
+top::Curve->(Double,Double)
+top cur=
+  let (x,y)=
+        foldl (\m@(_,yy) (s,t)->
+                let m'@(_,yy')=evalCurve cur $ interval $ (s+t)/2 in
+                if iup yy>iup yy' then m else m') (-1/0,-1/0) $
+        (t0 cur,t0 cur):(t1 cur,t1 cur):(snd $ derivRoots 1e-5 cur)
+  in
+   (rnd x,rnd y)
+\end{code}
diff --git a/Graphics/Typography/Geometry/Outlines.lhs b/Graphics/Typography/Geometry/Outlines.lhs
new file mode 100644
--- /dev/null
+++ b/Graphics/Typography/Geometry/Outlines.lhs
@@ -0,0 +1,291 @@
+\begin{code}
+{-# OPTIONS -XUnboxedTuples -cpp -XRecordWildCards -XNamedFieldPuns -XBangPatterns -XMagicHash -XScopedTypeVariables #-}
+-- | This module contains the necessary calls to the other modules of Metafont'
+-- to compute the outlines of a given number of pen strokes. The normal way of
+-- using it is by calling 'outlines'. One other possible way would be :
+--
+-- @
+-- let curves=cutAll curvesList in
+-- remerge $ contour curves $ intersections curves
+-- @
+
+module Graphics.Typography.Geometry.Outlines (cutAll, intersections, contour, remerge, outlines) where
+
+import Algebra.Polynomials.Bernstein
+import Algebra.Polynomials.Numerical 
+import Graphics.Typography.Geometry.Bezier
+import Graphics.Typography.Geometry
+import Data.List (sort)
+import qualified Data.Map as M
+import qualified Data.Vector as V
+
+import Control.Parallel
+
+(!)::V.Vector a->Int->a
+(!)=(V.!)
+
+-- | Cuts a curve into a list of consecutive non-selfintersecting curves.
+cutNoSelf::Curve->[Curve]
+cutNoSelf c@(Circle{})=[c]
+cutNoSelf bez@(Bezier{..})=
+  let ix=intervalize cx
+      dx=derivate ix
+      solutions=
+        sort $ filter (\(s,t)->(ilow $ eval ix (Interval s s))*
+                               (iup $ eval ix (Interval t t)) <= 0) $
+        solve 1e-10 (V.singleton dx) (t0,t1)
+      roots lastU []=
+        if lastU>=t1 then
+          []
+        else
+          [bez { t0=lastU }]
+      roots lastU (u:s)
+        | u<=lastU = roots lastU s -- on ne coupe pas au debut
+        | otherwise =
+          (bez { t0=lastU, t1=u }):
+          (roots u s)
+  in
+   roots t0 $ map (\(s,t)->(s+t)/2) solutions
+  
+cutNoSelf off@(Offset{..})= -- offset
+  let thr=1e-2
+      ix=intervalize cx
+      iy=intervalize cy
+      x'=derivate ix
+      y'=derivate iy
+      (Matrix2 a b c d)=intervalize matrix
+      (Matrix2 a_ b_ c_ d_)=inverse $ intervalize matrix
+                            
+      xx'=(scale a_ x')+(scale b_ y')
+      yy'=(scale c_ x')+(scale d_ y')
+          
+      xx''=derivate xx'
+      yy''=derivate yy'
+      
+      evalC (t::Interval)=
+        let norm=sqrt $ (eval xx' t)*(eval xx' t)+(eval yy' t)*(eval yy' t)
+            derx=(eval yy'' t)/norm - 
+                 ((eval yy' t)*((eval xx' t)*(eval xx'' t)+
+                                (eval yy' t)*(eval yy'' t)))/(norm*norm*norm)
+            dery=(eval xx'' t)/norm - 
+                 ((eval xx' t)*((eval xx' t)*(eval xx'' t)+
+                                (eval yy' t)*(eval yy'' t)))/(norm*norm*norm)
+        in
+         ((eval x' t)+(a*derx-b*dery), (eval y' t)+(c*derx-d*dery))
+      
+      zerosx=
+        let verif t lastxx
+              | t>=t1 = []
+              | otherwise =
+                let (xx,_)=evalC (Interval t t) in
+                if (iup $ xx*lastxx)<=0 then
+                  t:verif (t+thr) xx
+                else
+                  verif (t+thr) xx
+                  
+                  
+                  
+        in
+         verif t0 $ fst $ evalC (Interval t0 t0)
+             
+      roots lastU []=
+        if lastU>=t1 then
+          []
+        else
+          [off { t0=lastU }]
+      roots lastU (u:s)
+        | u<=lastU = roots lastU s -- on ne coupe pas au debut
+        | otherwise =
+          (off { t0=lastU, t1=u }):
+          (roots u s)
+  in
+   roots t0 zerosx
+
+-- | @'cutAll' curves@ is the array of all the curves, cut such that
+-- each part does not intersect itself.
+cutAll::[[Curve]]->V.Vector (V.Vector Curve)
+cutAll l=V.fromList $ map (\c->V.fromList $ concatMap cutNoSelf c) l
+
+
+data Topology=Dedans | SurLaLigne | Dehors deriving (Eq, Ord, Show)
+
+minsert::Ord a=>a->b->M.Map a [b]->M.Map a [b]
+minsert x y m=M.insertWith' (++) x [y] m
+
+munion::Ord a=>M.Map a [b]->M.Map a [b]->M.Map a [b]
+munion=M.unionWith (++)
+  
+  
+mdeleteFindMin::Ord a=>M.Map a [b]->(Maybe (a,b),M.Map a [b])
+mdeleteFindMin m=
+  if M.null m then
+    (Nothing, m)
+  else
+    let ((a,b),m')=M.deleteFindMin m in
+    case b of
+      []->mdeleteFindMin m'
+      (h:s)->(Just (a,h), if null s then m' else M.insert a s m')
+
+
+-- | Computes the intersections between any pair of curves given
+-- as input, in parallel in GHC using @+RTS -N@.
+intersections::V.Vector (V.Vector Curve)->
+               M.Map (Int,Int,Double) [(Int,Int,Double,Double)]
+intersections curves=
+  let interAll ci cj
+        | ci>=V.length curves = M.empty
+        | cj>=V.length curves = interAll (ci+1) (ci+1)
+        | otherwise = 
+            -- traceShow (ci,i,cj,j) $
+          let next=interAll ci (cj+1)
+              inters
+                | ci==cj =
+                  V.ifoldl'
+                  (\s0 i curvei->
+                    V.ifoldl' 
+                    (\s1 j curvej->
+                      foldl (\s2 (ti,ti',tj,tj')->
+                              minsert (ci,i,ti) (cj,j+i+1,tj,tj') $
+                              minsert (cj,j+i+1,tj) (ci,i,ti,ti') $ s2) s1 $
+                      inter curvei curvej
+                    )
+                    s0 $ V.drop (i+1) (curves!cj)
+                  ) M.empty $ V.take (V.length (curves!ci)-1) (curves!ci)
+                | otherwise = 
+                  V.ifoldl'
+                  (\s0 i curvei->
+                    V.ifoldl'
+                    (\s1 j curvej->
+                      foldl (\s2 (ti,ti',tj,tj')->
+                              minsert (ci,i,ti) (cj,j,tj,tj') $
+                              minsert (ci,i,ti') (cj,j,tj,tj') $
+                              minsert (cj,j,tj) (ci,i,ti,ti') $
+                              minsert (cj,j,tj') (ci,i,ti,ti') $ s2) s1 $
+                      inter curvei curvej
+                    )
+                    s0 (curves!cj)
+                  ) M.empty $ V.take (V.length (curves!ci)-1) (curves!ci)
+          in
+           (next`par`inters)`seq`
+           (next`munion`inters)
+  in
+   interAll 0 0
+   
+-- | 'contour' takes the curves and the intersections computed as in 'intersections',
+-- and outputs a list of all simple closed paths defined by the curves in the input.
+contour::V.Vector (V.Vector Curve)->
+         M.Map (Int,Int,Double) [(Int,Int,Double,Double)]->
+         [[(Int,Int,Double,Double)]]
+contour curves inters0=
+  
+  let allPaths inters1 passages1=
+        let (first,inters2)=mdeleteFindMin inters1 in
+        case first of
+          Nothing->[]
+          Just ((ci0,i0,ti0),(cj0,j0,tj0a,tj0b))->
+            --traceShow ("new path",pi0,pj0) $
+            let walk ci i tia tib inters passages=
+                  --traceShow ("point",ci,i,tia,tib) $ traceShow (inters) $
+                  let (a,b)=M.split (ci,i,tib) inters
+                      (next,b')=mdeleteFindMin b
+                  in
+                   case next of
+                     Nothing-> -- traceShow ("echec 1") $
+                              ([],a,passages)
+                     Just ((ci',i',ti'),(cj,j,tja,tjb))
+                       | ci==ci0 && i==i0 && (ci',i',ti')>=(ci,i,ti0)->
+                         -- fin du chemin
+                         ([(ci,i,tia,ti0)],a`munion`b',passages)
+                     
+                       | ci==ci' && i==i' ->
+                         let isVisible=
+                               let tt=(tia+ti')/2
+                                   (xi,yi)=evalCurve (curves!ci!i) (Interval tt tt) 
+                               in
+                                V.foldl (\vis cur->
+                                          vis && 
+                                          iup (distance xi yi $ (cur!0) {t0=0,t1=1})>=1)
+                                True curves
+                         in
+                          if (not isVisible) then
+                            --traceShow ("invisible",pi') $
+                            ([],a`munion`b',passages)
+                          else
+                           let alreadyPassed=
+                                 let (_,p1)=M.split (ci,i,ti') passages in
+                                 (not $ M.null p1) &&
+                                 (let ((ci_,i_,_),ti'_)=M.findMin p1 in
+                                   ci_==ci && i_==i && ti'_<=ti')
+                           in
+                            if alreadyPassed then
+                              --traceShow ("already passed",pi') $
+                              ([],a`munion`b',passages)
+                            else
+                              --traceShow ("trying",pi') $
+                              let (nextPath,nextInters,nextPassages)=
+                                    walk cj j tja tjb (a`munion`b') $
+                                    M.insert (ci,i,ti') tia passages
+                              in
+                               if null nextPath then
+                                 walk ci i tia tib (a`munion`b') passages
+                               else
+                                 ((ci,i,tia,ti'):nextPath,
+                                  nextInters,
+                                  M.insert (ci,i,ti') tia nextPassages)
+                       | otherwise -> --traceShow ("echec 2",ci',i',ti') $
+                           ([],inters,passages)
+                                      
+                (path,inters3,passages1')=walk cj0 j0 tj0a tj0b inters2 passages1
+            in                       
+             if null path then
+               --traceShow ("abandon") $
+               allPaths inters3 passages1'
+             else
+               --traceShow ("reussi") $
+               path:(allPaths inters3 passages1')
+  in
+   allPaths inters0 M.empty
+   
+-- | 'remerge' takes the curves, the output of 'contour', and outputs
+-- the list of "remerged" curves, i.e. where the parts free of self-intersections
+-- are glued back to each other.
+remerge::V.Vector (V.Vector Curve)->[(Int,Int,Double,Double)]->[Curve]
+remerge _ []=[]
+remerge curves [(ci,i,ti0,ti1)]=[(curves!ci!i) { t0=ti0,t1=ti1 }]
+remerge curves (l@((ci,i,ti0,_):s))=
+  
+  let (cj,j,_,tj1)=last s in
+  if ci==cj && j+1==i && tj1==ti0 then
+    -- dans ce cas, le dernier est colle au premier
+    let takeFirsts []=(# [],[] #)
+        takeFirsts ((h@(ci',_,_,_)):ss)
+          | ci'==ci = 
+            let (# u,v #)=takeFirsts ss in
+            (# h:u, v #)
+          | otherwise = (# [],h:ss #)
+        (# uu,vv #)=takeFirsts l
+    in
+     remerge_ $ vv++uu
+  else
+    remerge_ l
+  
+  where
+    remerge_ []=[]
+    remerge_ [(cj,j,tj0,tj1)]=[(curves!cj!j) { t0=tj0,t1=tj1 }]
+    remerge_ ((cj,j,tj0,tj1):(cck@(ck,k,tk0,_)):ss)
+      | cj==ck && k==j+1 && tj1==tk0 =
+        let (h':s')=remerge_ $ cck:ss in
+        (h' { t0=tj0 }) : s'
+        
+      | otherwise = 
+          ((curves!cj!j) { t0=tj0,t1=tj1 }) : (remerge_ $ cck:ss)
+    
+
+-- | Takes a list of curves, potentially offset, and outputs the relevants part
+-- of the outlines.
+outlines::[[Curve]]->[[Curve]]
+outlines curves=
+  let curves'=cutAll curves in
+  map (remerge curves') $ contour curves' $ intersections curves'
+
+\end{code}
diff --git a/Graphics/Typography/Outlines.lhs b/Graphics/Typography/Outlines.lhs
deleted file mode 100644
--- a/Graphics/Typography/Outlines.lhs
+++ /dev/null
@@ -1,291 +0,0 @@
-\begin{code}
-{-# OPTIONS -XUnboxedTuples -cpp -XRecordWildCards -XNamedFieldPuns -XBangPatterns -XMagicHash -XScopedTypeVariables #-}
--- | This module contains the necessary calls to the other modules of Metafont'
--- to compute the outlines of a given number of pen strokes. The normal way of
--- using it is by calling 'outlines'. One other possible way would be :
---
--- @
--- let curves=cutAll curvesList in
--- remerge $ contour curves $ intersections curves
--- @
-
-module Graphics.Typography.Outlines (cutAll, intersections, contour, remerge, outlines) where
-
-import Algebra.Polynomials.Bernstein
-import Algebra.Polynomials.Numerical 
-import Graphics.Typography.Bezier
-import Graphics.Typography
-import Data.List (sort)
-import qualified Data.Map as M
-import qualified Data.Vector as V
-
-import Control.Parallel
-
-(!)::V.Vector a->Int->a
-(!)=(V.!)
-
--- | Cuts a curve into a list of consecutive non-selfintersecting curves.
-cutNoSelf::Curve->[Curve]
-cutNoSelf c@(Circle{})=[c]
-cutNoSelf bez@(Bezier{..})=
-  let ix=intervalize cx
-      dx=derivate ix
-      solutions=
-        sort $ filter (\(s,t)->(ilow $ eval ix (Interval s s))*
-                               (iup $ eval ix (Interval t t)) <= 0) $
-        solve 1e-10 (V.singleton dx) (t0,t1)
-      roots lastU []=
-        if lastU>=t1 then
-          []
-        else
-          [bez { t0=lastU }]
-      roots lastU (u:s)
-        | u<=lastU = roots lastU s -- on ne coupe pas au debut
-        | otherwise =
-          (bez { t0=lastU, t1=u }):
-          (roots u s)
-  in
-   roots t0 $ map (\(s,t)->(s+t)/2) solutions
-  
-cutNoSelf off@(Offset{..})= -- offset
-  let thr=1e-2
-      ix=intervalize cx
-      iy=intervalize cy
-      x'=derivate ix
-      y'=derivate iy
-      (Matrix2 a b c d)=intervalize matrix
-      (Matrix2 a_ b_ c_ d_)=inverse $ intervalize matrix
-                            
-      xx'=(scale a_ x')+(scale b_ y')
-      yy'=(scale c_ x')+(scale d_ y')
-          
-      xx''=derivate xx'
-      yy''=derivate yy'
-      
-      evalC (t::Interval)=
-        let norm=sqrt $ (eval xx' t)*(eval xx' t)+(eval yy' t)*(eval yy' t)
-            derx=(eval yy'' t)/norm - 
-                 ((eval yy' t)*((eval xx' t)*(eval xx'' t)+
-                                (eval yy' t)*(eval yy'' t)))/(norm*norm*norm)
-            dery=(eval xx'' t)/norm - 
-                 ((eval xx' t)*((eval xx' t)*(eval xx'' t)+
-                                (eval yy' t)*(eval yy'' t)))/(norm*norm*norm)
-        in
-         ((eval x' t)+(a*derx-b*dery), (eval y' t)+(c*derx-d*dery))
-      
-      zerosx=
-        let verif t lastxx
-              | t>=t1 = []
-              | otherwise =
-                let (xx,_)=evalC (Interval t t) in
-                if (iup $ xx*lastxx)<=0 then
-                  t:verif (t+thr) xx
-                else
-                  verif (t+thr) xx
-                  
-                  
-                  
-        in
-         verif t0 $ fst $ evalC (Interval t0 t0)
-             
-      roots lastU []=
-        if lastU>=t1 then
-          []
-        else
-          [off { t0=lastU }]
-      roots lastU (u:s)
-        | u<=lastU = roots lastU s -- on ne coupe pas au debut
-        | otherwise =
-          (off { t0=lastU, t1=u }):
-          (roots u s)
-  in
-   roots t0 zerosx
-
--- | @'cutAll' curves@ is the array of all the curves, cut such that
--- each part does not intersect itself.
-cutAll::[[Curve]]->V.Vector (V.Vector Curve)
-cutAll l=V.fromList $ map (\c->V.fromList $ concatMap cutNoSelf c) l
-
-
-data Topology=Dedans | SurLaLigne | Dehors deriving (Eq, Ord, Show)
-
-minsert::Ord a=>a->b->M.Map a [b]->M.Map a [b]
-minsert x y m=M.insertWith' (++) x [y] m
-
-munion::Ord a=>M.Map a [b]->M.Map a [b]->M.Map a [b]
-munion=M.unionWith (++)
-  
-  
-mdeleteFindMin::Ord a=>M.Map a [b]->(Maybe (a,b),M.Map a [b])
-mdeleteFindMin m=
-  if M.null m then
-    (Nothing, m)
-  else
-    let ((a,b),m')=M.deleteFindMin m in
-    case b of
-      []->mdeleteFindMin m'
-      (h:s)->(Just (a,h), if null s then m' else M.insert a s m')
-
-
--- | Computes the intersections between any pair of curves given
--- as input, in parallel in GHC using @+RTS -N@.
-intersections::V.Vector (V.Vector Curve)->
-               M.Map (Int,Int,Double) [(Int,Int,Double,Double)]
-intersections curves=
-  let interAll ci cj
-        | ci>=V.length curves = M.empty
-        | cj>=V.length curves = interAll (ci+1) (ci+1)
-        | otherwise = 
-            -- traceShow (ci,i,cj,j) $
-          let next=interAll ci (cj+1)
-              inters
-                | ci==cj =
-                  V.ifoldl'
-                  (\s0 i curvei->
-                    V.ifoldl' 
-                    (\s1 j curvej->
-                      foldl (\s2 (ti,ti',tj,tj')->
-                              minsert (ci,i,ti) (cj,j+i+1,tj,tj') $
-                              minsert (cj,j+i+1,tj) (ci,i,ti,ti') $ s2) s1 $
-                      inter curvei curvej
-                    )
-                    s0 $ V.drop (i+1) (curves!cj)
-                  ) M.empty $ V.take (V.length (curves!ci)-1) (curves!ci)
-                | otherwise = 
-                  V.ifoldl'
-                  (\s0 i curvei->
-                    V.ifoldl'
-                    (\s1 j curvej->
-                      foldl (\s2 (ti,ti',tj,tj')->
-                              minsert (ci,i,ti) (cj,j,tj,tj') $
-                              minsert (ci,i,ti') (cj,j,tj,tj') $
-                              minsert (cj,j,tj) (ci,i,ti,ti') $
-                              minsert (cj,j,tj') (ci,i,ti,ti') $ s2) s1 $
-                      inter curvei curvej
-                    )
-                    s0 (curves!cj)
-                  ) M.empty $ V.take (V.length (curves!ci)-1) (curves!ci)
-          in
-           (next`par`inters)`seq`
-           (next`munion`inters)
-  in
-   interAll 0 0
-   
--- | 'contour' takes the curves and the intersections computed as in 'intersections',
--- and outputs a list of all simple closed paths defined by the curves in the input.
-contour::V.Vector (V.Vector Curve)->
-         M.Map (Int,Int,Double) [(Int,Int,Double,Double)]->
-         [[(Int,Int,Double,Double)]]
-contour curves inters0=
-  
-  let allPaths inters1 passages1=
-        let (first,inters2)=mdeleteFindMin inters1 in
-        case first of
-          Nothing->[]
-          Just ((ci0,i0,ti0),(cj0,j0,tj0a,tj0b))->
-            --traceShow ("new path",pi0,pj0) $
-            let walk ci i tia tib inters passages=
-                  --traceShow ("point",ci,i,tia,tib) $ traceShow (inters) $
-                  let (a,b)=M.split (ci,i,tib) inters
-                      (next,b')=mdeleteFindMin b
-                  in
-                   case next of
-                     Nothing-> -- traceShow ("echec 1") $
-                              ([],a,passages)
-                     Just ((ci',i',ti'),(cj,j,tja,tjb))
-                       | ci==ci0 && i==i0 && (ci',i',ti')>=(ci,i,ti0)->
-                         -- fin du chemin
-                         ([(ci,i,tia,ti0)],a`munion`b',passages)
-                     
-                       | ci==ci' && i==i' ->
-                         let isVisible=
-                               let tt=(tia+ti')/2
-                                   (xi,yi)=evalCurve (curves!ci!i) (Interval tt tt) 
-                               in
-                                V.foldl (\vis cur->
-                                          vis && 
-                                          iup (distance xi yi $ (cur!0) {t0=0,t1=1})>=1)
-                                True curves
-                         in
-                          if (not isVisible) then
-                            --traceShow ("invisible",pi') $
-                            ([],a`munion`b',passages)
-                          else
-                           let alreadyPassed=
-                                 let (_,p1)=M.split (ci,i,ti') passages in
-                                 (not $ M.null p1) &&
-                                 (let ((ci_,i_,_),ti'_)=M.findMin p1 in
-                                   ci_==ci && i_==i && ti'_<=ti')
-                           in
-                            if alreadyPassed then
-                              --traceShow ("already passed",pi') $
-                              ([],a`munion`b',passages)
-                            else
-                              --traceShow ("trying",pi') $
-                              let (nextPath,nextInters,nextPassages)=
-                                    walk cj j tja tjb (a`munion`b') $
-                                    M.insert (ci,i,ti') tia passages
-                              in
-                               if null nextPath then
-                                 walk ci i tia tib (a`munion`b') passages
-                               else
-                                 ((ci,i,tia,ti'):nextPath,
-                                  nextInters,
-                                  M.insert (ci,i,ti') tia nextPassages)
-                       | otherwise -> --traceShow ("echec 2",ci',i',ti') $
-                           ([],inters,passages)
-                                      
-                (path,inters3,passages1')=walk cj0 j0 tj0a tj0b inters2 passages1
-            in                       
-             if null path then
-               --traceShow ("abandon") $
-               allPaths inters3 passages1'
-             else
-               --traceShow ("reussi") $
-               path:(allPaths inters3 passages1')
-  in
-   allPaths inters0 M.empty
-   
--- | 'remerge' takes the curves, the output of 'contour', and outputs
--- the list of "remerged" curves, i.e. where the parts free of self-intersections
--- are glued back to each other.
-remerge::V.Vector (V.Vector Curve)->[(Int,Int,Double,Double)]->[Curve]
-remerge _ []=[]
-remerge curves [(ci,i,ti0,ti1)]=[(curves!ci!i) { t0=ti0,t1=ti1 }]
-remerge curves (l@((ci,i,ti0,_):s))=
-  
-  let (cj,j,_,tj1)=last s in
-  if ci==cj && j+1==i && tj1==ti0 then
-    -- dans ce cas, le dernier est colle au premier
-    let takeFirsts []=(# [],[] #)
-        takeFirsts ((h@(ci',_,_,_)):ss)
-          | ci'==ci = 
-            let (# u,v #)=takeFirsts ss in
-            (# h:u, v #)
-          | otherwise = (# [],h:ss #)
-        (# uu,vv #)=takeFirsts l
-    in
-     remerge_ $ vv++uu
-  else
-    remerge_ l
-  
-  where
-    remerge_ []=[]
-    remerge_ [(cj,j,tj0,tj1)]=[(curves!cj!j) { t0=tj0,t1=tj1 }]
-    remerge_ ((cj,j,tj0,tj1):(cck@(ck,k,tk0,_)):ss)
-      | cj==ck && k==j+1 && tj1==tk0 =
-        let (h':s')=remerge_ $ cck:ss in
-        (h' { t0=tj0 }) : s'
-        
-      | otherwise = 
-          ((curves!cj!j) { t0=tj0,t1=tj1 }) : (remerge_ $ cck:ss)
-    
-
--- | Takes a list of curves, potentially offset, and outputs the relevants part
--- of the outlines.
-outlines::[[Curve]]->[[Curve]]
-outlines curves=
-  let curves'=cutAll curves in
-  map (remerge curves') $ contour curves' $ intersections curves'
-
-\end{code}
diff --git a/typography-geometry.cabal b/typography-geometry.cabal
--- a/typography-geometry.cabal
+++ b/typography-geometry.cabal
@@ -1,5 +1,5 @@
 Name:		typography-geometry
-Version: 	1.0
+Version: 	1.0.0
 Synopsis:	Drawings for printed text documents
 Description: 	Drawings for printed text documents
 Category:	Typography
@@ -14,5 +14,6 @@
         tag: 1.0
 Library
         Build-Depends:	base<5, vector,polynomials-bernstein,containers,parallel
-        Exposed-modules: Graphics.Typography, Graphics.Typography.Bezier,
-                         Graphics.Typography.Approximation, Graphics.Typography.Outlines
+        Exposed-modules: Graphics.Typography.Geometry, Graphics.Typography.Geometry.Bezier,
+                         Graphics.Typography.Geometry.Approximation,
+                         Graphics.Typography.Geometry.Outlines