diff --git a/Changelog.md b/Changelog.md
new file mode 100644
--- /dev/null
+++ b/Changelog.md
@@ -0,0 +1,3 @@
+0.1
+
+First public version
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright Author name here (c) 2021
+
+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/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,8 @@
+# ad-delcont
+
+Reverse-mode Automatic Differentiation using delimited continuations (`shift` and `reset`) , as shown in [1]
+
+
+## References
+
+1) Wang et al, Demystifying Differentiable Programming : Shift/Reset the Penultimate Backpropagator,  https://www.cs.purdue.edu/homes/rompf/papers/wang-icfp19.pdf 
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/ad-delcont.cabal b/ad-delcont.cabal
new file mode 100644
--- /dev/null
+++ b/ad-delcont.cabal
@@ -0,0 +1,50 @@
+name:                ad-delcont
+version:             0.1.0.0
+synopsis:            Reverse-mode automatic differentiation with delimited continuations
+description:         Reverse-mode automatic differentiation using delimited continuations (@shift@/@reset@), inspired by the papers
+                     .
+                     * F. Wang et al, Backpropagation with Continuation Callbacks : Foundations for Efficient and Expressive Differentiable Programming, NeurIPS 2018 - https://papers.nips.cc/paper/2018/file/34e157766f31db3d2099831d348a7933-Paper.pdf
+                     .
+                     * F. Wang et al, Demystifying Differentiable Programming : Shift\/Reset the Penultimate Backpropagator, ICFP 2019 - https://www.cs.purdue.edu/homes/rompf/papers/wang-icfp19.pdf
+                     .
+                     The package exposes a small and easily extensible user interface to automatic differentiation combinators. It also manages to be very lightweight, as it only depends on @base@ and @transformers@.
+                     .
+                     To use the library, import "Numeric.AD.DelCont", which also contains all documentation.
+homepage:            https://github.com/ocramz/ad-delcont
+license:             BSD3
+license-file:        LICENSE
+author:              Marco Zocca
+maintainer:          ocramz
+copyright:           (c) 2021 Marco Zocca
+category:            Math, Numeric, Machine Learning, Optimization, Optimisation
+build-type:          Simple
+extra-source-files:  README.md
+cabal-version:       >=1.10
+tested-with:         GHC == 8.10.4
+extra-source-files:  README.md
+                     Changelog.md
+
+library
+  default-language:    Haskell2010
+  ghc-options:         -Wall
+  hs-source-dirs:      src
+  exposed-modules:     Numeric.AD.DelCont
+  other-modules:       Numeric.AD.DelCont.Internal
+  build-depends:       base >= 4.7 && < 5
+                     , transformers >= 0.5
+
+-- -- TODO
+-- test-suite spec
+--   default-language:    Haskell2010
+--   ghc-options:         -Wall
+--   type:                exitcode-stdio-1.0
+--   hs-source-dirs:      test
+--   main-is:             Spec.hs
+--   build-depends:       base
+--                      , ad-delcont
+--                      , hspec
+--                      , QuickCheck
+
+source-repository head
+  type:     git
+  location: https://github.com/ocramz/ad-delcont
diff --git a/src/Numeric/AD/DelCont.hs b/src/Numeric/AD/DelCont.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/DelCont.hs
@@ -0,0 +1,68 @@
+{-| Reverse-mode automatic differentiation using delimited continuations.
+
+== Quickstart
+
+Most users will only need to import 'rad1', 'rad2' and leverage the 'Num', 'Fractional', 'Floating' instances of the 'AD' type.
+
+Similarly to @ad@, a user supplies a /polymorphic/ function to be differentiated, e.g.
+
+@
+f :: Num a => a -> a
+f x = x + (x * x)
+@
+
+and the library takes care of the rest :
+
+@
+>>> 'rad1' f 1.2
+(2.6399999999999997,3.4000000000000004)
+@
+
+It's important to emphasize that the library cannot differentiate functions of concrete types, e.g. @Double -> Double@. On the other hand, it's easy to experiment with other numerical interfaces that support one, zero and plus.
+
+== Advanced usage
+
+The library is small and easily extensible.
+
+For example, a user might want to supply their own numerical typeclass other than 'Num', and build up a library of 'AD' combinators based on that, specializing 'op1' and 'op2' with custom implementations of @zero@, @one@ and @plus@. This insight first appeared in the user interface of @backprop@, as the Backprop typeclass.
+
+Exposing unconstrained AD combinators lets users specialize this library to e.g. exotic number-like types or discrete data structures such as dictionaries, automata etc.
+
+== Implementation details and design choices
+
+This is the first (known) Haskell implementation of the ideas presented in Wang et al. Here the role of variable mutation and delimited continuations is made explicit by the use of 'ST' and 'ContT', as compared to the reference Scala implementation.
+
+@ad-delcont@ relies on non-standard interpretation of the user-provided function; in order to compute the adjoint values (the /sensitivities/) of the function parameters, the function is first evaluated ("forwards"), while keeping track of continuation points, and all the intermediate adjoints are accumulated upon returning from the respective continuations ("backwards") via safe mutation in the ST monad.
+
+As a result of this design, the main 'AD' type cannot be given 'Eq' and 'Ord' instances (since it's unclear how equality and ordering predicates would apply to continuations and state threads).
+
+The user interface is inspired by that of @ad@ and @backprop@, however the internals are completely different in that this library doesn't reify the function to be differentiated into a "tape" data structure.
+
+Another point in common with @backprop@ is that users can differentiate heterogeneous functions: the input and output types can be different. This makes it possible to differentiate functions of statically-typed vectors and matrices.
+
+
+== References
+
+* @backprop@ - https://hackage.haskell.org/package/backprop
+
+* @ad@ - https://hackage.haskell.org/package/ad
+
+* F. Wang et al, Backpropagation with Continuation Callbacks : Foundations for Efficient and Expressive Differentiable Programming, NeurIPS 2018 - https://papers.nips.cc/paper/2018/file/34e157766f31db3d2099831d348a7933-Paper.pdf
+
+* F. Wang et al, Demystifying Differentiable Programming : Shift\/Reset the Penultimate Backpropagator, ICFP 2019 - https://doi.org/10.1145/3341700 - https://www.cs.purdue.edu/homes/rompf/papers/wang-icfp19.pdf
+
+* M. Innes, Don't unroll adjoint: Differentiating SSA-Form Programs  https://arxiv.org/abs/1810.07951
+-}
+module Numeric.AD.DelCont (-- * Quickstart
+                            rad1, rad2
+                          , auto
+                          -- * Advanced usage
+                          , rad1g, rad2g,
+                            -- ** Lift functions into AD
+                            op1, op2
+                            -- *** Num instances
+                          -- * Types
+                          , AD, AD') where
+
+import Numeric.AD.DelCont.Internal (rad1, rad2, auto, rad1g, rad2g, op1, op2, AD, AD')
+
diff --git a/src/Numeric/AD/DelCont/Internal.hs b/src/Numeric/AD/DelCont/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/AD/DelCont/Internal.hs
@@ -0,0 +1,282 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE DeriveFunctor #-}
+{-# LANGUAGE LambdaCase #-}
+module Numeric.AD.DelCont.Internal
+  (rad1, rad2,
+   auto,
+   rad1g, rad2g,
+   op1Num, op2Num,
+   op1, op2,
+   AD, AD')
+  where
+
+import Control.Monad.ST (ST, runST)
+import Data.Bifunctor (Bifunctor(..))
+import Data.STRef (STRef, newSTRef, readSTRef, modifySTRef')
+
+-- transformers
+import Control.Monad.Trans.Class (MonadTrans(..))
+import Control.Monad.Trans.Cont (ContT, shiftT, resetT, evalContT)
+
+import Prelude hiding (read)
+
+-- | Dual numbers
+data D a da = D a da deriving (Show, Functor)
+instance Eq a => Eq (D a da) where
+  D x _ == D y _ = x == y
+instance Ord a => Ord (D a db) where
+  compare (D x _) (D y _) = compare x y
+instance Bifunctor D where
+  bimap f g (D a b) = D (f a) (g b)
+
+-- | Modify the adjoint part of a 'D'
+withD :: (da -> db) -> D a da -> D a db
+withD = second
+
+-- | Differentiable variable
+--
+-- A (safely) mutable reference to a dual number
+type DVar s a da = STRef s (D a da)
+-- | Introduce a fresh DVar
+var :: a -> da -> ST s (DVar s a da)
+var x dx = newSTRef (D x dx)
+
+-- | Lift a constant into 'AD'
+auto :: a -- ^ primal
+     -> da -- ^ adjoint (in most cases this can be set to (@0 :: a@))
+     -> AD s a da
+auto x dx = AD $ lift $ var x dx
+
+-- | Mutable references to dual numbers in the continuation monad
+--
+-- Here the @a@ and @da@ type parameters are respectively the /primal/ and /dual/ quantities tracked by the AD computation.
+newtype AD s a da = AD { unAD :: forall x dx . ContT (DVar s x dx) (ST s) (DVar s a da) }
+-- | Like 'AD' but the types of primal and dual coincide
+type AD' s a = AD s a a
+
+-- runAD :: (forall s . AD s a da) -> D a da
+-- runAD go = runST (evalContT  (unAD go) >>= readSTRef)
+
+-- | Lift a unary function
+--
+-- This is a polymorphic combinator for tracking how primal and adjoint values are transformed by a function.
+--
+-- How does this work :
+--
+-- 1) Compute the function result and bind the function inputs to the adjoint updating function (the "pullback")
+--
+-- 2) Allocate a fresh STRef @rb@ with the function result and @zero@ adjoint part
+--
+-- 3) @rb@ is passed downstream as an argument to the continuation @k@, with the expectation that the STRef will be mutated
+--
+-- 4) Upon returning from the @k@ (bouncing from the boundary of @resetT@), the mutated STRef is read back in
+--
+-- 5) The adjoint part of the input variable is updated using @rb@ and the result of the continuation is returned.
+op1_ :: db -- ^ zero
+     -> (da -> da -> da) -- ^ plus
+     -> (a -> (b, db -> da)) -- ^ returns : (function result, pullback)
+     -> ContT x (ST s) (DVar s a da)
+     -> ContT x (ST s) (DVar s b db)
+op1_ zero plusa f ioa = do
+  ra <- ioa
+  (D xa _) <- lift $ readSTRef ra
+  let (xb, g) = f xa -- 1)
+  shiftT $ \ k -> lift $ do
+    rb <- var xb zero -- 2)
+    ry <- k rb -- 3)
+    (D _ yd) <- readSTRef rb -- 4)
+    modifySTRef' ra (withD (\rda0 -> rda0 `plusa` g yd)) -- 5)
+    pure ry
+
+-- | Lift a unary function
+--
+-- The first two arguments constrain the types of the adjoint values of the output and input variable respectively, see 'op1Num' for an example.
+--
+-- The third argument is the most interesting: it specifies at once how to compute the function value and how to compute the sensitivity with respect to the function parameter.
+--
+-- Note : the type parameters are completely unconstrained.
+op1 :: db -- ^ zero
+    -> (da -> da -> da) -- ^ plus
+    -> (a -> (b, db -> da)) -- ^ returns : (function result, pullback)
+    -> AD s a da
+    -> AD s b db
+op1 z plusa f (AD ioa) = AD $ op1_ z plusa f ioa
+
+-- | Helper for constructing Num instances (= 'op1' specialized to Num)
+op1Num :: (Num da, Num db) =>
+          (a -> (b, db -> da))
+       -> AD s a da
+       -> AD s b db
+op1Num = op1 0 (+)
+
+-- | Lift a binary function
+op2_ :: dc -- ^ zero
+     -> (da -> da -> da) -- ^ plus
+     -> (db -> db -> db) -- ^ plus
+     -> (a -> b -> (c, dc -> da, dc -> db)) -- ^ returns : (function result, pullbacks)
+     -> ContT x (ST s) (DVar s a da)
+     -> ContT x (ST s) (DVar s b db)
+     -> ContT x (ST s) (DVar s c dc)
+op2_ zero plusa plusb f ioa iob = do
+  ra <- ioa
+  rb <- iob
+  (D xa _) <- lift $ readSTRef ra
+  (D xb _) <- lift $ readSTRef rb
+  let (xc, g, h) = f xa xb
+  shiftT $ \ k -> lift $ do
+    rc <- var xc zero
+    ry <- k rc
+    (D _ yd) <- readSTRef rc
+    modifySTRef' ra (withD (\rda0 -> rda0 `plusa` g yd))
+    modifySTRef' rb (withD (\rdb0 -> rdb0 `plusb` h yd))
+    pure ry
+
+-- | Lift a binary function
+--
+-- See 'op1' for more information.
+op2 :: dc -- ^ zero
+    -> (da -> da -> da) -- ^ plus
+    -> (db -> db -> db) -- ^ plus
+    -> (a -> b -> (c, dc -> da, dc -> db)) -- ^ returns : (function result, pullbacks)
+    -> (AD s a da -> AD s b db -> AD s c dc)
+op2 z plusa plusb f (AD ioa) (AD iob) = AD $ op2_ z plusa plusb f ioa iob
+
+-- | Helper for constructing Num instances (= 'op2' specialized to Num)
+op2Num :: (Num da, Num db, Num dc) =>
+          (a -> b -> (c, dc -> da, dc -> db))
+       -> AD s a da
+       -> AD s b db
+       -> AD s c dc
+op2Num = op2 0 (+) (+)
+
+-- | The numerical methods of (Num, Fractional, Floating etc.) can be read off their @backprop@ counterparts : https://hackage.haskell.org/package/backprop-0.2.6.4/docs/src/Numeric.Backprop.Op.html#%2A.
+instance (Num a) => Num (AD s a a) where
+  (+) = op2Num $ \x y -> (x + y, id, id)
+  (-) = op2Num $ \x y -> (x - y, id, negate)
+  (*) = op2Num $ \x y -> (x*y, (y *), (x *))
+  fromInteger x = auto (fromInteger x) 0
+  abs = op1Num $ \x -> (abs x, (* signum x))
+  signum = op1Num $ \x -> (signum x, const 0)
+
+instance (Fractional a) => Fractional (AD s a a) where
+  (/) = op2Num $ \x y -> (x / y, (/ y), (\g -> -g*x/(y*y) ))
+  fromRational x = auto (fromRational x) 0
+  recip = op1Num $ \x -> (recip x, (/(x*x)) . negate)
+
+instance Floating a => Floating (AD s a a) where
+  pi = auto pi 0
+  exp = op1Num $ \x -> (exp x, (exp x *))
+  log = op1Num $ \x -> (log x, (/x))
+  sqrt = op1Num $ \x -> (sqrt x, (/ (2 * sqrt x)))
+  logBase = op2Num $ \x y ->
+                       let
+                         dx = - logBase x y / (log x * x)
+                       in ( logBase x y
+                          , ( * dx)
+                          , (/(y * log x))
+                          )
+  sin = op1Num $ \x -> (sin x, (* cos x))
+  cos = op1Num $ \x -> (cos x, (* (-sin x)))
+  tan = op1Num $ \x -> (tan x, (/ cos x^(2::Int)))
+  asin = op1Num $ \x -> (asin x, (/ sqrt(1 - x*x)))
+  acos = op1Num $ \x -> (acos x, (/ sqrt (1 - x*x)) . negate)
+  atan = op1Num $ \x -> (atan x, (/ (x*x + 1)))
+  sinh = op1Num $ \x -> (sinh x, (* cosh x))
+  cosh = op1Num $ \x -> (cosh x, (* sinh x))
+  tanh = op1Num $ \x -> (tanh x, (/ cosh x^(2::Int)))
+  asinh = op1Num $ \x -> (asinh x, (/ sqrt (x*x + 1)))
+  acosh = op1Num $ \x -> (acosh x, (/ sqrt (x*x - 1)))
+  atanh = op1Num $ \x -> (atanh x, (/ (1 - x*x)))
+
+-- instance Eq a => Eq (AD s a da) where -- ??? likely impossible
+-- instance Ord (AD s a da) where -- ??? see above
+
+-- | Evaluate (forward mode) and differentiate (reverse mode) a unary function, without committing to a specific numeric typeclass
+rad1g :: da -- ^ zero
+      -> db -- ^ one
+      -> (forall s . AD s a da -> AD s b db)
+      -> a -- ^ function argument
+      -> (b, da) -- ^ (result, adjoint)
+rad1g zero one f x = runST $ do
+  xr <- var x zero
+  zr' <- evalContT $
+    resetT $ do
+      let
+        z = f (AD (pure xr))
+      zr <- unAD z
+      lift $ modifySTRef' zr (withD (const one))
+      pure zr
+  (D z _) <- readSTRef zr'
+  (D _ x_bar) <- readSTRef xr
+  pure (z, x_bar)
+
+
+
+-- | Evaluate (forward mode) and differentiate (reverse mode) a binary function, without committing to a specific numeric typeclass
+rad2g :: da -- ^ zero
+      -> db -- ^ zero
+      -> dc -- ^ one
+      -> (forall s . AD s a da -> AD s b db -> AD s c dc)
+      -> a -> b
+      -> (c, (da, db)) -- ^ (result, adjoints)
+rad2g zeroa zerob one f x y = runST $ do
+  xr <- var x zeroa
+  yr <- var y zerob
+  zr' <- evalContT $
+    resetT $ do
+      let
+        z = f (AD (pure xr)) (AD (pure yr))
+      zr <- unAD z
+      lift $ modifySTRef' zr (withD (const one))
+      pure zr
+  (D z _) <- readSTRef zr'
+  (D _ x_bar) <- readSTRef xr
+  (D _ y_bar) <- readSTRef yr
+  pure (z, (x_bar, y_bar))
+
+
+-- | Evaluate (forward mode) and differentiate (reverse mode) a unary function
+--
+-- >>> rad1 (\x -> x * x) 1
+-- (1, 2)
+rad1 :: (Num a, Num b) =>
+        (forall s . AD' s a -> AD' s b) -- ^ function to be differentiated
+     -> a -- ^ function argument
+     -> (b, a) -- ^ (result, adjoint)
+rad1 = rad1g 0 1
+
+-- | Evaluate (forward mode) and differentiate (reverse mode) a binary function
+--
+-- >>> rad2 (\x y -> x + y + y) 1 1
+-- (1,2)
+--
+-- >>> rad2 (\x y -> (x + y) * x) 3 2
+-- (15,(8,3))
+rad2 :: (Num a, Num b, Num c) =>
+        (forall s . AD' s a -> AD' s b -> AD' s c) -- ^ function to be differentiated
+     -> a
+     -> b
+     -> (c, (a, b)) -- ^ (result, adjoints)
+rad2 = rad2g 0 0 1
+
+
+
+
+-- -- playground
+
+-- -- | Dual numbers DD (alternative take, using a type family for the first variation)
+-- data DD a = Dd a (Adj a)
+-- class Diff a where type Adj a :: *
+-- instance Diff Double where type Adj Double = Double
+
+
+-- -- product type (simplified version of vinyl's Rec)
+-- data Rec :: [*] -> * where
+--   RNil :: Rec '[]
+--   (:*) :: !a -> !(Rec as) -> Rec (a ': as)
+
+-- data SDRec s as where
+--   SDNil :: SDRec s '[]
+--   (:&) :: DVar s a a -> !(SDRec s as) -> SDRec s (a ': as)