diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 337ce0f..3da6feb 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,8 +18,9 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.0'
           - '1.7'
+          - '1.8'
+          - '1.9'
           - '1.10'
           - 'nightly'
         os:
@@ -48,3 +49,4 @@ jobs:
       - uses: codecov/codecov-action@v4
         with:
           files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.gitignore b/.gitignore
index 0f84bed..db41b91 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
 *.jl.cov
 *.jl.mem
 /Manifest.toml
+benchmark/tune.json
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 270fb85..c131f97 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,20 +1,18 @@
 name = "UniformIsingModels"
 uuid = "91393fb9-e38b-4a31-9409-aa579b4163ad"
 authors = ["Stefano Crotti <stefano.crotti@polito.it> and contributors"]
-version = "0.3.0"
+version = "0.4.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
-OffsetArrays = "1"
-UnPack = "1"
 LogExpFunctions = "0.3"
-julia = "1"
+OffsetArrays = "1"
+julia = "1.7"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index 91a1b76..880238d 100644
--- a/README.md
+++ b/README.md
@@ -3,21 +3,24 @@
 [![Build Status](https://github.com/stecrotti/UniformIsingModels.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/stecrotti/UniformIsingModels.jl/actions/workflows/CI.yml?query=branch%3Amain)
 [![Coverage](https://codecov.io/gh/stecrotti/UniformIsingModels.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/stecrotti/UniformIsingModels.jl)
 
-A fully-connected ferromagnetic Ising model with uniform coupling strength, described by a Boltzmann distribution
+A fully-connected ferromagnetic [Ising model](https://en.wikipedia.org/wiki/Ising_model) with uniform coupling strength, described by a Boltzmann distribution
 
->![equation](https://latex.codecogs.com/svg.image?p(\boldsymbol\sigma|J,&space;\boldsymbol{h},&space;\beta)&space;=&space;\frac{1}{Z_{J,&space;\boldsymbol{h},&space;\beta}}\exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j&space;&plus;\sum_{i=1}^Nh_i\sigma_i\right)\right],\quad\boldsymbol\sigma\in\\{-1,1\\}^N)
+$p(\boldsymbol{\sigma}) = \frac{1}{Z} \exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right],\quad \boldsymbol{\sigma}\in\{-1,1\}^N $
 
 is exactly solvable in polynomial time.
 
 
-| Quantity | Cost          |
-| ------------- | ----------- |
-| Normalization, Free energy      |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2)) |
-| Sample a configuration      |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2)) |
-| Average energy, Entropy |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2))  |
-| Distribution of the sum of the N spins | ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^2))     |
-| Site magnetizations     | ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^3))     |
-| Pair magnetizations, Correlations |  ![equation](https://latex.codecogs.com/svg.image?\mathcal{O}(N^5))  |
+| Quantity | Expression | Cost          |
+| ------------- | ----------| ----------- |
+| Normalization | $Z=\sum\limits_{\boldsymbol{\sigma}}\exp\left[\beta\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right]$ | $\mathcal O (N^2)$ |
+| Free energy | $F = -\frac{1}{\beta}\log Z$ | $\mathcal O (N^2)$ |
+| Sample a configuration | $\boldsymbol{\sigma} \sim p(\boldsymbol{\sigma})$ | $\mathcal O (N^2)$ |
+| Average energy | $U = \sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\left[-\left(\frac{J}{N}\sum_{i<j}\sigma_i\sigma_j+\sum_{i=1}^Nh_i\sigma_i\right)\right]$ | $\mathcal O (N^2)$ |
+| Entropy | $S = -\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\log p(\boldsymbol{\sigma})$ | $\mathcal O (N^2)$ |
+| Distribution of the sum of the N spins | $p_S(s)=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\delta\left(s-\sum_{i=1}^N\sigma_i\right)$ | $\mathcal O (N^2)$ |
+| Site magnetizations     | $m_i=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\sigma_i,\quad\forall i\in\{1,2,\ldots,N\}$ | $\mathcal O (N^3)$ |
+| Correlations     | $r_{ij}=\sum\limits_{\boldsymbol{\sigma}}p(\boldsymbol{\sigma})\sigma_i\sigma_j,\quad\forall j\in\{1,2,\ldots,N\},i<j$ | $\mathcal O (N^5)$ |
+
 
 ## Example
 ```
@@ -37,7 +40,7 @@ x = UniformIsing(N, J, h, β)
 Compute stuff
 ```
 # normalization and free energy
-Z = exp(x.logZ)
+Z = normalization(x)
 F = free_energy(x)
 
 # energy and probability of a configuration
@@ -60,10 +63,9 @@ U = avg_energy(x)
 # entropy
 S = entropy(x)
 
-# pairwise magnetizations <σᵢσⱼ> and correlations
-p = pair_magnetizations(x)
-c = correlations(x)
-
+# correlations <σᵢσⱼ> and covariances <σᵢσⱼ>-<σᵢ><σⱼ>
+p = correlations(x)
+c = covariances(x)
 ```
 
 ## Notes
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
new file mode 100644
index 0000000..cfbea13
--- /dev/null
+++ b/benchmark/benchmarks.jl
@@ -0,0 +1,20 @@
+using BenchmarkTools
+using UniformIsingModels
+
+SUITE = BenchmarkGroup()
+
+N = 100
+J = 0.5
+h = 1.2 * randn(N)
+β = 2.3
+x = UniformIsing(N, J, h, β)
+
+SUITE["constructor"] = BenchmarkGroup()
+SUITE["constructor"]["constructor"] = @benchmarkable UniformIsing($N, $J, $h, $β)
+
+SUITE["observables"] = BenchmarkGroup()
+SUITE["observables"]["normalization"] = @benchmarkable normalization(x)
+SUITE["observables"]["entropy"] = @benchmarkable entropy(x)
+SUITE["observables"]["site_magnetizations"] = @benchmarkable site_magnetizations(x)
+# SUITE["observables"]["pair_magnetizations"] = @benchmarkable pair_magnetizations(x)
+SUITE["observables"]["sum_distribution"] = @benchmarkable sum_distribution(x)
\ No newline at end of file
diff --git a/src/UniformIsingModels.jl b/src/UniformIsingModels.jl
index 460735e..86e05d0 100644
--- a/src/UniformIsingModels.jl
+++ b/src/UniformIsingModels.jl
@@ -1,193 +1,20 @@
 module UniformIsingModels
 
-import OffsetArrays: OffsetVector, fill
-import LinearAlgebra: dot
-import Base: show
-import UnPack: @unpack
-import Random: GLOBAL_RNG, AbstractRNG
-import LogExpFunctions: logsumexp, logaddexp, logsubexp
-
-export UniformIsing, energy, normalization, pdf, 
+using OffsetArrays: OffsetVector, fill
+using LinearAlgebra: dot
+using Random: default_rng, AbstractRNG
+using LogExpFunctions: logsumexp, logaddexp
+
+export UniformIsing, nvariables, variables, recompute_partials!,
+        energy, lognormalization, normalization, pdf,
+        avg_energy, entropy, free_energy, 
         site_magnetizations!, site_magnetizations,
-        pair_magnetizations!, pair_magnetizations,
         correlations!, correlations,
+        covariances, covariances,
         sum_distribution!, sum_distribution,
-        sample!, sample,
-        avg_energy, entropy, free_energy
+        sample!, sample
 
 include("accumulate.jl")
-
-struct UniformIsing{T<:Real, U<:OffsetVector}
-    N    :: Int                # number of spins
-    J    :: T                  # uniform coupling strength
-    h    :: Vector{T}          # external fields
-    β    :: T                  # inverse temperature
-    L    :: OffsetVector{U, Vector{U}} # partial sums from the left
-    R    :: OffsetVector{U, Vector{U}} # partial sums from the right
-    dLdB ::  OffsetVector{U, Vector{U}} # Derivative of L wrt β
-    logZ :: T                  # normalization
-    function UniformIsing(N::Int, J::T, h::Vector{T}, β::T=1.0) where T
-        @assert length(h) == N
-        @assert β ≥ 0
-        # L = accumulate_left(h, β)
-        R = accumulate_right(h, β)
-        L, dLdB = accumulate_d_left(h, β)
-        logZ = logsumexp( β*J/2*(s^2/N-1) + L[N][s] for s in -N:N )
-        new{T, eltype(L)}(N, J, h, β, L, R, dLdB, logZ)
-    end
-end
-
-function show(io::IO, x::UniformIsing)
-    @unpack N, J, h, β = x
-    println(io, "UniformIsing with N = $N variables at temperature β = $β")
-end
-
-function energy(x::UniformIsing, σ)
-    s = sum(σ)
-    f = dot(σ, x.h)
-    -( x.J/2*(s^2/x.N-1) + f ) 
-end
-
-pdf(x::UniformIsing, σ) = exp(-x.β*energy(x, σ) - x.logZ)
-
-free_energy(x::UniformIsing) = -1/x.β*x.logZ
-
-function sample_spin(rng::AbstractRNG, p::Real)
-    @assert 0 ≤ p ≤ 1
-    r = rand(rng)
-    r < p ? 1 : -1
-end
-
-# return a sample along with its probability
-function sample!(rng::AbstractRNG, σ, x::UniformIsing)
-    @unpack N, J, h, β, L, R, logZ = x
-    a = 0.0; b = 0
-    f(s) = β*J/2*(s^2/N-1)
-    for i in 1:N
-        tmp_plus = tmp_minus = 0.0
-        for s in -N:N
-            tmp_plus += exp(f(b+1+s) + R[i+1][s])
-            tmp_minus += exp(f(b-1+s) + R[i+1][s])
-        end
-        p_plus = exp(β*h[i]) * tmp_plus
-        p_minus = exp(-β*h[i]) * tmp_minus
-        p_i = p_plus / (p_plus + p_minus)
-        σi = sample_spin(rng, p_i)
-        σ[i] = σi
-        a += h[i]*σi
-        b += σi
-    end
-    p = exp(f(b) + β*a - logZ)
-    @assert a == dot(h, σ); @assert b == sum(σ)
-    σ, p
-end
-sample!(σ, x::UniformIsing) = sample!(GLOBAL_RNG, σ, x)
-sample(rng::AbstractRNG, x::UniformIsing) = sample!(rng, zeros(Int, x.N), x)
-sample(x::UniformIsing) = sample(GLOBAL_RNG, x)
-
-# first store in `p[i]` the quantity log(p(σᵢ=+1)), then transform at the end 
-function site_magnetizations!(p, x::UniformIsing)
-    @unpack N, J, h, β, L, R, logZ = x
-    f(s) = β*J/2*(s^2/N-1)
-    for i in eachindex(p)
-        p[i] = -Inf
-        for sL in -N:N
-            for sR in -N:N
-                s = sL + sR
-                p[i] = logaddexp(p[i], f(s+1) + β*h[i] + L[i-1][sL] + R[i+1][sR])
-            end
-        end
-        # include normalization
-        p[i] = exp(p[i] - logZ)
-        # transform form p(+) to m=2p(+)-1
-        p[i] = 2*p[i] - 1
-    end
-    p
-end
-function site_magnetizations(x::UniformIsing{T,U}) where {T,U} 
-    site_magnetizations!(zeros(T,x.N), x)
-end
-
-# distribution of the sum of all variables as an OffsetVector
-function sum_distribution!(p, x::UniformIsing)
-    @unpack N, J, β, L, logZ = x
-    for s in -N:N
-        p[s] = exp( β*J/2*(s^2/N-1) + L[N][s] - logZ ) 
-    end
-    p
-end
-function sum_distribution(x::UniformIsing{T,U}) where {T,U} 
-    p = fill(zero(T), -x.N:x.N)
-    sum_distribution!(p, x)
-end
-
-function avg_energy(x::UniformIsing{T}) where T
-    @unpack N, J, h, β, L, dLdB, logZ = x
-    Zt = sum( exp( β*J/2*(s^2/N-1) + L[N][s]) * (J/2*(s^2/N-1)+dLdB[N][s]) for s in -N:N)
-    -exp(log(Zt) - logZ)
-end
-
-entropy(x::UniformIsing; kw...) = x.β * (avg_energy(x; kw...) - free_energy(x)) 
-
-function pair_magnetizations!(m, x::UniformIsing{T,U};
-        M = accumulate_middle(x.h, x.β)) where {T,U}
-    @unpack N, J, h, β, L, R, logZ = x
-    Z = exp(logZ)
-    f(s) = β*J/2*(s^2/N-1)
-    for i in 1:N
-        # j = i
-        m[i,i] = 1
-        # j = i+1
-        j = i + 1
-        j > N && break
-        m[i,j] = 0
-        for sL in -N:N
-            for sR in -N:N
-                s = sL + sR
-                m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
-                            exp( f(s-2) - β*(h[i]+h[j]) ) -
-                            exp( f(s)   + β*(h[i]-h[j]) ) - 
-                            exp( f(s)   - β*(h[i]-h[j]) )   ) * 
-                            exp(L[i-1][sL] + R[j+1][sR]) 
-            end
-        end
-        m[i,j] /= Z; m[j,i] = m[i,j]
-        # j > i + 1
-        for j in i+2:N
-            m[i,j] = 0
-            for sM in -N:N
-                for sL in -N:N
-                    for sR in -N:N
-                        s = sL + sM + sR 
-                        m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
-                                  exp( f(s-2) - β*(h[i]+h[j]) ) -
-                                  exp( f(s)   + β*(h[i]-h[j]) ) - 
-                                  exp( f(s)   - β*(h[i]-h[j]) )   ) *
-                                  exp(L[i-1][sL] + M[i+1,j-1][sM] + R[j+1][sR])
-                    end 
-                end
-            end
-            m[i,j] /= Z; m[j,i] = m[i,j]
-        end
-    end
-    m
-end
-function pair_magnetizations(x::UniformIsing{T,U}; kw...) where {T,U} 
-    pair_magnetizations!(zeros(T,x.N,x.N), x; kw...)
-end
-
-function correlations!(c, x::UniformIsing;
-        m = site_magnetizations(x), p = pair_magnetizations(x))
-    N = x.N
-    for i in 1:N
-        for j in 1:N
-            c[i,j] = p[i,j] - m[i]*m[j]
-        end
-    end
-    c
-end
-function correlations(x::UniformIsing{T,U}; kw...) where {T,U} 
-    correlations!(zeros(T,x.N,x.N), x; kw...)
-end
+include("uniform_ising.jl")
 
 end # end module
\ No newline at end of file
diff --git a/src/accumulate.jl b/src/accumulate.jl
index 6b24416..e1efeed 100644
--- a/src/accumulate.jl
+++ b/src/accumulate.jl
@@ -17,30 +17,30 @@ function accumulate_left(h, β)
     accumulate_left!(L, h, β)
 end
 
-function accumulate_d_left!(L, dLdB, h, β)
+function accumulate_d_left!(L, dLdβ, h, β)
     N = length(h)
-    L[0] .= -Inf; L[0][0] = 0.0; dLdB[0] .= 0.0
+    L[0] .= -Inf; L[0][0] = 0.0; dLdβ[0] .= 0.0
     for k in 1:N
         for s in -(k-1):k-1
             L[k][s] = logaddexp( β*h[k] + L[k-1][s-1],
                                 -β*h[k] + L[k-1][s+1] )
-            dLdB[k][s] = (+h[k] + dLdB[k-1][s-1])*exp(+β*h[k]+L[k-1][s-1]) + 
-                       (-h[k] + dLdB[k-1][s+1])*exp(-β*h[k]+L[k-1][s+1])
-            dLdB[k][s] = dLdB[k][s] / exp(L[k][s])
-            isnan(dLdB[k][s]) && (dLdB[k][s] = 0.0)
+            dLdβ[k][s] = (+h[k] + dLdβ[k-1][s-1])*exp(+β*h[k]+L[k-1][s-1]) + 
+                       (-h[k] + dLdβ[k-1][s+1])*exp(-β*h[k]+L[k-1][s+1])
+            dLdβ[k][s] = dLdβ[k][s] / exp(L[k][s])
+            isnan(dLdβ[k][s]) && (dLdβ[k][s] = 0.0)
         end
         L[k][k] = β*h[k] + L[k-1][k-1]
         L[k][-k] = -β*h[k] + L[k-1][-k+1]
-        dLdB[k][k] = h[k] + dLdB[k-1][k-1]
-        dLdB[k][-k] = -h[k] + dLdB[k-1][-k+1]
+        dLdβ[k][k] = h[k] + dLdβ[k-1][k-1]
+        dLdβ[k][-k] = -h[k] + dLdβ[k-1][-k+1]
     end
-    L, dLdB
+    L, dLdβ
 end
 function accumulate_d_left(h, β)
     N = length(h)
     L = OffsetVector([fill(-Inf, -N:N) for i in 0:N], 0:N)
-    dLdB = OffsetVector([fill(0.0, -N:N) for i in 0:N], 0:N)
-    accumulate_d_left!(L, dLdB, h, β)
+    dLdβ = OffsetVector([fill(0.0, -N:N) for i in 0:N], 0:N)
+    accumulate_d_left!(L, dLdβ, h, β)
 end
 
 function accumulate_right!(R, h, β)
diff --git a/src/uniform_ising.jl b/src/uniform_ising.jl
new file mode 100644
index 0000000..1c71c7b
--- /dev/null
+++ b/src/uniform_ising.jl
@@ -0,0 +1,195 @@
+mutable struct UniformIsing{T<:Real, U<:OffsetVector}
+    J    :: T                  # uniform coupling strength
+    h    :: Vector{T}          # external fields
+    β    :: T                  # inverse temperature
+    L    :: OffsetVector{U, Vector{U}} # partial sums from the left
+    R    :: OffsetVector{U, Vector{U}} # partial sums from the right
+    dLdβ ::  OffsetVector{U, Vector{U}} # Derivative of L wrt β
+
+    function UniformIsing(N::Integer, J::T, h::Vector{T}, β::T=1.0) where T
+        @assert length(h) == N
+        @assert β ≥ 0
+        R = accumulate_right(h, β)
+        L, dLdβ = accumulate_d_left(h, β)
+        U = eltype(L)
+        return new{T, U}(J, h, β, L, R, dLdβ)
+    end
+end
+function UniformIsing(N::Integer, J::T, β::T=1.0)  where {T<:Real}
+    h = zeros(T, N)
+    return UniformIsing(N, J, h, β)
+end
+
+# re-compute the partial quantities needed to compute observables, in case some parameter (`J,h,β`) was modified
+function recompute_partials!(x::UniformIsing)
+    (; h, β, L, R, dLdβ) = x
+    accumulate_left!(L, h, β)
+    accumulate_right!(R, h, β)
+    accumulate_d_left!(L, dLdβ, h, β)
+end
+
+nvariables(x::UniformIsing) = length(x.h)
+variables(x::UniformIsing) = 1:nvariables(x)
+
+function energy(x::UniformIsing, σ)
+    s = sum(σ)
+    f = dot(σ, x.h)
+    N = nvariables(x)
+    return -( x.J/2*(s^2/N-1) + f ) 
+end
+
+function lognormalization(x::UniformIsing)
+    (; β, J, L) = x
+    N = nvariables(x)
+    return logsumexp( β*J/2*(s^2/N-1) + L[N][s] for s in -N:N )
+end
+function normalization(x::UniformIsing; logZ = lognormalization(x))
+    return exp(logZ)
+end
+
+pdf(x::UniformIsing, σ) = exp(-x.β*energy(x, σ) - lognormalization(x))
+
+free_energy(x::UniformIsing; logZ = lognormalization(x)) = -logZ / x.β
+
+function sample_spin(rng::AbstractRNG, p::Real)
+    @assert 0 ≤ p ≤ 1
+    r = rand(rng)
+    return r < p ? 1 : -1
+end
+
+# return a sample along with its probability
+function sample!(rng::AbstractRNG, σ, x::UniformIsing; logZ = lognormalization(x))
+    (; J, h, β, R) = x
+    N = nvariables(x)
+    a = 0.0; b = 0
+    f(s) = β*J/2*(s^2/N-1)
+    for i in 1:N
+        tmp_plus = tmp_minus = 0.0
+        for s in -N:N
+            tmp_plus += exp(f(b+1+s) + R[i+1][s])
+            tmp_minus += exp(f(b-1+s) + R[i+1][s])
+        end
+        p_plus = exp(β*h[i]) * tmp_plus
+        p_minus = exp(-β*h[i]) * tmp_minus
+        p_i = p_plus / (p_plus + p_minus)
+        σi = sample_spin(rng, p_i)
+        σ[i] = σi
+        a += h[i]*σi
+        b += σi
+    end
+    p = exp(f(b) + β*a - logZ)
+    @assert a == dot(h, σ); @assert b == sum(σ)
+    return σ, p
+end
+sample!(σ, x::UniformIsing; kw...) = sample!(default_rng(), σ, x; kw...)
+sample(rng::AbstractRNG, x::UniformIsing; kw...) = sample!(rng, zeros(Int, nvariables(x)), x; kw...)
+sample(x::UniformIsing; kw...) = sample(default_rng(), x; kw...)
+
+# first store in `p[i]` the quantity log(p(σᵢ=+1)), then transform at the end 
+function site_magnetizations!(p, x::UniformIsing; logZ = lognormalization(x))
+    (; J, h, β, L, R) = x
+    N = nvariables(x)
+    f(s) = β*J/2*(s^2/N-1)
+    for i in eachindex(p)
+        p[i] = -Inf
+        for sL in -N:N
+            for sR in -N:N
+                s = sL + sR
+                p[i] = logaddexp(p[i], f(s+1) + β*h[i] + L[i-1][sL] + R[i+1][sR])
+            end
+        end
+        # include normalization
+        p[i] = exp(p[i] - logZ)
+        # transform form p(+) to m=2p(+)-1
+        p[i] = 2*p[i] - 1
+    end
+    return p
+end
+function site_magnetizations(x::UniformIsing{T,U}; kw...) where {T,U} 
+    return site_magnetizations!(zeros(T, nvariables(x)), x; kw...)
+end
+
+# distribution of the sum of all variables as an OffsetVector
+function sum_distribution!(p, x::UniformIsing; logZ = lognormalization(x))
+    (; J, β, L) = x
+    N = nvariables(x)
+    for s in -N:N
+        p[s] = exp( β*J/2*(s^2/N-1) + L[N][s] - logZ ) 
+    end
+    return p
+end
+function sum_distribution(x::UniformIsing{T,U}; kw...) where {T,U} 
+    p = fill(zero(T), -nvariables(x):nvariables(x))
+    return sum_distribution!(p, x; kw...)
+end
+
+function avg_energy(x::UniformIsing{T}; logZ = lognormalization(x)) where T
+    (; J, β, L, dLdβ) = x
+    N = nvariables(x)
+    Zt = sum( exp( β*J/2*(s^2/N-1) + L[N][s]) * (J/2*(s^2/N-1)+dLdβ[N][s]) for s in -N:N)
+    return -exp(log(Zt) - logZ)
+end
+
+entropy(x::UniformIsing; kw...) = x.β * (avg_energy(x; kw...) - free_energy(x; kw...)) 
+
+function correlations!(m, x::UniformIsing{T,U};
+        M = accumulate_middle(x.h, x.β), logZ = lognormalization(x)) where {T,U}
+    (; J, h, β, L, R) = x
+    N = nvariables(x)
+    Z = exp(logZ)
+    f(s) = β*J/2*(s^2/N-1)
+    for i in 1:N
+        # j = i
+        m[i,i] = 1
+        # j = i+1
+        j = i + 1
+        j > N && break
+        m[i,j] = 0
+        for sL in -N:N
+            for sR in -N:N
+                s = sL + sR
+                m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
+                            exp( f(s-2) - β*(h[i]+h[j]) ) -
+                            exp( f(s)   + β*(h[i]-h[j]) ) - 
+                            exp( f(s)   - β*(h[i]-h[j]) )   ) * 
+                            exp(L[i-1][sL] + R[j+1][sR]) 
+            end
+        end
+        m[i,j] /= Z; m[j,i] = m[i,j]
+        # j > i + 1
+        for j in i+2:N
+            m[i,j] = 0
+            for sM in -N:N
+                for sL in -N:N
+                    for sR in -N:N
+                        s = sL + sM + sR 
+                        m[i,j] += ( exp( f(s+2) + β*(h[i]+h[j]) ) + 
+                                  exp( f(s-2) - β*(h[i]+h[j]) ) -
+                                  exp( f(s)   + β*(h[i]-h[j]) ) - 
+                                  exp( f(s)   - β*(h[i]-h[j]) )   ) *
+                                  exp(L[i-1][sL] + M[i+1,j-1][sM] + R[j+1][sR])
+                    end 
+                end
+            end
+            m[i,j] /= Z; m[j,i] = m[i,j]
+        end
+    end
+    return m
+end
+function correlations(x::UniformIsing{T,U}; kw...) where {T,U} 
+    correlations!(zeros(T,nvariables(x),nvariables(x)), x; kw...)
+end
+
+function covariances!(c, x::UniformIsing; logZ = lognormalization(x),
+        m = site_magnetizations(x; logZ), p = correlations(x; logZ))
+    N = nvariables(x)
+    for i in 1:N
+        for j in 1:N
+            c[i,j] = p[i,j] - m[i]*m[j]
+        end
+    end
+    return c
+end
+function covariances(x::UniformIsing{T,U}; kw...) where {T,U} 
+    covariances!(zeros(T,nvariables(x),nvariables(x)), x; kw...)
+end
\ No newline at end of file
diff --git a/test/basics.jl b/test/basics.jl
new file mode 100644
index 0000000..a8141ca
--- /dev/null
+++ b/test/basics.jl
@@ -0,0 +1,30 @@
+N = 10
+J = 0.5
+β = 2.3
+
+x = UniformIsing(N, J, β)
+
+@testset "basics" begin
+
+    @testset "outer constructor" begin
+        @test all(isequal(0), x.h)
+    end
+
+    @testset "left accumulator" begin
+        L = UniformIsingModels.accumulate_left(x.h, x.β)
+        @test L == x.L
+    end
+
+    @testset "mutate and recompute partials" begin
+        hnew = ones(N)
+        Jnew = -1.1
+        βnew = 0.1
+        x.h = hnew
+        x.J = Jnew
+        x.β = βnew
+        recompute_partials!(x)
+        xnew = UniformIsing(N, Jnew, hnew, βnew)
+        @test lognormalization(x) == lognormalization(xnew)
+    end 
+
+end
\ No newline at end of file
diff --git a/test/observables.jl b/test/observables.jl
index 5d56056..db0f507 100644
--- a/test/observables.jl
+++ b/test/observables.jl
@@ -7,10 +7,11 @@ function Obs(f::Function)
 end
 function observables_bruteforce(x::UniformIsing, 
         observables::Vector{<:Function})
-    if x.N > 10
+    N = nvariables(x)
+    if N > 10
         @warn "Exponential scaling alert"
     end
-    for s in Iterators.product(fill((-1,1),x.N)...)
+    for s in Iterators.product(fill((-1,1), N)...)
         for f! in observables  
             f!(x, s)
         end
@@ -20,59 +21,70 @@ end
 
 N = 10
 J = 0.5
-h = 1.2*randn(N)
+h = 1.2 * randn(N)
 β = 2.3
 
 x = UniformIsing(N, J, h, β)
 
-@testset "normalization" begin
-    _normaliz = (x, s) -> exp(-x.β*energy(x, s))
-    Z_bruteforce = observables_bruteforce(x, [Obs(_normaliz)])[1]
-    @test x.logZ ≈ log(Z_bruteforce)
-end
+@testset "observables" begin
 
-@testset "magnetizations" begin
-    m = site_magnetizations(x)
-    _magnetiz = [Obs((x, s) -> pdf(x, s)*s[i]) for i in 1:x.N]
-    magnetiz_bruteforce = observables_bruteforce(x, _magnetiz)
-    @test all(1:x.N) do i 
-        m[i] ≈ magnetiz_bruteforce[i]
+    @testset "normalization" begin
+        _normaliz = (x, s) -> exp(-x.β*energy(x, s))
+        Z_bruteforce = observables_bruteforce(x, [Obs(_normaliz)])[1]
+        @test normalization(x) ≈ Z_bruteforce
     end
-end
 
-@testset "pair magnetizations" begin
-    p = pair_magnetizations(x)
-    _pair_magnetiz = [Obs((x, s) -> pdf(x, s)*s[i]*s[j]) 
-                                        for i in 1:x.N for j in 1:x.N]
-    pair_magnetiz_bruteforce = observables_bruteforce(x, _pair_magnetiz)
-    @test all(Iterators.product(1:x.N,1:x.N)) do (i,j) 
-        k = Int( (j-1)*N + i )
-        p[i,j] ≈ pair_magnetiz_bruteforce[k]
+    @testset "magnetizations" begin
+        m = site_magnetizations(x)
+        _magnetiz = [Obs((x, s) -> pdf(x, s)*s[i]) for i in variables(x)]
+        magnetiz_bruteforce = observables_bruteforce(x, _magnetiz)
+        @test all(variables(x)) do i 
+            m[i] ≈ magnetiz_bruteforce[i]
+        end
     end
-end
 
-@testset "correlations" begin
-    m = site_magnetizations(x)
-    c = correlations(x)
-    _correl = [Obs((x, s) -> pdf(x, s)*(s[i]*s[j]-m[i]*m[j])) 
-                                        for i in 1:x.N for j in 1:x.N]
-    correl_bruteforce = observables_bruteforce(x, _correl)
-    @test all(Iterators.product(1:x.N,1:x.N)) do (i,j) 
-        k = Int( (j-1)*N + i )
-        isapprox( c[i,j], correl_bruteforce[k], atol=1e-4 )
+    @testset "correlations" begin
+        p = correlations(x)
+        _pair_magnetiz = [Obs((x, s) -> pdf(x, s)*s[i]*s[j]) 
+                                            for i in variables(x) for j in variables(x)]
+        pair_magnetiz_bruteforce = observables_bruteforce(x, _pair_magnetiz)
+        @test all(Iterators.product(variables(x), variables(x))) do (i,j) 
+            k = Int( (j-1)*N + i )
+            p[i,j] ≈ pair_magnetiz_bruteforce[k]
+        end
     end
-end
 
-@testset "average energy" begin
-    U = avg_energy(x)
-    _energy = Obs((x,s) -> pdf(x,s)*energy(x,s))
-    avg_energy_bruteforce = observables_bruteforce(x, [_energy])[1]
-    @test U ≈ avg_energy_bruteforce
-end
+    @testset "covariances" begin
+        m = site_magnetizations(x)
+        c = covariances(x)
+        _correl = [Obs((x, s) -> pdf(x, s)*(s[i]*s[j]-m[i]*m[j])) 
+                                            for i in variables(x) for j in variables(x)]
+        correl_bruteforce = observables_bruteforce(x, _correl)
+        @test all(Iterators.product(variables(x),variables(x))) do (i,j) 
+            k = Int( (j-1)*N + i )
+            isapprox( c[i,j], correl_bruteforce[k], atol=1e-4 )
+        end
+    end
+
+    @testset "distribution of sum" begin
+        p = sum_distribution(x)
+        _sum_distr = [Obs((x, s) -> pdf(x, s) * (sum(s) == σ)) for σ in -N:N]
+        sum_distr_bruteforce = observables_bruteforce(x, _sum_distr)
+        @test p.parent ≈ sum_distr_bruteforce
+    end
+
+    @testset "average energy" begin
+        U = avg_energy(x)
+        _energy = Obs((x,s) -> pdf(x,s)*energy(x,s))
+        avg_energy_bruteforce = observables_bruteforce(x, [_energy])[1]
+        @test U ≈ avg_energy_bruteforce
+    end
+
+    @testset "entropy" begin
+        S = entropy(x)
+        _entropy = Obs((x,s) -> -pdf(x,s)*log(pdf(x,s)))
+        entropy_bruteforce = observables_bruteforce(x, [_entropy])[1]
+        @test S ≈ entropy_bruteforce
+    end
 
-@testset "entropy" begin
-    S = entropy(x)
-    _entropy = Obs((x,s) -> -pdf(x,s)*log(pdf(x,s)))
-    entropy_bruteforce = observables_bruteforce(x, [_entropy])[1]
-    @test S ≈ entropy_bruteforce
 end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 8d9e04b..b7437ca 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using UniformIsingModels
 using Test
 
+include("basics.jl")
 include("observables.jl")
 include("sampling.jl")
\ No newline at end of file