diff --git a/helper.jl b/helper.jl
new file mode 100644
index 0000000..420c819
--- /dev/null
+++ b/helper.jl
@@ -0,0 +1,10 @@
+"Sum over array while minimizing catastrophic cancellation as much as possible"
+function better_sum(arr::Array{Float64})
+    pos_arr = arr[arr .> 0]
+    neg_arr = arr[arr .< 0]
+
+    sort!(pos_arr)
+    sort!(neg_arr, rev=true)
+
+    return sum(pos_arr) + sum(neg_arr)
+end
\ No newline at end of file
diff --git a/ho_basis.jl b/ho_basis.jl
index 94f8547..ce9968e 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,12 +1,13 @@
 using NuclearToolkit
 using SpecialFunctions
+include("helper.jl")
 
 # Gaussian potentials in HO space
 inv_factorial(n) = Iterators.prod(inv.(1:n))
 sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1))
 N_lnk(l, n, k) = (-1)^(n-k) * inv_factorial(k) * binomial(n+l+1/2, n-k) * sqrt_factorial(n) / sqrt(gamma(n+l+3/2))
 Talmi(l, R, k1, k2) = (1 + 1/R^2)^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2) = sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+V_Gaussian(R, l, n1, n2) = better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
 
 function get_sp_basis(E_max)
     Es = Int[]