diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2010, Edward Kmett
+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 Edward Kmett 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/Numeric/RAD.hs b/Numeric/RAD.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/RAD.hs
@@ -0,0 +1,446 @@
+{-# LANGUAGE Rank2Types, TypeFamilies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.RAD
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Reverse Mode Automatic Differentiation via overloading to perform
+-- nonstandard interpretation that replaces original numeric type with
+-- a bundle that contains a value of the original type and the tape that
+-- will be used to recover the value of the sensitivity.
+-- 
+-- This package uses StableNames internally to recover sharing information from 
+-- the tape to avoid combinatorial explosion, and thus runs asymptotically faster
+-- than it could without such sharing information, but the use of side-effects
+-- contained herein is benign.
+--
+-- The API has been built to be close to the design of 'Numeric.FAD' from the 'fad' package
+-- by Barak Pearlmutter and Jeffrey Mark Siskind and contains portions of that code, with minor liberties taken.
+-- 
+-----------------------------------------------------------------------------
+
+module Numeric.RAD 
+    ( 
+    -- * First-Order Reverse Mode Automatic Differentiation
+      RAD
+    , lift
+    -- * First-Order Differentiation Operators
+    , diffUU
+    , diffUF
+    , diff2UU
+    , diff2UF
+    -- * Common access patterns 
+    , diff
+    , diff2
+    , jacobian
+    , jacobian2
+    , grad
+    , grad2
+    -- * Optimization Routines 
+    , zeroNewton
+    , inverseNewton
+    , fixedPointNewton
+    , extremumNewton
+    , argminNaiveGradient
+    ) where
+
+import Prelude hiding (mapM)
+import Control.Applicative (Applicative(..),(<$>))
+import Control.Monad.ST
+import Control.Monad (forM_)
+import Data.List (foldl')
+import Data.Array.ST
+import Data.Array
+import Data.Ix
+import Text.Show
+import Data.Graph (graphFromEdges', topSort, Vertex)
+import Data.Reify (reifyGraph, MuRef(..))
+import qualified Data.Reify.Graph as Reified
+import Data.Traversable (Traversable, mapM)
+import System.IO.Unsafe (unsafePerformIO)
+
+newtype RAD s a = RAD (Tape a (RAD s a))
+
+data Tape a t
+    = C a 
+    | V a Int
+    | B a a a t t
+    | U a a t 
+
+instance Show a => Show (RAD s a) where
+    showsPrec d = disc1 (showsPrec d)
+
+-- | The 'lift' function injects a primal number into the RAD data type with a 0 derivative.
+-- If reverse-mode AD numbers formed a monad, then 'lift' would be 'return'.
+lift :: a -> RAD s a 
+lift = RAD . C 
+{-# INLINE lift #-}
+
+primal :: RAD s a -> a
+primal (RAD (C y)) = y
+primal (RAD (V y _)) = y
+primal (RAD (B y _ _ _ _)) = y
+primal (RAD (U y _ _)) = y
+{-# INLINE primal #-}
+
+var :: a -> Int -> RAD s a 
+var a v = RAD (V a v)
+
+-- TODO: A higher-order data-reify
+-- mapDeRef :: (Applicative f) => (forall a . Num a => RAD s a -> f (u a)) -> a -> f (Tape a (u a))
+instance MuRef (RAD s a) where
+    type DeRef (RAD s a) = Tape a
+    mapDeRef f (RAD (C a)) = pure (C a)
+    mapDeRef f (RAD (V a v)) = pure (V a v)
+    mapDeRef f (RAD (B a jb jc x1 x2)) = B a jb jc <$> f x1 <*> f x2
+    mapDeRef f (RAD (U a j x)) = U a j <$> f x
+
+on :: (a -> a -> c) -> (b -> a) -> b -> b -> c
+on f g a b = f (g a) (g b)
+
+instance Eq a =>  Eq (RAD s a) where
+    (==) = (==) `on` primal
+
+instance Ord a => Ord (RAD s a) where
+    compare = compare `on` primal
+
+instance Bounded a => Bounded (RAD s a) where
+    maxBound = lift maxBound
+    minBound = lift minBound
+
+unary_ :: (a -> a) -> a -> RAD s a -> RAD s a
+unary_ f _ (RAD (C b)) = RAD (C (f b))
+unary_ f g b = RAD (U (disc1 f b) g b)
+{-# INLINE unary_ #-}
+
+unary :: (a -> a) -> (a -> a) -> RAD s a -> RAD s a
+unary f _ (RAD (C b)) = RAD (C (f b))
+unary f g b = RAD (U (disc1 f b) (disc1 g b) b)
+{-# INLINE unary #-}
+
+binary_ :: (a -> a -> a) -> a -> a -> RAD s a -> RAD s a -> RAD s a
+binary_ f _ _ (RAD (C b)) (RAD (C c)) = RAD (C (f b c))
+binary_ f gb gc b c = RAD (B (f vb vc) gb gc b c)
+    where vb = primal b; vc = primal c
+{-# INLINE binary_ #-}
+
+-- binary_ with partials
+binary :: (a -> a -> a) -> (a -> a) -> (a -> a) -> RAD s a -> RAD s a -> RAD s a
+binary f _ _ (RAD (C b)) (RAD (C c)) = RAD (C (f b c))
+binary f gb gc b c = RAD (B (f vb vc) (gb vc) (gc vb) b c)
+    where vb = primal b; vc = primal c
+{-# INLINE binary #-}
+
+disc1 :: (a -> b) -> RAD s a -> b
+disc1 f x = f (primal x)
+{-# INLINE disc1 #-}
+
+disc2 :: (a -> b -> c) -> RAD s a -> RAD s b -> c
+disc2 f x y = f (primal x) (primal y)
+{-# INLINE disc2 #-}
+
+disc3 :: (a -> b -> c -> d) -> RAD s a -> RAD s b -> RAD s c -> d
+disc3 f x y z = f (primal x) (primal y) (primal z)
+{-# INLINE disc3 #-}
+
+from :: Num a => RAD s a -> a -> RAD s a
+from (RAD (C a)) x = RAD (C x)
+from a x = RAD (U x 1 a)
+
+fromBy :: Num a => RAD s a -> RAD s a -> Int -> a -> RAD s a 
+fromBy (RAD (C a)) _ _ x = RAD (C x)
+fromBy a delta n x = RAD (B x 1 (fromIntegral n) a delta)
+
+instance (Num a, Enum a) => Enum (RAD s a) where
+    succ = unary_ succ 1
+    pred = unary_ pred 1
+    toEnum   = lift . toEnum
+    fromEnum = disc1 fromEnum
+    -- the enumerated results vary with the lower bound and so their derivatives reflect that
+    enumFrom a           = from a <$> disc1 enumFrom a
+    enumFromTo a b       = from a <$> disc2 enumFromTo a b
+    -- these results vary with respect to both the lower bound and the delta between that and the second argument
+    enumFromThen a b     = zipWith (fromBy a delta) [0..] $ disc2 enumFromThen a b where delta = b - a
+    enumFromThenTo a b c = zipWith (fromBy a delta) [0..] $ disc3 enumFromThenTo a b c where delta = b - a
+
+instance Num a => Num (RAD s a) where
+    fromInteger = lift . fromInteger
+    (+) = binary_ (+) 1 1 
+    (-) = binary_ (-) 1 (-1)
+    negate = unary_ negate (-1)
+    (*) = binary (*) id id
+    -- incorrect if the argument is complex
+    abs = unary abs signum
+    signum = lift . signum . primal
+
+-- notComplex :: Num a => a -> Bool
+-- notComplex x = s == 0 || s == 1 || s == -1
+--     where s = signum x 
+
+instance Real a => Real (RAD s a) where
+    toRational = disc1 toRational
+
+instance RealFloat a => RealFloat (RAD s a) where
+    floatRadix = disc1 floatRadix
+    floatDigits = disc1 floatDigits
+    floatRange = disc1 floatRange
+
+    decodeFloat = disc1 decodeFloat
+    encodeFloat m e = lift (encodeFloat m e)
+
+    scaleFloat n = unary_ (scaleFloat n) (scaleFloat n 1) 
+    isNaN = disc1 isNaN
+    isInfinite = disc1 isInfinite
+    isDenormalized = disc1 isDenormalized
+    isNegativeZero = disc1 isNegativeZero
+    isIEEE = disc1 isIEEE
+
+    exponent x
+        | m == 0 = 0 
+        | otherwise = n + floatDigits x
+        where (m,n) = decodeFloat x 
+
+    significand x =  unary_ significand (scaleFloat (- floatDigits x) 1) x
+
+    atan2 (RAD (C x)) (RAD (C y)) = RAD (C (atan2 x y))
+    atan2 x y = RAD (B (atan2 vx vy) (vy*r) (-vx*r) x y)
+        where vx = primal x
+              vy = primal y
+              r = recip (vx^2 + vy^2)
+
+instance RealFrac a => RealFrac (RAD s a) where
+    properFraction (RAD (C a)) = (w, RAD (C p))
+        where (w, p) = properFraction a
+    properFraction a = (w, RAD (U p 1 a))
+        where (w, p) = properFraction (primal a)
+    truncate = disc1 truncate
+    round = disc1 truncate
+    ceiling = disc1 truncate
+    floor = disc1 truncate
+
+instance Fractional a => Fractional (RAD s a) where
+    (/) = binary (/) recip id
+--   recip = unary recip  (const . negate . (^2))
+    fromRational r = lift $ fromRational r
+
+instance Floating a => Floating (RAD s a) where
+    pi      = lift pi
+    exp     = unary exp id
+    log     = unary log recip
+    sqrt    = unary sqrt (recip . (2*) . sqrt)
+    RAD (C x) ** RAD (C y) = lift (x ** y)
+    x ** y  = RAD (B (vx ** vy) (vy/vx) (log (vx)) x y)
+        where vx = primal x
+              vy = primal y
+    sin     = unary sin cos
+    cos     = unary cos (negate . sin)
+    asin    = unary asin (recip . sqrt . (1-) . (^2))
+    acos    = unary acos (negate . recip . sqrt . (1-) . (^2))
+    atan    = unary atan (recip . (1+) . (^2))
+    sinh    = unary sinh cosh
+    cosh    = unary cosh sinh
+    asinh   = unary asinh (recip . sqrt . (1+) . (^2))
+    acosh   = unary acosh (recip . sqrt . (-1+) . (^2))
+    atanh   = unary atanh (recip . (1-) . (^2))
+
+-- back propagate sensitivities along the tape.
+backprop :: (Ix t, Ord t, Num a) => (Vertex -> (Tape a t, t, [t])) -> STArray s t a -> Vertex -> ST s ()
+backprop vmap ss v = do
+        case node of
+            U _ g b -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + g*da)
+            B _ gb gc b c -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + gb*da)
+                dc <- readArray ss c
+                writeArray ss c (dc + gc*da)
+            _ -> return ()
+    where 
+        (node, i, _) = vmap v
+
+
+runTape :: Num a => (Int, Int) -> RAD s a -> Array Int a 
+runTape vbounds tape = accumArray (+) 0 vbounds [ (id, sensitivities ! ix) | (ix, V _ id) <- xs ]
+    where
+        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape
+        (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)
+        sensitivities = runSTArray $ do
+            ss <- newArray (sbounds xs) 0
+            writeArray ss start 1
+            forM_ (topSort g) $ 
+                backprop vmap ss
+            return ss
+        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> (min lo b, max hi b)) (a,a) as
+        edgeSet (i, t) = (t, i, successors t)
+        nonConst (_, C{}) = False
+        nonConst _ = True
+        successors (U _ _ b)   = [b]
+        successors (B _ _ _ b c) = [b,c]
+        successors _ = []
+
+d :: Num a => RAD s a -> a
+d r = runTape (0,0) r ! 0 
+
+d2 :: Num a => RAD s a -> (a,a)
+d2 r = (primal r, d r)
+
+
+-- | The 'diffUU' function calculates the first derivative of a
+-- scalar-to-scalar function.
+diffUU :: Num a => (forall s. RAD s a -> RAD s a) -> a -> a
+diffUU f a = d $ f (var a 0)
+
+
+-- | The 'diffUF' function calculates the first derivative of
+-- scalar-to-nonscalar function.
+diffUF :: (Functor f, Num a) => (forall s. RAD s a -> f (RAD s a)) -> a -> f a
+diffUF f a = d <$> f (var a 0)
+
+-- diffMU :: Num a => (forall s. [RAD s a] -> RAD s a) -> [a] -> [a] -> a
+-- TODO: finish up diffMU and their ilk
+
+-- avoid dependency on MTL
+newtype S a = S { runS :: Int -> (a,Int) } 
+
+instance Monad S where
+    return a = S (\s -> (a,s))
+    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
+    
+bind :: Traversable f => f a -> (f (RAD s a), (Int,Int))
+bind xs = (r,(0,s)) 
+    where 
+        (r,s) = runS (mapM freshVar xs) 0
+        freshVar a = S (\s -> let s' = s + 1 in s' `seq` (RAD (V a s), s'))
+
+unbind :: Functor f => f (RAD s b) -> Array Int a -> f a 
+unbind xs ys = fmap (\(RAD (V _ i)) -> ys ! i) xs
+
+-- | The 'diff2UU' function calculates the value and derivative, as a
+-- pair, of a scalar-to-scalar function.
+diff2UU :: Num a => (forall s. RAD s a -> RAD s a) -> a -> (a, a)
+diff2UU f a = d2 $ f (var a 0)
+
+-- | Note that the signature differs from that used in Numeric.FAD, because while you can always
+-- 'unzip' an arbitrary functor, not all functors can be zipped.
+diff2UF :: (Functor f, Num a) => (forall s. RAD s a -> f (RAD s a)) -> a -> f (a, a)
+diff2UF f a = d2 <$> f (var a 0)
+
+-- | The 'diff' function is a synonym for 'diffUU'.
+diff :: Num a => (forall s. RAD s a -> RAD s a) -> a -> a
+diff = diffUU 
+
+-- | The 'diff2' function is a synonym for 'diff2UU'.
+diff2 :: Num a => (forall s. RAD s a -> RAD s a) -> a -> (a, a)
+diff2 = diff2UU
+
+-- requires the input list to be finite in length
+grad :: (Traversable f, Num a) => (forall s. f (RAD s a) -> RAD s a) -> f a -> f a
+grad f as = unbind s (runTape bounds $ f s)
+    where (s,bounds) = bind as
+
+-- compute the primal and gradient
+grad2 :: (Traversable f, Num a) => (forall s. f (RAD s a) -> RAD s a) -> f a -> (a, f a)
+grad2 f as = (primal r, unbind s (runTape bounds r))
+    where (s,bounds) = bind as
+          r = f s
+
+-- | The 'jacobian' function calcualtes the Jacobian of a
+-- nonscalar-to-nonscalar function, using m invocations of reverse AD,
+-- where m is the output dimensionality. When the output dimensionality is
+-- significantly greater than the input dimensionality you should use 'Numeric.FAD.jacobian' instead.
+jacobian :: (Traversable f, Functor g, Num a) => (forall s. f (RAD s a) -> g (RAD s a)) -> f a -> g (f a)
+jacobian f as = unbind s . runTape bounds <$> f s
+    where (s,bounds) = bind as
+
+-- | The 'jacobian2' function calcualtes both the result and the Jacobian of a
+-- nonscalar-to-nonscalar function, using m invocations of reverse AD,
+-- where m is the output dimensionality. 
+-- 'fmap snd' on the result will recover the result of 'jacobian'
+jacobian2 :: (Traversable f, Functor g, Num a) => (forall s. f (RAD s a) -> g (RAD s a)) -> f a -> g (a, f a)
+jacobian2 f as = row <$> f s
+    where (s,bounds) = bind as
+          row a = (primal a, unbind s (runTape bounds a))
+
+-- | The 'zeroNewton' function finds a zero of a scalar function using
+-- Newton's method; its output is a stream of increasingly accurate
+-- results.  (Modulo the usual caveats.)
+--
+-- TEST CASE:
+--  @take 10 $ zeroNewton (\\x->x^2-4) 1  -- converge to 2.0@
+--
+-- TEST CASE
+--  :module Data.Complex Numeric.RAD
+--  @take 10 $ zeroNewton ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@
+--
+zeroNewton :: Fractional a => (forall s. RAD s a -> RAD s a) -> a -> [a]
+zeroNewton f x0 = iterate (\x -> let (y,y') = diff2UU f x in x - y/y') x0
+
+-- | The 'inverseNewton' function inverts a scalar function using
+-- Newton's method; its output is a stream of increasingly accurate
+-- results.  (Modulo the usual caveats.)
+--
+-- TEST CASE:
+--   @take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converge to 10@
+--
+inverseNewton :: Fractional a => (forall s. RAD s a -> RAD s a) -> a -> a -> [a]
+inverseNewton f x0 y = zeroNewton (\x -> f x - lift y) x0
+
+-- | The 'fixedPointNewton' function find a fixedpoint of a scalar
+-- function using Newton's method; its output is a stream of
+-- increasingly accurate results.  (Modulo the usual caveats.)
+fixedPointNewton :: Fractional a => (forall s. RAD s a -> RAD s a) -> a -> [a]
+fixedPointNewton f = zeroNewton (\x -> f x - x)
+
+-- | The 'extremumNewton' function finds an extremum of a scalar
+-- function using Newton's method; produces a stream of increasingly
+-- accurate results.  (Modulo the usual caveats.)
+extremumNewton :: Fractional a => (forall s t. RAD t (RAD s a) -> RAD t (RAD s a)) -> a -> [a]
+extremumNewton f x0 = zeroNewton (diffUU f) x0
+
+-- | The 'argminNaiveGradient' function performs a multivariate
+-- optimization, based on the naive-gradient-descent in the file
+-- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the
+-- VLAD compiler Stalingrad sources.  Its output is a stream of
+-- increasingly accurate results.  (Modulo the usual caveats.)  
+-- This is /O(n)/ faster than 'Numeric.FAD.argminNaiveGradient'
+argminNaiveGradient :: (Fractional a, Ord a) => (forall s. [RAD s a] -> RAD s a) -> [a] -> [[a]]
+argminNaiveGradient f x0 =
+    let
+        gf = grad f
+        loop x fx gx eta i =
+            -- should check gx = 0 here
+            let
+                x1 = zipWith (+) x (map ((-eta)*) gx)
+                fx1 = lowerFU f x1
+                gx1 = gf x1
+            in
+              if eta == 0 then []
+              else if (fx1 > fx) then loop x fx gx (eta/2) 0
+                   else if all (==0) gx then []
+                        -- else if fx1 == fx then loop x1 fx1 gx1 eta (i+1)
+                        else x1:(if (i==10)
+                                 then loop x1 fx1 gx1 (eta*2) 0
+                                 else loop x1 fx1 gx1 eta (i+1))
+    in
+      loop x0 (lowerFU f x0) (gf x0) 0.1 0
+
+{-
+lowerUU :: (forall s. RAD s a -> RAD s b) -> a -> b
+lowerUU f = primal . f . lift
+
+lowerUF :: Functor f => (forall s. RAD s a -> f (RAD s b)) -> a -> f b
+lowerUF f = fmap primal . f . lift
+
+lowerFF :: (Functor f, Functor g) => (forall s. f (RAD s a) -> g (RAD s b)) -> f a -> g b
+lowerFF f = fmap primal . f . fmap lift
+-}
+
+lowerFU :: Functor f => (forall s. f (RAD s a) -> RAD s b) -> f a -> b
+lowerFU f = primal . f . fmap lift
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,15 @@
+* Implement the diffMF, etc. functionality from Numeric.FAD
+
+* Allow the type to vary within our RAD data type container, in the same fashion as Numeric.FAD.
+
+    Although, while Pearlmutter and Siskind provided the functionality to permit it in derivative combinator, they provided
+    no combinators to convert, for instance, a Dual tag Float to a Dual tag Double, so that extra functionality cannot
+    currently be leveraged.
+
+* Implement a reverse-mode cotangent tower.
+
+    Can we just play back the tape in a forward-mode tower?
+
+* Provide the ability to use reverse mode locally on FAD inputs, i.e.
+
+    reverseCheckpoint :: (forall s. RAD s a -> RAD s a) -> FAD t a -> FAD t a 
diff --git a/rad.cabal b/rad.cabal
new file mode 100644
--- /dev/null
+++ b/rad.cabal
@@ -0,0 +1,21 @@
+Name:                rad
+Version:             0.1
+License:             BSD3
+License-File:        LICENSE
+Copyright:           Edward Kmett 2010
+Author:              Edward Kmett 2010
+Maintainer:          ekmett@gmail.com
+Stability:           Experimental
+Homepage:            http://comonad.com/reader/
+Synopsis:            Reverse Automatic Differentiation.
+Description:
+    Reverse-Mode Automatic Differentiation via overloading.
+    Existential type \"branding\" is used to prevent sensitivity confusion.
+Category:            Math
+Build-Type:          Simple
+Build-Depends:       base >= 4 && < 6,
+                     data-reify >= 0.5 && < 0.6, 
+                     containers >= 0.2 && < 0.3,
+                     array >= 0.2 && < 0.3
+Exposed-Modules:     Numeric.RAD
+Extra-Source-Files:  TODO