diff --git a/Project.toml b/Project.toml
index be7c1a9..a47afec 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,6 +6,7 @@ version = "0.1.0"
 [deps]
 AngularMomentumAlgebra = "0bb1b1ac-0216-11e9-361b-bb110c659ce1"
 AtomicLevels = "10933b4c-d60f-11e8-1fc6-bd9035a249a1"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 EnergyExpressions = "f4b57a2e-27c7-11e9-0cb0-edd185b706f6"
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -15,6 +16,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SpecialPolynomials = "a25cea48-d430-424a-8ee7-0d3ad3742e9e"
 
 [compat]
+Combinatorics = "1.0"
 HypergeometricFunctions = "0.2, 0.3"
 Polynomials = "1"
 SpecialFunctions = "0.10"
diff --git a/src/Hydrogen.jl b/src/Hydrogen.jl
index fa79ac9..5bb6bcd 100644
--- a/src/Hydrogen.jl
+++ b/src/Hydrogen.jl
@@ -2,6 +2,7 @@ module Hydrogen
 
 using LinearAlgebra
 using SparseArrays
+using Combinatorics
 
 using AtomicLevels
 using AngularMomentumAlgebra
@@ -13,7 +14,6 @@ using SpecialPolynomials: Laguerre
 using SpecialFunctions: gamma
 
 include("sparse_builder.jl")
-include("factorial_ratio.jl")
 
 include("constants.jl")
 include("energies.jl")
diff --git a/src/dipoles.jl b/src/dipoles.jl
index e04c1eb..53ca16f 100644
--- a/src/dipoles.jl
+++ b/src/dipoles.jl
@@ -50,7 +50,7 @@ function radial_dipole_moment(a::Orbital, b::Orbital)
     n == n′ && return 3/2*n*√(n^2-ℓ^2)
 
     A = powneg1(n′-ℓ)/(4*factorial(2ℓ-1))
-    B = √(factorial_ratio(n+ℓ,n-ℓ-1)*factorial_ratio(n′+ℓ-1,n′-ℓ))
+    B = √(factorial(n+ℓ,n-ℓ-1)*factorial(n′+ℓ-1,n′-ℓ))
     C = (4*n*n′)^(ℓ+1)*(n-n′)^(n+n′-2ℓ-2)/((n+n′)^(n+n′))
 
     nᵣ = n-ℓ-1
diff --git a/src/factorial_ratio.jl b/src/factorial_ratio.jl
deleted file mode 100644
index ac07033..0000000
--- a/src/factorial_ratio.jl
+++ /dev/null
@@ -1,25 +0,0 @@
-@doc raw"""
-    factorial_ratio(a,b)
-
-Compute
-
-```math
-\frac{a!}{b!}
-```
-
-via simple looping. This is efficient if ``\abs{a-b}`` is not too
-large.
-"""
-function factorial_ratio(a,b)
-    a < b && return inv(factorial_ratio(b,a))
-    d = a-b
-    isinteger(d) || throw(DomainError("factorial_ratio only valid for integer distances"))
-    a == b && return one(a)
-
-    f = a
-    while a > b+1
-        a -= 1
-        f *= a
-    end
-    f
-end
diff --git a/src/orbitals.jl b/src/orbitals.jl
index 0e1bc81..b0fca4a 100644
--- a/src/orbitals.jl
+++ b/src/orbitals.jl
@@ -43,4 +43,39 @@ function non_relativistic_orbital(n::Integer, ℓ::Integer, Z = 1)
     r -> C(r) * exp(q*r)
 end
 
-export non_relativistic_orbital
+"""
+    radial_moment(k, n, ℓ[, Z = 1])
+
+Compute the radial moment ``⟨rᵏ⟩`` for a hydrogenic orbital ``(n,ℓ)``
+in a Coulomb potential of charge ``Z``. The formulas are taken from Table 2⁵ of
+
+- Condon, E. U., & Shortley, G. H. (1951). The theory of atomic
+  spectra. Cambridge: Cambridge University Press.
+
+and are available for ``k ∈ -4:4``.
+"""
+function radial_moment(k::Integer, n::Integer, ℓ::Integer, Z = 1)
+    if k == 1
+        (3n^2 - ℓ*(ℓ+1))/2Z
+    elseif k == 2
+        n^2*(5n^2 + 1 - 3ℓ*(ℓ+1))/2Z^2
+    elseif k == 3
+        n^2*(35n^2*(n^2-1) - 30n^2*(ℓ+2)*(ℓ-1) + 3*(ℓ+2)*(ℓ+1)*ℓ*(ℓ-1))/8Z^3
+    elseif k == 4
+        n^4*(63n^4 - 35n^2*(2ℓ^2+2ℓ-3) + 5ℓ*(ℓ+1)*(3ℓ^2+3ℓ-10) + 12)/8Z^4
+    elseif k == -1
+        Z/n^2
+    elseif k == -2
+        Z^2/(n^3*(ℓ+1//2))
+    elseif k == -3
+        Z^3/(n^3*(ℓ+1)*(ℓ+1//2)*ℓ)
+    elseif k == -4
+        Z^4/2*(3n^2 - ℓ*(ℓ+1))/(n^5*(ℓ+3//2)*(ℓ+1)*(ℓ+1//2)*ℓ*(ℓ-1//2))
+    elseif k == 0
+        one(Z)
+    else
+        throw(ArgumentError("Radial moment not available for k = $(k)"))
+    end
+end
+
+export non_relativistic_orbital, radial_moment
diff --git a/test/runtests.jl b/test/runtests.jl
index 3487c44..aee41e3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -31,14 +31,6 @@ function lyman_balmer²(o, n, ℓ)
 end
 
 @testset "Hydrogen.jl" begin
-    @testset "Factorial ratios" begin
-        for a = 1:6
-            for b = 1:6
-                @test Hydrogen.factorial_ratio(a,b) == factorial(a)/factorial(b)
-            end
-        end
-    end
-
     @testset "Radial dipole moments" begin
         # We cannot compare with Table 13 of Bethe and Salpeter 1977,
         # since it provides so few figures and some of the values