diff --git a/aern2-mfun.cabal b/aern2-mfun.cabal
new file mode 100644
--- /dev/null
+++ b/aern2-mfun.cabal
@@ -0,0 +1,91 @@
+cabal-version: 1.12
+
+-- This file has been generated from package.yaml by hpack version 0.34.4.
+--
+-- see: https://github.com/sol/hpack
+
+name:           aern2-mfun
+version:        0.2.9.0
+synopsis:       Multi-variate real function optimisation and proving
+description:    Please see the README on GitHub at <https://github.com/michalkonecny/aern2/#readme>
+category:       Math
+homepage:       https://github.com/michalkonecny/aern2#readme
+bug-reports:    https://github.com/michalkonecny/aern2/issues
+author:         Eike Neumann, Junaid Rasheed, Michal Konecny
+maintainer:     mikkonecny@gmail.com
+copyright:      2019-2022 Eike Neumann, Junaid Rasheed, Michal Konecny
+license:        BSD3
+build-type:     Simple
+
+source-repository head
+  type: git
+  location: https://github.com/michalkonecny/aern2
+
+library
+  exposed-modules:
+      AERN2.AD.Differential
+      AERN2.AD.GenericOperations
+      AERN2.AD.MPBallOperations
+      AERN2.AD.Type
+      AERN2.BoxFun.Box
+      AERN2.BoxFun.Optimisation
+      AERN2.BoxFun.TestFunctions
+      AERN2.BoxFun.Type
+      AERN2.Linear.Matrix.Inverse
+      AERN2.Linear.Matrix.Type
+      AERN2.Linear.Vector.Type
+      AERN2.Util.Util
+  other-modules:
+      Paths_aern2_mfun
+  hs-source-dirs:
+      src
+  default-extensions:
+      RebindableSyntax,
+      ScopedTypeVariables,
+      TypeFamilies,
+      TypeOperators,
+      ConstraintKinds,
+      DefaultSignatures,
+      MultiParamTypeClasses,
+      FlexibleContexts,
+      FlexibleInstances,
+      UndecidableInstances
+  other-extensions:
+      TemplateHaskell
+  ghc-options: -Wall
+  build-depends:
+      aern2-fun >=0.2.9
+    , aern2-mp >=0.2.9
+    , base ==4.*
+    , collect-errors >=0.1.5
+    , mixed-types-num >=0.5.9
+    , vector
+  default-language: Haskell2010
+
+executable aern2-mfun-benchmark-optimisation
+  main-is: bench/OptimisationBenchmark.hs
+  other-modules:
+      Paths_aern2_mfun
+  default-extensions:
+      RebindableSyntax,
+      ScopedTypeVariables,
+      TypeFamilies,
+      TypeOperators,
+      ConstraintKinds,
+      DefaultSignatures,
+      MultiParamTypeClasses,
+      FlexibleContexts,
+      FlexibleInstances,
+      UndecidableInstances
+  other-extensions:
+      TemplateHaskell
+  ghc-options: -threaded -rtsopts -with-rtsopts=-N -Wall -O2
+  build-depends:
+      aern2-fun >=0.2.9
+    , aern2-mfun
+    , aern2-mp >=0.2.9
+    , base ==4.*
+    , collect-errors >=0.1.5
+    , mixed-types-num >=0.5.9
+    , vector
+  default-language: Haskell2010
diff --git a/bench/OptimisationBenchmark.hs b/bench/OptimisationBenchmark.hs
new file mode 100644
--- /dev/null
+++ b/bench/OptimisationBenchmark.hs
@@ -0,0 +1,30 @@
+module Main where
+
+import MixedTypesNumPrelude
+import AERN2.BoxFun.TestFunctions
+import AERN2.BoxFun.Optimisation
+import AERN2.MP.Ball
+-- import AERN2.Real
+-- import AERN2.MP.Dyadic
+-- import AERN2.MP.Accuracy
+-- import AERN2.AD.Differential
+
+-- import qualified Data.List as List
+
+-- import AERN2.BoxFun.Type
+-- import AERN2.Linear.Vector.Type as V
+
+main :: IO ()
+main = 
+    let
+        ac = bits 100
+    in
+    do
+    putStrLn $ "ratz: "       ++ (show $ minFun ratz4        ac (prec 53))
+    putStrLn $ "griewank 2: " ++ (show $ minFun (griewank 2) ac (prec 53))
+    putStrLn $ "shekel 2: "   ++ (show $ minFun shekel       ac (prec 53))
+    putStrLn $ "himmelblau: " ++ (show $ minFun himmelblau   ac (prec 53))
+    putStrLn $ "rosenbrock: " ++ (show $ minFun rosenbrock   ac (prec 53))
+    putStrLn $ "trefethen: "  ++ (show $ minFun siam4        ac (prec 53))
+    putStrLn $ "griewank 5: " ++ (show $ minFun (griewank 5) ac (prec 53))
+    putStrLn $ "griewank 7: " ++ (show $ minFun (griewank 7) ac (prec 53))
diff --git a/src/AERN2/AD/Differential.hs b/src/AERN2/AD/Differential.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/AD/Differential.hs
@@ -0,0 +1,11 @@
+module AERN2.AD.Differential 
+(
+      module AERN2.AD.Type
+    , module AERN2.AD.GenericOperations
+    , module AERN2.AD.MPBallOperations
+)
+where
+
+import AERN2.AD.Type
+import AERN2.AD.GenericOperations
+import AERN2.AD.MPBallOperations
diff --git a/src/AERN2/AD/GenericOperations.hs b/src/AERN2/AD/GenericOperations.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/AD/GenericOperations.hs
@@ -0,0 +1,220 @@
+{-# OPTIONS_GHC -Wno-orphans #-}
+module AERN2.AD.GenericOperations where
+
+import MixedTypesNumPrelude
+import AERN2.AD.Type
+
+instance 
+    (CanDiv a a, CanSubSameType a, CanMulSameType a, CanMulBy a Integer, CanAddSameType a,
+    CanSubSameType (DivType a a)) 
+    =>
+    CanDiv (Differential a) (Differential a)
+    where
+    type DivType     (Differential a) (Differential a) = Differential (DivType a a)
+    divide a b =
+        case min (order a) (order b) of
+            2 -> 
+                let
+                    dtDiff    = (a_dxt * b_x - a_x * b_dxt)
+                    ySqrd     = b_x * b_x
+                    exp0      = a_d2x * b_x + a_dxt * b_dx - b_dxt * a_dx - a_x * b_d2x
+                in
+                OrderTwo  (a_x / b_x) ((a_dx * b_x - a_x * b_dx)/ySqrd) (dtDiff/ySqrd) 
+                          (exp0/ySqrd - (2*b_dx * dtDiff) / (ySqrd * b_x))
+            1 -> OrderOne  (a_x / b_x) ((a_dx * b_x - a_x * b_dx)/(b_x * b_x)) 
+            0 -> OrderZero (a_x / b_x)
+            _ -> error "illegal Differential order"
+        where
+        a_x = diff_x a
+        b_x = diff_x b
+        a_dx = diff_dx a
+        b_dx = diff_dx b
+        a_dxt = diff_dxt a
+        b_dxt = diff_dxt b
+        a_d2x = diff_d2x a
+        b_d2x = diff_d2x b
+
+instance 
+    (CanExpSameType a, CanMulSameType a, CanAddSameType a) =>
+    CanExp (Differential a)
+    where
+    type ExpType (Differential a) = Differential a
+    exp (OrderZero x)           = OrderZero (exp x)
+    exp (OrderOne x dx)         = OrderOne  (exp x) (dx * exp x)
+    exp (OrderTwo x dx dxt d2x) = OrderTwo  (exp x) (dx * exp x) (dxt * exp x) ((dxt * dx + d2x) * exp x)
+
+
+{- TODO: fix this in AERN2.MP -}
+clampedCos :: (CanSinCosSameType a, CanMinMaxSameType a, HasIntegers a) => a -> a
+clampedCos (x :: a) = max (convertExactly $ -1 :: a) $ min ((convertExactly 1) :: a) (cos x)
+
+clampedSin :: (CanSinCosSameType a, CanMinMaxSameType a, HasIntegers a) => a -> a
+clampedSin (x :: a) = max (convertExactly $ -1 :: a) $ min ((convertExactly 1) :: a) (sin x)
+
+instance 
+    (CanSinCosSameType a, CanMulSameType a, CanNegSameType a, CanSubSameType a, CanAddSameType a
+    , HasIntegers a, CanMinMaxSameType a) 
+    =>
+    CanSinCos (Differential a)
+    where
+    type SinCosType (Differential a) = Differential a
+    cos (OrderZero x)            = OrderZero (clampedCos x)
+    cos (OrderOne  x dx)         = OrderOne  (clampedCos x) (-dx * clampedSin x)
+    cos (OrderTwo  x dx dxt d2x) = OrderTwo  (clampedCos x) (-dx * clampedSin x) (-dxt * clampedSin x) 
+                                                (-dxt * dx * clampedCos x - d2x * clampedSin x)
+    
+    sin (OrderZero x)            = OrderZero (clampedSin x)
+    sin (OrderOne x dx)          = OrderOne  (clampedSin x) (dx * clampedCos x)
+    sin (OrderTwo x dx dxt d2x)  = OrderTwo  (clampedSin x) (dx * clampedCos x) (dxt * clampedCos x)
+                                                (d2x * clampedCos x - dxt * dx * clampedSin x)
+
+instance 
+    (CanAddSameType a) => 
+    CanAddAsymmetric (Differential a) (Differential a)
+    where
+    type AddType (Differential a) (Differential a) = Differential a
+    add a b =
+        case min (order a) (order b) of
+            2 -> OrderTwo  (a_x + b_x) (a_dx + b_dx) (a_dxt + b_dxt) (a_d2x + b_d2x)
+            1 -> OrderOne  (a_x + b_x) (a_dx + b_dx)
+            0 -> OrderZero (a_x + b_x)
+            _ -> error "illegal Differential order"
+        where
+        a_x = diff_x a
+        b_x = diff_x b
+        a_dx = diff_dx a
+        b_dx = diff_dx b
+        a_dxt = diff_dxt a
+        b_dxt = diff_dxt b
+        a_d2x = diff_d2x a
+        b_d2x = diff_d2x b
+        
+instance 
+    (CanMulSameType a, CanAddSameType a) =>
+    CanMulAsymmetric (Differential a) (Differential a)
+    where
+    type MulType (Differential a) (Differential a) = Differential a
+    mul a b =
+        case min (order a) (order b) of
+            2 -> OrderTwo  (a_x * b_x) (a_dx * b_x + a_x * b_dx) (a_x * b_dxt + b_x * a_dxt) 
+                           (a_x * b_d2x + a_dx * b_dxt + a_dxt * b_dx + a_d2x * b_x)
+            1 -> OrderOne  (a_x * b_x) (a_dx * b_x + a_x * b_dx)
+            0 -> OrderZero (a_x * b_x)
+            _ -> error "illegal Differential order"
+        where
+        a_x = diff_x a
+        b_x = diff_x b
+        a_dx = diff_dx a
+        b_dx = diff_dx b
+        a_dxt = diff_dxt a
+        b_dxt = diff_dxt b
+        a_d2x = diff_d2x a
+        b_d2x = diff_d2x b
+
+instance 
+    (CanMulSameType a, CanAddSameType a, CanPowBy a a, CanSubThis a Integer, CanLogSameType a, CanDivSameType a) =>
+    CanPow (Differential a) (Differential a)
+    where
+    type PowType        (Differential a) (Differential a) = (Differential a)
+    pow a b =
+        case min (order a) (order b) of
+            2 -> OrderOne  (a_x ^ b_x) ((a_x ^ (b_x - 1)) * (b_x * a_dx + a_x * (log a_x) * b_dx)) --FIXME: Add real OrderTwo definition here
+            -- 2 -> OrderTwo   (x a ^ x b) 
+            --                 ((x a ^ (x b - (convertExactly 1 :: a))) * (x b * dx a + x a * log (x a) * dx b)) 
+            --                 ((x a ^ (x b - (convertExactly 1 :: a))) * (x b * dxt a + x a * log (x a) * dxt b))
+            --                 (x a)
+                            
+            --                 where
+            --                     ta = convertExactly 2 :: a
+            1 -> OrderOne  (a_x ^ b_x) ((a_x ^ (b_x - 1)) * (b_x * a_dx + a_x * (log a_x) * b_dx))
+            0 -> OrderZero (a_x ^ b_x)
+            _ -> undefined
+        where
+        a_x = diff_x a
+        b_x = diff_x b
+        a_dx = diff_dx a
+        b_dx = diff_dx b
+        -- a_dxt = diff_dxt a
+        -- b_dxt = diff_dxt b
+        -- a_d2x = diff_d2x a
+        -- b_d2x = diff_d2x b
+
+instance 
+    (CanSubSameType a) =>
+    CanSub (Differential a) (Differential a)
+    where
+    type SubType (Differential a) (Differential a) = Differential a
+    sub a b =
+        case min (order a) (order b) of
+            2 -> OrderTwo  (a_x - b_x) (a_dx - b_dx) (a_dxt - b_dxt) (a_d2x - b_d2x)
+            1 -> OrderOne  (a_x - b_x) (a_dx - b_dx)
+            0 -> OrderZero (a_x - b_x)
+            _ -> error "illegal Differential order"
+        where
+        a_x = diff_x a
+        b_x = diff_x b
+        a_dx = diff_dx a
+        b_dx = diff_dx b
+        a_dxt = diff_dxt a
+        b_dxt = diff_dxt b
+        a_d2x = diff_d2x a
+        b_d2x = diff_d2x b
+
+instance 
+    (CanMulBy (Differential a) Integer) =>
+    CanNeg (Differential a)
+    where
+    type NegType (Differential a) = Differential a
+    negate x = (-1) * x
+
+instance 
+    (CanSqrtSameType a, CanMulSameType a, CanNegSameType a, CanAddSameType a, CanMulBy a Integer, CanRecipSameType a) 
+    =>
+    CanSqrt (Differential a)
+    where
+    type SqrtType (Differential a) = Differential a
+    sqrt (OrderZero x)             = OrderZero (sqrt x)
+    sqrt (OrderOne x dx)           = OrderOne  (sqrt x) (dx * sqrtx')
+                                        where sqrtx' = recip (2 * sqrt x)
+    sqrt (OrderTwo x dx dxt d2x)   = OrderTwo  (sqrt x) (dx * sqrtx') (dxt * sqrtx') 
+                                               ((d2x * sqrtx') + (dx * dxt * sqrtx'')) 
+                                        where 
+                                        sqrtx'  = recip (2 * sqrt x)
+                                        sqrtx'' = negate $ recip (4 * x * sqrt x)
+                                            -- sqrtx'  == 1 / (2 * sqrt(x))
+                                            -- sqrtx'' == -1 / (4 * x * sqrt(x)) == -1 / (4 * x^(3/2))
+
+instance
+    CanMinMaxSameType a =>
+    CanMinMaxAsymmetric (Differential a) (Differential a)
+    where
+    type MinMaxType (Differential a) (Differential a) = Differential a
+    min a b = OrderZero (min (diff_x a) (diff_x b))
+    max a b = OrderZero (max (diff_x a) (diff_x b))
+
+
+-- instance
+--     (CanMinMaxSameType a, HasIntegers a) =>
+--     CanMinMaxAsymmetric (Differential a) (Differential a)
+--     where
+--     type MinMaxType (Differential a) (Differential a) = Differential a
+--     min a b = 
+--         case min (order a) (order b) of
+--             2 -> OrderTwo  (min (a_x) (b_x)) (min (a_dx) (b_dx)) (min (a_dxt) (b_dxt)) 
+--                            (min (a_d2x) (b_d2x))
+--             1 -> OrderOne  (min (a_x) (b_x)) (min (a_dx) (b_dx))
+--             0 -> OrderZero (min (a_x) (b_x))
+--     max a b =
+--         case min (order a) (order b) of
+--             2 -> OrderTwo  (max (a_x) (b_x)) (max (a_dx) (b_dx)) (max (a_dxt) (b_dxt)) 
+--                            (max (a_d2x) (b_d2x))
+--             1 -> OrderOne  (max (a_x) (b_x)) (max (a_dx) (b_dx))
+--             0 -> OrderZero (max (a_x) (b_x))
+
+-- instance
+--     (CanAbsSameType a) =>
+--     CanAbs (Differential a)
+--     where
+--     type AbsType (Differential a) = Differential a
+--     abs = fmap abs
+
diff --git a/src/AERN2/AD/MPBallOperations.hs b/src/AERN2/AD/MPBallOperations.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/AD/MPBallOperations.hs
@@ -0,0 +1,223 @@
+{-# OPTIONS_GHC -Wno-orphans #-}
+module AERN2.AD.MPBallOperations where
+
+import MixedTypesNumPrelude
+
+import AERN2.AD.Type
+import AERN2.AD.GenericOperations ()
+import AERN2.MP.Ball
+import Numeric.CollectErrors.Type (noValueNumErrorCertain, NumError (NumError))
+
+{-- addition --}
+
+instance
+    CanAddAsymmetric MPBall (Differential (CN MPBall))
+    where
+    type AddType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)
+    add a b = add (differential (order b) (cn a)) b
+
+instance 
+    CanAddAsymmetric (Differential (CN MPBall)) MPBall
+    where
+    type AddType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)
+    add b a = add b (differential (order b) (cn a))
+
+instance 
+    CanAddAsymmetric Integer (Differential (CN MPBall))
+    where
+    type AddType Integer (Differential (CN MPBall)) = Differential (CN MPBall)
+    add a b = add (differential (order b) (cn $ mpBall a)) b
+
+instance 
+    CanAddAsymmetric (Differential (CN MPBall)) Integer
+    where
+    type AddType (Differential (CN MPBall)) Integer = Differential (CN MPBall)
+    add b a = add b (differential (order b) (cn $ mpBall a))
+
+instance
+    (CanBeMPBall a) =>
+    CanAddAsymmetric (CN a) (Differential (CN MPBall))
+    where
+    type AddType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)
+    add a b = add (differential (order b) (fmap mpBall a)) b
+
+instance 
+    (CanBeMPBall a) =>
+    CanAddAsymmetric (Differential (CN MPBall)) (CN a)
+    where
+    type AddType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)
+    add b a = add b (differential (order b) (fmap mpBall a))
+
+{-- subtraction --}
+
+instance 
+    CanSub Integer (Differential (CN MPBall))
+    where
+    type SubType Integer (Differential (CN MPBall)) = Differential (CN MPBall)
+    sub a b = sub (differential (order b) (cn $ mpBall a)) b
+
+instance 
+    CanSub (Differential (CN MPBall)) Integer
+    where
+    type SubType (Differential (CN MPBall)) Integer = Differential (CN MPBall)
+    sub b a = sub b (differential (order b) (cn $ mpBall a))
+
+instance 
+    CanSub MPBall (Differential (CN MPBall))
+    where
+    type SubType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)
+    sub a b = sub (differential (order b) (cn a)) b
+
+instance 
+    CanSub (Differential (CN MPBall)) MPBall
+    where
+    type SubType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)
+    sub b a = sub b (differential (order b) (cn a))
+
+instance 
+    (CanBeMPBall a) =>
+    CanSub (CN a) (Differential (CN MPBall))
+    where
+    type SubType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)
+    sub a b = sub (differential (order b) (fmap mpBall a)) b
+
+instance 
+    (CanBeMPBall a) =>
+    CanSub (Differential (CN MPBall)) (CN a) 
+    where
+    type SubType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)
+    sub b a = sub b (differential (order b) (fmap mpBall a))
+
+{-- multiplication --}
+
+instance 
+    CanMulAsymmetric Integer (Differential (CN MPBall))
+    where
+    type MulType Integer (Differential (CN MPBall)) = Differential (CN MPBall)
+    mul a b = mul (differential (order b) (cn $ mpBall a)) b
+
+instance 
+    CanMulAsymmetric (Differential (CN MPBall)) Integer
+        where
+        type MulType (Differential (CN MPBall)) Integer = Differential (CN MPBall)
+        mul b a = mul b (differential (order b) (cn $ mpBall a))
+
+instance 
+    CanMulAsymmetric MPBall (Differential (CN MPBall))
+    where
+    type MulType MPBall (Differential (CN MPBall)) = Differential (CN MPBall)
+    mul a b = mul (differential (order b) (cn a)) b
+
+instance 
+    CanMulAsymmetric (Differential (CN MPBall)) MPBall
+        where
+        type MulType (Differential (CN MPBall)) MPBall = Differential (CN MPBall)
+        mul b a = mul b (differential (order b) (cn a))
+
+{-- division --}
+
+instance 
+    CanDiv Integer (Differential (CN MPBall))
+    where
+    type DivType     Integer (Differential (CN MPBall)) = Differential (CN MPBall)
+    divide     a b = divide     (differential (order b) (cn $ mpBall a)) b
+    
+instance 
+    CanDiv (Differential (CN MPBall)) Integer
+    where
+    type DivType     (Differential (CN MPBall)) Integer = Differential (CN MPBall)
+    divide     b a = divide     b (differential (order b) (cn $ mpBall a))        
+
+instance 
+    CanDiv MPBall (Differential (CN MPBall))
+    where
+    type DivType     MPBall (Differential (CN MPBall)) = Differential (CN MPBall)
+    divide     a b = divide     (differential (order b) (cn a)) b
+    
+instance 
+    CanDiv (Differential (CN MPBall)) MPBall
+    where
+    type DivType     (Differential (CN MPBall)) MPBall = Differential (CN MPBall)
+    divide     b a = divide     b (differential (order b) (cn a))        
+
+instance 
+    (CanBeMPBall a) =>
+    CanMulAsymmetric (CN a) (Differential (CN MPBall))
+    where
+    type MulType (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)
+    mul a b = mul (differential (order b) (fmap mpBall a)) b
+
+instance 
+    (CanBeMPBall a) =>
+    CanMulAsymmetric (Differential (CN MPBall)) (CN a)
+        where
+        type MulType (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)
+        mul b a = mul b (differential (order b) (fmap mpBall a))
+
+instance 
+    (CanBeMPBall a) =>
+    CanDiv (CN a) (Differential (CN MPBall))
+    where
+    type DivType     (CN a) (Differential (CN MPBall)) = Differential (CN MPBall)
+    divide     a b = divide     (differential (order b) (fmap mpBall a)) b
+
+instance 
+    (CanBeMPBall a) =>
+    CanDiv (Differential (CN MPBall)) (CN a) 
+    where
+    type DivType     (Differential (CN MPBall)) (CN a) = Differential (CN MPBall)
+    divide     b a = divide     b (differential (order b) (fmap mpBall a))
+
+
+instance 
+    CanPow (Differential (CN MPBall)) Integer
+    where
+        type PowType     (Differential (CN MPBall)) Integer = (Differential (CN MPBall))
+        pow (OrderZero x) n =
+            OrderZero p
+            where
+            rawP = pow x n
+            p = if even n then abs rawP else rawP -- TODO: this should go in MPBall
+        pow (OrderOne x dx) n = 
+            OrderOne p dp
+            where
+            nEven = even n
+            rawP' = pow x (n - 1)
+            rawP  = x * rawP'
+            p  = if nEven then abs rawP else rawP  -- TODO: this should go in MPBall
+            dp = n * dx * if nEven then rawP' else abs rawP'
+        pow (OrderTwo x dx dxt d2x) n =
+            OrderTwo p dp dpt d2p
+            where
+            nEven  = even n
+            rawP'' = pow x (n - 2)
+            rawP'  = x * rawP''
+            rawP   = x * rawP'
+            p    = if nEven then abs rawP else rawP -- TODO: this should go in MPBall
+            xnm1 = if nEven then rawP' else abs rawP'
+            dp   = n * dx  * xnm1
+            dpt  = n * dxt * xnm1
+            xnm2 = if nEven then abs rawP'' else rawP''
+            d2p  = n*((n - 1)*dx*dxt*xnm2 + d2x*xnm1)
+
+instance
+    CanAbs (Differential (CN MPBall))
+    where
+    type AbsType (Differential (CN MPBall)) = Differential (CN MPBall)
+    abs (OrderZero x)   = OrderZero $ abs x 
+    abs (OrderOne x dx) = OrderOne (abs x) newDx
+        where
+            newDx = do 
+                dx_ <- dx
+                pure $ (hullMPBall dx_ (-dx_))
+    abs (OrderTwo _ _ _ _) = error "Abs for differential order two undefined"
+
+instance 
+    CanDivIMod (Differential (CN MPBall)) (Differential (CN MPBall))
+    where
+    type DivIType (Differential (CN MPBall)) (Differential (CN MPBall)) = (Differential (CN MPBall))
+    divIMod a b = (error "Integer division for Differential (CN MPBall) undefined", OrderTwo (mod ax bx) err err err)
+        where   
+            ax = diff_x a
+            bx = diff_x b
+            err = noValueNumErrorCertain $ NumError "No derivatives after modulus"
diff --git a/src/AERN2/AD/Type.hs b/src/AERN2/AD/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/AD/Type.hs
@@ -0,0 +1,46 @@
+module AERN2.AD.Type where
+
+import MixedTypesNumPrelude
+import AERN2.MP.Precision
+
+data Differential a =
+    OrderZero  {diff_x :: a}
+    | OrderOne {diff_x :: a, diff_dx :: a}
+    | OrderTwo {diff_x :: a, diff_dx :: a, diff_dxt :: a, diff_d2x :: a}
+    deriving (Show)
+
+order :: Differential a -> Integer
+order (OrderZero _)     = 0
+order (OrderOne  _ _)   = 1
+order (OrderTwo  _ _ _ _) = 2
+
+class CanBeDifferential a where
+    differential :: Integer -> a -> Differential a
+
+instance 
+    (HasIntegers a) =>
+    CanBeDifferential a
+    where 
+    differential 0 a = OrderZero a
+    differential 1 a = OrderOne  a (convertExactly 0)
+    differential _ a = OrderTwo  a (convertExactly 0) (convertExactly 0) (convertExactly 0)
+
+instance Functor Differential where
+    fmap f (OrderZero x)        = OrderZero (f x)
+    fmap f (OrderOne x dx)      = OrderOne  (f x) (f dx)
+    fmap f (OrderTwo x dx dxt d2x)  = OrderTwo  (f x) (f dx) (f dxt) (f d2x)
+
+instance 
+    (HasPrecision a) => (HasPrecision (Differential a))
+    where
+    getPrecision a = getPrecision (diff_x a) -- TODO: safe?
+
+instance 
+    (CanSetPrecision a) => (CanSetPrecision (Differential a))
+    where
+    setPrecision p = fmap (setPrecision p)
+
+setValue :: Differential a -> a -> Differential a
+setValue (OrderZero _x)           v = OrderZero v
+setValue (OrderOne _x dx)         v = OrderOne  v dx
+setValue (OrderTwo _x dx dxt d2x) v = OrderTwo  v dx dxt d2x
diff --git a/src/AERN2/BoxFun/Box.hs b/src/AERN2/BoxFun/Box.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/BoxFun/Box.hs
@@ -0,0 +1,150 @@
+{-# OPTIONS_GHC -Wno-orphans #-}
+module AERN2.BoxFun.Box where
+
+import qualified Prelude as Prelude
+import MixedTypesNumPrelude
+import AERN2.MP.Dyadic
+import AERN2.MP.Ball
+import AERN2.Linear.Vector.Type (Vector, (!))
+import qualified AERN2.Linear.Vector.Type as V
+import AERN2.Util.Util
+
+type Box = Vector (CN MPBall)
+
+-- instance HasEqAsymmetric Box Box where
+--     type EqCompareType Box Box = (CN Kleenean)
+--     equalTo box0 box1 = 
+--         V.foldl' (&&) (cn CertainTrue) $ V.zipWith equalTo box0 box1
+
+-- intersection :: Box -> Box -> Maybe Box
+-- intersection = undefined
+
+intersectionCertainlyEmpty :: Box -> Box -> Bool
+intersectionCertainlyEmpty vx vy =
+    V.foldl' (||) False $ V.zipWith ballIntersectionCertainlyEmpty vx vy
+    where 
+    ballIntersectionCertainlyEmpty x y = 
+        (lx !>! ry || rx !<! ly)
+        where
+        (lx, rx) = endpointsAsIntervals x
+        (ly, ry) = endpointsAsIntervals y
+
+nonEmptyIntersection :: Box -> Box -> Box
+nonEmptyIntersection vx vy = 
+    V.zipWith intersectCN vx vy
+
+instance IsBall Box where
+    type CentreType (Box) = Box
+    centre v = V.map (\x -> (cnMPBallP (getPrecision x) . centre) x) v
+
+inftyNorm :: Box -> CN MPBall
+inftyNorm box = 
+    V.foldl' (\n x -> max n (abs x)) (cn $ mpBall 0) box
+
+ellOneNorm :: Box -> CN MPBall
+ellOneNorm box = 
+    V.foldl' (\n x -> n + (abs x)) (cn $ mpBall 0) box
+
+width :: Box -> CN MPBall
+width box = 
+    V.foldl' (\n x -> max n (2 * ((fmap mpBall) $ (fmap radius) x) ) ) (cn $ mpBall 0) box
+
+widestDirection :: Box -> Integer
+widestDirection box =
+    aux (V.length box - 1) (V.length box - 1) (dyadic 0)
+    where
+    aux k i lth =
+        if k < 0 then 
+            i
+        else
+            let 
+                x    = box ! k
+                lth' = 
+                    2 * (dyadic $ radius $ unCN x) -- TODO: unsafe
+            in
+                if (lth' > lth) then 
+                    aux (k - 1) k lth'
+                else 
+                    aux (k - 1) i lth
+
+
+bisect :: Integer -> Box -> (Box, Box)
+bisect k box =
+    if exponent_ > 0
+    && (abs $ mc - lbc) < (dyadic 0.5)^exponent_ then
+        (setPrecision (increasePrecision p) lb, setPrecision (increasePrecision p) rb)
+    else 
+        (lb, rb)
+    where
+    exponent_      = integer p - ilog
+    NormBits ilog = 2 + (getNormLog $  mpBall $ 1 + abs lbc )
+    lbc = unCN (centre $ lb ! k)
+    -- rbc = unCN (centre $ rb ! k)
+    mc  = unCN (centre $ m)
+    increasePrecision p1 =
+        p1 + (prec $ (integer p1) `Prelude.div` 2)
+    lb = setPrecision p leftBox
+    rb = setPrecision p rightBox
+    p   = getPrecision box
+    interval = box ! k
+    m = (fmap $ mpBallP p) $ (fmap centre) interval
+    (l, r) = endpointsAsIntervals interval
+    leftBox  = V.imap (\i x -> if i == k then fromEndpointsAsIntervals l m else x) box
+    rightBox = V.imap (\i x -> if i == k then fromEndpointsAsIntervals m r else x) box
+
+-- Bisects a into 2^d boxes of the same size, where d is the dimension of the given box
+fullBisect :: Box -> [Box]
+fullBisect b =
+    case V.length b of
+        0 -> [b]
+        l ->
+            -- y is the dimension bisected in the current iteration
+            -- x is a bisection of the previous dimension (tail recursion)
+            concatMap (\x -> map (\y -> x V.+++ V.singleton y) (bisectDimension (l-1))) (fullBisect (V.take (fromIntegral (l-1)) b))
+
+            where
+                bisectDimension n = [fst bn ! n, snd bn ! n]
+                    where bn = bisect n b
+
+-- Get the endpoints of a box as a list containing a pair of MPBalls for each dimension
+getEndpoints :: Box -> [(MPBall, MPBall)]
+getEndpoints b  = 
+    case V.length b of
+        0 -> []
+        _ -> endpointsAsIntervals (unCN (b ! 0)) : getEndpoints (V.drop (int 1) b)
+
+lowerBounds :: Box -> Box
+lowerBounds = V.map lowerBound
+
+upperBounds :: Box -> Box
+upperBounds = V.map upperBound
+
+instance 
+    CanNeg Box
+    where
+    type NegType Box = Box
+    negate box = V.map (\x -> -x) box
+    
+createEnclosingBox :: Box -> Box -> Box
+createEnclosingBox box1 box2 =
+    enclosingBox
+    where
+        indexedBox1 = V.zip (V.fromList [0 .. integer (V.length box1) - 1]) box1
+        enclosingBox = 
+            V.map 
+            (\(i, x) -> 
+                let
+                    y = box2 V.! i
+                    (l, r) = endpointsAsIntervals x
+                    (l', r') = endpointsAsIntervals y
+                    newL = min l l'
+                    newR = max r r'
+                in
+                    fromEndpointsAsIntervals newL newR)
+            indexedBox1
+
+intersectList :: [Box] -> Box
+intersectList []            = V.empty
+intersectList [b]           = b
+intersectList [b1,b2]       = if intersectionCertainlyEmpty b1 b2 then V.empty else nonEmptyIntersection b1 b2
+intersectList (b1:b2:bs)    = if intersectionCertainlyEmpty b1 b2 then V.empty else intersectList $ nonEmptyIntersection b1 b2 : bs
diff --git a/src/AERN2/BoxFun/Optimisation.hs b/src/AERN2/BoxFun/Optimisation.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/BoxFun/Optimisation.hs
@@ -0,0 +1,319 @@
+module AERN2.BoxFun.Optimisation where
+
+import qualified Prelude
+import MixedTypesNumPrelude
+import qualified Numeric.CollectErrors as CN
+import AERN2.MP.Dyadic
+import AERN2.MP.Ball
+import AERN2.BoxFun.Box (Box)
+import qualified AERN2.BoxFun.Box as Box
+import AERN2.BoxFun.Type
+import AERN2.Kleenean
+import AERN2.Linear.Vector.Type as V
+import AERN2.Linear.Matrix.Type
+import AERN2.Linear.Matrix.Inverse
+
+import qualified Data.List as List
+
+import qualified AERN2.PQueue as Q
+
+import AERN2.Util.Util
+
+import Debug.Trace (trace)
+
+globalMinimumGreaterThanN :: BoxFun -> Accuracy -> CN Rational -> Precision -> Bool
+globalMinimumGreaterThanN f ac n initialPrecision =
+    trace (show x)
+    x !>! n
+    where x = globalMinimum f ac initialPrecision 
+
+minFun :: BoxFun -> Accuracy -> Precision -> (Integer, CN MPBall)
+minFun f ac initialPrecision = 
+    bestLocalMinimum f (domain f) ac initialPrecision
+
+data SearchBox =
+    SearchBox
+    {
+            extents :: Box
+        ,   minimum :: CN MPBall 
+    } deriving (Show)
+
+instance 
+    HasPrecision SearchBox
+    where
+    getPrecision (SearchBox b _) = getPrecision b
+
+instance 
+    CanSetPrecision SearchBox
+    where
+    setPrecision p (SearchBox b m) = SearchBox (setPrecision p b) m
+
+instance Prelude.Eq SearchBox where
+    (==) (SearchBox _ _) (SearchBox _ _) =
+        False -- TODO: safe?
+
+instance Prelude.Ord SearchBox where
+    (<=) (SearchBox _ min0) (SearchBox _ min1) = 
+        case (CN.toEither $ (lowerBound min0 :: CN MPBall), CN.toEither $ (lowerBound min1 :: CN MPBall)) of
+            (Left _, Left _) -> True
+            (Left _, Right _ ) -> True
+            (Right _ , Left _) -> False
+            (Right m0, Right m1) -> 
+                centre m0 - (dyadic $ radius m0) <= centre m1 - (dyadic $ radius m1) -- TODO: radius should be 0
+
+---
+
+globalMinimumWithCutoff :: BoxFun -> Accuracy -> CN MPBall -> Precision -> CN MPBall
+globalMinimumWithCutoff f ac cutoff initialPrecision =
+    if dimension f == 1 then
+        let
+            fl       = apply f (V.map lowerBound $ domain f)
+            fr       = apply f (V.map upperBound $ domain f)
+            localMin = snd $ bestLocalMinimumWithCutoff f (domain f) ac cutoff initialPrecision
+        in
+            min fl $ min localMin fr
+    else 
+        let
+            localMin       = snd $ bestLocalMinimumWithCutoff f (domain f) ac cutoff initialPrecision
+            boundaryFuns   = boundaryRestrictions f
+            boundaryMinima = List.map (\g -> globalMinimumWithCutoff g ac (min cutoff ((upperBound localMin :: CN MPBall))) initialPrecision) boundaryFuns
+        in
+            List.foldl' min localMin boundaryMinima
+
+
+globalMinimum :: BoxFun -> Accuracy -> Precision -> CN MPBall
+globalMinimum f ac initialPrecision =
+    globalMinimumWithCutoff f ac (apply f (centre boxp)) initialPrecision
+    where
+    boxp = setPrecision initialPrecision (domain f)
+
+bestLocalMinimum :: BoxFun -> Box -> Accuracy -> Precision -> (Integer, CN MPBall)
+bestLocalMinimum f box ac initialPrecision =
+    bestLocalMinimumWithCutoff f box ac (apply f (centre boxp)) initialPrecision
+    where
+    boxp = setPrecision initialPrecision box
+
+bestLocalMinimumWithCutoff :: BoxFun -> Box -> Accuracy -> CN MPBall -> Precision -> (Integer, CN MPBall)
+bestLocalMinimumWithCutoff f box ac initialCutoff initialPrecision =
+    aux initialQueue initialCutoff 0 dummyBox
+    where
+    boxp             = setPrecision initialPrecision box
+    initialRange     = apply f boxp
+    initialSearchBox = SearchBox boxp initialRange
+    initialQueue     = Q.singleton initialSearchBox
+    dummyBox         = SearchBox (V.fromList [cn $ mpBall $ 10^6]) initialRange -- TODO: hack...
+
+    aux q cutoff steps (SearchBox _lastBox rng) =  
+        case Q.minView q of
+            Nothing -> trace ("no local minimum.") $ (steps, rng)
+            Just (minBox, q') ->
+                --trace ("value: "++ (show $ val)) $
+                trace ("min box: "++ (show $ minBox)) $
+                --trace ("box acc: "++ (show $ getAccuracy $ ext)) $
+                --trace (show $ Box.width (extents minBox)) $
+                --trace ("lower bound "++ (show $ Box.lowerBound $ val)) $
+                --trace ("val' "++ (show $ val')) $
+                trace ("cutoff: "++ (show $ cutoff)) $
+                trace ("queue size: "++ (show $ Q.size q)) $
+                --trace ("cutoff == 0? "++(show $ cutoff == (mpBall 0))) $
+                --trace ("precision: "++ (show $ precision)) $
+                --trace ("dist to last "++ (show $ distToLast)) $
+                --trace ("accuracy: "++ (show $ getAccuracy val')) $
+                --trace ("precision centre: "++ (show $ fmap (getPrecision . centre) val)) $
+                if getAccuracy val' >= ac then
+                    (steps, val')
+                else 
+                    aux q'' newCutoff (steps + 1) (SearchBox ext rng)
+                where
+                val' = fromEndpointsAsIntervals (lowerBound val) (cutoff)
+                SearchBox ext val = minBox
+
+                (newCutoff, newBoxes) = 
+                    processBox f ac cutoff minBox
+
+                q'' = foldr (Q.insert) q' newBoxes
+
+lipschitzContraction :: BoxFun -> Box -> SearchBox -> SearchBox
+lipschitzContraction f g (SearchBox box m) =
+    {-trace("fa: "++(show $ getAccuracy (apply f box))) $
+    trace("la: "++(show $ getAccuracy $ dotProduct)) $
+    trace("ba: "++(show $ getAccuracy $ box ! int 0)) $-}
+    {-if (radius $ (~!) $ newRange) < (radius $ (~!) $ m) then
+        trace ("Lipschitz better.")
+        SearchBox box m'
+    else -}
+    SearchBox box m'
+    where
+    boxCentre      = centre box
+    centreValue    = apply f boxCentre
+    difference     = box - boxCentre
+    dotProduct     = g * difference 
+    newRange       = centreValue + dotProduct
+    m'             = intersectCN m newRange
+
+lipschitzRange :: BoxFun -> CN MPBall -> Box -> Box -> Box -> CN MPBall -> CN MPBall
+lipschitzRange _f fc c g box m =
+    m'
+    where
+    difference     = box - c
+    normG          = Box.ellOneNorm g
+    normDiff       = Box.inftyNorm  difference
+    dotProduct     = normG * normDiff
+    newRange       = fc + (fromEndpointsAsIntervals (-dotProduct) dotProduct :: CN MPBall)
+    m'             = intersectCN m newRange
+
+applyLipschitz :: BoxFun -> Box -> CN MPBall
+applyLipschitz f box = 
+    lipschitzRange f fbc bc dfb' box fb
+    where
+    (fb, dfb') = valueGradient f box
+    bc  = centre box
+    fbc  = apply f bc
+
+increasePrecision :: Precision -> Precision
+increasePrecision p =
+    p + (prec $ (integer p) `Prelude.div` 2)
+
+newtonStep :: BoxFun -> Accuracy -> Vector (CN MPBall) -> Vector (CN MPBall) -> Matrix (CN MPBall) -> SearchBox -> Bool -> Maybe (Bool, SearchBox)
+newtonStep f ac c dfc hInv b@(SearchBox box m) newtonSuccesful = 
+    --Just $ SearchBox box' m'
+    {-trace ("precision m "++(show $ (fmap getPrecision) m)) $
+    trace ("precision m' "++(show $ (fmap getPrecision) m')) $
+    trace ("precision box centre "++(show $ getPrecision c)) $
+    trace ("precision box "++(show $ getPrecision box)) $
+    trace ("precision newton box "++(show $ getPrecision newtonBox)) $
+    trace ("precision box' "++(show $ getPrecision box')) $
+    trace ("precision hInv "++(show $ getPrecision (entries hInv ! int 0))) $-}
+    if getAccuracy m >= ac then
+        Just (newtonSuccesful, b)
+    --else if not hInvDefined then
+    --    Just (newtonSuccesful, b)
+    else if Box.intersectionCertainlyEmpty box newtonBox then
+        Nothing
+    else if Box.width box' !<=! (dyadic $ 0.75) * Box.width box then
+        if getAccuracy m' > getAccuracy m then
+            newtonStep f ac c dfc hInv (SearchBox box' m') True
+        else 
+            Just (True, SearchBox (setPrecision (increasePrecision $ getPrecision box') box') m')
+    else 
+        Just (newtonSuccesful, SearchBox box' m')
+    where
+    {-c           = centre box
+    dfc         = gradient f c-}
+    -- hInvDefined = V.foldl' (&&) (True) $ V.map (isJust . fst . ensureNoCN) (entries hInv)
+    newtonBox   = c - hInv * (dfc)
+    box'        = Box.nonEmptyIntersection box newtonBox
+    m'          = apply f box'
+
+processBox :: BoxFun -> Accuracy -> CN MPBall -> SearchBox -> (CN MPBall, [SearchBox])
+processBox f ac cutoff box =
+    if getAccuracy ext < bits 10 then 
+        split f (gradient f ext) cutoff ext  
+    else 
+        result
+    where
+    ext            = extents box
+    (_fb, dfb, hfb) = valueGradientHessian f ext
+    c              = centre ext
+    dfc            = gradient f c
+    maybeHinv      = inverse hfb
+    -- p              = getPrecision box
+    box'           = --Just (False, box)
+        case maybeHinv of 
+            Nothing   -> Just (False, box)
+            Just hInv -> newtonStep f ac c dfc hInv box False
+    result =
+        case box' of 
+            Nothing -> (cutoff, [])
+            Just (newtonSuccesful, bx@(SearchBox bxe m)) ->
+                let
+                    c' = min (upperBound $ apply f $ centre bxe :: CN MPBall) cutoff    
+                in
+                if newtonSuccesful then
+                    if getAccuracy m >= ac then
+                        (c', [bx])
+                    else
+                        processBox f ac c' bx
+                else
+                    split f dfb c' bxe
+
+split :: BoxFun -> Vector (CN MPBall) -> CN MPBall -> Box -> (CN MPBall, [SearchBox])
+split f dfb cutoff bxe = 
+    let
+    diff    = bxe - centre bxe
+    dir i   = (fmap dyadic) $ (fmap radius) $ (dfb ! i) * (diff ! i) :: CN Dyadic
+    dirs    = V.map dir $ V.enumFromTo 0 (V.length bxe - 1)
+    dirsDefined = V.foldl' (&&) True $ V.map (not . CN.hasError) dirs
+    aux k j d = 
+        if k == V.length bxe then 
+            j 
+        else 
+            let
+                d' = unCN $ dirs ! k
+            in
+            if d' > d then 
+                aux (k + 1) k d'
+            else
+                aux (k + 1) j d
+    splittingIndex = 
+        if dirsDefined then (aux 1 0 (unCN $ dirs ! 0)) else Box.widestDirection bxe
+    (a , b)    = Box.bisect splittingIndex bxe
+    (fa, dfa') = valueGradient f a
+    (fb, dfb') = valueGradient f b
+    ac  = centre a
+    bc  = centre b
+    fac  = apply f ac
+    fbc  = apply f bc
+    fa'  = lipschitzRange f fac ac dfa' a fa
+    fb'  = lipschitzRange f fbc bc dfb' b fb
+    cutoff'       = min (upperBound fac :: CN MPBall) $ min (upperBound fbc :: CN MPBall) cutoff
+    leftMonotone  = V.foldl' (||) False $ V.map (!/=! 0) dfa'
+    rightMonotone = V.foldl' (||) False $ V.map (!/=! 0) dfb'
+    boxes = 
+        case (leftMonotone || fa' !>! cutoff', rightMonotone || fb' !>! cutoff') of
+            (True,  True)  -> []
+            (True,  False) -> [SearchBox b fb']
+            (False, True)  -> [SearchBox a fa']
+            (False, False) -> [SearchBox a fa', SearchBox b fb']
+    in
+        (cutoff', boxes)
+
+-- Precondition: f and g must have the same domain
+maxBoxFunGreaterThanN :: BoxFun -> BoxFun -> CN Rational -> Precision -> Bool
+maxBoxFunGreaterThanN f g n initialPrecision =
+    case Box.getEndpoints fbox == Box.getEndpoints gbox of
+        CertainTrue ->
+            checkMaxAboveN f g ||
+                (Box.width fboxp !>! cutoff && Box.width gboxp !>! cutoff) &&
+                    let
+                        newBoxes = Box.fullBisect fboxp
+
+                        updateDomain z = BoxFun (dimension z) (bf_eval z)
+
+                        checkBoxes [] = True
+                        checkBoxes (box : boxes) = 
+                            if checkMaxAboveN (updateDomain f box) (updateDomain g box) 
+                                then checkBoxes boxes
+                                else maxBoxFunGreaterThanN f' g' n initialPrecision && checkBoxes boxes
+                                where
+                                    f' = updateDomain f box
+                                    g' = updateDomain g box
+
+                    in
+                        checkBoxes newBoxes
+
+        _ ->
+            trace "Domain of f not equal to domain of g"
+            False
+    where
+        cutoff = 1/2^10
+
+        fbox                = domain f
+        fboxp               = setPrecision initialPrecision fbox 
+
+        gbox                = domain g
+        gboxp               = setPrecision initialPrecision gbox
+
+        checkMaxAboveN h i = applyMinimum h !>! n || applyMinimum i !>! n
+
+
diff --git a/src/AERN2/BoxFun/TestFunctions.hs b/src/AERN2/BoxFun/TestFunctions.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/BoxFun/TestFunctions.hs
@@ -0,0 +1,425 @@
+module AERN2.BoxFun.TestFunctions where
+
+import MixedTypesNumPrelude
+import AERN2.MP.Ball
+import qualified Data.List as List
+
+import AERN2.AD.Differential
+import AERN2.BoxFun.Type
+import AERN2.Linear.Vector.Type ((!), Vector)
+import qualified AERN2.Linear.Vector.Type as V
+
+fromListDomain :: [(Rational, Rational)] -> Vector (CN MPBall)
+fromListDomain [] = V.empty
+fromListDomain [x] = V.singleton $ fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ (fst x)) (cn $  mpBallP (prec 53) (snd x))
+fromListDomain (x:xs) = V.cons (fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ (fst x)) (cn $  mpBallP (prec 53) (snd x))) (fromListDomain xs)
+
+symmetricDomain :: Integer -> Rational -> Rational -> Vector (CN MPBall)
+symmetricDomain n l r = 
+    V.map (\_ -> fromEndpointsAsIntervals (cn $ mpBallP (prec 53) $ l) (cn $  mpBallP (prec 53) r)) $ V.enumFromTo 1 n
+
+griewank :: Integer -> BoxFun
+griewank n =
+    BoxFun
+    n
+    (\v ->
+        let
+            ord = order (v ! 0)
+            sm  = List.foldl' (+) (differential ord (cn $ mpBall 0)) [let x = (v ! k) in x^2 | k <- [0 .. V.length v - 1]]
+            prd = List.foldl' (*) (differential ord (cn $ mpBall 1)) [cos $ (v ! k) / (sqrt $ 1 + mpBallP p k) | k <- [0 .. V.length v - 1]]
+            p   = getPrecision v
+        in
+        1 + ((setPrecision p $ mpBall 1)/(setPrecision p $ mpBall 4000)) * sm - prd
+    )
+    (symmetricDomain n (-600.0) 600.0)
+
+rosenbrock :: BoxFun
+rosenbrock =
+    BoxFun
+    2
+    (\v ->
+        let
+            p = getPrecision v
+            x = v ! 0
+            y = v ! 1
+            a = mpBallP p 1
+            b = mpBallP p 100
+            amx  = a - x
+            ymxs = y - x*x
+        in
+            amx^2 + b*ymxs^2
+    )
+    (symmetricDomain 2 (-1.2) 1.2)
+
+himmelblau :: BoxFun
+himmelblau =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v ! 0
+            y = v ! 1
+            a = (x^2 + y - 11)
+            b = (x + y^2 - 7)
+        in
+            a^2 + b^2
+    )
+    (symmetricDomain 2 (-600.0) 600.0)
+
+shekel :: BoxFun
+shekel = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x   = v ! 0
+            y   = v ! 1
+            c0  = 1;  a00 = 43; a01 = 23
+            c1  = 17; a10 = 6;  a11 = 9
+        in
+            - 1/(c0 + (x - a00)^2 + (y - a01)^2)
+            - 1/(c1 + (x - a10)^2 + (y - a11)^2)
+    )
+    (symmetricDomain 2 (-600.0) 600.0)
+
+siam4 :: BoxFun
+siam4 = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x   = v ! 0
+            y   = v ! 1
+        in
+        exp(sin(50 * x)) + sin(60 * exp y) 
+        + sin(70 * sin(x)) + sin(sin(80 * y)) 
+        - sin(10 * (x + y)) 
+        + (x^2 + y^2) / 4
+    )
+    (symmetricDomain 2 (-10.0) 10.0)
+
+ratz4 :: BoxFun
+ratz4 =
+    BoxFun
+    2
+    (\v ->
+        let 
+            x  = v ! 0
+            y  = v ! 1
+            xs = x^2
+            ys = y^2
+        in 
+            sin(xs + 2 * ys) * exp (-xs - ys)
+    )
+    (symmetricDomain 2 (-3.0) 3.0)
+
+-- global minimum at bukin(-10, 1) ~ 0
+-- bukin :: BoxFun
+-- bukin =
+--     BoxFun
+--     2
+--     (\v ->
+--         let
+--             x = v!0
+--             y = v!1
+--         in
+--             100 * sqrt (abs (y - x^2 / 100)) + abs(x + 10) / 100
+--     )
+--     (fromListDomain [(-15.0, 5.0), (-3.0, 3.0)])
+
+ackley :: BoxFun
+ackley =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            pi = piBallP p
+        in
+            -20 * exp(sqrt((x^2 + y^2) / 2) / (-5)) - exp((cos (2 * pi * x) + cos(2 * pi * y)) / 2) + exp(mpBallP p 1) + 20
+    )
+    (symmetricDomain 2 (-5.0) 5.0)
+
+-- eggholder :: BoxFun
+-- eggholder =
+--     BoxFun
+--     2
+--     (\v ->
+--         let
+--             x = v!0
+--             y = v!1
+--         in
+--             -(y + 47) * sin (sqrt (abs (x / 2 + (y + 47)))) - x * sin (sqrt (abs (x - (y + 47))))
+--     )
+--     (symmetricDomain 2 (-512.0) 512.0)
+    
+-- heron :: BoxFun
+-- heron = 
+--     BoxFun
+--     2
+--     (\v ->
+--         let
+--             x = v!0
+--             y = v!1
+--             p = getPrecision v
+--             eps = 1/2^(23)
+--             i = 2
+--         in
+--             max
+--                 min
+--                     max
+--                         (y - sqrt x)
+--                         ((sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - 6 * eps * (i-1))
+--                     max
+--                         (sqrt x - y)
+--                         (- (sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - 6 * eps * (i-1))
+--                 min
+--                     max
+--                         ((y + x/y)/2 - sqrt x)
+--                         (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + 6 * eps * (i-1))
+--                     max
+--                         (sqrt x - (y+x/y)/2)
+--                         ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + 6 * eps * (i-1))
+
+--     )
+--     (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+-- max (min (max 1p 1q) (max 2p 2q)) (min (max 3p 3q) (max 4p 4q))
+
+i :: Integer
+i = 3
+
+heron1p :: BoxFun
+heron1p = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            _p = getPrecision v
+            _eps = 1/2^(23)
+        in
+            (y - sqrt x)
+     )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron1q :: BoxFun
+heron1q = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            eps = 1/2^(23)
+        in
+            ((sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - (mpBallP p 6) * eps * (i-1))
+
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron2p :: BoxFun
+heron2p = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            _p = getPrecision v
+            _eps = 1/2^(23)
+        in
+            (sqrt x - y)
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron2q :: BoxFun
+heron2q = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            eps = 1/2^(23)
+        in
+            (- (sqrt x - y) - (mpBallP p 1/2)^(2^(i-1)) - (mpBallP p 6) * eps * (i-1))
+
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron3p :: BoxFun
+heron3p = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            _p = getPrecision v
+            _eps = 1/2^(23)
+        in
+            ((y + x/y)/2 - sqrt x)
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+    
+heron3q :: BoxFun
+heron3q = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            eps = 1/2^(23)
+        in
+            (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron3p2 :: BoxFun
+heron3p2 = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            _p = getPrecision v
+            _eps = 1/2^(23)
+            _i = 4
+        in
+            ((y + x/y)/2 - sqrt x)
+    )
+    (fromListDomain [(0.633758544921875, 0.63385009765625), (0.79999999999999982236431605997495353221893310546875, 0.80006103515624982239142111428709114306911942549049854278564453125)])
+    
+heron3q2 :: BoxFun
+heron3q2 = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            eps = 1/2^(23)
+            _i = 4
+        in
+            (- (sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))
+    )
+    (fromListDomain [(0.633758544921875, 0.63385009765625), (0.79999999999999982236431605997495353221893310546875, 0.80006103515624982239142111428709114306911942549049854278564453125)])
+
+heron4p :: BoxFun
+heron4p = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            _p = getPrecision v
+            _eps = 1/2^(23)
+        in
+            (sqrt x - (y+x/y)/2)
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+    
+heron4q :: BoxFun
+heron4q = 
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            y = v!1
+            p = getPrecision v
+            eps = 1/2^(23)
+        in
+            ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))
+
+    )
+    (fromListDomain [(0.5, 2.0), (0.8, 1.8)])
+
+-- heronFull :: BoxFun
+-- heronFull =
+--     BoxFun
+--     3
+--     (\v ->
+--         let
+--             x = v!0
+--             y = v!1
+--             i = v!2
+--             p = getPrecision v
+--             eps = 1/2^(23)
+--         in
+--             ((sqrt x - (y+x/y)/2) + (mpBallP p 1/2)^(2^i) + (mpBallP p 6) * eps * (i-1))
+
+--     )
+--     (fromListDomain [(0.5, 2.0), (0.8, 1.8), (1.0, 5.0)])
+
+
+mxp1 :: BoxFun
+mxp1 =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            -- y = v!1
+        in
+            -x+1
+
+    )
+    (fromListDomain [(0.0, 2.0), (0.0, 2.0)])
+
+xm1 :: BoxFun
+xm1 =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            -- y = v!1
+        in
+            x-1
+
+    )
+    (fromListDomain [(0.0, 2.0), (0.0, 2.0)])
+
+xe2p1 :: BoxFun
+xe2p1 =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            -- y = v!1
+        in
+            x^2+1
+    )
+    (fromListDomain [(0.0, 2.0), (0.0, 2.0)])
+
+xe2m1 :: BoxFun
+xe2m1 =
+    BoxFun
+    2
+    (\v ->
+        let
+            x = v!0
+            -- y = v!1
+        in
+            x^2-1
+    )
+    (fromListDomain [(0.0, 2.0), (0.0, 2.0)])
diff --git a/src/AERN2/BoxFun/Type.hs b/src/AERN2/BoxFun/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/BoxFun/Type.hs
@@ -0,0 +1,155 @@
+module AERN2.BoxFun.Type where
+
+import MixedTypesNumPrelude
+
+import AERN2.Linear.Vector.Type (Vector, (!))
+import qualified AERN2.Linear.Vector.Type as V
+import AERN2.MP.Ball
+import AERN2.AD.Differential
+import qualified AERN2.Linear.Matrix.Type as M
+import AERN2.Util.Util
+import AERN2.BoxFun.Box
+import Numeric.CollectErrors
+
+data BoxFun =
+    BoxFun
+    {
+            dimension :: Integer
+        ,   bf_eval   :: Vector (Differential (CN MPBall)) -> Differential (CN MPBall)
+        ,   domain    :: Vector (CN MPBall)
+    }
+
+instance Show BoxFun where
+    show (BoxFun d _f b) = show d ++ " dimensional BoxFun, domain: " ++ show b
+
+boundaryRestrictions :: BoxFun -> [BoxFun]
+boundaryRestrictions (BoxFun d ev dom) =
+    concat
+    [
+        [
+            BoxFun 
+                (d - 1)
+                (
+                \v ->
+                    ev $ V.map (\j -> if j == i then (setPrecision $ getPrecision v) $ differential 2 $ upperBound (dom ! i) else if j < i then v ! j else v ! (j - 1)) $ V.enumFromTo 0 (d - 1)
+                ) 
+                (V.map (\j -> if j >= i then dom ! (j + 1) else dom ! j) $ V.enumFromTo 0 (d - 2))
+            ,
+            BoxFun 
+                (d - 1)
+                (
+                    \v ->
+                        ev $ V.map (\j -> if j == i then (setPrecision $ getPrecision v) $ differential 2 $ lowerBound (dom ! i) else if j < i then v ! j else v ! (j - 1)) $ V.enumFromTo 0 (d - 1)
+                )
+                (V.map (\j -> if j >= i then dom ! (j + 1) else dom ! j) $ V.enumFromTo 0 (d - 2))
+        ]
+        |
+        i <- [0 .. d - 1]
+    ]
+
+valueGradientHessian :: BoxFun -> Vector (CN MPBall) -> (CN MPBall, Vector (CN MPBall), M.Matrix (CN MPBall))
+valueGradientHessian (BoxFun d e _) v =
+    (value, grad, hess)
+    where
+    vgh i j = e (w i j)
+
+    triangle =  
+        V.map (\i-> V.map (\j -> vgh i j) $ V.enumFromTo 0 i) $ V.enumFromTo 0 (d - 1)
+
+    value = diff_x $ (triangle ! 0) ! 0
+    grad  = V.map (\i -> diff_dx $ (triangle ! i) ! 0) $ V.enumFromTo 0 (d - 1)
+    hess  = M.create d d (\i j -> diff_d2x $ if i > j then (triangle ! i) ! j else (triangle ! j) ! i)
+
+    w i j = V.imap (\k x -> OrderTwo x (delta i k) (delta j k) (pure $ mpBall 0)) v
+    delta :: Integer -> Integer -> CN MPBall
+    delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)
+
+valueGradient :: BoxFun -> Vector (CN MPBall) -> (CN MPBall, Vector (CN MPBall))
+valueGradient (BoxFun d e _) v =
+    aux (d - 1) [] (pure $ mpBall 0)
+    where
+    tangent k = 
+        V.imap (\i x -> OrderOne x (delta i k)) v
+    valgrad k =
+        let
+            etk = e (tangent k)
+        in 
+            (diff_x etk, diff_dx etk)
+    aux k ret val =
+        if k < 0 then
+            (val, V.fromList ret)
+        else 
+            let
+                (val2,g) = valgrad k
+            in
+            aux (k - 1) (g : ret) val2
+    delta :: Integer -> Integer -> CN MPBall
+    delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)
+
+apply :: BoxFun -> Vector (CN MPBall) -> CN MPBall
+apply (BoxFun _d e _) v = 
+    diff_x (e v')
+    where
+    v' = V.map (\x -> differential 0 x) v
+
+applyMinimum :: BoxFun -> CN MPBall
+applyMinimum h = fst $ endpointsAsIntervals (apply h (domain h))
+
+applyMinimumOnBox :: BoxFun -> Vector (CN MPBall) -> CN MPBall
+applyMinimumOnBox h hbox = fst $ endpointsAsIntervals (apply h hbox)
+
+applyMaximum :: BoxFun -> CN MPBall
+applyMaximum h = snd $ endpointsAsIntervals (apply h (domain h))
+
+applyMaximumOnBox :: BoxFun -> Vector (CN MPBall) -> CN MPBall
+applyMaximumOnBox h hbox = snd $ endpointsAsIntervals (apply h hbox)
+
+gradient :: BoxFun -> Vector (CN MPBall) -> Vector (CN MPBall)
+gradient (BoxFun d e _) v =
+    aux (d - 1) []
+    where
+    tangent k = 
+        V.imap (\i x -> OrderOne x (delta i k)) v
+    grad k =
+        -- diff_dx $ e (tangent k)
+        case e (tangent k) of
+            OrderZero _ -> noValueNumErrorCertain (NumError "details")
+            val -> diff_dx $ val
+    aux k ret =
+        if k < 0 then
+            V.fromList ret
+        else 
+            aux (k - 1) (grad k : ret)
+    delta :: Integer -> Integer -> CN MPBall
+    delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)
+
+hessian :: BoxFun -> Vector (CN MPBall) -> M.Matrix (CN MPBall)
+hessian (BoxFun d e _) v = 
+    M.create d d a
+    where
+    a i j = diff_d2x $ e (w i j)
+    w i j = V.imap (\k x -> OrderTwo x (delta i k) (delta j k) (pure $ mpBall 0)) v
+    delta :: Integer -> Integer -> CN MPBall
+    delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)
+
+jacobian :: [BoxFun] -> Vector (CN MPBall) -> M.Matrix (CN MPBall)
+jacobian fs v =
+    M.create (length fs) highestDiminseionInFs a
+    where
+    highestDiminseionInFs = maximum (map dimension fs)
+
+    a i j = diff_dx $ bf_eval (fs!!i) (w j)
+    w j = V.imap (\k x -> OrderOne x (delta j k)) v
+    delta :: Integer -> Integer -> CN MPBall
+    delta i k = if i == k then (cn $ mpBall 1) else (cn $ mpBall 0)
+
+gradientUsingGradient :: BoxFun -> Box -> Box
+gradientUsingGradient f v =
+    V.zipWith fromEndpointsAsIntervals lowerBounds upperBounds
+    where
+        lowerBounds = firstDerivatives - secondDerivatives * V.map (fmap (mpBall . radius)) v
+        upperBounds = firstDerivatives + secondDerivatives * V.map (fmap (mpBall . radius)) v
+        firstDerivatives  = head $ M.rows $ jacobian [f] (centre v)
+        secondDerivatives = fmap abs $ hessian f v
+        p = getPrecision v
+        
diff --git a/src/AERN2/Linear/Matrix/Inverse.hs b/src/AERN2/Linear/Matrix/Inverse.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/Linear/Matrix/Inverse.hs
@@ -0,0 +1,126 @@
+module AERN2.Linear.Matrix.Inverse 
+(inverse) 
+where
+
+import MixedTypesNumPrelude
+import qualified Numeric.CollectErrors as CN
+
+import AERN2.Linear.Matrix.Type
+import qualified AERN2.Linear.Vector.Type as V
+import Data.Maybe
+import AERN2.MP.Float
+import AERN2.MP.Dyadic
+import AERN2.MP
+
+inverse :: Matrix (CN MPBall) -> Maybe (Matrix (CN MPBall))
+inverse m = 
+    if  CN.hasError m
+    ||  (isNothing maybeY)
+    ||  (not $ nerr !<! 1) then
+        Nothing
+    else
+        Just $ aux x0
+    where
+    cm     = Matrix (width m) (V.map (\x -> mpFloat $ centre $ unCN x) $ entries m)
+    idM    = identity (width m) (width m)                    :: Matrix (CN MPBall)
+    maybeY = inverseGauss cm
+    y      = fmap (cn . mpBall . dyadic) $ fromJust maybeY
+    yid    = y * idM
+    nye    = V.map V.inftyNorm (columns yid)
+    err    = idM - y * m
+    nerr   = inftyNorm err
+    ui     = cn $ fromEndpointsAsIntervals (mpBall $ -1) (mpBall 1)
+    x0     = create (width m) (width m) (\_ j -> (ui *  (nye V.! j)) / (1 - nerr))
+    aux x  = let x' = it x in if getAccuracy x' <= getAccuracy x then x' else aux x'
+    it x   = intersect x $ yid + err * x
+
+inverseGauss :: Matrix MPFloat -> Maybe (Matrix MPFloat)
+inverseGauss m = 
+    aux 0 m inv0
+    where
+    inv0        = identity (width m) (width m) :: Matrix MPFloat
+    aux j n inv = 
+        if j == width m then
+            Just inv
+        else
+            let
+                n'    = pivot n j n
+                inv'  = pivot n j inv
+                a     = get n' j j 
+                ia    = ((mpFloat 1)/.a)
+                n''   = multiplyRow n'   j ia 
+                inv'' = multiplyRow inv' j ia 
+                elim i l = 
+                    if i /= j then
+                        let
+                            b = get l i j
+                        in
+                            addRows l (-b) j i
+                    else 
+                        l
+                n'''   = 
+                    foldr elim n'' [0 .. (height m) - 1]
+                inv''' = 
+                    foldr elim inv'' [0 .. (height m) - 1]
+            in 
+            if a /= mpFloat 0 then
+                aux (j + 1) n''' inv'''
+            else 
+                Nothing
+
+
+largestRowIndex :: Matrix MPFloat -> Integer -> Integer
+largestRowIndex m j =
+    integer $ aux (j + 1) j (get m j j)
+    where
+    h = height m
+    aux k i x =
+        if k == h then
+            i
+        else 
+            let
+                x' = get m i j
+            in 
+            if x' > x then
+                aux (k + 1) k x'
+            else 
+                aux (k + 1) i x
+    
+pivot :: Matrix MPFloat -> Integer -> Matrix MPFloat -> Matrix MPFloat
+pivot m j n = 
+    swapRows n j i
+    where
+    i = largestRowIndex m j
+
+multiplyRow :: Matrix MPFloat -> Integer -> MPFloat -> Matrix MPFloat
+multiplyRow m i a = 
+    imap h m
+    where
+    h k _ x = 
+        if k == i then 
+            a *. x
+        else 
+            x
+
+{- add a*(row i0) to row i1 -}
+addRows :: Matrix MPFloat -> MPFloat -> Integer -> Integer -> Matrix MPFloat
+addRows m a i0 i1 = 
+    imap h m
+    where
+    h k l x = 
+        if k == i1 then 
+            x +. a *. get m i0 l
+        else 
+            x
+
+swapRows :: Matrix MPFloat -> Integer -> Integer -> Matrix MPFloat
+swapRows m i0 i1 =
+    imap h m
+    where
+    h k l x = 
+        if k == i0 then 
+            get m i1 l
+        else if k == i1 then
+            get m i0 l
+        else 
+            x
diff --git a/src/AERN2/Linear/Matrix/Type.hs b/src/AERN2/Linear/Matrix/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/Linear/Matrix/Type.hs
@@ -0,0 +1,132 @@
+module AERN2.Linear.Matrix.Type where
+
+import qualified Prelude as P
+import MixedTypesNumPrelude
+import qualified Numeric.CollectErrors as CN
+import AERN2.Linear.Vector.Type (Vector, (!))
+import qualified AERN2.Linear.Vector.Type as V
+import qualified Data.List as List
+
+import AERN2.MP.Ball
+
+data (Matrix a) = 
+    Matrix
+    {
+        width     :: Integer,
+        entries   :: Vector a
+    } deriving (Show)
+
+height :: Matrix a -> Integer
+height (Matrix w e) = 
+    (V.length e) `P.div` w
+
+get :: Matrix a -> Integer -> Integer -> a
+get m i j =
+    entries m ! (i * (width m) + j)
+
+identity :: (HasIntegers a) => Integer -> Integer -> Matrix a
+identity m n = 
+    diag m n (convertExactly 1)
+
+diag :: (HasIntegers a) => Integer -> Integer -> a -> Matrix a
+diag m n x = 
+    create m n (\i j -> if i == j then x else (convertExactly 0))
+
+rows :: Matrix a -> [Vector a]
+rows m@(Matrix w e) = 
+    [V.slice (i*w) w e| i <- [0 .. height m - 1]]
+
+columns :: Matrix a -> Vector (Vector a)
+columns m = 
+    V.map (\j -> V.map (\i -> get m i j) $ V.enumFromTo 0 (height m - 1)) $ V.enumFromTo 0 (width m - 1)
+
+create :: Integer -> Integer -> (Integer -> Integer -> a) -> Matrix a
+create m n f =
+    Matrix n $ V.map (\x -> f (i x) (j x)) $ V.enumFromTo 0 (m*n - 1)
+    where
+    j x = x `mod` n
+    i x = (x - j x) `P.div` n
+
+imap :: (Integer -> Integer -> a -> a) -> Matrix a -> Matrix a
+imap f (Matrix w ents) =
+    Matrix w (V.imap g ents)
+    where
+    j x = x `mod` w
+    i x = (x - j x) `P.div` w
+    g k x = f (i k) (j k) x
+
+instance CanIntersectAsymmetric (Matrix (CN MPBall)) (Matrix (CN MPBall)) where
+    type IntersectionType (Matrix (CN MPBall)) (Matrix (CN MPBall)) = Matrix (CN MPBall)
+    intersect (Matrix w0 v0) (Matrix _w1 v1) =
+        Matrix w0 $ V.zipWith intersect v0 v1
+
+inftyNorm :: (CanAddSameType a, CanSubSameType a, CanAbsSameType a, HasIntegers a, CanMinMaxSameType a) => Matrix a -> a
+inftyNorm (m :: Matrix a) = 
+    -- TODO: could be optimised.
+    List.foldl' max (convertExactly 0 :: a)
+    [
+        V.foldl' (+) (convertExactly 0 :: a) $ V.map abs r 
+        |
+        r <- rows m
+    ]
+
+instance Functor Matrix where
+    fmap h m = 
+        Matrix (width m) (V.map h (entries m))
+
+instance 
+    (CanAddSameType a, CanMulSameType a, HasIntegers a) =>
+    CanMulAsymmetric (Matrix a) (Matrix a)
+    where
+    type MulType (Matrix a) (Matrix a) = Matrix a
+    mul m0 m1 = 
+        create (height m0) (width m1) (aux 0 (convertExactly 0))
+        where
+        aux k sm i j = 
+            if k == width m0 then 
+                sm
+            else 
+                aux (k + 1) (sm + (get m0 i k) * (get m1 k j)) i j
+
+
+instance 
+    (CanAddSameType a) =>
+    CanAddAsymmetric (Matrix a) (Matrix a)
+    where
+    type AddType (Matrix a) (Matrix a) = Matrix a
+    add (Matrix w e) (Matrix _ e') =
+        Matrix w (e + e')
+
+instance 
+    (CanSubSameType a) =>
+    CanSub (Matrix a) (Matrix a)
+    where
+    type SubType (Matrix a) (Matrix a) = Matrix a
+    sub (Matrix w e) (Matrix _ e') =
+        Matrix w (e - e')
+    
+
+instance 
+    (CanAddSameType a, CanMulSameType a, HasIntegers a) =>
+    CanMulAsymmetric (Matrix a) (Vector a)
+    where
+    type MulType (Matrix a) (Vector a) = Vector a
+    mul m@(Matrix _w _e) v =
+        V.fromList [r * v| r <- rows m]
+
+instance 
+    (HasAccuracy a, HasPrecision a) => HasAccuracy (Matrix a)
+    where
+    getAccuracy m = 
+        V.foldl' max NoInformation $ V.map getAccuracy (entries m)
+
+instance 
+    (HasPrecision a) => HasPrecision (Matrix a)
+    where
+    getPrecision m = 
+        V.foldl' max (prec 2) $ V.map getPrecision (entries m)
+
+instance 
+    (CN.CanTestErrorsPresent a) => CN.CanTestErrorsPresent (Matrix a)
+    where
+    hasError m = V.foldl' (||) False $ V.map (CN.hasError) (entries m)
diff --git a/src/AERN2/Linear/Vector/Type.hs b/src/AERN2/Linear/Vector/Type.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/Linear/Vector/Type.hs
@@ -0,0 +1,176 @@
+module AERN2.Linear.Vector.Type where
+
+import Control.Monad.ST
+import Data.STRef
+import MixedTypesNumPrelude hiding (length)
+import qualified Data.Vector as V
+import qualified Data.Vector.Generic.Mutable as M
+import AERN2.MP.Precision
+import AERN2.MP.Ball
+import qualified Prelude as P
+
+type (Vector a) = V.Vector a
+
+(+++) :: Vector a -> Vector a -> Vector a
+(+++) = (V.++)
+
+drop :: Int -> Vector a -> Vector a 
+drop = V.drop
+
+take :: Int -> Vector a -> Vector a 
+take = V.take
+
+empty :: Vector a
+empty = V.empty
+
+singleton :: a -> Vector a
+singleton = V.singleton
+
+cons :: a -> Vector a -> Vector a
+cons = V.cons
+
+fromList :: [a] -> Vector a
+fromList = V.fromList
+
+map :: (a -> b) -> Vector a -> Vector b
+map = V.map
+
+imap :: (Integer -> a -> b) -> Vector a -> Vector b
+imap h = V.imap (\i x -> h (integer i) x)
+
+enumFromTo :: Enum a => a -> a -> Vector a
+enumFromTo = V.enumFromTo
+
+slice :: Integer -> Integer -> Vector a -> Vector a
+slice i j = V.slice (int i) (int j)
+
+foldl' :: (b -> a -> b) -> b -> Vector a -> b
+foldl' = V.foldl'
+
+zipWith :: (a -> b -> c) -> Vector a -> Vector b -> Vector c
+zipWith = V.zipWith
+
+(!) :: Vector a -> Integer -> a
+(!) v i = (V.!) v (int i)
+
+length :: Vector a -> Integer
+length = integer . V.length
+
+intLength :: Vector a -> Int 
+intLength = V.length
+
+inftyNorm :: (HasIntegers a, CanMinMaxSameType a) => Vector a -> a
+inftyNorm (v :: Vector a) =
+    V.foldl' max (convertExactly 0 :: a) v
+
+find :: (a -> Bool) -> Vector a -> Maybe a 
+find = V.find
+
+elem :: P.Eq a => a -> Vector a -> Bool
+elem = V.elem
+
+toList :: Vector a -> [a] 
+toList = V.toList
+
+zip :: Vector a -> Vector b -> Vector (a, b) 
+zip = V.zip
+
+null :: Vector a -> Bool
+null = V.null
+
+any :: (a -> Bool) -> V.Vector a -> Bool
+any = V.any
+
+all :: (a -> Bool) -> V.Vector a -> Bool
+all = V.all
+instance 
+    (HasAccuracy a, HasPrecision a) => HasAccuracy (Vector a)
+    where
+    getAccuracy v = 
+        V.foldl' max NoInformation $ V.map getAccuracy v
+
+instance 
+    (HasPrecision a) => HasPrecision (Vector a)
+    where
+    getPrecision v = 
+        if V.null v then 
+            (prec 2) 
+        else 
+            getPrecision $ v ! 0 -- TODO: safe? Alternative: V.foldl' max (prec 2) $ V.map getPrecision v 
+
+instance 
+    (CanSetPrecision a) => CanSetPrecision (Vector a)
+    where
+    setPrecision p = V.map (setPrecision p)
+
+instance 
+    (CanAddSameType a) =>
+    CanAddAsymmetric (Vector a) (Vector a)
+    where
+    type AddType (Vector a) (Vector a) = Vector a
+    add v w =
+        runST $
+        do
+        mv <- M.new (intLength v)
+        aux mv 0
+        V.freeze mv
+        where
+        lth = length v
+        aux :: (V.MVector s a) -> Integer -> (ST s ())
+        aux mv k = 
+            if k == lth then
+                return ()
+            else
+                do                
+                M.write mv (int k) (v ! k + w ! k)
+                aux mv (k + 1)
+
+instance 
+    (CanSubSameType a) =>
+    CanSub (Vector a) (Vector a)
+    where
+    type SubType (Vector a) (Vector a) = Vector a
+    sub v w =
+        runST $
+        do
+        mv <- M.new (intLength v)
+        aux mv 0
+        V.freeze mv
+        where
+        lth = length v
+        aux :: (V.MVector s a) -> Integer -> (ST s ())
+        aux mv k = 
+            if k == lth then
+                return ()
+            else
+                do                
+                M.write mv (int k) (v ! k - w ! k)
+                aux mv (k + 1)
+
+
+instance 
+    (CanAddSameType a, CanMulSameType a, HasIntegers a) =>
+    CanMulAsymmetric (Vector a) (Vector a) 
+    where
+    type MulType (Vector a) (Vector a) = a
+    mul v w =
+        runST $
+        do
+        sum <- newSTRef (convertExactly 0)
+        aux sum 0
+        readSTRef sum
+        where
+        lth = length v
+        aux :: (STRef s a) -> Integer -> (ST s ())
+        aux sum k = 
+            if k == lth then
+                return ()
+            else
+                do                
+                modifySTRef sum (\x -> x + (v ! k) * (w ! k))
+                aux sum (k + 1)
+
+instance 
+    CanMulAsymmetric (CN MPBall) (Vector (CN MPBall)) where
+    type MulType (CN MPBall) (Vector (CN MPBall)) = Vector (CN MPBall)
+    mul x v = V.map (\y -> x * y) v
diff --git a/src/AERN2/Util/Util.hs b/src/AERN2/Util/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/AERN2/Util/Util.hs
@@ -0,0 +1,14 @@
+module AERN2.Util.Util where
+
+import MixedTypesNumPrelude
+
+import AERN2.MP.Ball
+
+upperBound :: (IsInterval i) => i -> i
+upperBound = endpointRAsInterval
+
+lowerBound :: (IsInterval i) => i -> i
+lowerBound = endpointLAsInterval
+
+intersectCN :: CN MPBall -> CN MPBall -> CN MPBall
+intersectCN = intersectCNMPBall