diff --git a/Control/Monad/MC.hs b/Control/Monad/MC.hs
new file mode 100644
--- /dev/null
+++ b/Control/Monad/MC.hs
@@ -0,0 +1,14 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Control.Monad.MC
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Control.Monad.MC (
+    module Control.Monad.MC.GSL
+    ) where
+
+import Control.Monad.MC.GSL
diff --git a/Control/Monad/MC/GSL.hs b/Control/Monad/MC/GSL.hs
new file mode 100644
--- /dev/null
+++ b/Control/Monad/MC/GSL.hs
@@ -0,0 +1,271 @@
+{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances, UndecidableInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module     : Control.Monad.MC.GSL
+-- Copyright  : Copyright (c) , Patrick Perry <patperry@stanford.edu>
+-- License    : BSD3
+-- Maintainer : Patrick Perry <patperry@stanford.edu>
+-- Stability  : experimental
+--
+
+module Control.Monad.MC.GSL (
+    -- * The Monte Carlo monad
+    MC,
+    runMC,
+    evalMC,
+    execMC,
+    unsafeInterleaveMC,
+    
+    -- * The Monte Carlo monad transformer
+    MCT,
+    runMCT,
+    evalMCT,
+    execMCT,
+    liftMCT,
+    unsafeInterleaveMCT,
+
+    -- | Pure random number generator creation
+    RNG,
+    mt19937,
+
+    -- | Random distributions
+    uniform,
+    uniformInt,
+    normal,
+    poisson,
+
+    ) where
+
+import Control.Monad            ( liftM, MonadPlus(..) )
+import Control.Monad.Cont       ( MonadCont(..) )
+import Control.Monad.Error      ( MonadError(..) )
+import Control.Monad.Reader     ( MonadReader(..) )
+import Control.Monad.State      ( MonadState(..) )
+import Control.Monad.Writer     ( MonadWriter(..) )
+import Control.Monad.Trans      ( MonadTrans(..), MonadIO(..) )
+import Data.Word
+import System.IO.Unsafe         ( unsafePerformIO )
+        
+import GSL.Random.Gen hiding ( mt19937 )
+import qualified GSL.Random.Gen as Gen
+import GSL.Random.Dist
+
+-- | A Monte Carlo monad with an internal random number generator.
+newtype MC a = MC (RNG -> (a,RNG))
+
+-- | Run this Monte Carlo monad with the given initial random number generator,
+-- getting the result and the new random number generator.
+runMC :: MC a -> RNG -> (a, RNG)
+runMC (MC g) r =
+    let r' = unsafePerformIO $ cloneRNG r
+    in r' `seq` g r'
+{-# NOINLINE runMC #-}
+    
+-- | Evaluate this Monte Carlo monad and throw away the final random number
+-- generator.  Very much like @fst@ composed with @runMC@.
+evalMC :: MC a -> RNG -> a
+evalMC g r = fst $ runMC g r
+
+-- | Exicute this Monte Carlo monad and return the final random number
+-- generator.  Very much like @snd@ composed with @runMC@.
+execMC :: MC a -> RNG -> RNG
+execMC g r = snd $ runMC g r
+
+-- | Get the baton from the Monte Carlo monad without performing any
+-- computations.  Useful but dangerous.
+unsafeInterleaveMC :: MC a -> MC a
+unsafeInterleaveMC (MC m) = MC $ \r -> let
+    (a,_) = m r
+    in (a,r)
+
+
+instance Functor MC where
+    fmap f (MC m) = MC $ \r -> let
+        (a,r') = m r
+        in (f a, r')
+
+instance Monad MC where
+    return a = MC $ \r -> (a,r)
+    (MC m) >>= k =
+        MC $ \r -> let
+            (a, r') = m r
+            (MC m') = k a
+            in m' r'
+
+instance MonadState RNG MC where
+    get = MC $ getHelp 
+    put r' = MC $ putHelp r' 
+
+getHelp :: RNG -> (RNG,RNG)
+getHelp r = unsafePerformIO $ do
+    r' <- cloneRNG r
+    r' `seq` return (r',r)
+{-# NOINLINE getHelp #-}
+
+putHelp :: RNG -> RNG -> ((),RNG)
+putHelp r' r = unsafePerformIO $ do
+    io <- copyRNG r r'
+    io `seq` return ((),r)
+{-# NOINLINE putHelp #-}
+
+-- | A parameterizable Monte Carlo monad for encapsulating an inner
+-- monad.
+newtype MCT m a = MCT (RNG -> m (a,RNG))
+
+-- | Similar to 'runMC'.
+runMCT :: (Monad m) => MCT m a -> RNG -> m (a,RNG)
+runMCT (MCT g) r =
+    let r' = unsafePerformIO $ cloneRNG r
+    in r' `seq` g r'
+{-# NOINLINE runMCT #-}
+
+-- | Similar to 'evalMC'.
+evalMCT :: (Monad m) => MCT m a -> RNG -> m a
+evalMCT g r = do
+    ~(a,_) <- runMCT g r
+    return a
+    
+-- | Similar to 'execMC'.    
+execMCT :: (Monad m) => MCT m a -> RNG -> m RNG
+execMCT g r = do
+    ~(_,r') <- runMCT g r
+    return r'
+
+-- | Take a Monte Carlo computations and lift it to an MCT computation.
+liftMCT :: (Monad m) => MC a -> MCT m a
+liftMCT (MC m) = MCT $ return . m
+
+-- | Similar to 'unsafeInterleaveMC'.
+unsafeInterleaveMCT :: (Monad m) => MCT m a -> MCT m a
+unsafeInterleaveMCT (MCT g) = MCT $ \r -> do
+    ~(a,_) <- g r
+    return (a,r)
+
+instance (Monad m) => Functor (MCT m) where
+    fmap f (MCT m) = MCT $ \r -> do
+        ~(x, r') <- m r
+        return (f x, r')    
+
+instance (Monad m) => Monad (MCT m) where
+    return a = MCT $ \r -> return (a,r)
+    
+    (MCT m) >>= k =
+        MCT $ \r -> do
+            ~(a,r') <- m r
+            let (MCT m') = k a
+            m' r'
+            
+    fail str = MCT $ \_ -> fail str
+
+instance (MonadPlus m) => MonadPlus (MCT m) where
+    mzero = MCT $ \_ -> mzero
+        
+    (MCT m) `mplus` (MCT n) = 
+        MCT $ \r ->
+            let r' = unsafePerformIO $ cloneRNG r
+            in r' `seq` (m r `mplus` n r')
+
+instance (Monad m) => MonadState RNG (MCT m) where
+    get    = MCT $ return . getHelp 
+    put r' = MCT $ return . (putHelp r')
+
+instance MonadTrans MCT where
+    lift m = MCT $ \r -> do
+        a <- m
+        return (a,r)
+
+instance (MonadCont m) => MonadCont (MCT m) where
+    callCC f = MCT $ \r ->
+        callCC $ \c ->
+        let (MCT m) = (f (\a -> MCT $ \r' -> c (a, r'))) 
+        in m r
+
+instance (MonadError e m) => MonadError e (MCT m) where
+    throwError              = lift . throwError
+    (MCT m) `catchError` h = MCT $ \r -> 
+        m r `catchError` \e -> let (MCT m') = h e in m' r
+
+instance (MonadIO m) => MonadIO (MCT m) where
+    liftIO = lift . liftIO
+
+instance (MonadReader r m) => MonadReader r (MCT m) where
+    ask              = lift ask
+    local f (MCT m) = MCT $ \r ->
+        local f (m r)
+
+instance (MonadState s m) => MonadState s (MCT m) where
+    get = lift get 
+    put = lift . put
+
+instance (MonadWriter w m) => MonadWriter w (MCT m) where
+    tell            = lift . tell
+    listen (MCT m) = MCT $ \r -> do
+        ~((a,r'),w) <- listen (m r)
+        return ((a,w),r')
+    pass (MCT m) = MCT $ \r -> pass $ do
+        ~((a,f),r') <- m r
+        return ((a,r'),f)
+
+-- | Get a Mersenne Twister random number generator seeded with the given
+-- value.
+mt19937 :: Word64 -> RNG
+mt19937 s = unsafePerformIO $ do
+    r <- newRNG Gen.mt19937
+    setSeed r s
+    return r
+{-# NOINLINE mt19937 #-}
+
+-- | @uniformInt n@ generates an integer uniformly in the range @[0,n-1]@.
+-- It is an error to call this function with a non-positive value.
+uniformInt :: Int -> MC Int
+uniformInt n = MC $ uniformIntHelp n
+
+uniformIntHelp :: Int -> RNG -> (Int,RNG)
+uniformIntHelp n r = unsafePerformIO $ do
+    x <- getUniformInt r n
+    x `seq` return (x,r)
+
+-- | @uniform  a b@ generates a value uniformly distributed in @[a,b)@.
+uniform :: Double -> Double -> MC Double
+uniform a b = MC $ uniformHelp a b
+
+uniformHelp :: Double -> Double -> RNG -> (Double,RNG)
+uniformHelp a b r = unsafePerformIO $ do
+    x <- getFlat r a b
+    x `seq` return (x,r)
+{-# NOINLINE uniformHelp #-}
+    
+-- | @normal mu sigma@ generates a Normal random variable with mean
+-- @mu@ and standard deviation @sigma@.
+normal :: Double -> Double -> MC Double
+normal mu sigma = MC $ normalHelp mu sigma
+
+normalHelp :: Double -> Double -> RNG -> (Double,RNG)
+normalHelp mu sigma r = unsafePerformIO $ do
+    x <- liftM (mu +) $ getGaussian r sigma
+    x `seq` return (x,r)
+{-# NOINLINE normalHelp #-}
+
+-- | @poisson mu@ generates a Poisson random variable with mean @mu@.
+poisson :: Double -> MC Int
+poisson mu = MC $ poissonHelp mu
+
+poissonHelp :: Double -> RNG -> (Int,RNG)
+poissonHelp mu r = unsafePerformIO $ do
+    x <- getPoisson r mu
+    x `seq` return (x,r)
+{-# NOINLINE poissonHelp #-}
+
+
+{-
+
+
+unifInt :: (Monad m) => Int -> MCT m Int
+unifInt n = MCT $ unifInt' n
+
+unifInt' :: (Monad m) => Int -> RNG -> m (Int,RNG)
+unifInt' n r =
+    unsafePerformIO $ do
+        i <- rngUnifInt r n
+        i `seq` (return . return) (i,r)
+-}
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) Patrick Perry <patperry@stanford.edu> 2008
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+1. Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the author nor the names of his contributors
+   may be used to endorse or promote products derived from this software
+   without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``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 AUTHORS OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGE.
diff --git a/Setup.lhs b/Setup.lhs
new file mode 100644
--- /dev/null
+++ b/Setup.lhs
@@ -0,0 +1,6 @@
+#!/usr/bin/env runhaskell
+> import Distribution.Simple
+> import System.Cmd
+>
+> main = defaultMain
+>
diff --git a/examples/Pi.hs b/examples/Pi.hs
new file mode 100644
--- /dev/null
+++ b/examples/Pi.hs
@@ -0,0 +1,77 @@
+
+import Control.Monad.MC
+import Control.Monad
+import System.Environment( getArgs )
+import Text.Printf( printf )
+
+-- | Generate a point in the box [-1,1) x [-1,1)
+unitBox :: MC (Double,Double)
+unitBox = liftM2 (,) (uniform (-1) 1) 
+                     (uniform (-1) 1)
+
+-- | Indicates whether or not a point is in the unit circle
+inUnitCircle :: (Double,Double) -> Bool
+inUnitCircle (x,y) = x*x + y*y <= 1
+
+-- | Generate @n@ points in the unit box and count how many are in the
+-- unit circle.
+countInBox :: Int -> MC Int
+countInBox n = do
+    xs <- replicateM n unitBox
+    return $ count inUnitCircle xs
+
+-- | Count how many times the predicate is true
+count :: (a -> Bool) -> [a] -> Int
+count f = length . filter f
+
+-- | Compute a Monte Carlo estimate of pi based on @n@ samples.  Return
+-- the estimate and the standard error of the estimate.
+computePi :: Int -> MC (Double,Double)
+computePi n = do
+    m  <- countInBox n
+    let p  = toDouble m / toDouble n
+        se = sqrt (p * (1 - p) / toDouble n)
+    return (4*p, 4*se)
+  where
+    toDouble = realToFrac . toInteger
+
+-- | Given an estimate and standard error, produce a 99% confidence
+-- interval based on the Central Limit Theorem
+interval :: Double -> Double -> (Double,Double)
+interval mu se = let
+    delta = 2.575*se
+    in (mu-delta, mu+delta)
+
+-- | Tests if the value is in the interval [a,b]
+inInterval :: Double -> (Double,Double) -> Bool
+x `inInterval` (a,b) = x >= a && x <= b
+
+-- | Compute an estimate of pi based on @n@ points and see if the true
+-- value is in the confidence interval
+covers :: Int -> MC Bool
+covers n = do
+    (mu,se) <- computePi n
+    return $ pi `inInterval` (interval mu se)
+
+-- | Compute @r@ estimates of pi based on @n@ samples each, and count
+-- how many times the true values is included in the 99% confidence
+-- inverval
+coverage :: Int -> Int -> MC Int
+coverage r n = do
+    liftM (count id) $ replicateM r (covers n)
+    
+    
+main = do
+    [n] <- map read `fmap` getArgs
+    main' n
+    
+main' n = let
+    seed   = 0
+    (mu,se) = evalMC (computePi n) $ mt19937 seed
+    (l,u)   = interval mu se
+    r       = 500
+    c       = evalMC (coverage r n) $ mt19937 seed
+    in do
+        printf "Estimate from one simulation: %g\n" mu
+        printf "99%% Confidence Interval:    (%g,%g)\n" l u
+        printf "\nOf %d intervals, %d contain the true value.\n" r c
diff --git a/monte-carlo.cabal b/monte-carlo.cabal
new file mode 100644
--- /dev/null
+++ b/monte-carlo.cabal
@@ -0,0 +1,32 @@
+name:            monte-carlo
+version:         0.1
+homepage:        http://stat.stanford.edu/~patperry/code/monte-carlo
+synopsis:        A monad and transformer for Monte Carlo calculations.
+description:
+    A monad and transformer for Monte Carlo calculations.  The monads
+    carry and provide access to a random number generator.  Importantly,
+    they only keep one copy of the generator state, and so are much more
+    efficient than MonadRandom.  Currently, only the generator that comes
+    with the GNU Scientific Library is supported.
+    .
+category:        Math
+license:         BSD3
+license-file:    LICENSE
+copyright:       (c) 2008. Patrick Perry <patperry@stanford.edu>
+author:          Patrick Perry
+maintainer:      Patrick Perry <patperry@stanford.edu>
+cabal-version: >= 1.2.0
+build-type:      Simple
+tested-with:     GHC ==6.8.3
+
+extra-source-files:     examples/Pi.hs
+
+library
+    exposed-modules:    Control.Monad.MC
+                        Control.Monad.MC.GSL
+                        
+    ghc-options:        -Wall
+    extensions:         MultiParamTypeClasses, FlexibleInstances, 
+                        UndecidableInstances 
+
+    build-depends:      base, mtl, gsl-random