diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index ece1117..7b67558 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -15,13 +15,17 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v2
+      - uses: julia-actions/setup-julia@latest
         with:
           version: '1.10'
       - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.instantiate()'
+        run: |
+          import Pkg
+          Pkg.develop("DACE")
+          Pkg.instantiate()
+        shell: julia --color=yes --project=docs/ {0}
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # If authenticating with GitHub Actions token
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # If authenticating with SSH deploy key
-        run: julia --project=docs/ docs/make.jl
+        run: julia --color=yes --project=docs/ docs/make.jl
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 6bb1d58..9c16181 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -9,7 +9,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: ['1.9', '1.10']
+        julia-version: ['1.10', '1.11']
         julia-arch: [x64]
         os: [ubuntu-latest, macOS-latest, windows-latest]
         exclude:
diff --git a/Project.toml b/Project.toml
index 9f81c5e..cae02ca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,19 +7,25 @@ version = "0.1.0"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
 DACE_jll = "40de70a5-cf80-53d9-bda4-3aa67fea2f4f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+libcxxwrap_julia_jll = "3eaa8342-bff7-56a5-9981-c04077f7cee7"
 
 [compat]
-julia = "1.6"
-CxxWrap = "~0.16"
-DACE_jll = "~0.5"
+CxxWrap = "^0.17"
+DACE_jll = "^0.7.1"
 DiffEqBase = "^6"
-SpecialFunctions = "^2"
 DocStringExtensions = "^0.9"
+SpecialFunctions = "^2"
+julia = ">=1.6"
 
 [extras]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DifferentiableEigen = "73a20539-4e65-4dcb-a56d-dc20f210a01b"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["DelimitedFiles", "DifferentiableEigen", "ForwardDiff", "GenericLinearAlgebra", "LinearAlgebra", "Test"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 26061b8..ad6d37b 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,16 +1,15 @@
 [deps]
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
+DACE = "f4439622-e757-4457-9ca0-f65a61d23f40"
 DACE_jll = "40de70a5-cf80-53d9-bda4-3aa67fea2f4f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-CxxWrap = "~0.16"
-DACE_jll = "~0.5"
 DiffEqBase = "^6"
-SpecialFunctions = "^2"
 DocStringExtensions = "^0.9"
+SpecialFunctions = "^2"
diff --git a/src/DACE.jl b/src/DACE.jl
index 8691bec..fbc0c26 100644
--- a/src/DACE.jl
+++ b/src/DACE.jl
@@ -16,7 +16,7 @@ module DACE
         m_ub::Float64
     end
 
-    # load the C++ interface and initialise it
+    # load the C++ interface and initialize it
     @wrapmodule(() -> libdace, :define_julia_module)
     function __init__()
         @initcxx
@@ -29,6 +29,9 @@ module DACE
     # extend DACE functionality #
     # ========================= #
 
+    # type alias for convenience
+    const DAOrF64 = Union{DA, Float64}
+
     # ------------------- #
     # custom constructors #
     # ------------------- #
@@ -40,15 +43,21 @@ module DACE
     DAAllocated() = DA(0.0)
     DAAllocated(x::Real) = DA(x)
 
-    # AlgebraicVector objects
-    # TODO is there a better way to do this?
-    AlgebraicVector(v::Vector{<:DA}) = AlgebraicVector{DA}(v)
-    AlgebraicVector(v::Vector{Float64}) = AlgebraicVector{Float64}(v)
-    AlgebraicVector{T}(v::Vector{<:T}) where T = begin
+    # AlgebraicVector and AlgebraicMatrix constructors from Julia vectors and matrices
+    AlgebraicVector(v::AbstractVector{<:DA}) = AlgebraicVector{DA}(v)
+    AlgebraicMatrix(m::AbstractMatrix{<:DA}) = AlgebraicMatrix{DA}(m)
+    AlgebraicVector(v::AbstractVector{Float64}) = AlgebraicVector{Float64}(v)
+    AlgebraicMatrix(m::AbstractMatrix{Float64}) = AlgebraicMatrix{Float64}(m)
+
+    AlgebraicVector{T}(v::AbstractVector{U}) where {T<:DAOrF64, U<:DAOrF64} = begin
         res = AlgebraicVector{T}(length(v))
-        for i in eachindex(v)
-            res[i] = v[i]
-        end
+        res .= v
+        return res
+    end
+
+    AlgebraicMatrix{T}(m::AbstractMatrix{U}) where {T<:DAOrF64, U<:DAOrF64} = begin
+        res = AlgebraicMatrix{T}(size(m)...)
+        res .= m
         return res
     end
 
@@ -56,20 +65,22 @@ module DACE
     # overload functions in Base #
     # -------------------------- #
 
-    # addittive and multiplicative identities
+    # additive and multiplicative identities
     Base.zero(::Type{DA}) = DA(0.0)
     Base.one(::Type{DA}) = DA(1.0)
 
-    # similar function for AlgebraicVector
-    function Base.similar(foo::AlgebraicVector{DA})
-        sz = size(foo)[1]
-        bar = AlgebraicVector{DA}(sz)
-        return bar
-    end
+    # similar function for AlgebraicVector and AlgebraicMatrix
+    Base.similar(v::AlgebraicVector{<:DA}) = AlgebraicVector{DA}(v)
+    Base.similar(m::AlgebraicMatrix{<:DA}) = AlgebraicMatrix{DA}(m)
+    Base.similar(v::AlgebraicVector{Float64}) = AlgebraicVector{Float64}(v)
+    Base.similar(m::AlgebraicMatrix{Float64}) = AlgebraicMatrix{Float64}(m)
 
     # promotions and conversions to DA
     @inline Base.promote_rule(::Type{T}, ::Type{R}) where {T<:DA, R<:Real} = T
     @inline Base.promote_rule(::Type{R}, ::Type{T}) where {T<:DA, R<:Real} = T
+    # FIXME: how to handle these type promotions properly?
+    @inline Base.promote_rule(::Type{T}, ::Type{DAAllocated}) where T<:DA = DAAllocated
+    @inline Base.promote_rule(::Type{DAAllocated}, ::Type{T}) where T<:DA = DAAllocated
 
     for R in (AbstractFloat, AbstractIrrational, Integer, Rational)
         @eval Base.convert(::Type{<:DA}, x::$R) = DA(convert(Float64, x))
@@ -88,31 +99,17 @@ module DACE
         @eval Base.:^(a::$R, b::DA) = a^DACE.cons(b)
     end
 
-    # comparison operators (considering only the constant part)
-    for op = (:(==), :(!=), :<, :(<=), :>, :(>=))
-
-        # both arguments are DA objects
-        @eval Base.$op(a::DA, b::DA) = $op(DACE.cons(a), DACE.cons(b))
-
-        # one argument is a number
-        for R in (AbstractFloat, AbstractIrrational, Integer, Rational)
-            @eval Base.$op(a::DA, b::$R) = $op(DACE.cons(a), b)
-            @eval Base.$op(a::$R, b::DA) = $op(a, DACE.cons(b))
-        end
-    end
-
-    Base.isfinite(a::DA) = isfinite(DACE.cons(a))
-    Base.isinf(a::DA) = isinf(DACE.cons(a))
-    Base.isnan(a::DA) = isnan(DACE.cons(a))
+    Base.isfinite(a::DA) = !(isinf(a) || isnan(a))
 
     # assignment of vector and matrix elements
     Base.setindex!(v::AlgebraicVector{<:DA}, x::Real, i::Integer) = v[i] = convert(DA, x)
     Base.setindex!(m::AlgebraicMatrix{<:DA}, x::Real, i::Integer, j::Integer) = m[i,j] = convert(DA, x)
 
-    # custom show functions for DA and AlgebraicVector objects
+    # custom show functions for DA-related objects
+    Base.show(io::IO, m::Monomial) = print(io, toString(m))
     Base.show(io::IO, da::DA) = print(io, toString(da))
-    Base.show(io::IO, vec::AlgebraicVector) = print(io, toString(vec))
-    Base.show(io::IO, mat::AlgebraicMatrix) = print(io, toString(mat))
+    Base.show(io::IO, vec::AlgebraicVector{T}) where T<:DAOrF64 = print(io, toString(vec))
+    Base.show(io::IO, mat::AlgebraicMatrix{T}) where T<:DAOrF64 = print(io, toString(mat))
 
     # -------------------------------------- #
     # overload functions in SpecialFunctions #
@@ -133,26 +130,29 @@ module DACE
     # ---------------------------------------------- #
 
     # map inversion
-    invert(v::Vector{<:DA}) = Vector{DA}(invert(AlgebraicVector(v)))
+    invert(v::AbstractVector{<:DA}) = invert(AlgebraicVector(v))
 
     # linear part, Jacobian and Hessian
-    linear(v::Vector{<:DA}) = linear(AlgebraicVector(v))
-    jacobian(v::Vector{<:DA}) = jacobian(AlgebraicVector(v))
-    hessian(v::Vector{<:DA}) = stack(hessian(AlgebraicVector(v)), dims=3)
+    linear(v::AbstractVector{<:DA}) = linear(AlgebraicVector(v))
+    jacobian(v::AbstractVector{<:DA}) = jacobian(AlgebraicVector(v))
+    hessian(v::AbstractVector{<:DA}) = hessian(AlgebraicVector(v)) # vector of Hessian matrices
+    hess_stack(a::AbstractVector{<:DA}) = stack(hessian(a), dims=3)
 
     # compilation and evaluation of DA objects
-    compile(v::Vector{<:DA}) = compile(StdVector{DA}(v))
-    for R in (DA, Float64)
-        @eval eval(cda::compiledDA, v::Vector{<:$R}) = Vector{$R}(eval(cda, AlgebraicVector(v)))
-        @eval eval(da::DA, v::Vector{<:$R}) = eval(da, AlgebraicVector(v))
-        @eval eval(a::Vector{<:DA}, v::Vector{<:$R}) = Vector{$R}(eval(AlgebraicVector(a), AlgebraicVector(v)))
-        @eval eval(a::AlgebraicVector{DA}, v::Vector{<:$R}) = Vector{$R}(eval(a, AlgebraicVector(v)))
-    end
+    compile(v::AbstractVector{<:DA}) = compile(AlgebraicVector(v))
+
+    eval(cda::compiledDA, v::AbstractVector{<:DA}) = eval(cda, AlgebraicVector(v))
+    eval(a::AbstractVector{<:DA}, v::AbstractVector{<:DA}) = eval(AlgebraicVector(a), AlgebraicVector(v))
+    eval(a::AlgebraicVector{<:DA}, v::AbstractVector{<:DA}) = eval(a, AlgebraicVector(v))
+
+    eval(cda::compiledDA, v::AbstractVector{Float64}) = eval(cda, AlgebraicVector(v))
+    eval(a::AbstractVector{<:DA}, v::AbstractVector{Float64}) = eval(AlgebraicVector(a), AlgebraicVector(v))
+    eval(a::AlgebraicVector{<:DA}, v::AbstractVector{Float64}) = eval(a, AlgebraicVector(v))
 
     # ------- #
     # exports #
     # ------- #
 
-    export DA, AlgebraicVector, AlgebraicMatrix, compiledDA
+    export DA, AlgebraicVector, AlgebraicMatrix, compiledDA, Monomial
 
 end  # module DACE
diff --git a/test/comparison_operators.jl b/test/comparison_operators.jl
index 2d6f454..aab0eb8 100644
--- a/test/comparison_operators.jl
+++ b/test/comparison_operators.jl
@@ -13,7 +13,7 @@
 
     @testset "Test >=" begin
         DACE.init(10, 4)
-        
+
         x = 1 + DA(1, 1)
         y = DA(2)
         @test y >= x
@@ -22,4 +22,36 @@
 
     end
 
+    @testset "Test isnan" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isnan(x) == false
+        @test isnan(y) == false
+        @test isnan(z) == true
+    end
+
+    @testset "Test isinf" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isinf(x) == false
+        @test isinf(y) == true
+        @test isinf(z) == false
+    end
+
+    @testset "Test isfinite" begin
+        DACE.init(2, 2)
+
+        x = 1.0 + DA(1, 1.0)
+        y = 1.0 + DA(2, Inf)
+        z = 1.0 + DA(2, NaN)
+        @test isfinite(x) == true
+        @test isfinite(y) == false
+        @test isfinite(z) == false
+    end
 end
diff --git a/test/data/double_moments_nonzero_6D4.txt b/test/data/double_moments_nonzero_6D4.txt
new file mode 100644
index 0000000..c989a3b
--- /dev/null
+++ b/test/data/double_moments_nonzero_6D4.txt
@@ -0,0 +1,215 @@
+6 4 0 0 0
+3.861346813186461002e+00 4.564531354848048311e-02 -2.003168404365469055e-04 1.568048483265416021e-03 -1.008404804051421001e-03
+-1.487123146365427973e+00 -2.003168404365469055e-04 5.204846143436810957e-01 1.003843436114855001e-05 -1.008403740486084072e-03
+-2.364168486116942969e+00 1.568048483265416021e-03 1.003843436114855001e-05 3.257040860848609692e-02 2.108634834860479981e-06
+1.061040254384064907e+00 -1.008404804051421001e-03 -1.008403740486084072e-03 2.108634834860479981e-06 2.903406840638488862e-01
+0 0 0 0 1.000000000000000000e+00
+1 0 0 0 3.861346813186461002e+00
+0 1 0 0 -1.487123146365427973e+00
+0 0 1 0 -2.364168486116942969e+00
+0 0 0 1 1.061040254384064907e+00
+2 0 0 0 1.495564452525371735e+01
+1 1 0 0 -5.742498538874404979e+00
+0 2 0 0 2.732019866799491048e+00
+1 0 1 0 -9.127306401220252496e+00
+0 1 1 0 3.515819716046580279e+00
+0 0 2 0 5.621863039356963299e+00
+1 0 0 1 4.096036000124408716e+00
+0 1 0 1 -1.578905925260490895e+00
+0 0 1 1 -2.508475823281475670e+00
+0 0 0 2 1.416147105487250091e+00
+3 0 0 0 5.810143509875263135e+01
+2 1 0 0 -2.224243212790509361e+01
+1 2 0 0 1.054987199784835639e+01
+0 3 0 0 -5.610899414882860370e+00
+2 0 1 0 -3.534555391814467384e+01
+1 1 1 0 1.357394095786189858e+01
+0 2 1 0 -6.458985129308948636e+00
+1 0 2 0 2.170054866957408279e+01
+0 1 2 0 -8.360450116624367212e+00
+0 0 3 0 -1.344503529813761666e+01
+2 0 0 1 1.586075327019980818e+01
+1 1 0 1 -6.095416284932143824e+00
+0 2 0 1 2.901782295537967737e+00
+1 0 1 1 -9.682047324762656615e+00
+0 1 1 1 3.732807146427444156e+00
+0 0 2 1 5.965013019055636079e+00
+1 0 0 2 5.466095196596799433e+00
+0 1 0 2 -2.108125053151148176e+00
+0 0 1 2 -3.348005883805800220e+00
+0 0 0 3 2.118715391605747733e+00
+4 0 0 0 2.263977565111529202e+02
+3 1 0 0 -8.641297657477772987e+01
+2 2 0 0 4.086371916018460837e+01
+1 3 0 0 -2.166727038353090151e+01
+0 4 0 0 1.261002131205176369e+01
+3 0 1 0 -1.372912283314975923e+02
+2 1 1 0 5.256699819039605615e+01
+1 2 1 0 -2.493750626166159989e+01
+0 3 1 0 1.326519385104437809e+01
+2 0 2 0 8.402133203289817232e+01
+1 1 2 0 -3.227269761692694061e+01
+0 2 2 0 1.535918268504924278e+01
+1 0 3 0 -5.188949814020720197e+01
+0 1 3 0 1.999459249966997731e+01
+0 0 4 0 3.233564787558528053e+01
+3 0 0 1 6.160271744589882559e+01
+2 1 0 1 -2.360361564476562179e+01
+1 2 0 1 1.120266540037057545e+01
+0 3 0 1 -5.961655079649003319e+00
+2 0 1 1 -3.748461594541477382e+01
+1 1 1 1 1.440814429851531209e+01
+0 2 1 1 -6.860328195337675794e+00
+1 0 2 1 2.301944807462959020e+01
+0 1 2 1 -8.876428399065940056e+00
+0 0 3 1 -1.426568810956993083e+01
+2 0 0 2 2.116286882185341867e+01
+1 1 0 2 -8.137301281257123264e+00
+0 2 0 2 3.875308691309187648e+00
+1 0 1 2 -1.292053214450011112e+01
+0 1 1 2 4.983970372681242544e+00
+0 0 2 2 7.961363912784209518e+00
+1 0 0 3 8.176810776792713398e+00
+0 1 0 3 -3.155074843532354123e+00
+0 0 1 3 -5.008991201473874888e+00
+0 0 0 4 3.481537676103419798e+00
+5 0 0 0 8.848084885077015542e+02
+4 1 0 0 -3.367278987765552642e+02
+3 2 0 0 1.587610072520565154e+02
+2 3 0 0 -8.392729666005463685e+01
+1 4 0 0 4.869616143807560604e+01
+0 5 0 0 -3.043420184161947262e+01
+4 0 1 0 -5.348780178025162968e+02
+3 1 1 0 2.041907876212068231e+02
+2 2 1 0 -9.657607820086181505e+01
+1 3 1 0 5.121659737119421152e+01
+0 4 1 0 -2.981244029379772087e+01
+3 0 2 0 3.263057122968832005e+02
+2 1 2 0 -1.249594712533143763e+02
+1 2 2 0 5.929022458829660991e+01
+0 3 2 0 -3.154409706577002126e+01
+2 0 3 0 -2.008749685961057310e+02
+1 1 3 0 7.716942052317870093e+01
+0 2 3 0 -3.673260709788132061e+01
+1 0 4 0 1.247748210078840572e+02
+0 1 4 0 -4.808763027692226899e+01
+0 0 5 0 -7.819856085929647804e+01
+4 0 0 1 2.399827740954719673e+02
+3 1 0 1 -9.167894822523282983e+01
+2 2 0 1 4.338163259310227460e+01
+1 3 0 1 -2.301610361277774786e+01
+0 4 0 1 1.340237242855771882e+01
+3 0 1 1 -1.455644694397240357e+02
+2 1 1 1 5.578392076521424769e+01
+1 2 1 1 -2.648053849578640850e+01
+0 3 1 1 1.409443262112051976e+01
+2 0 2 1 8.910610057708380793e+01
+1 1 2 1 -3.425602624054199197e+01
+0 2 2 1 1.631354528233337930e+01
+1 0 3 1 -5.504315099276296763e+01
+0 1 3 1 2.122857264863208826e+01
+0 0 4 1 3.430931064490537352e+01
+3 0 0 2 8.218419122182882575e+01
+2 1 0 2 -3.150597005616125657e+01
+1 2 0 2 1.495890310878606400e+01
+0 3 0 2 -7.963406900954058010e+00
+2 0 1 2 -5.001537845081414702e+01
+1 1 1 2 1.923467477391860925e+01
+0 2 1 2 -9.161912768919760808e+00
+1 0 2 2 3.071905720658160277e+01
+0 1 2 2 -1.185162605169537287e+01
+0 0 3 2 -1.904007235236298001e+01
+2 0 0 3 3.165367555340006689e+01
+1 1 0 3 -1.217688507651806873e+01
+0 2 0 3 5.801121115403214290e+00
+1 0 1 3 -1.932800152898223089e+01
+0 1 1 3 7.459136449207157682e+00
+0 0 2 3 1.191108539262575761e+01
+1 0 0 4 1.343487831949328815e+01
+0 1 0 4 -5.186021345180377828e+00
+0 0 1 4 -8.230923786684197552e+00
+0 0 0 5 6.154648726042156071e+00
+6 0 0 0 3.468222420292661809e+03
+5 1 0 0 -1.316045939776740170e+03
+4 2 0 0 6.186616411813527066e+02
+3 3 0 0 -3.260749752717806018e+02
+2 4 0 0 1.886257174365447611e+02
+1 5 0 0 -1.175296382911461706e+02
+0 6 0 0 7.807601639717404396e+01
+5 0 1 0 -2.090061331485109349e+03
+4 1 1 0 7.955417604155155686e+02
+3 2 1 0 -3.751472761929320541e+02
+2 3 1 0 1.983515498510538180e+02
+1 4 1 0 -1.151070271638552640e+02
+0 5 1 0 7.195221381843514052e+01
+4 0 2 0 1.271054503833614035e+03
+3 1 2 0 -4.853100265892080074e+02
+2 2 2 0 2.295759175685338960e+02
+1 3 2 0 -1.217703278544002359e+02
+0 4 2 0 7.089287803064581794e+01
+3 0 3 0 -7.799896960956992871e+02
+2 1 3 0 2.987491343332151814e+02
+1 2 3 0 -1.417730940352845437e+02
+0 3 3 0 7.544012832034754013e+01
+2 0 4 0 4.829493672553531951e+02
+1 1 4 0 -1.855640853365281373e+02
+0 2 4 0 8.834323811036863106e+01
+1 0 5 0 -3.016982444517786917e+02
+0 1 5 0 1.162925128627191071e+02
+0 0 6 0 1.901404995628859922e+02
+5 0 0 1 9.376759208009291342e+02
+4 1 0 1 -3.571615986777495664e+02
+3 2 0 1 1.685024763464636806e+02
+2 3 0 1 -8.913016282061434481e+01
+1 4 0 1 5.174326893971637276e+01
+0 5 0 1 -3.235549322730126676e+01
+4 0 1 1 -5.669728501466158832e+02
+3 1 1 1 2.166338815747311060e+02
+2 2 1 1 -1.025267437279678830e+02
+1 3 1 1 5.440489035362881509e+01
+0 4 1 1 -3.168567932760130645e+01
+3 0 2 1 3.459687344437305683e+02
+2 1 2 1 -1.326064470159769826e+02
+1 2 2 1 6.295880935584480653e+01
+0 3 2 1 -3.351596560389176460e+01
+2 0 3 1 -2.130312450292713891e+02
+1 1 3 1 8.191182034946592694e+01
+0 2 3 1 -3.901500286233520143e+01
+1 0 4 1 1.323580627362582049e+02
+0 1 4 1 -5.105535020484413877e+01
+0 0 5 1 -8.297147998624799925e+01
+4 0 0 2 3.201153812651620001e+02
+3 1 0 2 -1.223549716897481119e+02
+2 2 0 2 5.791906896193700049e+01
+1 3 0 2 -3.073978120463948827e+01
+0 4 0 2 1.790570588587489453e+01
+3 0 1 2 -1.941974619353107414e+02
+2 1 1 2 7.446001506471382925e+01
+1 2 1 2 -3.535940677162169266e+01
+0 3 1 2 1.882692720155041854e+01
+2 0 2 2 1.188933867150845316e+02
+1 1 2 2 -4.573130100060818393e+01
+0 2 2 2 2.178659695905734850e+01
+1 0 3 2 -7.345410010996637595e+01
+0 1 3 2 2.834394320905924403e+01
+0 0 4 2 4.579179349374999930e+01
+3 0 0 3 1.229082631914641297e+02
+2 1 0 3 -4.714021163704033057e+01
+1 2 0 3 2.238968092141435662e+01
+0 3 0 3 -1.192304093924768793e+01
+2 0 1 3 -7.480884506735876016e+01
+1 1 1 3 2.878329125323609716e+01
+0 2 1 3 -1.371486655437458424e+01
+1 0 2 3 4.595303810744565709e+01
+0 1 2 3 -1.773743615801475570e+01
+0 0 3 3 -2.848605213816248494e+01
+2 0 0 4 5.200265832230332563e+01
+1 1 0 4 -2.001299803443895797e+01
+0 2 0 4 9.537065531727536793e+00
+1 0 1 4 -3.175678775025720668e+01
+0 1 1 4 1.226063657018980102e+01
+0 0 2 4 1.957264348427232292e+01
+1 0 0 5 2.374767924799464680e+01
+0 1 0 5 -9.170274556321912129e+00
+0 0 1 5 -1.455058985477054279e+01
+0 0 0 6 1.158449020229396353e+01
diff --git a/test/data/double_moments_zero_6D4.txt b/test/data/double_moments_zero_6D4.txt
new file mode 100644
index 0000000..4867f98
--- /dev/null
+++ b/test/data/double_moments_zero_6D4.txt
@@ -0,0 +1,215 @@
+6 4 0 0 0
+0.000000000000000000e+00 4.564531354848048311e-02 -2.003168404365469055e-04 1.568048483265416021e-03 -1.008404804051421001e-03
+0.000000000000000000e+00 -2.003168404365469055e-04 5.204846143436810957e-01 1.003843436114855001e-05 -1.008403740486084072e-03
+0.000000000000000000e+00 1.568048483265416021e-03 1.003843436114855001e-05 3.257040860848609692e-02 2.108634834860479981e-06
+0.000000000000000000e+00 -1.008404804051421001e-03 -1.008403740486084072e-03 2.108634834860479981e-06 2.903406840638488862e-01
+0 0 0 0 1.000000000000000000e+00
+1 0 0 0 0.000000000000000000e+00
+0 1 0 0 0.000000000000000000e+00
+0 0 1 0 0.000000000000000000e+00
+0 0 0 1 0.000000000000000000e+00
+2 0 0 0 4.564531354848048311e-02
+1 1 0 0 -2.003168404365318622e-04
+0 2 0 0 5.204846143436810957e-01
+1 0 1 0 1.568048483265416237e-03
+0 1 1 0 1.003843436117590748e-05
+0 0 2 0 3.257040860848609692e-02
+1 0 0 1 -1.008404804051413195e-03
+0 1 0 1 -1.008403740486080169e-03
+0 0 1 1 2.108634834852018969e-06
+0 0 0 2 2.903406840638488862e-01
+3 0 0 0 0.000000000000000000e+00
+2 1 0 0 0.000000000000000000e+00
+1 2 0 0 0.000000000000000000e+00
+0 3 0 0 0.000000000000000000e+00
+2 0 1 0 0.000000000000000000e+00
+1 1 1 0 0.000000000000000000e+00
+0 2 1 0 0.000000000000000000e+00
+1 0 2 0 0.000000000000000000e+00
+0 1 2 0 0.000000000000000000e+00
+0 0 3 0 0.000000000000000000e+00
+2 0 0 1 0.000000000000000000e+00
+1 1 0 1 0.000000000000000000e+00
+0 2 0 1 0.000000000000000000e+00
+1 0 1 1 0.000000000000000000e+00
+0 1 1 1 0.000000000000000000e+00
+0 0 2 1 0.000000000000000000e+00
+1 0 0 2 0.000000000000000000e+00
+0 1 0 2 0.000000000000000000e+00
+0 0 1 2 0.000000000000000000e+00
+0 0 0 3 0.000000000000000000e+00
+4 0 0 0 6.250483946817288021e-03
+3 1 0 0 -2.743057497230430575e-05
+2 2 0 0 2.375776367255038810e-02
+1 3 0 0 -3.127855003233204201e-04
+0 4 0 0 8.127127013054713522e-01
+3 0 1 0 2.147221940336076305e-04
+2 1 1 0 -1.700055516883069995e-07
+1 2 1 0 8.161410883496880708e-04
+0 3 1 0 1.567455191114897062e-05
+2 0 2 0 1.491604065428217638e-03
+1 1 2 0 -6.492919840630616973e-06
+0 2 2 0 1.695239676514431401e-02
+1 0 3 0 1.532159394536145676e-04
+0 1 3 0 9.808677267925808296e-07
+0 0 4 0 3.182494550771236032e-03
+3 0 0 1 -1.380868603941660076e-04
+2 1 0 1 -4.562490398948848558e-05
+1 2 0 1 -5.244551850366249712e-04
+0 3 0 1 -1.574575895908758572e-03
+2 0 1 1 -3.066205948819983318e-06
+1 1 1 1 -1.591771156288571132e-06
+0 2 1 1 1.077266399342191594e-06
+1 0 2 1 -3.283754362740694960e-05
+0 1 2 1 -3.284407953517247403e-05
+0 0 3 1 2.060372945320368410e-07
+2 0 0 2 1.325472532047236238e-02
+1 1 0 2 -5.612637012920515533e-05
+0 2 0 2 1.511198927294606087e-01
+1 0 1 2 4.552640165615701867e-04
+0 1 1 2 2.910313188848184796e-06
+0 0 2 2 9.456514724519604967e-03
+1 0 0 3 -8.783428218646985308e-04
+0 1 0 3 -8.783418954758787089e-04
+0 0 1 3 1.836667441176831339e-06
+0 0 0 4 2.528931384679911076e-01
+5 0 0 0 0.000000000000000000e+00
+4 1 0 0 0.000000000000000000e+00
+3 2 0 0 0.000000000000000000e+00
+2 3 0 0 0.000000000000000000e+00
+1 4 0 0 0.000000000000000000e+00
+0 5 0 0 0.000000000000000000e+00
+4 0 1 0 0.000000000000000000e+00
+3 1 1 0 0.000000000000000000e+00
+2 2 1 0 0.000000000000000000e+00
+1 3 1 0 0.000000000000000000e+00
+0 4 1 0 0.000000000000000000e+00
+3 0 2 0 0.000000000000000000e+00
+2 1 2 0 0.000000000000000000e+00
+1 2 2 0 0.000000000000000000e+00
+0 3 2 0 0.000000000000000000e+00
+2 0 3 0 0.000000000000000000e+00
+1 1 3 0 0.000000000000000000e+00
+0 2 3 0 0.000000000000000000e+00
+1 0 4 0 0.000000000000000000e+00
+0 1 4 0 0.000000000000000000e+00
+0 0 5 0 0.000000000000000000e+00
+4 0 0 1 0.000000000000000000e+00
+3 1 0 1 0.000000000000000000e+00
+2 2 0 1 0.000000000000000000e+00
+1 3 0 1 0.000000000000000000e+00
+0 4 0 1 0.000000000000000000e+00
+3 0 1 1 0.000000000000000000e+00
+2 1 1 1 0.000000000000000000e+00
+1 2 1 1 0.000000000000000000e+00
+0 3 1 1 0.000000000000000000e+00
+2 0 2 1 0.000000000000000000e+00
+1 1 2 1 0.000000000000000000e+00
+0 2 2 1 0.000000000000000000e+00
+1 0 3 1 0.000000000000000000e+00
+0 1 3 1 0.000000000000000000e+00
+0 0 4 1 0.000000000000000000e+00
+3 0 0 2 0.000000000000000000e+00
+2 1 0 2 0.000000000000000000e+00
+1 2 0 2 0.000000000000000000e+00
+0 3 0 2 0.000000000000000000e+00
+2 0 1 2 0.000000000000000000e+00
+1 1 1 2 0.000000000000000000e+00
+0 2 1 2 0.000000000000000000e+00
+1 0 2 2 0.000000000000000000e+00
+0 1 2 2 0.000000000000000000e+00
+0 0 3 2 0.000000000000000000e+00
+2 0 0 3 0.000000000000000000e+00
+1 1 0 3 0.000000000000000000e+00
+0 2 0 3 0.000000000000000000e+00
+1 0 1 3 0.000000000000000000e+00
+0 1 1 3 0.000000000000000000e+00
+0 0 2 3 0.000000000000000000e+00
+1 0 0 4 0.000000000000000000e+00
+0 1 0 4 0.000000000000000000e+00
+0 0 1 4 0.000000000000000000e+00
+0 0 0 5 0.000000000000000000e+00
+6 0 0 0 1.426526497911094525e-03
+5 1 0 0 -6.260385977129067150e-06
+4 2 0 0 3.253302705745017565e-03
+3 3 0 0 -4.283162493568665186e-05
+2 4 0 0 3.709677670073341371e-02
+1 5 0 0 -8.140002025409653896e-04
+0 6 0 0 2.115022284555948318e+00
+5 0 1 0 4.900530936240848900e-05
+4 1 1 0 -1.093048130977039220e-07
+3 2 1 0 1.117594251575015560e-04
+2 3 1 0 -2.654534051599322098e-07
+1 4 1 0 1.274360359105469742e-03
+0 5 1 0 4.079181553295896379e-05
+4 0 2 0 2.049275953913331713e-04
+3 1 2 0 -8.920692914056796430e-07
+2 2 2 0 7.763595646170691478e-04
+1 3 2 0 -1.013839475891981090e-05
+0 4 2 0 2.647038539221760467e-02
+3 0 3 0 2.100390167125216383e-05
+2 1 3 0 -1.646335728118161756e-08
+1 2 3 0 7.974614713724542669e-05
+0 3 3 0 1.531579687706619008e-06
+2 0 4 0 1.462269617221728787e-04
+1 1 4 0 -6.313550605117684497e-07
+0 2 4 0 1.656439488294534149e-03
+1 0 5 0 2.495152876668646795e-05
+0 1 5 0 1.597363132629751881e-07
+0 0 6 0 5.182757395644976492e-04
+5 0 0 1 -3.151509019808367441e-05
+4 1 0 1 -6.192366897499277517e-06
+3 2 0 1 -7.181700687397020444e-05
+2 3 0 1 -7.124142445665698764e-05
+1 4 0 1 -8.182817360357919068e-04
+0 5 0 1 -4.097712639684392144e-03
+4 0 1 1 -8.529275798153370877e-07
+3 1 1 1 -2.160702014494526917e-07
+2 2 1 1 -1.595561872002275120e-06
+1 3 1 1 -2.485465022616492313e-06
+0 4 1 1 1.650489206000971265e-06
+3 0 2 1 -4.511516574564927891e-06
+2 1 2 1 -1.491044852742074409e-06
+1 2 2 1 -1.707834148179174843e-05
+0 3 2 1 -5.128451482083046242e-05
+2 0 3 1 -2.995716339280910055e-07
+1 1 3 1 -1.555337117656591592e-07
+0 2 3 1 1.052610216802085663e-07
+1 0 4 1 -3.207950487996169293e-06
+0 1 4 1 -3.209231135906315605e-06
+0 0 5 1 3.355359435783245026e-08
+4 0 0 2 1.815326774662629969e-03
+3 1 0 2 -7.686939478087576744e-06
+2 2 0 2 6.898995099419431669e-03
+1 3 0 2 -8.763995852045791679e-05
+0 4 0 2 2.359699128773015409e-01
+3 0 1 2 6.235157345507074192e-05
+2 1 1 2 -4.315346152320175036e-08
+1 2 1 2 2.369599799749549636e-04
+0 3 1 2 4.544380960369576172e-06
+2 0 2 2 4.331395588513227771e-04
+1 1 2 2 -1.818931871618711994e-06
+0 2 2 2 4.922036718042210185e-03
+1 0 3 2 4.448440517428498171e-05
+0 1 3 2 2.843702694823105337e-07
+0 0 4 2 9.240076466382207438e-04
+3 0 0 3 -1.202828530827790982e-04
+2 1 0 3 -3.974645005902499584e-05
+1 2 0 3 -4.568181841017177058e-04
+0 3 0 3 -1.371496480727846246e-03
+2 0 1 3 -2.670720132615621301e-06
+1 1 1 3 -1.386454913809398438e-06
+0 2 1 3 9.383356551692645375e-07
+1 0 2 3 -2.860222466616237774e-05
+0 1 2 3 -2.860791758595160328e-05
+0 0 3 3 1.794630271698821611e-07
+2 0 0 4 1.154692950011542568e-02
+1 1 0 4 -4.711585771801501954e-05
+0 2 0 4 1.316305305386868130e-01
+1 0 1 4 3.965412937858708853e-04
+0 1 1 4 2.531242761609486802e-06
+0 0 2 4 8.236832869676392888e-03
+1 0 0 5 -1.275093278713779199e-03
+0 1 0 5 -1.275091933871791428e-03
+0 0 1 5 2.666296406356183368e-06
+0 0 0 6 3.671258340892509797e-01
diff --git a/test/data/integer_moments_nonzero_6D4.txt b/test/data/integer_moments_nonzero_6D4.txt
new file mode 100644
index 0000000..d2e5133
--- /dev/null
+++ b/test/data/integer_moments_nonzero_6D4.txt
@@ -0,0 +1,215 @@
+6 4 0 0 0
+3.0 4.0 2.0 1.0 1.0
+2.0 2.0 3.0 1.0 2.0
+-1.0 1.0 1.0 2.0 1.0
+2.0 1.0 2.0 1.0 3.0
+0 0 0 0 1.0
+1 0 0 0 3.0
+0 1 0 0 2.0
+0 0 1 0 -1.0
+0 0 0 1 2.0
+2 0 0 0 13.0
+1 1 0 0 8.0
+0 2 0 0 7.0
+1 0 1 0 -2.0
+0 1 1 0 -1.0
+0 0 2 0 3.0
+1 0 0 1 7.0
+0 1 0 1 6.0
+0 0 1 1 -1.0
+0 0 0 2 7.0
+3 0 0 0 63.0
+2 1 0 0 38.0
+1 2 0 0 29.0
+0 3 0 0 26.0
+2 0 1 0 -7.0
+1 1 1 0 -3.0
+0 2 1 0 -3.0
+1 0 2 0 7.0
+0 1 2 0 4.0
+0 0 3 0 -7.0
+2 0 0 1 32.0
+1 1 0 1 24.0
+0 2 0 1 22.0
+1 0 1 1 -2.0
+0 1 1 1 -2.0
+0 0 2 1 4.0
+1 0 0 2 25.0
+0 1 0 2 22.0
+0 0 1 2 -3.0
+0 0 0 3 26.0
+4 0 0 0 345.0
+3 1 0 0 204.0
+2 2 0 0 147.0
+1 3 0 0 120.0
+0 4 0 0 115.0
+3 0 1 0 -24.0
+2 1 1 0 -9.0
+1 2 1 0 -6.0
+0 3 1 0 -5.0
+2 0 2 0 29.0
+1 1 2 0 16.0
+0 2 2 0 15.0
+1 0 3 0 -12.0
+0 1 3 0 -5.0
+0 0 4 0 25.0
+3 0 0 1 165.0
+2 1 0 1 118.0
+1 2 0 1 97.0
+0 3 0 1 94.0
+2 0 1 1 -5.0
+1 1 1 1 -3.0
+0 2 1 1 -3.0
+1 0 2 1 13.0
+0 1 2 1 12.0
+0 0 3 1 -5.0
+2 0 0 2 117.0
+1 1 0 2 92.0
+0 2 0 2 89.0
+1 0 1 2 -4.0
+0 1 1 2 -3.0
+0 0 2 2 15.0
+1 0 0 3 99.0
+0 1 0 3 94.0
+0 0 1 3 -5.0
+0 0 0 4 115.0
+5 0 0 0 2043.0
+4 1 0 0 1194.0
+3 2 0 0 825.0
+2 3 0 0 638.0
+1 4 0 0 553.0
+0 5 0 0 542.0
+4 0 1 0 -93.0
+3 1 1 0 -27.0
+2 2 1 0 -13.0
+1 3 1 0 -7.0
+0 4 1 0 -11.0
+3 0 2 0 129.0
+2 1 2 0 72.0
+1 2 2 0 55.0
+0 3 2 0 48.0
+2 0 3 0 -43.0
+1 1 3 0 -17.0
+0 2 3 0 -19.0
+1 0 4 0 47.0
+0 1 4 0 22.0
+0 0 5 0 -81.0
+4 0 0 1 942.0
+3 1 0 1 648.0
+2 2 0 1 504.0
+1 3 0 1 440.0
+0 4 0 1 438.0
+3 0 1 1 -6.0
+2 1 1 1 -0.0
+1 2 1 1 2.0
+0 3 1 1 -2.0
+2 0 2 1 58.0
+1 1 2 1 44.0
+0 2 2 1 40.0
+1 0 3 1 -10.0
+0 1 3 1 -12.0
+0 0 4 1 22.0
+3 0 0 2 615.0
+2 1 0 2 462.0
+1 2 0 2 399.0
+0 3 0 2 398.0
+2 0 1 2 -3.0
+1 1 1 2 3.0
+0 2 1 2 -1.0
+1 0 2 2 47.0
+0 1 2 2 40.0
+0 0 3 2 -19.0
+2 0 0 3 476.0
+1 1 0 3 400.0
+0 2 0 3 398.0
+1 0 1 3 2.0
+0 1 1 3 -2.0
+0 0 2 3 48.0
+1 0 0 4 449.0
+0 1 0 4 438.0
+0 0 1 4 -11.0
+0 0 0 5 542.0
+6 0 0 0 13029.0
+5 1 0 0 7536.0
+4 2 0 0 5055.0
+3 3 0 0 3756.0
+2 4 0 0 3079.0
+1 5 0 0 2776.0
+0 6 0 0 2809.0
+5 0 1 0 -318.0
+4 1 1 0 -33.0
+3 2 1 0 24.0
+2 3 1 0 43.0
+1 4 1 0 42.0
+0 5 1 0 33.0
+4 0 2 0 687.0
+3 1 2 0 384.0
+2 2 2 0 277.0
+1 3 2 0 224.0
+0 4 2 0 221.0
+3 0 3 0 -138.0
+2 1 3 0 -47.0
+1 2 3 0 -32.0
+0 3 3 0 -23.0
+2 0 4 0 193.0
+1 1 4 0 96.0
+0 2 4 0 99.0
+1 0 5 0 -118.0
+0 1 5 0 -37.0
+0 0 6 0 331.0
+5 0 0 1 5811.0
+4 1 0 1 3894.0
+3 2 0 1 2907.0
+2 3 0 1 2398.0
+1 4 0 1 2181.0
+0 5 0 1 2234.0
+4 0 1 1 63.0
+3 1 1 1 75.0
+2 2 1 1 73.0
+1 3 1 1 65.0
+0 4 1 1 53.0
+3 0 2 1 297.0
+2 1 2 1 216.0
+1 2 2 1 177.0
+0 3 2 1 176.0
+2 0 3 1 -23.0
+1 1 3 1 -15.0
+0 2 3 1 -13.0
+1 0 4 1 71.0
+0 1 4 1 74.0
+0 0 5 1 -37.0
+4 0 0 2 3579.0
+3 1 0 2 2592.0
+2 2 0 2 2115.0
+1 3 0 2 1916.0
+0 4 0 2 1973.0
+3 0 1 2 66.0
+2 1 1 2 75.0
+1 2 1 2 68.0
+0 3 1 2 57.0
+2 0 2 2 219.0
+1 1 2 2 168.0
+0 2 2 2 167.0
+1 0 3 2 -22.0
+0 1 3 2 -13.0
+0 0 4 2 99.0
+3 0 0 3 2571.0
+2 1 0 3 2050.0
+1 2 0 3 1837.0
+0 3 0 3 1894.0
+2 0 1 3 73.0
+1 1 1 3 69.0
+0 2 1 3 57.0
+1 0 2 3 179.0
+0 1 2 3 176.0
+0 0 3 3 -23.0
+2 0 0 4 2203.0
+1 1 0 4 1920.0
+0 2 0 4 1973.0
+1 0 1 4 62.0
+0 1 1 4 53.0
+0 0 2 4 221.0
+1 0 0 5 2201.0
+0 1 0 5 2234.0
+0 0 1 5 33.0
+0 0 0 6 2809.0
diff --git a/test/data/integer_moments_zero_6D4.txt b/test/data/integer_moments_zero_6D4.txt
new file mode 100644
index 0000000..66c029f
--- /dev/null
+++ b/test/data/integer_moments_zero_6D4.txt
@@ -0,0 +1,215 @@
+6 4 0 0 0
+0.0 4.0 2.0 1.0 1.0
+0.0 2.0 3.0 1.0 2.0
+0.0 1.0 1.0 2.0 1.0
+0.0 1.0 2.0 1.0 3.0
+0 0 0 0 1.0
+1 0 0 0 0.0
+0 1 0 0 0.0
+0 0 1 0 0.0
+0 0 0 1 0.0
+2 0 0 0 4.0
+1 1 0 0 2.0
+0 2 0 0 3.0
+1 0 1 0 1.0
+0 1 1 0 1.0
+0 0 2 0 2.0
+1 0 0 1 1.0
+0 1 0 1 2.0
+0 0 1 1 1.0
+0 0 0 2 3.0
+3 0 0 0 0.0
+2 1 0 0 0.0
+1 2 0 0 0.0
+0 3 0 0 0.0
+2 0 1 0 0.0
+1 1 1 0 0.0
+0 2 1 0 0.0
+1 0 2 0 0.0
+0 1 2 0 0.0
+0 0 3 0 0.0
+2 0 0 1 0.0
+1 1 0 1 0.0
+0 2 0 1 0.0
+1 0 1 1 0.0
+0 1 1 1 0.0
+0 0 2 1 0.0
+1 0 0 2 0.0
+0 1 0 2 0.0
+0 0 1 2 0.0
+0 0 0 3 0.0
+4 0 0 0 48.0
+3 1 0 0 24.0
+2 2 0 0 20.0
+1 3 0 0 18.0
+0 4 0 0 27.0
+3 0 1 0 12.0
+2 1 1 0 8.0
+1 2 1 0 7.0
+0 3 1 0 9.0
+2 0 2 0 10.0
+1 1 2 0 6.0
+0 2 2 0 8.0
+1 0 3 0 6.0
+0 1 3 0 6.0
+0 0 4 0 12.0
+3 0 0 1 12.0
+2 1 0 1 12.0
+1 2 0 1 11.0
+0 3 0 1 18.0
+2 0 1 1 6.0
+1 1 1 1 5.0
+0 2 1 1 7.0
+1 0 2 1 4.0
+0 1 2 1 6.0
+0 0 3 1 6.0
+2 0 0 2 14.0
+1 1 0 2 10.0
+0 2 0 2 17.0
+1 0 1 2 5.0
+0 1 1 2 7.0
+0 0 2 2 8.0
+1 0 0 3 9.0
+0 1 0 3 18.0
+0 0 1 3 9.0
+0 0 0 4 27.0
+5 0 0 0 0.0
+4 1 0 0 0.0
+3 2 0 0 0.0
+2 3 0 0 0.0
+1 4 0 0 0.0
+0 5 0 0 0.0
+4 0 1 0 0.0
+3 1 1 0 0.0
+2 2 1 0 0.0
+1 3 1 0 0.0
+0 4 1 0 0.0
+3 0 2 0 0.0
+2 1 2 0 0.0
+1 2 2 0 0.0
+0 3 2 0 0.0
+2 0 3 0 0.0
+1 1 3 0 0.0
+0 2 3 0 0.0
+1 0 4 0 0.0
+0 1 4 0 0.0
+0 0 5 0 0.0
+4 0 0 1 0.0
+3 1 0 1 0.0
+2 2 0 1 0.0
+1 3 0 1 0.0
+0 4 0 1 0.0
+3 0 1 1 0.0
+2 1 1 1 0.0
+1 2 1 1 0.0
+0 3 1 1 0.0
+2 0 2 1 0.0
+1 1 2 1 0.0
+0 2 2 1 0.0
+1 0 3 1 0.0
+0 1 3 1 0.0
+0 0 4 1 0.0
+3 0 0 2 0.0
+2 1 0 2 0.0
+1 2 0 2 0.0
+0 3 0 2 0.0
+2 0 1 2 0.0
+1 1 1 2 0.0
+0 2 1 2 0.0
+1 0 2 2 0.0
+0 1 2 2 0.0
+0 0 3 2 0.0
+2 0 0 3 0.0
+1 1 0 3 0.0
+0 2 0 3 0.0
+1 0 1 3 0.0
+0 1 1 3 0.0
+0 0 2 3 0.0
+1 0 0 4 0.0
+0 1 0 4 0.0
+0 0 1 4 0.0
+0 0 0 5 0.0
+6 0 0 0 960.0
+5 1 0 0 480.0
+4 2 0 0 336.0
+3 3 0 0 264.0
+2 4 0 0 252.0
+1 5 0 0 270.0
+0 6 0 0 405.0
+5 0 1 0 240.0
+4 1 1 0 144.0
+3 2 1 0 108.0
+2 3 1 0 96.0
+1 4 1 0 99.0
+0 5 1 0 135.0
+4 0 2 0 144.0
+3 1 2 0 84.0
+2 2 2 0 70.0
+1 3 2 0 66.0
+0 4 2 0 90.0
+3 0 3 0 78.0
+2 1 3 0 54.0
+1 2 3 0 48.0
+0 3 3 0 60.0
+2 0 4 0 72.0
+1 1 4 0 48.0
+0 2 4 0 60.0
+1 0 5 0 60.0
+0 1 5 0 60.0
+0 0 6 0 120.0
+5 0 0 1 240.0
+4 1 0 1 192.0
+3 2 0 1 156.0
+2 3 0 1 156.0
+1 4 0 1 171.0
+0 5 0 1 270.0
+4 0 1 1 96.0
+3 1 1 1 72.0
+2 2 1 1 66.0
+1 3 1 1 69.0
+0 4 1 1 99.0
+3 0 2 1 54.0
+2 1 2 1 48.0
+1 2 2 1 46.0
+0 3 2 1 66.0
+2 0 3 1 42.0
+1 1 3 1 36.0
+0 2 3 1 48.0
+1 0 4 1 36.0
+0 1 4 1 48.0
+0 0 5 1 60.0
+4 0 0 2 192.0
+3 1 0 2 132.0
+2 2 0 2 130.0
+1 3 0 2 138.0
+0 4 0 2 225.0
+3 0 1 2 66.0
+2 1 1 2 58.0
+1 2 1 2 59.0
+0 3 1 2 87.0
+2 0 2 2 50.0
+1 1 2 2 42.0
+0 2 2 2 62.0
+1 0 3 2 36.0
+0 1 3 2 48.0
+0 0 4 2 60.0
+3 0 0 3 114.0
+2 1 0 3 120.0
+1 2 0 3 123.0
+0 3 0 3 210.0
+2 0 1 3 60.0
+1 1 1 3 57.0
+0 2 1 3 87.0
+1 0 2 3 42.0
+0 1 2 3 66.0
+0 0 3 3 60.0
+2 0 0 4 144.0
+1 1 0 4 126.0
+0 2 0 4 225.0
+1 0 1 4 63.0
+0 1 1 4 99.0
+0 0 2 4 90.0
+1 0 0 5 135.0
+0 1 0 5 270.0
+0 0 1 5 135.0
+0 0 0 6 405.0
diff --git a/test/data/map_partial_orders_6D4.txt b/test/data/map_partial_orders_6D4.txt
new file mode 100644
index 0000000..8d5ff3c
--- /dev/null
+++ b/test/data/map_partial_orders_6D4.txt
@@ -0,0 +1,211 @@
+# map prtial orders 6D4
+0 0 0 0
+1 0 0 0
+0 1 0 0
+0 0 1 0
+0 0 0 1
+2 0 0 0
+1 1 0 0
+0 2 0 0
+1 0 1 0
+0 1 1 0
+0 0 2 0
+1 0 0 1
+0 1 0 1
+0 0 1 1
+0 0 0 2
+3 0 0 0
+2 1 0 0
+1 2 0 0
+0 3 0 0
+2 0 1 0
+1 1 1 0
+0 2 1 0
+1 0 2 0
+0 1 2 0
+0 0 3 0
+2 0 0 1
+1 1 0 1
+0 2 0 1
+1 0 1 1
+0 1 1 1
+0 0 2 1
+1 0 0 2
+0 1 0 2
+0 0 1 2
+0 0 0 3
+4 0 0 0
+3 1 0 0
+2 2 0 0
+1 3 0 0
+0 4 0 0
+3 0 1 0
+2 1 1 0
+1 2 1 0
+0 3 1 0
+2 0 2 0
+1 1 2 0
+0 2 2 0
+1 0 3 0
+0 1 3 0
+0 0 4 0
+3 0 0 1
+2 1 0 1
+1 2 0 1
+0 3 0 1
+2 0 1 1
+1 1 1 1
+0 2 1 1
+1 0 2 1
+0 1 2 1
+0 0 3 1
+2 0 0 2
+1 1 0 2
+0 2 0 2
+1 0 1 2
+0 1 1 2
+0 0 2 2
+1 0 0 3
+0 1 0 3
+0 0 1 3
+0 0 0 4
+5 0 0 0
+4 1 0 0
+3 2 0 0
+2 3 0 0
+1 4 0 0
+0 5 0 0
+4 0 1 0
+3 1 1 0
+2 2 1 0
+1 3 1 0
+0 4 1 0
+3 0 2 0
+2 1 2 0
+1 2 2 0
+0 3 2 0
+2 0 3 0
+1 1 3 0
+0 2 3 0
+1 0 4 0
+0 1 4 0
+0 0 5 0
+4 0 0 1
+3 1 0 1
+2 2 0 1
+1 3 0 1
+0 4 0 1
+3 0 1 1
+2 1 1 1
+1 2 1 1
+0 3 1 1
+2 0 2 1
+1 1 2 1
+0 2 2 1
+1 0 3 1
+0 1 3 1
+0 0 4 1
+3 0 0 2
+2 1 0 2
+1 2 0 2
+0 3 0 2
+2 0 1 2
+1 1 1 2
+0 2 1 2
+1 0 2 2
+0 1 2 2
+0 0 3 2
+2 0 0 3
+1 1 0 3
+0 2 0 3
+1 0 1 3
+0 1 1 3
+0 0 2 3
+1 0 0 4
+0 1 0 4
+0 0 1 4
+0 0 0 5
+6 0 0 0
+5 1 0 0
+4 2 0 0
+3 3 0 0
+2 4 0 0
+1 5 0 0
+0 6 0 0
+5 0 1 0
+4 1 1 0
+3 2 1 0
+2 3 1 0
+1 4 1 0
+0 5 1 0
+4 0 2 0
+3 1 2 0
+2 2 2 0
+1 3 2 0
+0 4 2 0
+3 0 3 0
+2 1 3 0
+1 2 3 0
+0 3 3 0
+2 0 4 0
+1 1 4 0
+0 2 4 0
+1 0 5 0
+0 1 5 0
+0 0 6 0
+5 0 0 1
+4 1 0 1
+3 2 0 1
+2 3 0 1
+1 4 0 1
+0 5 0 1
+4 0 1 1
+3 1 1 1
+2 2 1 1
+1 3 1 1
+0 4 1 1
+3 0 2 1
+2 1 2 1
+1 2 2 1
+0 3 2 1
+2 0 3 1
+1 1 3 1
+0 2 3 1
+1 0 4 1
+0 1 4 1
+0 0 5 1
+4 0 0 2
+3 1 0 2
+2 2 0 2
+1 3 0 2
+0 4 0 2
+3 0 1 2
+2 1 1 2
+1 2 1 2
+0 3 1 2
+2 0 2 2
+1 1 2 2
+0 2 2 2
+1 0 3 2
+0 1 3 2
+0 0 4 2
+3 0 0 3
+2 1 0 3
+1 2 0 3
+0 3 0 3
+2 0 1 3
+1 1 1 3
+0 2 1 3
+1 0 2 3
+0 1 2 3
+0 0 3 3
+2 0 0 4
+1 1 0 4
+0 2 0 4
+1 0 1 4
+0 1 1 4
+0 0 2 4
+1 0 0 5
+0 1 0 5
+0 0 1 5
+0 0 0 6
\ No newline at end of file
diff --git a/test/extraction.jl b/test/extraction.jl
new file mode 100644
index 0000000..e6b4c53
--- /dev/null
+++ b/test/extraction.jl
@@ -0,0 +1,17 @@
+@testset "Hessian" begin
+
+    DACE.init(4,4)
+
+    v = [DACE.random(-1) for i in 1:4]
+    a = AlgebraicVector(v)
+
+    H1 = DACE.hessian(a)
+    H2 = DACE.hessian(v)
+
+    H3 = DACE.hess_stack(a)
+    H4 = DACE.hess_stack(v)
+
+    @test all(stack(H1,dims=3) .== stack(H2,dims=3))
+    @test all(H3 .== H4)
+    @test all(stack(H1,dims=3) .== H3)
+end
\ No newline at end of file
diff --git a/test/factory.jl b/test/factory.jl
new file mode 100644
index 0000000..bf07f30
--- /dev/null
+++ b/test/factory.jl
@@ -0,0 +1,33 @@
+@testset verbose = true "Identity" begin
+
+    @testset verbose = true "First N" begin
+
+        DACE.init(4,4)
+        av = DACE.identity(2)
+        ae = [DA(1, 1.0), DA(2, 1.0)]
+
+        @test length(av) == 2
+        @test all(av .== ae)
+    end
+
+    @testset verbose = true "Unsorted" begin
+
+        DACE.init(4,4)
+        av = DACE.identity(Vector{UInt32}([3, 4, 1, 1, 3]), false)
+        ae = [DA(3, 1.0), DA(4, 1.0), DA(1, 1.0), DA(1, 1.0), DA(3, 1.0)]
+
+        @test length(av) == 5
+        @test all(av .== ae)
+    end
+
+    @testset verbose = true "Sorted & Unique" begin
+
+        DACE.init(4,4)
+        av = DACE.identity(Vector{UInt32}([3, 1, 1, 4, 3, 3]), true)
+        ae = [DA(1, 1.0), 0.0, DA(3, 1.0), DA(4, 1.0)]
+
+        @test length(av) == 4 # max number of variables
+        @test all(av .== ae)
+    end
+
+end
\ No newline at end of file
diff --git a/test/linear_algebra.jl b/test/linear_algebra.jl
new file mode 100644
index 0000000..488c804
--- /dev/null
+++ b/test/linear_algebra.jl
@@ -0,0 +1,95 @@
+using LinearAlgebra, GenericLinearAlgebra
+using ForwardDiff, DifferentiableEigen
+import ForwardDiff: Dual
+
+@testset verbose = true "Eigendecomposition" begin
+
+    DACE.init(1,6)
+
+    A = [4.0 2.0 1.0; 2.0 3.0 1.0; 1.0 1.0 2.0]
+    Ada = Matrix{DA}(undef,3,3)
+    Afd = Matrix{Dual{Nothing,Float64,6}}(undef,3,3)
+
+    k = 1
+    p = zeros(6)
+    for i = 1:3
+        for j = 1:i
+            p[k] = 1.0
+            Ada[i,j] = A[i,j] + DA(k,1.0)
+            Afd[i,j] = Dual(A[i,j], p...)
+            p[k] = 0.0
+            k += 1
+        end
+    end
+
+    for i = 1:3
+        for j = 1:i-1
+            Ada[j,i] = Ada[i,j]
+            Afd[j,i] = Afd[i,j]
+        end
+    end
+
+    Acc = AlgebraicMatrix{DA}(3,3)
+    Acc .= Ada
+
+    Fda = LinearAlgebra.eigen(Hermitian(Ada))
+    Ffd = LinearAlgebra.eigen(Hermitian(Afd))
+
+    λda, Vda = Fda.values, Fda.vectors          # GenericLinearAlgebra + DACE
+    λcc, Vcc = DACE.eigh(Acc)                   # Eigen3 + DACE
+    λfd, Vfd = Ffd.values, Ffd.vectors          # GenericLinearAlgebra + ForwardDiff
+    λde, Vde = DifferentiableEigen.eigen(Afd)   # DifferentiableEigen + ForwardDiff
+    λde, Vde = λde[1:2:end], Vde[1:2:end]
+
+    cλda = DACE.cons.(λda)
+    cλcc = DACE.cons.(λcc)
+    cλfd = [λ.value for λ in λfd]
+    cλde = [λ.value for λ in λde]
+
+    @testset verbose = true "Eigenvalues constant part" begin
+        @test all(isapprox.(cλda, cλcc, atol=1e-14, rtol=1e-14))
+        @test all(isapprox.(cλda, cλfd, atol=1e-14, rtol=1e-14))
+        @test all(isapprox.(cλda, cλde, atol=1e-14, rtol=1e-14))
+    end
+
+    Jλda = vcat([DACE.linear.(λda)...]'...)
+    Jλcc = vcat([DACE.linear.(λcc)...]'...)
+    Jλfd = vcat([[λ.partials...] for λ in λfd]'...)
+    Jλde = vcat([[λ.partials...] for λ in λde]'...)
+
+    @testset verbose = true "Eigenvalues partials" begin
+        @test all(isapprox.(Jλda, Jλcc, atol=1e-13, rtol=1e-13))
+        @test all(isapprox.(Jλda, Jλfd, atol=1e-14, rtol=1e-14))
+        @test all(isapprox.(Jλda, Jλde, atol=1e-14, rtol=1e-14))
+    end
+
+    cVda = DACE.cons.(Vda)
+    cVcc = DACE.cons.(Vcc)
+    cVfd = [V.value for V in Vfd]
+    cVde = reshape([V.value for V in Vde],3,3)
+
+    scc = [sign.(cVda[1,i]) == sign.(cVcc[1,i]) ? 1 : -1 for i in 1:3]'
+    sfd = [sign.(cVda[1,i]) == sign.(cVfd[1,i]) ? 1 : -1 for i in 1:3]'
+    sde = [sign.(cVda[1,i]) == sign.(cVde[1,i]) ? 1 : -1 for i in 1:3]'
+
+    @testset verbose = true "Eigenvectors constant part" begin
+        @test all(isapprox.(cVda, cVcc.*scc, atol=1e-14, rtol=1e-14))
+        @test all(isapprox.(cVda, cVfd.*sfd, atol=1e-14, rtol=1e-14))
+        @test all(isapprox.(cVda, cVde.*sde, atol=1e-14, rtol=1e-14))
+    end
+
+    JVda = [[DACE.linear(V)...] for V in Vda]
+    JVcc = [[DACE.linear(V)...] for V in Vcc]
+    JVfd = [[V.partials...] for V in Vfd]
+    JVde = reshape([[V.partials...] for V in Vde],3,3)
+
+    @testset verbose = true "Eigenvectors partials" begin
+        for i in 1:3
+            for j in 1:3
+                @test all(isapprox.(JVda[i,j], JVcc[i,j].*scc[j], atol=1e-13, rtol=1e-13))
+                @test all(isapprox.(JVda[i,j], JVfd[i,j].*sfd[j], atol=1e-13, rtol=1e-13))
+                @test all(isapprox.(JVda[i,j], JVde[i,j].*sde[j], atol=1e-13, rtol=1e-13))
+            end
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 69ebb5e..b2a50b8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,5 +17,23 @@ include("utils.jl")
         include("comparison_operators.jl")
     end
 
-    include("special_functions.jl")
+    @testset verbose = true "Special Functions" begin
+        include("special_functions.jl")
+    end
+
+    @testset verbose = true "Linear Algebra" begin
+        include("linear_algebra.jl")
+    end
+
+    @testset verbose = true "Access & extraction" begin
+        include("extraction.jl")
+    end
+
+    @testset verbose = true "Statistics" begin
+        include("statistics.jl")
+    end
+
+    @testset verbose = true "Factories" begin
+        include("factory.jl")
+    end
 end
diff --git a/test/special_functions.jl b/test/special_functions.jl
index 247f437..75b22c4 100644
--- a/test/special_functions.jl
+++ b/test/special_functions.jl
@@ -1,20 +1,17 @@
 using SpecialFunctions
 
-@testset verbose = true "Special functions" begin
-    @testset "Test erf" begin
-        DACE.init(1,3)
+@testset "Test erf" begin
+    DACE.init(1,3)
 
-        x = DACE.random(-1)
-        @test isapprox(DACE.cons(erf(x)), erf(DACE.cons(x)), atol=1e-15, rtol=1e-15)
+    x = DACE.random(-1)
+    @test isapprox(DACE.cons(erf(x)), erf(DACE.cons(x)), atol=1e-15, rtol=1e-15)
 
-    end
-
-    @testset "Test erfc" begin
-        DACE.init(1,3)
+end
 
-        x = DACE.random(-1)
-        @test isapprox(DACE.cons(erfc(x)), erfc(DACE.cons(x)), atol=1e-15, rtol=1e-15)
+@testset "Test erfc" begin
+    DACE.init(1,3)
 
-    end
+    x = DACE.random(-1)
+    @test isapprox(DACE.cons(erfc(x)), erfc(DACE.cons(x)), atol=1e-15, rtol=1e-15)
 
 end
diff --git a/test/statistics.jl b/test/statistics.jl
new file mode 100644
index 0000000..51b90b0
--- /dev/null
+++ b/test/statistics.jl
@@ -0,0 +1,76 @@
+using DelimitedFiles
+
+@testset verbose = true "Moments" begin
+    @testset verbose = true "Multi indices" begin
+
+        file = joinpath(@__DIR__,"data/map_partial_orders_6D4.txt")
+        data = readdlm(file,' ',Int64,comments=true,comment_char='#')
+        @test size(data,1) == binomial(6+4,4)
+
+        idx = DACE.getMultiIndices(6,4)
+        ida = [data[k,:] for k in axes(data,1)]
+        iin = indexin(idx,ida)
+
+        for k in eachindex(idx)
+            @test all(idx[k] .== ida[iin[k]])
+        end
+    end
+    @testset verbose = true "Raw Moments" begin
+
+        s1 = ("integer", "double")
+        s2 = ("", "non")
+        atol = [1e-15 1e-14; 1e-15 5e-12]
+        for i in 1:2
+            for j in 1:2
+
+                file = joinpath(@__DIR__,"data/"*s1[i]*"_moments_"*s2[j]*"zero_6D4.txt")
+                data = readdlm(file,' ',Float64)
+
+                ord = convert(Int64,data[1,1])
+                dim = convert(Int64,data[1,2])
+
+                μ = data[2:dim+1,1]
+                P = data[2:dim+1,2:end]
+
+                M(t) = exp(μ'*t + 0.5*t'*P*t)
+
+                DACE.init(ord,dim)
+                Mt = M([DA(i,1.0) for i in 1:dim])
+                ix, m0 = DACE.getRawMoments(Mt,ord)
+
+                ia = [UInt32.(data[dim+1+k,1:end-1]) for k in eachindex(ix)]
+                ii = indexin(ix,ia)
+                ma = data[dim+2:end,end][ii]
+
+                @test all(isapprox.(m0, ma; atol=atol[i,j], rtol=1e-15))
+            end
+        end
+    end
+    @testset verbose = true "Central Moments" begin
+        for s in ("integer", "double")
+
+            file_raw = joinpath(@__DIR__,"data/"*s*"_moments_nonzero_6D4.txt")
+            file_ctr = joinpath(@__DIR__,"data/"*s*"_moments_zero_6D4.txt")
+            data_raw = readdlm(file_raw,' ',Float64)
+            data_ctr = readdlm(file_ctr,' ',Float64)
+
+            ord = convert(Int64,data_raw[1,1])
+            dim = convert(Int64,data_raw[1,2])
+
+            μ_raw = data_raw[2:dim+1,1]
+            P_raw = data_raw[2:dim+1,2:end]
+
+            M(t) = exp(μ_raw'*t + 0.5*t'*P_raw*t)
+
+            DACE.init(ord,dim)
+            Mt = M([DA(i,1.0) for i in 1:dim])
+            ix, m0 = DACE.getCentralMoments(Mt,ord)
+
+            ia = [UInt32.(data_ctr[dim+1+k,1:end-1]) for k in eachindex(ix)]
+            ii = indexin(ix,ia)
+            ma = data_ctr[dim+2:end,end][ii]
+
+            @test all(isapprox.(m0, ma; atol=1e-11, rtol=1e-15))
+        end
+    end
+end