diff --git a/math.jl b/math.jl
index fb55ec9..4a3fec5 100644
--- a/math.jl
+++ b/math.jl
@@ -22,13 +22,8 @@ gauss_int(a, n) = double_factorial(n - 1) / (2 * a)^((n + 1)/2) * (iseven(n) ? s
 
 "Gives ∫dp p u' u where u' and u' may have different l (Ref: worked out in Mathematica)"
 function integral(np, lp, n, l, μω)
-    s = 0.0
-    for m in 0 : n
-        for mp in 0 : np
-            s += (-1)^(m + mp) * gauss_int(1, 2 * m + 2 * mp + l + lp + 3) / (factorial(n - m) * gamma(l + 3/2 + m) * factorial(m) * factorial(np - mp) * gamma(lp + 3/2 + mp) * factorial(mp))
-        end
-    end
-    return 2 * sqrt(μω * factorial(n) * factorial(np) * gamma(n + l + 3/2) * gamma(np + lp + 3/2)) * s
+    s = [(-1)^(m + mp) * gauss_int(1, 2 * m + 2 * mp + l + lp + 3) * N_lnk(l, n, m) * N_lnk(lp, np, mp) for (m, mp) in Iterators.product(0:n, 0:np)]
+    return 2 * sqrt(μω) * better_sum(s)
 end
 
 "Gives <n' l'|| p ||n l> for the HO basis"
diff --git a/test/racahs_reduction_formula.jl b/test/racahs_reduction_formula.jl
index 3bb87dd..1118833 100644
--- a/test/racahs_reduction_formula.jl
+++ b/test/racahs_reduction_formula.jl
@@ -9,7 +9,7 @@ function explicit_integral(np, lp, n, l, μω)
     return integral
 end
 
-n_list = 0 : 5
+n_list = 0 : 10
 l_list = 0 : 5
 μω_list = [1, 2.66]