diff --git a/test/gauss.jl b/test/gauss.jl
index 65c9620..edef5a5 100644
--- a/test/gauss.jl
+++ b/test/gauss.jl
@@ -1,11 +1,11 @@
 function gauss(p::Vector{Float64}; x::Float64=0.0)
   λ, μ, σ = p
-  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2/2σ^2)
+  return λ / sqrt(2π * σ^2) * exp(-(x - μ)^2 / 2σ^2)
 end
 
 function J_gauss(p::Vector{Float64}; x::Float64=0.0)
   λ, μ, σ = p
-  c = [(1/λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2)/σ^3]
+  c = [(1 / λ) ((x - μ) / σ^2) ((x - μ)^2 - σ^2) / σ^3]
   return c .* gauss(p; x=x)
 end
 
@@ -14,13 +14,13 @@ function H_gauss(p::Vector{Float64}; x::Float64=0.0)
   c = zeros(3, 3)
 
   # c[1,1] = 0
-  c[1, 2] = c[2, 1] = (x-μ) / (λ * σ^2)
-  c[1, 3] = c[3, 1] = -((μ + σ - x)*(-μ + σ + x))/(λ*σ^3)
+  c[1, 2] = c[2, 1] = (x - μ) / (λ * σ^2)
+  c[1, 3] = c[3, 1] = -((μ + σ - x) * (-μ + σ + x)) / (λ * σ^3)
 
   c[2, 2] = ((x - μ)^2 - σ^2) / σ^4
-  c[2, 3] = c[3, 2] = -((μ - x)*((μ - x)^2 - 3σ^2))/σ^5
+  c[2, 3] = c[3, 2] = -((μ - x) * ((μ - x)^2 - 3σ^2)) / σ^5
 
-  c[3, 3] = (-5*σ^2*(μ - x)^2 + (μ - x)^4 + 2*σ^4)/σ^6
+  c[3, 3] = (-5 * σ^2 * (μ - x)^2 + (μ - x)^4 + 2 * σ^4) / σ^6
   return c .* gauss(p; x=x)
 end
 
@@ -38,20 +38,54 @@ end
   @test TaylorTest.check(J_gauss, H_gauss, p; x=x)
 end
 
+function check_symmetry(a)
+  R = true
+  if ndims(a) == 2
+    R = norm(a - a') < 1e-10
+  else
+    for i in axes(a, 1)
+      r = a[i, :, :]' == a[i, :, :]
+      R &= r
+      if !r
+        n = maximum(abs, a[i, :, :]' - a[i, :, :])
+        @warn "Tensor not symmetric at i=$i: $n"
+      end
+    end
+  end
+  return R
+end
+
 @testset "Composed Gauss function" begin
-  p = 1 .+ 5 * rand(6)
+  p = 1 .+ 5rand(6)
 
-  x = p[3] + p[4] + 2p[5] + sqrt(p[6]) * (2rand() - 1)
-
-  p_aux = _p -> [_p[1]/_p[2], _p[3] + _p[4] + 2_p[5], _p[6]]
+  p_aux = _p -> [_p[1] / _p[2], _p[3] + _p[4] + 2_p[5], _p[6]]
 
   J_p_aux = _p -> [1/_p[2] -_p[1]/_p[2]^2 0 0 0 0;
-    0            0 1 1 2 0;
-    0            0 0 0 0 1]
+    0 0 1 1 2 0;
+    0 0 0 0 0 1]
 
-  g = (_p; x=x) -> gauss(p_aux(_p); x=x)
-  Jg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
-  Hg = (_p; x=x) -> J_gauss(p_aux(_p); x=x) * J_p_aux(p)
+  H_p_aux = _p -> begin
+    h = zeros(3, 6, 6)
+    h[1, 1, 2] = h[1, 2, 1] = -1 / _p[2]^2
+    h[1, 2, 2] = 2_p[1] / _p[2]^3
+    h
+  end
+
+  g = (_p; x = 0.0) -> gauss(p_aux(_p); x=x)
+  Jg = (_p; x = 0.0) -> J_gauss(p_aux(_p); x=x) * J_p_aux(_p)
+  Hg = (_p; x = 0.0) -> begin
+    p = p_aux(_p)
+    _J_gauss = vec(J_gauss(p; x=x))
+    _H_p = H_p_aux(_p)
+    _J_p = J_p_aux(_p)
+    r = _J_p' * H_gauss(p; x=x) * _J_p
+    r .+= TO.tensorcontract([2, 3], _J_gauss, [1], _H_p, [1, 2, 3])
+    r
+  end
+
+  _p = p_aux(p)
+  x = _p[2] + _p[2] * 2 * (rand() - 0.5)
 
   @test TaylorTest.check(g, Jg, p; x=x)
+  @test TaylorTest.check(Jg, Hg, p; x=x)
 end