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 f4a66d9..79a447d 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,10 +1,11 @@
 using NuclearToolkit
 using SpecialFunctions
+include("helper.jl")
 
 # Gaussian potentials in HO space
 log_N(l, n, k) = -logfactorial(n)/2 + logabsbinomial(n, k)[1] + loggamma(n + l + 3/2)/2 - loggamma(k + l + 3/2)
 log_Talmi(l, R, k1, k2) = -(3/2 + l + k1 + k2) * log(1 + 1/R^2) + loggamma(3/2 + l + k1 + k2)
-V_Gaussian(R, l, n1, n2) = (-1)^(n1 + n2) * sum([(-1)^(k1 + k2) * exp(log_N(l, n1, k1) + log_N(l, n2, k2) + log_Talmi(l, R, k1, k2)) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
+V_Gaussian(R, l, n1, n2) = (-1)^(n1 + n2) * better_sum([(-1)^(k1 + k2) * exp(log_N(l, n1, k1) + log_N(l, n2, k2) + log_Talmi(l, R, k1, k2)) for (k1, k2) in Iterators.product(0:n1, 0:n2)])
 
 function get_sp_basis(E_max)
     Es = Int[]