diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Bogdan Penkovsky (c) 2017
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Author name here nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/MEP.hs b/MEP.hs
new file mode 100644
--- /dev/null
+++ b/MEP.hs
@@ -0,0 +1,39 @@
+{- |
+Copyright Bogdan Penkovsky (c) 2017
+
+= Multiple Expression Programming
+
+-}
+
+module MEP (
+  Chromosome (..)
+  , Gene
+  , Population
+  , Phenotype
+  , Config (..)
+  , defaultConfig
+  , LossFunction
+
+  -- * Genetic algorithm
+  , initialize
+  , evaluateGeneration
+  , evolve
+  , binaryTournament
+  , crossover
+  , mutation3
+  , smoothMutation
+  , newChromosome
+
+  -- * Random
+  , Rand
+  , newPureMT
+  , runRandom
+  , evalRandom
+  ) where
+
+import System.Random.Mersenne.Pure64 ( newPureMT )
+
+import MEP.Types
+import MEP.Operators
+import MEP.Run
+import MEP.Random
diff --git a/MEP/Operators.hs b/MEP/Operators.hs
new file mode 100644
--- /dev/null
+++ b/MEP/Operators.hs
@@ -0,0 +1,292 @@
+-- |
+-- = Genetic operators
+
+module MEP.Operators (
+  Config (..)
+  , defaultConfig
+  , LossFunction
+  , Phenotype
+  -- * Genetic operators
+  , initialize
+  , evaluateGeneration
+  , evolve
+  , phenotype
+  , binaryTournament
+  , crossover
+  , mutation3
+  , smoothMutation
+  , newChromosome
+  ) where
+
+import           Data.Vector ( Vector )
+import qualified Data.Vector as V
+import           Data.List
+                 ( nub
+                 , sortBy
+                 )
+import           Data.Ord ( comparing )
+import qualified Control.Monad as CM
+
+import           MEP.Random
+import           MEP.Types
+import           MEP.Run ( evaluate )
+
+data Config a = Config
+  {
+    p'const :: Double        -- ^ Probability of constant generation
+    , p'var :: Double        -- ^ Probability of variable generation.
+                             -- The probability of operator generation is inferred
+                             -- automatically as @1 - p'const - p'var@.
+    , p'mutation :: Double   -- ^ Mutation probability
+    , p'crossover :: Double  -- ^ Crossover probability
+
+    , c'length :: Int        -- ^ The chromosome length
+    , c'popSize :: Int       -- ^ A (sub)population size
+    , c'popN :: Int          -- ^ Number of subpopulations (1 or more)  [not implemented]
+    , c'ops :: Vector (F a)  -- ^ Functions pool with their symbolic
+                             -- representations
+    , c'vars :: Int          -- ^ The input dimensionality
+  }
+
+-- |
+-- @
+-- defaultConfig = Config
+--   {
+--     p'const = 0.1
+--   , p'var = 0.4
+--   , p'mutation = 0.1
+--   , p'crossover = 0.9
+--
+--   , c'length = 50
+--   , c'popSize = 100
+--   , c'popN = 1
+--   , c'ops = V.empty  -- <-- To be overridden
+--   , c'vars = 1
+--   }
+-- @
+defaultConfig :: Config Double
+defaultConfig = Config
+  {
+    p'const = 0.1
+  , p'var = 0.4
+  , p'mutation = 0.1
+  , p'crossover = 0.9
+
+  , c'length = 50
+  , c'popSize = 100
+  , c'popN = 1
+  , c'ops = V.empty
+  , c'vars = 1
+  }
+
+-- | A function to minimize.
+--
+-- The argument is a vector evaluation function whose input
+-- is a vector (length @c'vars@) and ouput is
+-- a vector with a different length @c'length@.
+--
+-- The result is a vector of the best indices
+-- and a scalar loss value.
+type LossFunction a =
+  ((V.Vector a -> V.Vector a) -> (V.Vector Int, Double))
+
+-- | Evaluates a chromosome according to the given
+-- loss function.
+phenotype
+  :: Num a =>
+     LossFunction a
+     -> Chromosome a
+     -> Phenotype a
+phenotype loss chr = let (is, val) = loss (evaluate chr)
+                     in (val, chr, is)
+
+type Phenotype a = (Double, Chromosome a, V.Vector Int)
+
+-- | Randomly generate a new population
+initialize :: Config Double -> Rand (Population Double)
+initialize c@Config { c'popSize = size } = mapM (\_ -> newChromosome c) [1..size]
+
+evaluateGeneration
+  :: Num a =>
+     LossFunction a
+     -> [Chromosome a]
+     -> [Phenotype a]
+evaluateGeneration loss pop = map (phenotype loss) pop
+
+-- | Selection operator that produces the next evaluated population.
+--
+-- Standard algorithm: the best offspring O replaces the worst
+-- individual W in the current population if O is better than W.
+evolve
+  ::
+     Config Double
+     -- ^ Common configuration
+     -> LossFunction Double
+     -- ^ Custom loss function
+     -> (Chromosome Double -> Rand (Chromosome Double))
+     -- ^ Mutation
+     -> (Chromosome Double -> Chromosome Double -> Rand (Chromosome Double, Chromosome Double))
+     -- ^ Crossover
+     -> ([Phenotype Double] -> Rand (Chromosome Double))
+     -- ^ A chromosome selection algorithm. Does not need to be random, but may be.
+     -> [Phenotype Double]
+     -- ^ Evaluated population
+     -> Rand [Phenotype Double]
+     -- ^ New generation
+evolve c loss mut cross select phenotypes = do
+  let pc = p'crossover c
+      pm = p'mutation c
+      -- Sort in decreasing @val@ order so that
+      -- the worst (with the biggest loss) is in the head
+      sort' = sortBy (comparing (\(val, _, _) -> negate val))
+
+      ev phen0 _ = do
+        chr1 <- select phen0
+        chr2 <- select phen0
+        (of1, of2) <- cross chr1 chr2
+        of1' <- withProbability pm mut of1
+        of2' <- withProbability pm mut of2
+        let r1@(val1, _, _) = phenotype loss of1'
+            r2@(val2, _, _) = phenotype loss of2'
+            (worstVal, _, _) = head phen0
+            phen' | val1 < worstVal = r1 : tail phen0
+                  | val2 < worstVal = r2 : tail phen0
+                  -- No change
+                  | otherwise = phen0
+        let phen1 = sort' phen'
+        return phen1
+
+  pop' <- CM.foldM ev (sort' phenotypes) [1..c'popSize c `div` 2]
+  return pop'
+
+-- | Binary tournament selection
+binaryTournament :: Ord a => [Phenotype a] -> Rand (Chromosome a)
+binaryTournament phen = do
+  (val1, cand1, _) <- draw $ V.fromList phen
+  (val2, cand2, _) <- draw $ V.fromList phen
+  if val1 < val2
+    then return cand1
+    else return cand2
+
+-- | Uniform crossover operator
+crossover ::
+  Chromosome a
+  -> Chromosome a
+  -> Rand (Chromosome a, Chromosome a)
+crossover ca cb = do
+  r <- V.zipWithM (curry (swap 0.5)) ca cb
+  return $ V.unzip r
+
+swap :: Double -> (t, t) -> Rand (t, t)
+swap p = withProbability p (\(a, b) -> return (b, a))
+
+replaceAt :: Int -> a -> Vector a -> Vector a
+replaceAt i gene chr0 =
+  let (c1, c2) = V.splitAt i chr0
+  in c1 V.++ V.singleton gene V.++ V.tail c2
+
+-- | Mutation operator with up to three mutations per chromosome
+mutation3 ::
+  Config Double
+  -- ^ Common configuration
+  -> Chromosome Double
+  -> Rand (Chromosome Double)
+mutation3 c chr = do
+                                      -- Subtract 1 to get a non-zero head to
+                                      -- replace
+  is <- nub <$> CM.replicateM k (getMaxInt (chrLen - 1))
+  genes <- mapM new' is
+  let chr' = foldr (uncurry replaceAt)
+                   chr
+                   (zip is genes)
+  return chr'
+    where chrLen = V.length chr
+          k = 3
+          new' = new (p'const c) (p'var c) (c'vars c) (c'ops c)
+
+-- | Mutation operator with a fixed mutation probability
+-- of each gene
+smoothMutation
+  ::
+     Double
+     -- ^ Probability of gene mutation
+     -> Config Double
+     -- ^ Common configuration
+     -> Chromosome Double
+     -> Rand (Chromosome Double)
+smoothMutation p c chr =
+  let new' = new (p'const c) (p'var c) (c'vars c) (c'ops c)
+      mutate i = withProbability p (\_ -> new' i)
+  in V.zipWithM mutate (V.enumFromN 0 (V.length chr)) chr
+
+-- | Randomly initialize a new chromosome.
+-- By definition, the first gene is terminal (a constant
+-- or a variable).
+newChromosome ::
+  Config Double          -- ^ Common configuration
+  -> Rand (Chromosome Double)
+newChromosome c = do
+  let pConst = p'const c
+      pVar = p'var c
+  V.mapM (new pConst pVar (c'vars c) (c'ops c)) $ V.enumFromN 0 (c'length c)
+
+-- | Produce a new random gene
+new ::
+  Double    -- ^ Probability to produce a constant
+  -> Double    -- ^ Probability to produce a variable
+  -> Int       -- ^ Number of input variables
+  -> Vector (F Double)   -- ^ Operations vector
+  -> Int                 -- ^ Maximal operation index
+  -> Rand (Gene Double Int)
+new p1 p2 vars ops maxIndex = if maxIndex == 0
+  -- The head must be a terminal
+  -- p1' = p1 + (1 - p1 - p2) / 2 = 1/2 + p1/2 - p2/2
+  then let p1' = 0.5 * (1 + p1 - p2)
+       in newTerminal p1' vars
+  else do
+    p' <- getDouble
+    let sel | p' < p1 = newC
+            | p' < (p1 + p2) = newVar vars
+            | otherwise = newOp ops maxIndex
+    sel
+
+newTerminal ::
+  Double        -- ^ Probability @p@ of a constant generation.
+                   -- @1-p@ will be the probability of a variable generation.
+  -> Int           -- ^ Number of input variables
+  -> Rand (Gene Double i)
+newTerminal p vars = do
+  p' <- getDouble
+  if p' < p
+    then newC
+    else newVar vars
+
+-- | A randomly generated variable identifier
+newVar :: Int -> Rand (Gene a i)
+newVar vars = do
+  var <- draw $ V.enumFromN 0 vars
+  return $ Var var
+
+-- | A random operation from the operations vector
+newOp
+  :: Vector (F a)
+  -> Int
+  -> Rand (Gene a Int)
+newOp ops maxIndex = do
+  op <- draw ops
+  i1 <- getMaxInt maxIndex
+  i2 <- getMaxInt maxIndex
+  return $ Op op i1 i2
+
+-- | Draw a constant from the normal distribution
+newCNormal
+  :: Double  -- ^ Mean
+  -> Double  -- ^ Std deviation
+  -> Rand (Gene Double i)
+newCNormal mu sigma = do
+  n <- getNormal
+  return $ C (mu + sigma*n)
+
+-- | Draw a constant from the uniform distribution
+newC :: Rand (Gene Double i)
+newC = C <$> getDouble
diff --git a/MEP/Random.hs b/MEP/Random.hs
new file mode 100644
--- /dev/null
+++ b/MEP/Random.hs
@@ -0,0 +1,55 @@
+module MEP.Random
+    (
+    -- * Utilities
+    draw
+    , getNormal
+    , getMaxInt
+    , withProbability
+
+    -- * Re-exports
+    , getBool, getInt, getWord, getDouble
+    , runRandom, evalRandom
+    , Rand, Random
+    ) where
+
+import Control.Monad.Mersenne.Random
+import Data.Complex (Complex (..))
+import System.Random
+import Data.Vector as V
+
+-- | Randomly draw an element from a vector
+draw :: Vector a -> Rand a
+draw xs =
+  Rand $ \g -> let (n, g') = randomR (0, V.length xs - 1) g
+                   r = xs V.! n
+               in R r g'
+
+-- | Modify value with probability @p@
+withProbability
+  :: Double         -- ^ The probability @p@
+  -> (a -> Rand a)  -- ^ Modification function
+  -> (a -> Rand a)
+withProbability p modify x = do
+  t <- getDouble
+  if t < p
+     then modify x
+     else return x
+
+-- | Randomly generate Int between 0 and @n@.
+-- Should be strictly less than n if n > 1
+-- or zero otherwise. Therefore, getMaxInt 1
+-- should be always 0.
+getMaxInt :: Int  -- ^ @n@
+  -> Rand Int
+getMaxInt n = do
+  r <- getDouble
+  return $ floor (r * fromIntegral n)
+
+getNormal :: Rand Double
+getNormal = do
+  -- Box-Muller method
+  u <- getDouble
+  v <- getDouble
+  let (c :+ s) = exp (0 :+ (2*pi*v))
+  let r = sqrt $ (-2) * log u
+  return $ r*c
diff --git a/MEP/Run.hs b/MEP/Run.hs
new file mode 100644
--- /dev/null
+++ b/MEP/Run.hs
@@ -0,0 +1,45 @@
+{-# LANGUAGE BangPatterns #-}
+
+module MEP.Run where
+
+import qualified Data.Vector as V
+import qualified Data.Vector.Mutable as VM
+import           System.IO.Unsafe ( unsafePerformIO )
+
+import           MEP.Types
+
+-- | Evaluate each subexpression in a chromosome
+evaluate :: Num a
+         => Chromosome a  -- ^ Chromosome to evaluate
+         -> V.Vector a    -- ^ Variable values
+         -> V.Vector a    -- ^ Resulting vector of multiple evaluations
+evaluate chr vmap = unsafePerformIO $ do
+  -- Use dynamic programming to evaluate the chromosome
+  v <- VM.new chrLen
+
+  let -- Gene evaluation function
+      _f (C c) _ = return c
+      _f (Var n) _ = return $ vmap V.! n
+      _f (Op (_, f) i1 i2) v' = do
+        !r1 <- v' `VM.read` i1
+        !r2 <- v' `VM.read` i2
+        let !r = f r1 r2
+        return r
+
+      -- Chromosome evaluation
+      go !v' !j =
+        if j == chrLen
+           then return ()
+           else do
+             val <- _f (chr V.! j) v'
+             VM.write v' j val
+             go v' (j + 1)
+
+  go v 0
+
+  V.unsafeFreeze v
+    where chrLen = V.length chr
+{-# SPECIALIZE
+  evaluate :: Chromosome Double
+           -> V.Vector Double
+           -> V.Vector Double #-}
diff --git a/MEP/Types.hs b/MEP/Types.hs
new file mode 100644
--- /dev/null
+++ b/MEP/Types.hs
@@ -0,0 +1,30 @@
+{- | Provide the basic MEP data structures
+ -}
+{-# LANGUAGE GADTs #-}
+module MEP.Types where
+
+import qualified Data.Vector as V
+
+
+type Population a = [Chromosome a]
+
+-- | A chromosome is a vector of genes
+type Chromosome a = V.Vector (Gene a Int)
+
+-- | Either a terminal symbol or a three-address code (a function
+-- and two pointers)
+data Gene a i where
+  -- Terminal symbol: constant
+  C :: a -> Gene a i
+  -- Terminal symbol: variable
+  Var :: Int -> Gene a i
+  -- Operation
+  Op :: F a -> i -> i -> Gene a i
+
+instance (Show a, Show i) => Show (Gene a i) where
+  show (C c) = show c
+  show (Var n) = "v" ++ show n
+  show (Op (s, _) i1 i2) = show s ++ " " ++ show i1 ++ " " ++ show i2
+
+-- | A function and its symbolic representation
+type F a = (Char, a -> a -> a)
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,34 @@
+# Multi Expression Programming
+
+You say, Haskell has not enough machine learning libraries?
+
+Here is yet another one!
+
+## History
+
+There exist many other Genetic Algorithm (GA) Haskell packages.
+Personally I have used
+[simple genetic algorithm](http://hackage.haskell.org/package/moo),
+[GA](http://hackage.haskell.org/package/moo),
+and [moo](http://hackage.haskell.org/package/moo) for quite a long time.
+The last package was the most preferred, but the other two are
+also great.
+
+However, when I came up with this
+[MEP paper](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.5.4352&rep=rep1&type=pdf),
+to my surprise there was no MEP realization in Haskell.
+Soon I realized that existing GA packages are limited,
+and it would be more efficient to implement MEP from scratch.
+
+That is how this package was started. I also wish to say thank you
+to the authors of the [moo](http://hackage.haskell.org/package/moo) 
+GA library, which inspired the present 
+[hmep](http://github.com/masterdezign/hmep) package.
+
+## About MEP
+
+Multi Expression Programming is a genetic programming variant encoding multiple
+solutions in the same chromosome. A chromosome is a computer program.
+Each gene is featuring [code reuse](https://en.wikipedia.org/wiki/Code_reuse).
+For more details, please check http://mepx.org/papers.html and
+https://en.wikipedia.org/wiki/Multi_expression_programming.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/TODO b/TODO
new file mode 100644
--- /dev/null
+++ b/TODO
@@ -0,0 +1,5 @@
+1. Provide a Show instance for AI.MEP.Types.Chromosome
+   such that a Haskell code is generated.
+
+2. Provide a Storable instance for AI.MEP.Types.Gene
+   and make the Chromosome a Data.Vector.Storable.
diff --git a/app/Main.hs b/app/Main.hs
new file mode 100644
--- /dev/null
+++ b/app/Main.hs
@@ -0,0 +1,78 @@
+module Main where
+
+import qualified Data.Vector as V
+import           Data.List ( foldl' )
+import           Control.Monad ( foldM )
+import           Numeric.LinearAlgebra
+                 ( randomVector
+                 , RandDist( Uniform )
+                 , toList
+                 )
+
+import MEP
+
+ops = V.fromList [('*', (*)), ('+', (+)), ('/', (/)), ('-', (-))]
+
+config = defaultConfig {
+  c'ops = ops
+  , c'length = 50
+  }
+
+seed :: Int
+seed = 3
+
+randDomain :: Int -> [Double]
+randDomain = map (subtract pi. (2*pi *)). toList. randomVector seed Uniform
+
+dataset1 :: V.Vector (Double, Double)
+dataset1 = V.map (\x -> (x, sin x)) $ V.fromList $ randDomain nSamples
+  where nSamples = 50
+
+dist x y = abs $ x - y
+
+loss :: LossFunction Double
+loss evalf = (V.singleton i', loss')
+  where
+    (xs, ys) = unzip $ V.toList dataset1
+    -- Distances resulting from multiple expression evaluation
+    dss = zipWith (\x y -> V.map (dist y). evalf. V.singleton $ x) xs ys
+    -- Cumulative distances for each index
+    dcumul = sum' dss
+    -- Select index minimizing cumulative distances
+    i' = V.minIndex dcumul
+    -- The loss value with respect to the index of the best expression
+    loss' = dcumul V.! i'
+
+-- Could be optimized
+sum' :: Num a => [V.Vector a] -> V.Vector a
+sum' xss = foldl' (V.zipWith (+)) base xss
+  where
+    len = V.length $ head xss
+    base = V.replicate len 0
+
+nextGeneration
+  :: [Phenotype Double] -> Rand [Phenotype Double]
+nextGeneration = evolve config loss (mutation3 config) crossover binaryTournament
+
+avgLoss :: [Phenotype Double] -> Double
+avgLoss xs =
+  let (r, len) = foldl' (\(c, i) (val, _, _) -> (c + val, i + 1)) (0, 0) xs
+  in r / (fromIntegral len)
+
+runIO (pop, g') i = do
+  let (newPop, g2) = foldr (\_ xg -> run xg) (pop, g') [1..generations]
+  putStrLn $ "Population " ++ show (i * generations) ++ ": average loss " ++ show (avgLoss newPop)
+  return (newPop, g2)
+    where
+      run (x, g) = runRandom (nextGeneration x) g
+      generations = 40
+
+main :: IO ()
+main = do
+  g <- newPureMT
+  let (pop, g') = runRandom (initialize config) g
+      popEvaluated = evaluateGeneration loss pop
+  putStrLn $ "Average loss in the initial population " ++ show (avgLoss popEvaluated)
+
+  (final, _) <- foldM runIO (popEvaluated, g') [1..100]
+  print $ last final
diff --git a/hmep.cabal b/hmep.cabal
new file mode 100644
--- /dev/null
+++ b/hmep.cabal
@@ -0,0 +1,61 @@
+name:                hmep
+version:             0.0.0
+synopsis:            HMEP Multi Expression Programming –
+                     a genetic programming variant
+description:         A multi expression programming implementation with
+                     focus on speed.
+                     .
+                     https://en.wikipedia.org/wiki/Multi_expression_programming
+homepage:            https://github.com/masterdezign/hmep#readme
+license:             BSD3
+license-file:        LICENSE
+author:              Bogdan Penkovsky
+maintainer:          dev at penkovsky dot com
+copyright:           2017 Bogdan Penkovsky
+category:            AI
+build-type:          Simple
+extra-source-files:  README.md TODO
+cabal-version:       >=1.22
+
+library
+  exposed-modules:     MEP
+                     , MEP.Run
+                     , MEP.Types
+  other-modules:       MEP.Random
+                     , MEP.Operators
+  build-depends:       base >= 4.7 && < 5
+                     , containers
+                     , monad-mersenne-random
+                     , mersenne-random-pure64
+                     , random
+                     , vector
+  default-language:    Haskell2010
+
+executable hmep-demo
+  hs-source-dirs:      app
+  main-is:             Main.hs
+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+  build-depends:       base
+                     , containers
+                     , hmatrix
+                     , mersenne-random-pure64
+                     , monad-mersenne-random
+                     , vector
+                     , hmep
+  default-language:    Haskell2010
+
+test-suite hmep-test
+  type:                exitcode-stdio-1.0
+  hs-source-dirs:      test
+  main-is:             Spec.hs
+  build-depends:       base
+                     , containers
+                     , HUnit
+                     , vector
+                     , hmep
+  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+  default-language:    Haskell2010
+
+source-repository head
+  type:     git
+  location: https://github.com/masterdezign/hmep
diff --git a/test/Spec.hs b/test/Spec.hs
new file mode 100644
--- /dev/null
+++ b/test/Spec.hs
@@ -0,0 +1,33 @@
+import           Test.HUnit
+import           System.Exit ( exitSuccess
+                             , exitFailure )
+import qualified Data.Vector as V
+
+import           MEP.Types
+import           MEP.Run
+
+o'mult :: Num a => F a
+o'mult = ('*', (*))
+{-# SPECIALIZE o'mult :: F Double #-}
+
+-- Encodes x, x^2, x^4, x^8
+pow8Int :: Chromosome Int
+pow8Int = V.fromList [Var 0, Op o'mult 0 0, Op o'mult 1 1, Op o'mult 2 2]
+
+pow8 :: Chromosome Double
+pow8 = V.fromList [Var 0, Op o'mult 0 0, Op o'mult 1 1, Op o'mult 2 2]
+
+testEvaluate = test [
+  "2^8" ~: V.fromList [2, 4, 16, 256] ~=? evaluate pow8Int (V.singleton (2 :: Int)),
+
+  "2.5^8" ~: V.fromList [2.5,6.25,39.0625,1525.87890625] ~=? evaluate pow8 (V.singleton (2.5 :: Double))
+  ]
+
+allTests = TestList [ testEvaluate ]
+
+main :: IO ()
+main = do
+  result <- runTestTT allTests
+  if (errors result + failures result) > 0
+    then exitFailure
+    else exitSuccess