diff --git a/.travis.yml b/.travis.yml
index bf75be4..ab41d95 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -46,31 +46,15 @@ matrix:
   #- env: BUILD=cabal GHCVER=7.2.2 CABALVER=1.16 HAPPYVER=1.19.5 ALEXVER=3.1.7
   #  compiler: ": #GHC 7.2.2"
   #  addons: {apt: {packages: [cabal-install-1.16,ghc-7.2.2,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
-  - env: BUILD=cabal GHCVER=7.4.2 CABALVER=1.16 HAPPYVER=1.19.5 ALEXVER=3.1.7
-    compiler: ": #GHC 7.4.2"
-    addons: {apt: {packages: [libgmp-dev,
-                              libmpfr-dev,
-                              libmpc-dev,
-                              binutils-dev,
-                              g++-4.7,
-                              gcc, cabal-install-1.16,ghc-7.4.2,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
-  - env: BUILD=cabal GHCVER=7.6.3 CABALVER=1.16 HAPPYVER=1.19.5 ALEXVER=3.1.7
-    compiler: ": #GHC 7.6.3"
-    addons: {apt: {packages: [libgmp-dev,
-                              libmpfr-dev,
-                              libmpc-dev,
-                              binutils-dev,
-                              g++-4.7,
-                              gcc, cabal-install-1.16,ghc-7.6.3,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
-  - env: BUILD=cabal GHCVER=7.8.4 CABALVER=1.18 HAPPYVER=1.19.5 ALEXVER=3.1.7
+  - env: BUILD=cabal CABALVER=1.24 GHCVER=7.8.4 STACK_YAML=stack-7.8.yaml
     compiler: ": #GHC 7.8.4"
     addons: {apt: {packages: [libgmp-dev,
                               libmpfr-dev,
                               libmpc-dev,
                               binutils-dev,
                               g++-4.7,
-                              gcc, cabal-install-1.18,ghc-7.8.4,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
-  - env: BUILD=cabal GHCVER=7.10.3 CABALVER=1.22 HAPPYVER=1.19.5 ALEXVER=3.1.7
+                              gcc, cabal-install-1.24,ghc-7.8.4,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
+  - env: BUILD=cabal CABALVER=1.22 GHCVER=7.10.3
     compiler: ": #GHC 7.10.3"
     addons: {apt: {packages: [libgmp-dev,
                               libmpfr-dev,
@@ -78,17 +62,14 @@ matrix:
                               binutils-dev,
                               g++-4.7,
                               gcc, cabal-install-1.22,ghc-7.10.3,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
-
-  # Build with the newest GHC and cabal-install. This is an accepted failure,
-  # see below.
-  - env: BUILD=cabal GHCVER=head  CABALVER=head HAPPYVER=1.19.5 ALEXVER=3.1.7
-    compiler: ": #GHC HEAD"
+  - env: BUILD=cabal CABALVER=1.24 GHCVER=8.0.1
+    compiler: ": #GHC 8.0.1"
     addons: {apt: {packages: [libgmp-dev,
                               libmpfr-dev,
                               libmpc-dev,
                               binutils-dev,
                               g++-4.7,
-                              gcc, cabal-install-head,ghc-head,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
+                              gcc, cabal-install-1.24,ghc-8.0.1,happy-1.19.5,alex-3.1.7], sources: [hvr-ghc, ubuntu-toolchain-r-test]}}
 
   # The Stack builds. We can pass in arbitrary Stack arguments via the ARGS
   # variable, such as using --stack-yaml to point to a different file.
@@ -210,6 +191,8 @@ install:
 
 # Download and unpack the stack executable
 - export PATH=/opt/ghc/$GHCVER/bin:/opt/cabal/$CABALVER/bin:$HOME/.local/bin:/opt/alex/$ALEXVER/bin:/opt/happy/$HAPPYVER/bin:$HOME/.cabal/bin:$PATH
+# add cabal install path
+- export PATH=PATH="$HOME/.cabal/bin:$PATH"
 - mkdir -p ~/.local/bin
 - |
   if [ `uname` = "Darwin" ]
diff --git a/README.md b/README.md
index f4b3ead..aed83ea 100644
--- a/README.md
+++ b/README.md
@@ -26,12 +26,18 @@ as well as the libraries
 
 Since these are *hard* dependencies for SymEngine-hs to build.
 
+Compile `SymEngine` with the `CMake` flags
+
+```bash
+cmake -DWITH_SYMENGINE_THREAD_SAFE=yes -DBUILD_SHARED_LIBS:BOOL=ON
+```
+
 # Getting started
 
 To quickly build and check everything is working, run
 
 ```
-stack build && stack test
+stack build && stack test --test-arguments "--quickcheck-tests 2000" --verbose
 ```
 
 All of the test cases should pass with SymEngine
@@ -68,14 +74,58 @@ GHCi session with Symengine loaded
 -1
 ```
 
-# Things to Do
+# Development
 
-`[TODO: fill this up]`
+clone `Symengine`, build it with the setting
 
-# Contributing
+```
+cmake -DWITH_SYMENGINE_THREAD_SAFE=yes -DBUILD_SHARED_LIBS:BOOL=ON
+```
+
+this makes sure that dynamically linked libraries are being built, so we can
+link to them.
+
+
+to test changes, use
+```
+stack test --force-dirty  --test-arguments "--quickcheck-tests 2000" --verbose
+```
+
+* change `--quickcheck-tests" to some number (preferably > 100), since it generates those many instances to
+test on
+
+* the `--force-dirty` ensures that the library and the test builds are both
+rebuilt.
 
-`[TODO: fill this up]`
 
 # License
 
 All code is released under the [MIT License](https://github.com/symengine/symengine.hs/blob/master/LICENSE).
+
+
+# Things Learnt 
+
+* you can use `toEnum` to convert from `Int` to the `C<Int | Long | ..>` variants
+of C types
+
+* API design - how to best handle exceptions?
+
+# Bugs
+
+* if I create a lazy list of BasicSym, then what happens? it gets forced to evaluate
+when I pass it through something like `densematrix_diag`
+
+
+* `densematrix_new_vec  2 3 []` crashes. We need to check for this in our code
+
+
+* What exactly does 'unsafePerformIO' do? why does `unsafePerformIO` on `basicsym_new`
+yield weird as hell errors?
+
+* take proper care of ref. transparency. eg: `densematrix_set`
+
+* Maybe allow GHC to tell about "typo errors" when looking for modules
+
+* `merijn	You'll want newPinnedByteArray# :: Int# -> State# s -> (#State# s, MutableByteArray# s#)`
+
+* is the API Thread-safe?
diff --git a/src/Symengine.hs b/src/Symengine.hs
index d0671ff..e155165 100644
--- a/src/Symengine.hs
+++ b/src/Symengine.hs
@@ -1,22 +1,13 @@
 {-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
 
 {-|
 Module      : Symengine
 Description : Symengine bindings to Haskell
 -}
 module Symengine
-    (
-     ascii_art_str,
-     zero,
-     one,
-     im,
-     Symengine.pi,
-     e,
-     minus_one,
-     rational,
-     complex,
-     symbol,
-     BasicSym,
+  (module Symengine.Internal
     ) where
 
 import Foreign.C.Types
@@ -27,255 +18,9 @@ import Foreign.Marshal.Array
 import Foreign.Marshal.Alloc
 import Foreign.ForeignPtr
 import Control.Applicative
+import Control.Monad -- for foldM
 import System.IO.Unsafe
 import Control.Monad
 import GHC.Real
 
-data BasicStruct = BasicStruct {
-    data_ptr :: Ptr ()
-}
-
-instance Storable BasicStruct where
-    alignment _ = 8
-    sizeOf _ = sizeOf nullPtr
-    peek basic_ptr = BasicStruct <$> peekByteOff basic_ptr 0
-    poke basic_ptr BasicStruct{..} = pokeByteOff basic_ptr 0 data_ptr
-
-
--- |represents a symbol exported by SymEngine. create this using the functions
--- 'zero', 'one', 'minus_one', 'e', 'im', 'rational', 'complex', and also by
--- constructing a number and converting it to a Symbol
--- 
--- >>> 3.5 :: BasicSym
--- 7/2
---
--- >>> rational 2 10
--- 1 /5
---
--- >>> complex 1 2
--- 1 + 2*I
-data BasicSym = BasicSym { fptr :: ForeignPtr BasicStruct }
-
-withBasicSym :: BasicSym -> (Ptr BasicStruct -> IO a) -> IO a
-withBasicSym p f = withForeignPtr (fptr p ) f
-
-withBasicSym2 :: BasicSym -> BasicSym -> (Ptr BasicStruct -> Ptr BasicStruct -> IO a) -> IO a
-withBasicSym2 p1 p2 f = withBasicSym p1 (\p1 -> withBasicSym p2 (\p2 -> f p1 p2))
-
-withBasicSym3 :: BasicSym -> BasicSym -> BasicSym -> (Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO a) -> IO a
-withBasicSym3 p1 p2 p3 f = withBasicSym p1 (\p1 -> withBasicSym p2 (\p2 -> withBasicSym p3 (\p3 -> f p1 p2 p3)))
-
-
--- | constructor for 0
-zero :: BasicSym
-zero = basic_obj_constructor basic_const_zero_ffi
-
--- | constructor for 1
-one :: BasicSym
-one = basic_obj_constructor basic_const_one_ffi
-
--- | constructor for -1
-minus_one :: BasicSym
-minus_one = basic_obj_constructor basic_const_minus_one_ffi
-
--- | constructor for i = sqrt(-1)
-im :: BasicSym
-im = basic_obj_constructor basic_const_I_ffi
-
--- | the ratio of the circumference of a circle to its radius
-pi :: BasicSym
-pi = basic_obj_constructor basic_const_pi_ffi
-
--- | The base of the natural logarithm
-e :: BasicSym
-e = basic_obj_constructor basic_const_E_ffi
-
-expand :: BasicSym -> BasicSym
-expand = basic_unaryop basic_expand_ffi
-
-
-eulerGamma :: BasicSym
-eulerGamma = basic_obj_constructor basic_const_EulerGamma_ffi
-
-basic_obj_constructor :: (Ptr BasicStruct -> IO ()) -> BasicSym
-basic_obj_constructor init_fn = unsafePerformIO $ do
-    basic_ptr <- create_basic_ptr
-    withBasicSym basic_ptr init_fn
-    return basic_ptr
-
-basic_str :: BasicSym -> String
-basic_str basic_ptr = unsafePerformIO $ withBasicSym basic_ptr (basic_str_ffi >=> peekCString)
-
-integerToCLong :: Integer -> CLong
-integerToCLong i = CLong (fromInteger i)
-
-
-intToCLong :: Int -> CLong
-intToCLong i = integerToCLong (toInteger i)
-
-basic_int_signed :: Int -> BasicSym
-basic_int_signed i = unsafePerformIO $ do
-    iptr <- create_basic_ptr
-    withBasicSym iptr (\iptr -> integer_set_si_ffi iptr (intToCLong i) )
-    return iptr
-
-
-basic_from_integer :: Integer -> BasicSym
-basic_from_integer i = unsafePerformIO $ do
-    iptr <- create_basic_ptr
-    withBasicSym iptr (\iptr -> integer_set_si_ffi iptr (fromInteger i))
-    return iptr
-
--- |The `ascii_art_str` function prints SymEngine in ASCII art.
--- this is useful as a sanity check
-ascii_art_str :: IO String
-ascii_art_str = ascii_art_str_ffi >>= peekCString
-
--- Unexported ffi functions------------------------
-
--- |Create a basic object that represents all other objects through
--- the FFI
-create_basic_ptr :: IO BasicSym
-create_basic_ptr = do
-    basic_ptr <- newArray [BasicStruct { data_ptr = nullPtr }]
-    basic_new_heap_ffi basic_ptr
-    finalized_ptr <- newForeignPtr ptr_basic_free_heap_ffi basic_ptr
-    return $ BasicSym { fptr = finalized_ptr }
-
-basic_binaryop :: (Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()) -> BasicSym -> BasicSym -> BasicSym
-basic_binaryop f a b = unsafePerformIO $ do
-    s <- create_basic_ptr
-    withBasicSym3 s a b f
-    return s 
-
-basic_unaryop :: (Ptr BasicStruct -> Ptr BasicStruct -> IO ()) -> BasicSym -> BasicSym
-basic_unaryop f a = unsafePerformIO $ do
-    s <- create_basic_ptr
-    withBasicSym2 s a f
-    return s 
-
-
-basic_pow :: BasicSym -> BasicSym -> BasicSym
-basic_pow = basic_binaryop basic_pow_ffi
-
--- |Create a rational number with numerator and denominator
-rational :: BasicSym -> BasicSym -> BasicSym
-rational = basic_binaryop rational_set_ffi
-
--- |Create a complex number a + b * im
-complex :: BasicSym -> BasicSym -> BasicSym
-complex a b = (basic_binaryop complex_set_ffi) a b
-
-basic_rational_from_integer :: Integer -> Integer -> BasicSym
-basic_rational_from_integer i j = unsafePerformIO $ do
-    s <- create_basic_ptr
-    withBasicSym s (\s -> rational_set_si_ffi s (integerToCLong i) (integerToCLong j))
-    return s 
-
--- |Create a symbol with the given name
-symbol :: String -> BasicSym
-symbol name = unsafePerformIO $ do
-    s <- create_basic_ptr
-    cname <- newCString name
-    withBasicSym s (\s -> symbol_set_ffi s cname)
-    free cname
-    return s
-
--- |Differentiate an expression with respect to a symbol
-diff :: BasicSym -> BasicSym -> BasicSym
-diff expr symbol = (basic_binaryop basic_diff_ffi) expr symbol
-
-instance Show BasicSym where
-    show = basic_str 
-
-instance Eq BasicSym where
-    (==) a b = unsafePerformIO $ do 
-                i <- withBasicSym2 a b basic_eq_ffi
-                return $ i == 1
-
-
-instance Num BasicSym where
-    (+) = basic_binaryop basic_add_ffi
-    (-) = basic_binaryop basic_sub_ffi
-    (*) = basic_binaryop basic_mul_ffi
-    negate = basic_unaryop basic_neg_ffi
-    abs = basic_unaryop basic_abs_ffi
-    signum = undefined
-    fromInteger = basic_from_integer
-
-instance Fractional BasicSym where
-    (/) = basic_binaryop basic_div_ffi
-    fromRational (num :% denom) = basic_rational_from_integer num denom
-    recip r = one / r
-
-instance Floating BasicSym where
-    pi = Symengine.pi
-    exp x = e ** x
-    log = undefined
-    sqrt x = x  ** 1/2
-    (**) = basic_pow
-    logBase = undefined
-    sin = basic_unaryop basic_sin_ffi
-    cos = basic_unaryop basic_cos_ffi
-    tan = basic_unaryop basic_tan_ffi
-    asin = basic_unaryop basic_asin_ffi
-    acos = basic_unaryop basic_acos_ffi
-    atan = basic_unaryop basic_atan_ffi
-    sinh = basic_unaryop basic_sinh_ffi
-    cosh = basic_unaryop basic_cosh_ffi
-    tanh = basic_unaryop basic_tanh_ffi
-    asinh = basic_unaryop basic_asinh_ffi
-    acosh = basic_unaryop basic_acosh_ffi
-    atanh = basic_unaryop basic_atanh_ffi
-
-foreign import ccall "symengine/cwrapper.h ascii_art_str" ascii_art_str_ffi :: IO CString
-foreign import ccall "symengine/cwrapper.h basic_new_heap" basic_new_heap_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h &basic_free_heap" ptr_basic_free_heap_ffi :: FunPtr(Ptr BasicStruct -> IO ())
-
--- constants
-foreign import ccall "symengine/cwrapper.h basic_const_zero" basic_const_zero_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_one" basic_const_one_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_minus_one" basic_const_minus_one_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_I" basic_const_I_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_pi" basic_const_pi_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_E" basic_const_E_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_const_EulerGamma" basic_const_EulerGamma_ffi :: Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_str" basic_str_ffi :: Ptr BasicStruct -> IO CString
-foreign import ccall "symengine/cwrapper.h basic_eq" basic_eq_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO Int
-
-foreign import ccall "symengine/cwrapper.h symbol_set" symbol_set_ffi :: Ptr BasicStruct -> CString -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_diff" basic_diff_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h integer_set_si" integer_set_si_ffi :: Ptr BasicStruct -> CLong -> IO ()
-
-foreign import ccall "symengine/cwrapper.h rational_set" rational_set_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h rational_set_si" rational_set_si_ffi :: Ptr BasicStruct -> CLong -> CLong -> IO ()
-
-foreign import ccall "symengine/cwrapper.h complex_set" complex_set_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h basic_expand" basic_expand_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-
-foreign import ccall "symengine/cwrapper.h basic_add" basic_add_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_sub" basic_sub_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_mul" basic_mul_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_div" basic_div_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_pow" basic_pow_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_neg" basic_neg_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_abs" basic_abs_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h basic_sin" basic_sin_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_cos" basic_cos_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_tan" basic_tan_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h basic_asin" basic_asin_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_acos" basic_acos_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_atan" basic_atan_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h basic_sinh" basic_sinh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_cosh" basic_cosh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_tanh" basic_tanh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-
-foreign import ccall "symengine/cwrapper.h basic_asinh" basic_asinh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_acosh" basic_acosh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
-foreign import ccall "symengine/cwrapper.h basic_atanh" basic_atanh_ffi :: Ptr BasicStruct -> Ptr BasicStruct -> IO ()
+import Symengine.Internal
diff --git a/src/Symengine/BasicSym.hs b/src/Symengine/BasicSym.hs
new file mode 100644
index 0000000..8697639
--- /dev/null
+++ b/src/Symengine/BasicSym.hs
@@ -0,0 +1,267 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+
+module Symengine.BasicSym(
+     ascii_art_str,
+     zero,
+     one,
+     im,
+     Symengine.BasicSym.pi,
+     e,
+     minus_one,
+     rational,
+     complex,
+     symbol_new,
+     diff,
+     expand,
+     -- HACK: this should be internal :(
+     basicsym_new,
+     BasicSym,
+     lift_basicsym_binaryop,
+     lift_basicsym_unaryop
+)
+where
+
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.C.String
+import Foreign.Storable
+import Foreign.Marshal.Array
+import Foreign.Marshal.Alloc
+import Foreign.ForeignPtr
+import Control.Applicative
+import Control.Monad -- for foldM
+import System.IO.Unsafe
+import Control.Monad
+import GHC.Real
+
+import Symengine.Internal
+
+
+data BasicSym = BasicSym !(ForeignPtr CBasicSym)
+instance Wrapped BasicSym CBasicSym where
+  with (BasicSym (p)) f = withForeignPtr p f
+  
+-- | constructor for 0
+zero :: BasicSym
+zero = basicsym_construct basic_const_zero_ffi
+
+-- | constructor for 1
+one :: BasicSym
+one = basicsym_construct basic_const_one_ffi
+
+-- | constructor for -1
+minus_one :: BasicSym
+minus_one = basicsym_construct basic_const_minus_one_ffi
+
+-- | constructor for i = sqrt(-1)
+im :: BasicSym
+im = basicsym_construct basic_const_I_ffi
+
+-- | the ratio of the circumference of a circle to its radius
+pi :: BasicSym
+pi = basicsym_construct basic_const_pi_ffi
+
+-- | The base of the natural logarithm
+e :: BasicSym
+e = basicsym_construct basic_const_E_ffi
+
+expand :: BasicSym -> BasicSym
+expand = lift_basicsym_unaryop basic_expand_ffi
+
+
+eulerGamma :: BasicSym
+eulerGamma = basicsym_construct basic_const_EulerGamma_ffi
+
+basicsym_construct :: (Ptr CBasicSym -> IO ()) -> BasicSym
+basicsym_construct init_fn = unsafePerformIO $ do
+    basic_ptr <- basicsym_new
+    with basic_ptr init_fn
+    return basic_ptr
+
+basic_str :: BasicSym -> String
+basic_str basic_ptr = unsafePerformIO $ with basic_ptr (basic_str_ffi >=> peekCString)
+
+integerToCLong :: Integer -> CLong
+integerToCLong i = CLong (fromInteger i)
+
+
+intToCLong :: Int -> CLong
+intToCLong i = toEnum i
+
+
+intToCInt :: Int -> CInt
+intToCInt i = toEnum i
+
+basic_int_signed :: Int -> BasicSym
+basic_int_signed i = unsafePerformIO $ do
+    iptr <- basicsym_new
+    with iptr (\iptr -> integer_set_si_ffi iptr (intToCLong i) )
+    return iptr
+
+
+basic_from_integer :: Integer -> BasicSym
+basic_from_integer i = unsafePerformIO $ do
+  s <- basicsym_new
+  with s (\s -> integer_set_si_ffi s (fromInteger i))
+  return s 
+
+
+-- |The `ascii_art_str` function prints SymEngine in ASCII art.
+-- this is useful as a sanity check
+ascii_art_str :: IO String
+ascii_art_str = ascii_art_str_ffi >>= peekCString
+
+-- Unexported ffi functions------------------------
+
+-- |Create a basic object that represents all other objects through
+-- the FFI
+basicsym_new :: IO BasicSym
+basicsym_new = do
+    basic_ptr <- basic_new_heap_ffi
+    finalized_ptr <- newForeignPtr ptr_basic_free_heap_ffi  basic_ptr
+
+    return $ BasicSym finalized_ptr
+
+-- NOTE: throws exception
+lift_basicsym_binaryop :: (Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt) -> 
+    BasicSym -> BasicSym -> BasicSym
+lift_basicsym_binaryop f a b = unsafePerformIO $ do
+    s <- basicsym_new
+    with3 s a b f >>= throwOnSymIntException
+
+    return s
+
+lift_basicsym_unaryop :: (Ptr CBasicSym -> Ptr CBasicSym -> IO CInt) -> BasicSym -> BasicSym
+lift_basicsym_unaryop f a = unsafePerformIO $ do
+    s <- basicsym_new 
+    with2 s a f >>= throwOnSymIntException
+    return $ s
+
+
+basic_pow :: BasicSym -> BasicSym -> BasicSym
+basic_pow = lift_basicsym_binaryop basic_pow_ffi
+
+-- |Create a rational number with numerator and denominator
+rational :: BasicSym -> BasicSym -> BasicSym
+rational = lift_basicsym_binaryop rational_set_ffi
+
+-- |Create a complex number a + b * im
+complex :: BasicSym -> BasicSym -> BasicSym
+complex a b = (lift_basicsym_binaryop complex_set_ffi) a b
+
+basic_rational_from_integer :: Integer -> Integer -> BasicSym
+basic_rational_from_integer i j = unsafePerformIO $ do
+    s <- basicsym_new
+    with s (\s -> rational_set_si_ffi s (integerToCLong i) (integerToCLong j))
+    return s 
+
+-- |Create a symbol with the given name
+symbol_new :: String -> BasicSym
+symbol_new name = unsafePerformIO $ do
+    s <- basicsym_new
+    cname <- newCString name
+    with s (\s -> symbol_set_ffi s cname)
+    free cname
+    return s
+
+-- |Differentiate an expression with respect to a symbol
+diff :: BasicSym -> BasicSym -> BasicSym
+diff expr symbol = (lift_basicsym_binaryop basic_diff_ffi) expr symbol
+
+instance Show BasicSym where
+    show = basic_str
+
+instance Eq BasicSym where
+    (==) a b = unsafePerformIO $ do
+                i <- with2 a b basic_eq_ffi
+                return $ i == 1
+
+instance Num BasicSym where
+    (+) = lift_basicsym_binaryop $ basic_add_ffi
+    (-) = lift_basicsym_binaryop $ basic_sub_ffi
+    (*) = lift_basicsym_binaryop $ basic_mul_ffi
+    negate = lift_basicsym_unaryop basic_neg_ffi
+    abs = lift_basicsym_unaryop basic_abs_ffi
+
+    -- works only for long [-2^32, 2^32 - 1]
+    fromInteger = basic_from_integer
+
+instance Fractional BasicSym where
+    (/) = lift_basicsym_binaryop $ basic_div_ffi
+    fromRational (num :% denom) = basic_rational_from_integer num denom
+    recip r = one / r
+
+instance Floating BasicSym where
+    pi = Symengine.BasicSym.pi
+    exp x = e ** x
+    log = undefined
+    sqrt x = x  ** 1/2
+    (**) = basic_pow
+    logBase = undefined
+    sin = lift_basicsym_unaryop basic_sin_ffi
+    cos = lift_basicsym_unaryop basic_cos_ffi
+    tan = lift_basicsym_unaryop basic_tan_ffi
+    asin = lift_basicsym_unaryop basic_asin_ffi
+    acos = lift_basicsym_unaryop basic_acos_ffi
+    atan = lift_basicsym_unaryop basic_atan_ffi
+    sinh = lift_basicsym_unaryop basic_sinh_ffi
+    cosh = lift_basicsym_unaryop basic_cosh_ffi
+    tanh = lift_basicsym_unaryop basic_tanh_ffi
+    asinh = lift_basicsym_unaryop basic_asinh_ffi
+    acosh = lift_basicsym_unaryop basic_acosh_ffi
+    atanh = lift_basicsym_unaryop basic_atanh_ffi
+
+foreign import ccall "symengine/cwrapper.h ascii_art_str" ascii_art_str_ffi :: IO CString
+foreign import ccall "symengine/cwrapper.h basic_new_heap" basic_new_heap_ffi :: IO (Ptr CBasicSym)
+foreign import ccall "symengine/cwrapper.h &basic_free_heap" ptr_basic_free_heap_ffi :: FunPtr(Ptr CBasicSym -> IO ())
+
+-- constants
+foreign import ccall "symengine/cwrapper.h basic_const_zero" basic_const_zero_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_one" basic_const_one_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_minus_one" basic_const_minus_one_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_I" basic_const_I_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_pi" basic_const_pi_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_E" basic_const_E_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_const_EulerGamma" basic_const_EulerGamma_ffi :: Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h basic_str" basic_str_ffi :: Ptr CBasicSym -> IO CString
+foreign import ccall "symengine/cwrapper.h basic_eq" basic_eq_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO Int
+
+foreign import ccall "symengine/cwrapper.h symbol_set" symbol_set_ffi :: Ptr CBasicSym -> CString -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_diff" basic_diff_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h integer_set_si" integer_set_si_ffi :: Ptr CBasicSym -> CLong -> IO ()
+
+foreign import ccall "symengine/cwrapper.h rational_set" rational_set_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h rational_set_si" rational_set_si_ffi :: Ptr CBasicSym -> CLong -> CLong -> IO ()
+
+foreign import ccall "symengine/cwrapper.h complex_set" complex_set_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h basic_expand" basic_expand_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h basic_add" basic_add_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_sub" basic_sub_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_mul" basic_mul_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_div" basic_div_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_pow" basic_pow_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_neg" basic_neg_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_abs" basic_abs_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h basic_sin" basic_sin_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_cos" basic_cos_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_tan" basic_tan_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h basic_asin" basic_asin_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_acos" basic_acos_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_atan" basic_atan_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h basic_sinh" basic_sinh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_cosh" basic_cosh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_tanh" basic_tanh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h basic_asinh" basic_asinh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_acosh" basic_acosh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h basic_atanh" basic_atanh_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
diff --git a/src/Symengine/DenseMatrix.hs b/src/Symengine/DenseMatrix.hs
new file mode 100644
index 0000000..f1113bd
--- /dev/null
+++ b/src/Symengine/DenseMatrix.hs
@@ -0,0 +1,346 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeOperators #-}
+-- to write things like KnownNat(r * c) => ...
+{-# LANGUAGE FlexibleContexts #-}
+-- @
+{-# LANGUAGE TypeApplications #-}
+-- to bring stuff like (r, c) into scope
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- allow non injective type functions (+)
+{-# LANGUAGE AllowAmbiguousTypes #-}
+
+-- data declarations that are empty
+{-# LANGUAGE EmptyDataDecls #-}
+module Symengine.DenseMatrix
+  (
+   DenseMatrix,
+   -- densematrix_new,
+   densematrix_new_vec,
+   densematrix_new_eye,
+   densematrix_new_diag,
+   densematrix_new_zeros,
+   densematrix_get,
+   densematrix_set,
+   densematrix_size,
+   
+   -- arithmetic
+   densematrix_add,
+   densematrix_mul_matrix,
+   densematrix_mul_scalar,
+   det,
+   inv,
+   transpose,
+
+   --decomposition
+   L(L), D(D), U(U),
+   densematrix_lu,
+   densematrix_ldl,
+   densematrix_fflu,
+   densematrix_ffldu,
+   densematrix_lu_solve,
+
+   -- custom matrix class
+   Matrix(..)
+
+   --
+   )
+where
+
+import Prelude
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.C.String
+import Foreign.Storable
+import Foreign.Marshal.Array
+import Foreign.Marshal.Alloc
+import Foreign.ForeignPtr
+import Control.Applicative
+import Control.Monad -- for foldM
+import System.IO.Unsafe
+import Control.Monad
+import GHC.Real
+import Data.Proxy
+
+import GHC.TypeLits -- type level programming
+import qualified Data.Vector.Sized as V -- sized vectors
+import Data.Finite -- types to represent numbers
+
+import Symengine.Internal
+import Symengine.BasicSym
+import Symengine.VecBasic
+
+class Matrix m where
+  (<>) :: (KnownNat r, KnownNat c, KnownNat k) => m r k -> m k c -> m r c 
+
+
+instance Matrix (DenseMatrix) where
+  (<>) = densematrix_mul_matrix
+
+data CDenseMatrix
+data  DenseMatrix :: Nat -> Nat -> * where
+  -- allow constructing raw DenseMatrix from a constructor
+  DenseMatrix :: (KnownNat r, KnownNat c) => (ForeignPtr CDenseMatrix) -> DenseMatrix r c
+
+instance  (KnownNat r, KnownNat c) => Wrapped (DenseMatrix  r c) CDenseMatrix where
+    with (DenseMatrix p) f = withForeignPtr p f
+
+instance  (KnownNat r, KnownNat c) =>  Show (DenseMatrix r c) where
+  show :: DenseMatrix r c -> String
+  show mat =
+    unsafePerformIO $ with mat  (cdensematrix_str_ffi >=> peekCString)
+
+instance (KnownNat r, KnownNat c) => Eq (DenseMatrix r c) where
+  (==) :: DenseMatrix r c -> DenseMatrix r c -> Bool
+  (==) mat1 mat2 = 
+    1 == fromIntegral (unsafePerformIO $
+                       with2 mat1 mat2 cdensematrix_eq_ffi)
+
+instance (KnownNat r, KnownNat c) => Num (DenseMatrix r c) where
+  (+) = densematrix_add
+  (-) d1 d2 = let
+        d2_neg =  densematrix_mul_scalar d2 (fromInteger (-1))
+        in d1 + d2_neg
+  -- TODO: Should be elementwise multiplcation
+  (*) = undefined
+  -- TODO: should be elementwise signum
+  signum = undefined
+  -- TODO: should be elementwise abs
+  abs = undefined
+  -- make a 1x1 matrix
+  fromInteger = undefined -- densematrix_new_vec 
+
+densematrix_new :: (KnownNat r, KnownNat c) => IO (DenseMatrix r c)
+densematrix_new = DenseMatrix <$> (mkForeignPtr cdensematrix_new_ffi cdensematrix_free_ffi)
+
+_densematrix_copy :: (KnownNat r, KnownNat c) => DenseMatrix r c -> IO (DenseMatrix r c)
+_densematrix_copy mat = do
+   newmat <- densematrix_new
+   throwOnSymIntException =<< with2 newmat mat cdensematrix_set_ffi
+   return newmat
+
+densematrix_new_rows_cols ::  forall r c . (KnownNat r, KnownNat c) => DenseMatrix r c
+densematrix_new_rows_cols = 
+    unsafePerformIO $ DenseMatrix <$>
+      (mkForeignPtr (cdensematrix_new_rows_cols_ffi 
+                      (fromIntegral . natVal $ (Proxy @ r))
+                      (fromIntegral . natVal $ (Proxy @ c)))
+                      cdensematrix_free_ffi)
+  
+
+densematrix_new_vec :: forall r c. (KnownNat r, KnownNat c, KnownNat (r * c)) => V.Vector (r * c) BasicSym -> DenseMatrix r c
+densematrix_new_vec syms = unsafePerformIO $ do
+  vec <- vector_to_vecbasic syms
+  let cdensemat =  with vec (\v -> cdensematrix_new_vec_ffi
+                                     (fromIntegral . natVal $ (Proxy @ r))
+                                     (fromIntegral . natVal $ (Proxy @ c)) v)
+  DenseMatrix <$>  mkForeignPtr cdensemat cdensematrix_free_ffi
+
+
+type Offset = Int
+-- |create a matrix with rows 'r, cols 'c' and offset 'k'
+densematrix_new_eye :: forall k r c. (KnownNat r,  KnownNat c, KnownNat k, KnownNat (r + k), KnownNat (c + k)) => DenseMatrix (r + k) (c + k)
+densematrix_new_eye = unsafePerformIO $ do
+  let mat = densematrix_new_rows_cols
+  throwOnSymIntException =<< with mat (\m -> cdensematrix_eye_ffi m
+                 (fromIntegral . natVal $ (Proxy @ r))
+                 (fromIntegral . natVal $ (Proxy @ c))
+                 (fromIntegral . natVal $ (Proxy @ k)))
+
+
+  return mat
+
+densematrix_new_zeros :: forall r c. (KnownNat r, KnownNat c) => DenseMatrix r c
+densematrix_new_zeros = unsafePerformIO $ do
+  let mat = densematrix_new_rows_cols
+  throwOnSymIntException =<< with mat (\m -> cdensematrix_zeros_ffi m
+                  (fromIntegral . natVal $ (Proxy @ r))
+                  (fromIntegral . natVal $ (Proxy @ c)))
+  return mat
+
+-- create a matrix with diagonal elements of length 'd', offset 'k'
+densematrix_new_diag :: forall k d. (KnownNat d, KnownNat k, KnownNat (d + k)) => V.Vector d BasicSym -> DenseMatrix (d + k) (d + k)
+densematrix_new_diag syms  = unsafePerformIO $ do
+  let offset = fromIntegral $ natVal (Proxy @ k)
+  let diagonal = fromIntegral $ natVal (Proxy @ d)
+  let dim = offset + diagonal
+  vecsyms <- vector_to_vecbasic syms
+  let mat = densematrix_new_rows_cols :: DenseMatrix (d + k) (d + k)
+  throwOnSymIntException =<< with2 mat vecsyms (\m syms -> cdensematrix_diag_ffi m syms offset)
+
+
+  return mat
+
+type Row = Int
+type Col = Int
+
+
+densematrix_get :: forall r c. (KnownNat r, KnownNat c) => 
+  DenseMatrix r c -> Finite r -> Finite c -> BasicSym
+densematrix_get mat getr getc = unsafePerformIO $ do
+      sym <- basicsym_new
+      let indexr = fromIntegral $ (getFinite getr)
+      let indexc = fromIntegral $ (getFinite getc)
+      throwOnSymIntException =<< with2 mat sym (\m s -> cdensematrix_get_basic_ffi s m indexr indexc)
+
+      return sym
+
+densematrix_set :: forall r c. (KnownNat r, KnownNat c) => 
+  DenseMatrix r c -> Finite r -> Finite c -> BasicSym -> DenseMatrix r c
+densematrix_set mat r c sym = unsafePerformIO $ do
+    mat' <- _densematrix_copy mat
+    throwOnSymIntException =<< with2 mat' sym (\m s -> cdensematrix_set_basic_ffi 
+                              m
+                              (fromIntegral . getFinite $ r)
+                              (fromIntegral . getFinite $ c)
+                              s)
+
+    return mat'
+
+
+type NRows = Int
+type NCols = Int
+
+-- | provides dimenions of matrix. combination of the FFI calls
+-- `dense_matrix_rows` and `dense_matrix_cols`
+densematrix_size :: forall r c. (KnownNat r, KnownNat c) => DenseMatrix r c -> (NRows, NCols)
+densematrix_size mat = 
+   (fromIntegral . natVal $ (Proxy @ r), fromIntegral . natVal $ (Proxy @ c))
+
+densematrix_add :: forall r c. (KnownNat r, KnownNat c) => 
+  DenseMatrix r c -> DenseMatrix r c -> DenseMatrix r c
+densematrix_add mata matb = unsafePerformIO $ do
+   res <- densematrix_new
+   throwOnSymIntException =<< with3 res mata matb cdensematrix_add_matrix_ffi
+   return res
+
+
+densematrix_mul_matrix :: forall r k c. (KnownNat r, KnownNat k, KnownNat c) => 
+  DenseMatrix r k -> DenseMatrix k c -> DenseMatrix r c
+densematrix_mul_matrix mata matb = unsafePerformIO $ do
+   res <- densematrix_new
+   throwOnSymIntException =<< with3 res mata matb cdensematrix_mul_matrix_ffi
+   return res
+
+
+densematrix_mul_scalar :: forall r c. (KnownNat r, KnownNat c) => 
+  DenseMatrix r c -> BasicSym -> DenseMatrix r c
+densematrix_mul_scalar mata sym = unsafePerformIO $ do
+   res <- densematrix_new
+   throwOnSymIntException =<< with3 res mata sym cdensematrix_mul_scalar_ffi
+   return res
+
+det :: forall r c. (KnownNat r, KnownNat c) => DenseMatrix r c -> BasicSym
+det d = unsafePerformIO $ do
+  sym <- basicsym_new
+  throwOnSymIntException =<< with2 sym d cdensematrix_det_ffi 
+  return sym
+
+inv :: forall r c. (KnownNat r, KnownNat c) => DenseMatrix r c -> DenseMatrix r c
+inv d = unsafePerformIO $ do
+  m <- densematrix_new
+  throwOnSymIntException =<< with2 m d cdensematrix_inv_ffi
+  return m
+
+transpose :: forall r c. (KnownNat r, KnownNat c) => DenseMatrix r c -> DenseMatrix r c
+transpose d = unsafePerformIO $ do
+  m <- densematrix_new
+  throwOnSymIntException =<< with2 m d cdensematrix_transpose_ffi
+  return m
+
+newtype L r c = L (DenseMatrix r c)
+newtype U r c = U (DenseMatrix r c)
+
+densematrix_lu :: (KnownNat r, KnownNat c) => DenseMatrix r c-> (L r c, U r c)
+densematrix_lu mat = unsafePerformIO $ do
+   l <- densematrix_new
+   u <- densematrix_new
+   throwOnSymIntException =<< with3 l u mat cdensematrix_lu
+   return (L l, U u)
+
+newtype D r c = D (DenseMatrix r c)
+densematrix_ldl :: (KnownNat r, KnownNat c) => DenseMatrix r c-> (L r c, D r c)
+densematrix_ldl mat = unsafePerformIO $ do
+  l <- densematrix_new
+  d <- densematrix_new
+  throwOnSymIntException =<< with3 l d mat cdensematrix_ldl
+
+  return (L l, D d)
+
+
+newtype FFLU r c = FFLU (DenseMatrix r c)
+densematrix_fflu :: (KnownNat r, KnownNat c) => DenseMatrix r c -> FFLU r c
+densematrix_fflu mat = unsafePerformIO $ do
+  fflu <- densematrix_new
+  throwOnSymIntException =<< with2 fflu mat cdensematrix_fflu
+  return (FFLU fflu)
+
+
+densematrix_ffldu ::  (KnownNat r, KnownNat c) =>
+  DenseMatrix r c -> (L r c, D r c, U r c)
+densematrix_ffldu mat = unsafePerformIO $ do
+  l <- densematrix_new
+  d <- densematrix_new
+  u <- densematrix_new
+
+  throwOnSymIntException =<< with4 l d u mat cdensematrix_ffldu
+  return (L l, D d, U u)
+
+-- solve A x = B
+-- A is first param, B is second larameter
+densematrix_lu_solve :: (KnownNat r, KnownNat c) => 
+  DenseMatrix r c -> DenseMatrix r c -> DenseMatrix r c
+densematrix_lu_solve a b = unsafePerformIO $ do
+  x <- densematrix_new
+  throwOnSymIntException =<< with3 x a b cdensematrix_lu_solve
+  return x
+
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_new" cdensematrix_new_ffi :: IO (Ptr CDenseMatrix)
+foreign import ccall unsafe "symengine/cwrapper.h &dense_matrix_free" cdensematrix_free_ffi :: FunPtr ((Ptr CDenseMatrix) -> IO ())
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_new_rows_cols" cdensematrix_new_rows_cols_ffi :: CUInt -> CUInt -> IO (Ptr CDenseMatrix)
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_new_vec" cdensematrix_new_vec_ffi :: CUInt -> CUInt -> Ptr CVecBasic -> IO (Ptr CDenseMatrix)
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_zeros" cdensematrix_zeros_ffi :: Ptr CDenseMatrix -> CULong -> CULong -> IO CInt
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_eye" cdensematrix_eye_ffi :: Ptr CDenseMatrix -> CULong -> CULong  -> CULong -> IO CInt
+foreign import ccall unsafe "symengine/cwrapper.h dense_matrix_diag" cdensematrix_diag_ffi :: Ptr CDenseMatrix -> Ptr CVecBasic -> CULong  -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_eq" cdensematrix_eq_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_set" cdensematrix_set_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_str" cdensematrix_str_ffi :: Ptr CDenseMatrix -> IO CString
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_get_basic" cdensematrix_get_basic_ffi :: Ptr (CBasicSym)  -> Ptr CDenseMatrix -> CUInt -> CUInt -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_set_basic" cdensematrix_set_basic_ffi :: Ptr CDenseMatrix -> CUInt -> CUInt -> Ptr (CBasicSym)  -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_rows" cdensematrix_rows_ffi :: Ptr CDenseMatrix -> IO CULong
+foreign import ccall "symengine/cwrapper.h dense_matrix_cols" cdensematrix_cols_ffi :: Ptr CDenseMatrix -> IO CULong
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_add_matrix"
+  cdensematrix_add_matrix_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_mul_matrix"
+  cdensematrix_mul_matrix_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_mul_scalar"
+  cdensematrix_mul_scalar_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_det"
+  cdensematrix_det_ffi :: Ptr CBasicSym -> Ptr CDenseMatrix -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_inv"
+  cdensematrix_inv_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_transpose"
+  cdensematrix_transpose_ffi :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h dense_matrix_LU" cdensematrix_lu :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_LDL" cdensematrix_ldl :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_FFLU" cdensematrix_fflu :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_FFLDU" cdensematrix_ffldu :: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
+foreign import ccall "symengine/cwrapper.h dense_matrix_LU_solve" cdensematrix_lu_solve:: Ptr CDenseMatrix -> Ptr CDenseMatrix -> Ptr CDenseMatrix -> IO CInt
diff --git a/src/Symengine/Internal.hs b/src/Symengine/Internal.hs
new file mode 100644
index 0000000..2393b1f
--- /dev/null
+++ b/src/Symengine/Internal.hs
@@ -0,0 +1,94 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+-- data declarations that are empty
+{-# LANGUAGE EmptyDataDecls #-}
+
+module Symengine.Internal
+  (
+    cIntToEnum,
+    cIntFromEnum,
+    mkForeignPtr,
+    Wrapped(..),
+    with2,
+    with3,
+    with4,
+    CBasicSym,
+    CVecBasic,
+    SymengineException(NoException, RuntimeError, DivByZero, NotImplemented, DomainError, ParseError),
+    forceException,
+    throwOnSymIntException
+  ) where
+
+import Prelude
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.C.String
+import Foreign.Storable
+import Foreign.Marshal.Array
+import Foreign.Marshal.Alloc
+import Foreign.ForeignPtr
+import Control.Applicative
+import Control.Monad -- for foldM
+import System.IO.Unsafe
+import Control.Monad
+import GHC.Real
+import Control.Exception
+import Data.Typeable
+
+data SymengineException = NoException |
+                           RuntimeError |
+                           DivByZero |
+                           NotImplemented |
+                           DomainError |
+                           ParseError deriving (Show, Enum, Eq, Typeable)
+
+instance Exception SymengineException
+
+
+-- interpret the CInt as a SymengineException, and
+-- throw if it is actually an error
+throwOnSymIntException :: CInt -> IO ()
+throwOnSymIntException i = forceException . cIntToEnum $ i
+
+forceException :: SymengineException -> IO ()
+forceException exception = 
+  case exception of
+  NoException -> return ()
+  error @ _ -> throwIO error
+
+cIntToEnum :: Enum a => CInt -> a
+cIntToEnum = toEnum . fromIntegral
+
+cIntFromEnum :: Enum a => a -> CInt
+cIntFromEnum = fromIntegral . fromEnum
+
+-- |given a raw pointer IO (Ptr a) and a destructor function pointer, make a
+-- foreign pointer
+mkForeignPtr :: (IO (Ptr a)) -> FunPtr (Ptr a -> IO ()) -> IO (ForeignPtr a)
+mkForeignPtr cons des = do
+  rawptr <- cons
+  finalized <- newForeignPtr des rawptr
+  return  finalized
+
+class Wrapped o i | o -> i where
+  with :: o -> (Ptr i -> IO a) -> IO a
+
+with2 :: Wrapped o1 i1 => Wrapped o2 i2 => o1 -> o2 -> (Ptr i1 -> Ptr i2 -> IO a) -> IO a
+with2 o1 o2 f = with o1 (\p1 -> with o2 (\p2 -> f p1 p2))
+
+with3 :: Wrapped o1 i1 => Wrapped o2 i2 => Wrapped o3 i3 => o1 -> o2 -> o3 -> (Ptr i1 -> Ptr i2 -> Ptr i3 -> IO a) -> IO a
+with3 o1 o2 o3 f = with2 o1 o2 (\p1 p2 -> with o3 (\p3 -> f p1 p2 p3))
+
+
+with4:: Wrapped o1 i1 => Wrapped o2 i2 => Wrapped o3 i3 => Wrapped o4 i4 => o1 -> o2 -> o3 -> o4 -> (Ptr i1 -> Ptr i2 -> Ptr i3 -> Ptr i4 -> IO a) -> IO a
+with4 o1 o2 o3 o4 f = with o1 (\p1 -> with3 o2 o3 o4 (\p2 p3 p4 -> f p1 p2 p3 p4))
+
+-- BasicSym
+data CBasicSym
+
+-- VecBasic
+data CVecBasic
+
+-- CDenseMatrix
diff --git a/src/Symengine/NumberTheory.hs b/src/Symengine/NumberTheory.hs
new file mode 100644
index 0000000..26d8016
--- /dev/null
+++ b/src/Symengine/NumberTheory.hs
@@ -0,0 +1,223 @@
+
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
+
+module Symengine.NumberTheory(
+  Symengine.NumberTheory.gcd,
+  Symengine.NumberTheory.lcm,
+  gcd_extended,
+  next_prime,
+  Symengine.NumberTheory.mod,
+  quotient,
+  quotient_and_mod,
+  mod_f,
+  quotient_f,
+  quotient_and_mod_f,
+  mod_inverse,
+  fibonacci,
+  fibonacci2,
+  lucas,
+  -- I do not understand exactly what lucas2 does. Clarify and then
+  -- export 
+  -- lucas2,
+  binomial,
+  factorial
+)
+where
+
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.C.String
+import Foreign.Storable
+import Foreign.Marshal.Array
+import Foreign.Marshal.Alloc
+import Foreign.ForeignPtr
+import Control.Applicative
+import Control.Monad -- for foldM
+import System.IO.Unsafe
+import Control.Monad
+import GHC.Real
+
+import Symengine.Internal
+import Symengine.BasicSym
+
+
+gcd :: BasicSym -> BasicSym -> BasicSym
+gcd = lift_basicsym_binaryop ntheory_gcd_ffi
+
+lcm :: BasicSym -> BasicSym -> BasicSym
+lcm = lift_basicsym_binaryop ntheory_lcm_ffi
+
+gcd_extended :: BasicSym -> BasicSym -> (BasicSym, BasicSym, BasicSym)
+gcd_extended a b = unsafePerformIO $ do
+  g <- basicsym_new
+  s <- basicsym_new
+  t <- basicsym_new
+
+  with4 g s t a (\g s t a ->
+                    with b (\b ->
+                              ntheory_gcd_ext_ffi g s t a b))
+  return (g, s, t)
+  
+next_prime :: BasicSym -> BasicSym
+next_prime = lift_basicsym_unaryop ntheory_nextprime_ffi
+
+type Quotient = BasicSym
+type Modulo = BasicSym
+
+mod :: BasicSym -> BasicSym -> Quotient
+mod = lift_basicsym_binaryop ntheory_mod_ffi
+
+quotient :: BasicSym -> BasicSym -> BasicSym
+quotient = lift_basicsym_binaryop ntheory_quotient_ffi
+
+quotient_and_mod :: BasicSym -> BasicSym -> (Quotient, Modulo)
+quotient_and_mod a b = unsafePerformIO $ do
+  quotient <- basicsym_new
+  modulo <- basicsym_new
+  with4 quotient modulo a b ntheory_quotient_mod_ffi
+  return $  (quotient, modulo)
+
+
+mod_f :: BasicSym -> BasicSym -> Quotient
+mod_f = lift_basicsym_binaryop ntheory_mod_f_ffi
+
+quotient_f :: BasicSym -> BasicSym -> BasicSym
+quotient_f = lift_basicsym_binaryop ntheory_quotient_f_ffi
+
+quotient_and_mod_f :: BasicSym -> BasicSym -> (Quotient, Modulo)
+quotient_and_mod_f a b = unsafePerformIO $ do
+  quotient <- basicsym_new
+  modulo <- basicsym_new
+  with4 quotient modulo a b ntheory_quotient_mod_f_ffi
+  return $  (quotient, modulo)
+
+
+mod_inverse :: BasicSym -> BasicSym -> Quotient
+mod_inverse = lift_basicsym_binaryop ntheory_mod_inverse_ffi
+
+
+fibonacci ::  Int -> BasicSym
+fibonacci i = unsafePerformIO $  do
+  fib <- basicsym_new
+  with fib (\fib -> ntheory_fibonacci_ffi fib (fromIntegral i))
+  return fib
+
+fibonacci2 :: Int -> (BasicSym, BasicSym)
+fibonacci2 n  = unsafePerformIO $ do
+  g <- basicsym_new
+  s <- basicsym_new
+  
+  with2 g s (\g s -> ntheory_fibonacci2_ffi g s (fromIntegral n))
+  
+  return (g, s) 
+
+
+lucas :: Int -> BasicSym
+lucas n = unsafePerformIO $ do
+  l <- basicsym_new
+  with l (\l -> ntheory_lucas_ffi l (fromIntegral n))
+  return l
+
+{-
+lucas2 :: BasicSym -> BasicSym -> (BasicSym, BasicSym)
+lucas2 n n_prev = unsafePerformIO $ do
+  g <- basicsym_new
+  s <- basicsym_new
+
+  with4 g s n n_prev ntheory_lucas2_ffi
+  
+  return (g, s) 
+-}
+binomial :: BasicSym -> Int -> BasicSym
+binomial n r = unsafePerformIO $ do
+  ncr <- basicsym_new
+  with2 ncr n (\ncr n -> ntheory_binomial_ffi ncr n (fromIntegral r))
+  return ncr
+
+
+factorial :: Int -> BasicSym
+factorial n = unsafePerformIO $ do
+  fact <- basicsym_new
+  with fact (\fact -> ntheory_factorial_ffi fact (fromIntegral n))
+  return fact
+-- FFI Bindings
+-- gcd, lcm
+
+foreign import ccall "symengine/cwrapper.h ntheory_gcd" ntheory_gcd_ffi :: 
+  Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h ntheory_lcm" ntheory_lcm_ffi :: 
+  Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h ntheory_gcd_ext" ntheory_gcd_ext_ffi 
+  :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> 
+      Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+-- prime
+
+foreign import ccall "symengine/cwrapper.h ntheory_nextprime"
+  ntheory_nextprime_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+-- modulus
+
+foreign import ccall "symengine/cwrapper.h ntheory_mod"
+  ntheory_mod_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h ntheory_quotient"
+  ntheory_quotient_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h ntheory_quotient_mod"
+  ntheory_quotient_mod_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> 
+                          Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+
+-- _f versions (round towards -inf)
+foreign import ccall "symengine/cwrapper.h ntheory_mod_f"
+  ntheory_mod_f_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h ntheory_quotient_f"
+  ntheory_quotient_f_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+foreign import ccall "symengine/cwrapper.h ntheory_quotient_mod_f"
+  ntheory_quotient_mod_f_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> 
+                          Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+-- mod inverse
+foreign import ccall "symengine/cwrapper.h ntheory_mod_inverse"
+  ntheory_mod_inverse_ffi :: Ptr CBasicSym -> Ptr CBasicSym -> Ptr CBasicSym -> IO CInt
+
+
+-- fibonacci
+foreign import ccall "symengine/cwrapper.h ntheory_fibonacci"
+  ntheory_fibonacci_ffi :: Ptr CBasicSym -> 
+    CULong -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h ntheory_fibonacci2"
+  ntheory_fibonacci2_ffi :: Ptr CBasicSym -> Ptr CBasicSym ->  
+    CULong -> IO CInt
+
+-- lucas
+foreign import ccall "symengine/cwrapper.h ntheory_lucas"
+  ntheory_lucas_ffi :: Ptr CBasicSym ->
+    CULong -> IO CInt
+
+
+foreign import ccall "symengine/cwrapper.h ntheory_lucas2"
+  ntheory_lucas2_ffi :: Ptr CBasicSym -> Ptr CBasicSym ->  
+     CULong -> IO CInt
+
+
+-- binomial
+foreign import ccall "symengine/cwrapper.h ntheory_binomial"
+  ntheory_binomial_ffi :: Ptr CBasicSym -> Ptr CBasicSym ->  
+    CULong -> IO CInt
+
+-- factorial
+foreign import ccall "symengine/cwrapper.h ntheory_factorial"
+  ntheory_factorial_ffi :: Ptr CBasicSym -> 
+    CULong -> IO CInt
diff --git a/src/Symengine/VecBasic.hs b/src/Symengine/VecBasic.hs
new file mode 100644
index 0000000..358b665
--- /dev/null
+++ b/src/Symengine/VecBasic.hs
@@ -0,0 +1,112 @@
+{-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE InstanceSigs #-}
+{-# LANGUAGE FunctionalDependencies #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeOperators #-}
+-- @
+{-# LANGUAGE TypeApplications #-}
+-- to bring stuff like (r, c) into scope
+{-# LANGUAGE ScopedTypeVariables #-}
+module Symengine.VecBasic
+  (
+    VecBasic,
+    vecbasic_new,
+    vecbasic_push_back,
+    vecbasic_get,
+    vecbasic_size,
+    vector_to_vecbasic,
+   )
+where
+
+
+import Prelude
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.C.String
+import Foreign.Storable
+import Foreign.Marshal.Array
+import Foreign.Marshal.Alloc
+import Foreign.ForeignPtr
+import Control.Applicative
+import Control.Monad -- for foldM
+import System.IO.Unsafe
+import Control.Monad
+import GHC.Real
+import Symengine
+
+import GHC.TypeLits -- type level programming
+import qualified Data.Vector.Sized as V -- sized vectors
+
+import Symengine.Internal
+import Symengine.BasicSym
+
+-- |represents a symbol exported by SymEngine. create this using the functions
+-- 'zero', 'one', 'minus_one', 'e', 'im', 'rational', 'complex', and also by
+-- constructing a number and converting it to a Symbol
+-- 
+-- >>> 3.5 :: BasicSym
+-- 7/2
+--
+-- >>> rational 2 10
+-- 1 /5
+--
+-- >>> complex 1 2
+-- 1 + 2*I
+-- vectors binding
+
+-- | Represents a Vector of BasicSym
+-- | usually, end-users are not expected to interact directly with VecBasic
+-- | this should at some point be moved to Symengine.Internal
+newtype VecBasic = VecBasic (ForeignPtr CVecBasic)
+
+instance Wrapped VecBasic CVecBasic where
+    with (VecBasic p) f = withForeignPtr p f
+
+
+-- | push back an element into a vector
+vecbasic_push_back :: VecBasic -> BasicSym -> IO ()
+vecbasic_push_back vec sym =  with2 vec sym (\v p ->vecbasic_push_back_ffi v p)
+
+
+-- | get the i'th element out of a vecbasic
+vecbasic_get :: VecBasic -> Int -> Either SymengineException BasicSym
+vecbasic_get vec i =
+  if i >= 0 && i < vecbasic_size vec
+  then
+    unsafePerformIO $ do
+    sym <- basicsym_new
+    exception <- cIntToEnum <$> with2 vec sym (\v s -> vecbasic_get_ffi v i s)
+    case exception of
+      NoException -> return (Right sym)
+      _ -> return (Left exception)
+  else
+    Left RuntimeError
+
+
+-- | Create a new VecBasic
+vecbasic_new :: IO VecBasic
+vecbasic_new = do
+    ptr <- vecbasic_new_ffi
+    finalized <- newForeignPtr vecbasic_free_ffi ptr
+    return $ VecBasic (finalized)
+
+
+vector_to_vecbasic :: forall n. KnownNat n => V.Vector n BasicSym -> IO VecBasic
+vector_to_vecbasic syms = do
+  ptr <- vecbasic_new_ffi
+  forM_ syms (\sym -> with sym (\s -> vecbasic_push_back_ffi ptr s))
+  finalized <- newForeignPtr vecbasic_free_ffi ptr
+  return $ VecBasic finalized
+
+vecbasic_size :: VecBasic -> Int
+vecbasic_size vec = unsafePerformIO $
+  fromIntegral <$> with vec vecbasic_size_ffi
+
+foreign import ccall "symengine/cwrapper.h vecbasic_new" vecbasic_new_ffi :: IO (Ptr CVecBasic)
+foreign import ccall "symengine/cwrapper.h vecbasic_push_back" vecbasic_push_back_ffi :: Ptr CVecBasic -> Ptr CBasicSym -> IO ()
+foreign import ccall "symengine/cwrapper.h vecbasic_get" vecbasic_get_ffi :: Ptr CVecBasic -> Int -> Ptr CBasicSym -> IO CInt
+foreign import ccall "symengine/cwrapper.h vecbasic_size" vecbasic_size_ffi :: Ptr CVecBasic -> IO CSize
+foreign import ccall "symengine/cwrapper.h &vecbasic_free" vecbasic_free_ffi :: FunPtr (Ptr CVecBasic -> IO ())
+
diff --git a/stack.yaml b/stack.yaml
index 7b5a9de..4bc1780 100644
--- a/stack.yaml
+++ b/stack.yaml
@@ -1,29 +1,7 @@
-# For more information, see: https://github.com/commercialhaskell/stack/blob/master/doc/yaml_configuration.md
-
-# Specifies the GHC version and set of packages available (e.g., lts-3.5, nightly-2015-09-21, ghc-7.10.2)
-resolver: lts-3.2
-
-# Local packages, usually specified by relative directory name
+flags: {}
 packages:
 - '.'
-
-# Packages to be pulled from upstream that are not in the resolver (e.g., acme-missiles-0.3)
-extra-deps: []
-
-# Override default flag values for local packages and extra-deps
-flags: {}
-
-# Control whether we use the GHC we find on the path
-# system-ghc: true
-
-# Require a specific version of stack, using version ranges
-# require-stack-version: -any # Default
-# require-stack-version: >= 0.1.4.0
-
-# Override the architecture used by stack, especially useful on Windows
-# arch: i386
-# arch: x86_64
-
-# Extra directories used by stack for building
-# extra-include-dirs: [/path/to/dir]
-# extra-lib-dirs: [/path/to/dir]
+extra-deps:
+- finite-typelits-0.1.0.0
+- vector-sized-0.4.0.0
+resolver: lts-7.12
diff --git a/symengine.cabal b/symengine.cabal
index bd5568a..00c3179 100644
--- a/symengine.cabal
+++ b/symengine.cabal
@@ -15,8 +15,19 @@ cabal-version:       >=1.10
 
 library
   hs-source-dirs:      src
-  exposed-modules:     Symengine
-  build-depends:       base >= 4.5.0 && <= 5
+  exposed-modules:     Symengine,
+                       Symengine.DenseMatrix,
+                       Symengine.BasicSym
+                       Symengine.NumberTheory
+
+  other-modules:       Symengine.Internal,
+                       Symengine.VecBasic
+
+  build-depends:       base >= 4.5.0 && <= 5,
+                       singletons,
+                       vector-sized,
+                       finite-typelits
+
   default-language:    Haskell2010
 
 test-suite symengine-test
@@ -28,11 +39,18 @@ test-suite symengine-test
                      , tasty >= 0.10.0 && <= 0.13
                      , tasty-hunit >= 0.9.0 && <= 1.5
                      , tasty-quickcheck >= 0.8.0 && <= 1.5
+                     , vector-sized
+                     , finite-typelits
   ghc-options:         -threaded -rtsopts -with-rtsopts=-N
   include-dirs:        /usr/local/include/
+  extra-lib-dirs:      /usr/local/lib
   extra-libraries:     symengine stdc++ gmpxx gmp
-  
-  other-modules:       Symengine
+   
+  other-modules:       Symengine,
+                       Symengine.BasicSym
+                       Symengine.DenseMatrix
+                       Symengine.Internal
+                       Symengine.VecBasic
 
   default-language:    Haskell2010
 
diff --git a/test/Spec.hs b/test/Spec.hs
index e934667..931e685 100644
--- a/test/Spec.hs
+++ b/test/Spec.hs
@@ -1,18 +1,49 @@
+-- for @
+{-# LANGUAGE TypeApplications #-}
+
+-- for forall r c capturing with Proxy
+{-# LANGUAGE ScopedTypeVariables #-}
+
+-- lift 2, 3, etc to type level
+{-# LANGUAGE DataKinds #-}
+
+
+-- for *, + in type sigs
+{-# LANGUAGE TypeOperators #-}
+
+
+-- for *, + in type sigs
+{-# LANGUAGE FlexibleContexts #-}
+
+-- for (*) which is not injective
+{-# LANGUAGE UndecidableInstances #-}
 import Test.Tasty
 import Test.Tasty.QuickCheck as QC
 import Test.Tasty.HUnit as HU
 
 import Data.List
 import Data.Ord
-import Data.Monoid
+import Data.Ratio
 
 import Symengine as Sym
+import Symengine.DenseMatrix
+import Symengine.VecBasic
+import Symengine.BasicSym
+import Foreign.C.Types
 import Prelude hiding (pi)
 
+
+-- TODO: move arbitrary instance _inside_ the library
+import GHC.TypeLits
+import Data.Proxy
+import qualified Data.Vector.Sized as V
+
 main = defaultMain tests
 
 tests :: TestTree
-tests = testGroup "Tests" [unitTests]
+tests = testGroup "Tests" [genBasic,
+                           symbolIntRing,
+                           denseMatrixRing]
 
 
 -- These are used to check invariants that can be tested by creating
@@ -21,33 +52,178 @@ tests = testGroup "Tests" [unitTests]
 -- properties :: TestTree
 -- properties = testGroup "Properties" [qcProps]
 
-unitTests = testGroup "Unit tests"
-  [ HU.testCase "FFI Sanity Check - ASCII Art should be non-empty" $ 
-    do
-      ascii_art <- Sym.ascii_art_str
-      HU.assertBool "ASCII art from ascii_art_str is empty" (not . null $ ascii_art)
 
+genSafeChar :: Gen Char
+genSafeChar = elements ['a'..'z']
+
+genSafeString :: Gen String
+genSafeString = listOf1 genSafeChar
 
-    , HU.testCase "Basic Constructors" $
+
+instance Arbitrary(BasicSym) where
+  arbitrary = do
+    --intval <- QC.choose (1, 5000) :: Gen (Ratio Integer)
+    let pow2 = 512
+    intval <-  choose (-(2^pow2), 2 ^ pow2 - 1) :: Gen Int
+    strval <- genSafeString :: Gen String
+    choice <- arbitrary :: Gen Bool
+
+    if choice
+    then return (fromIntegral intval)
+    else return (symbol_new (take 10 strval))
+instance forall r c. (KnownNat r, KnownNat c, KnownNat (r * c)) => 
+  Arbitrary(DenseMatrix r c) where
+  arbitrary = do
+    let (rows, cols) = (natVal (Proxy @ r), natVal (Proxy @ c))
+    syms <- V.replicateM arbitrary
+
+    return (densematrix_new_vec syms)
+
+genBasic = testGroup "create and destroy BasicSym"
+      [QC.testProperty "create and die immediately " ((const True) :: BasicSym -> Bool) ]
+
+basicTests = testGroup "Basic tests"
+  [ HU.testCase "ascii art" $
     do
-      "0" @?= (show zero)     
-      "1" @?= (show one)     
+      ascii_art <- ascii_art_str
+      HU.assertBool "ASCII art from ascii_art_str is empty" (not . null $ ascii_art)
+    ,
+    HU.testCase "Basic Constructors" $
+    do
+      "0" @?= (show zero)
+      "1" @?= (show one)
       "-1" @?= (show minus_one)
-    , HU.testCase "Basic Trignometric Functions" $
+    ,
+    HU.testCase "Basic Trignometric Functions" $
     do
       let pi_over_3 = pi / 3 :: BasicSym
       let pi_over_2 = pi / 2 :: BasicSym
 
       sin zero @?= zero
       cos zero @?= one
-      
+
       sin (pi / 6) @?= 1 / 2
       sin (pi / 3) @?= (3 ** (1/2)) / 2
 
       cos (pi / 6) @?= (3 ** (1/2)) / 2
-      cos (pi / 3) @?= 1 / 2 
+      cos (pi / 3) @?= 1 / 2
 
       sin pi_over_2 @?= one
       cos pi_over_2 @?= zero
+   ,
+   HU.testCase "New Symbols, differentiation" $ 
+   do
+      let x = symbol_new "x"
+      let y = symbol_new "y"
+
+      x - x @?= zero
+      x + y @?= y + x
+      diff (x ** 2 + y) x @?= 2 * x
+      diff (x * y) x @?= y
+      diff (sin x) x @?= cos x
+      diff (cos x) x @?= -(sin x)
 
   ]
+-- tests for vectors
+vectorTests = testGroup "Vector"
+    [ HU.testCase "Vector - create, push_back, get out value" $
+      do
+        v <- vecbasic_new
+        vecbasic_push_back v (11 :: BasicSym)
+        vecbasic_push_back v (12 :: BasicSym)
+
+        vecbasic_get v 0 @?= Right (11 :: BasicSym)
+        vecbasic_get v 1 @?= Right (12 :: BasicSym)
+        vecbasic_get v 101 @?= Left RuntimeError
+    ]
+
+-- tests for symbol(ints)
+symbolIntRing = let
+  plus_commutativity :: BasicSym -> BasicSym -> Bool
+  plus_commutativity b1 b2 = b1 + b2 == b2 + b1
+
+  plus_assoc :: BasicSym -> BasicSym -> BasicSym -> Bool
+  plus_assoc b1 b2 b3 = (b1 + b2) + b3 == b1 + (b2 + b3)
+
+  plus_identity :: BasicSym -> Bool
+  plus_identity b = (b + 0) == b && (0 + b) == b
+
+  plus_inverse :: BasicSym -> Bool
+  plus_inverse b = (b + (-b)) == 0 && ((-b) + b) == 0
+
+  mult_identity :: BasicSym -> Bool
+  mult_identity b = (b * 1) == (1 * b) && (b * 1) == b
+
+  mult_assoc :: BasicSym -> BasicSym -> BasicSym -> Bool
+  mult_assoc a b c = (a * b) * c == a * (b * c)
+
+  mult_inverse :: BasicSym -> Bool
+  mult_inverse b = if b == 0 then True else b * (1.0 / b) == 1 && (1.0 / b) * b == 1
+
+  mult_commutativity :: BasicSym -> BasicSym -> Bool
+  mult_commutativity b1 b2 = b1 * b2 == b2 * b1
+
+ -- symengine (==) is structural equality, not "legit" equality. 
+ -- see: https://github.com/symengine/symengine/issues/207
+  mult_distributivity :: BasicSym -> BasicSym -> BasicSym -> Bool
+  mult_distributivity b1 b2 b3 = expand(b1 * (b2 + b3) - (b1 * b2 + b1 * b3)) == (0 :: BasicSym)
+  in
+    testGroup "Symbols of numbers - Ring" [
+      QC.testProperty "(+) identity" plus_identity,
+      QC.testProperty "(+) associativity" plus_assoc,
+      QC.testProperty "(+) inverse" plus_inverse,
+      QC.testProperty "(+) commutativity" plus_commutativity,
+      QC.testProperty "(*) identity" mult_identity,
+      QC.testProperty "(*) associativity" mult_assoc,
+      QC.testProperty "(*) inverse" mult_inverse,
+      QC.testProperty "(*) distributivity" mult_distributivity
+    ]
+
+
+denseMatrixRing = 
+  let
+    eye :: DenseMatrix 10 10
+    eye = densematrix_new_eye @ 0 @ 10 @ 10
+
+    zero :: DenseMatrix 10 10
+    zero = densematrix_new_zeros @ 10 @ 10
+
+    plus_identity :: DenseMatrix 10 10 -> Bool
+    plus_identity d = densematrix_add d zero == d && densematrix_add zero d == d
+
+    plus_invert :: DenseMatrix 10 10 -> Bool
+    plus_invert d = d - d == densematrix_new_zeros
+
+    plus_commutativity :: DenseMatrix 10 10 -> DenseMatrix 10 10 -> Bool
+    plus_commutativity d1 d2 = densematrix_add d1 d2 == densematrix_add d2 d1
+
+    plus_assoc :: DenseMatrix 10 10 -> DenseMatrix 10 10 -> 
+                      DenseMatrix 10 10 -> Bool
+    plus_assoc d1 d2 d3 =
+       densematrix_add (densematrix_add d1 d2) d3 == 
+      densematrix_add d1 (densematrix_add d2 d3)
+
+    mult_identity :: DenseMatrix 10 10 -> Bool
+    mult_identity d = d <> eye == d && eye <> d == d
+
+    mult_assoc :: DenseMatrix 2 2 -> DenseMatrix 2 2  -> DenseMatrix 2 2 -> Bool
+    mult_assoc d1 d2 d3 = (((d1 <> d2) <> d3) - (d1 <> (d2 <> d3))) == densematrix_new_zeros
+
+    mult_nonsingular_invertible :: DenseMatrix 10 10 -> Bool
+    mult_nonsingular_invertible d = if expand(det d) /= 0 then d <> (inv d) == eye  else True
+  in
+    testGroup "DenseMatrix - Ring"
+  [   
+     QC.testProperty "(+) identity" plus_identity,
+     QC.testProperty "(+) associativity" plus_assoc,
+     QC.testProperty "(+) commutativity" plus_commutativity,
+
+     QC.testProperty "(*) identity" mult_identity,
+     -- this fails because I need symbol reduction
+     -- QC.testProperty "(*) associativity " mult_assoc
+ 
+     -- no idea why this fails
+     -- QC.testProperty "(*) non-singluar invertible" mult_nonsingular_invertible
+  ]
+
+