# Clebsch-Gordan Coefficients Implementation

In scientific programming I often come into the position of needing to implement a rather obscure and complicated formula. Clebsch-Gordan coefficients are such a thing.

I'm mainly seeking advice on how to implement such complicated formulas in a more comprehensive way but please let that not stop you from giving me any other advice you can think of.

# Code

constexpr double clebsch_gordan(double j1, double m1, double j2, double m2, double J, double M){
if( !is_angular_momentum(j1) || !is_angular_momentum(m1) || !is_angular_momentum(j2) || !is_angular_momentum(m2) || !is_angular_momentum(J) ){
throw std::domain_error("clebsch_gordan: all parameters must be multiples of 0.5.");
}
if( (is_integer(j1) && is_half_integer(m1)) || (is_half_integer(j1) && is_integer(m1)) ){
throw std::domain_error("clebsch_gordan: j1 and m1 must both be integral or half integral.");
}
if( (is_integer(j2) && is_half_integer(m2)) || (is_half_integer(j2) && is_integer(m2)) ){
throw std::domain_error("clebsch_gordan: j2 and m2 must both be integral or half integral.");
}
if( (is_integer(j1 + j2) && !is_integer(J)) || (is_half_integer(j1 + j2) && !is_half_integer(J)) ){
throw std::domain_error("clebsch_gordan: J is no valid value for j1 and j2 combination.");
}
if( j1 < 0 || j2 < 0 ){
throw std::domain_error("clebsch_gordan: j1 and j2 must be non-negative.");
}
if( J > j1 + j2 || J < abs(j1 - j2) ){
throw std::domain_error("clebsch_gordan: J must be |j1-j2| <= J <= j1+j2.");
}
if( abs(m1) > j1 || abs(m2) > j2 || abs(M) > J ){
throw std::domain_error("clebsch_gordan: m1 and m2 must be |m1| <= j1 and |m2| <= j2 and |m1+m2| <= J");
}
if( m1 + m2 != M ){
throw std::domain_error("clebsch_gordan: m1 + m2 != M.");
}

auto const & f = constexpr_factorial;
double numerator = (2 * J + 1) * f(J + j1 - j2) * f(J - j1 + j2) * f(j1 + j2 - J);
numerator *= f(J + M) * f(J - M) * f(j1 - m1) * f(j1 + m1) * f(j2 - m2) * f(j2 + m2);
double denominator = f(j1 + j2 + J + 1);
int const min = std::max(0., std::max(j2 - J - m1, j1 + m2 - J));
int const max = std::min(j2 + m2, std::min(j1 - m1, j1 + j2 - J));
double sum = min > max ? 1 : 0;
for( int k = min; k <= max; ++k ){
sum += constexpr_pow(-1., k) / (f(k) * f(j1 + j2 - J - k) * f(j1 - m1 - k) * f(j2 + m2 - k) * f(J - j2 + m1 + k) * f(J - j1 - m2 + k));
}
return constexpr_sqrt(numerator / denominator) * sum;
}


# Helper functions

constexpr bool is_angular_momentum(double x)
{
return static_cast<int>(2*x) - 2*x == 0;
}

constexpr bool is_integer(double x)
{
return static_cast<int>(x) - x == 0;
}

constexpr bool is_half_integer(double x)
{
return !is_integer(x) && is_integer(2*x);
}

template<typename T>
constexpr bool is_almost_equal(T a, T b, T rel_epsilon = 1.e-15)
{
if( std::isnan(a) || std::isnan(b) )
{
return false;
}
if( a == b )
{
return true;
}
return abs(a-b) < std::abs(max(a, b)) * rel_epsilon;
}

constexpr double constexpr_sqrt(double x, double curr, double prev)
{
return curr == prev
? curr
: constexpr_sqrt(x, 0.5 * (curr + x / curr), curr);
}

constexpr double constexpr_sqrt(double x)
{
return x >= 0 && x < std::numeric_limits<double>::infinity()
? constexpr_sqrt(x, x, 0)
: std::numeric_limits<double>::quiet_NaN();
}

constexpr std::size_t constexpr_binomial(std::size_t n, std::size_t k)
{
if (k == 0 || n == k) {
return 1;
}
else {
return constexpr_binomial(n - 1, k - 1) + constexpr_binomial(n - 1, k);
}
}
constexpr std::size_t constexpr_factorial(std::size_t n)
{
return n == 0? 1 : n * constexpr_factorial(n-1);
}


I realise the is_almost_equal function is flawed as it is posted right now.

# Test cases (Framework: catch2)

TEST_CASE("clebsch_gordan", "[math]"){
REQUIRE(is_almost_equal(clebsch_gordan(0,0,0,0,0,0), 1.));
REQUIRE(is_almost_equal(clebsch_gordan(2,0,3,0,5,0), constexpr_sqrt(10./21.)));
REQUIRE(is_almost_equal(clebsch_gordan(2,1,3,0,5,1), 2*constexpr_sqrt(2./21.)));
REQUIRE(is_almost_equal(clebsch_gordan(2,0,3,1,5,1), constexpr_sqrt(3./7.)));
REQUIRE(is_almost_equal(clebsch_gordan(2,1,3,1,5,2), constexpr_sqrt(1./2.)));
REQUIRE(is_almost_equal(clebsch_gordan(0.5, 0.5, 0.5, -0.5, 1, 0), constexpr_sqrt(1./2.)));
REQUIRE(is_almost_equal(clebsch_gordan(0.5, 0.5, 0.5, -0.5, 0, 0), constexpr_sqrt(1./2.)));
REQUIRE(is_almost_equal(clebsch_gordan(2.5,0.5,1.5,0.5,2,1), -5./(2*constexpr_sqrt(21.))));
// exceptions
REQUIRE([](){
try{ clebsch_gordan(-1,0,0,0,0,0);return false; }
catch( std::domain_error& e ){ return true; }}()
);
REQUIRE([](){
try{ clebsch_gordan(1,0,1,0,3,0);return false; }
catch( std::domain_error& e ){ return true; }}()
);
REQUIRE([](){
try{ clebsch_gordan(1,0,2,0,3,1);return false; }
catch( std::domain_error& e ){ return true; }}()
);
}

• If all parameters must be exact multiples of 0.5, is it possible to just accept integers with an implicit factor? Jan 31, 2022 at 9:38
• @Reinderien in principle yes, that would be possible, but that would mean to refactor existing code in many places and might obscure the meaning of the spins (i.e. a 1 could be mistaken to be spin 1 instead of 2 * 1/2. Jan 31, 2022 at 10:17
• Is abs() equal to std::abs() here? Jan 31, 2022 at 10:50
• @TobySpeight basically yes. Jan 31, 2022 at 11:03

My primary concern is around accuracy and overflow. Given that:

• your factorial is implemented in integral math
• even though your figures are on the order of 1 through ~5 if Wikipedia is to be believed, you still use them in a product of nine (!!) factorials in your numerator

I would immediately be concerned about overflow. It only takes a factorial exceeding 20-bang to overflow a 64-bit size_t:

$$\log_2 21! = \sum _{i=1} ^{21} \log_2 i \approx 65$$

If you're lucky your compiler will yell at you. If you're unlucky it will silently produce a nonsense value, since integer overflow is undefined behaviour.

Typically when mathematicians write an expression with factorials in the numerator and denominator, that's not actually how they should be implemented verbatim. Instead, you should endeavour to identify overlaps, do the cancellation, and rewrite the formula in terms of pi-expressions.

As an example let's assume for the moment that

$$j_1 + j_2 = 3$$ $$J = 1$$

Then your denominator is

$$(j_1 + j_2 + J + 1)! = 5!$$

and one of the terms in the numerator will be

$$(j_1 + j_2 - J)! = 2!$$

but you shouldn't compute them verbatim, instead (showing the quotient of only those two terms) you would compute

$$\prod _{i=3} ^5 \frac 1 i$$

Another strategy to reduce the possibility of overflow: since your numerator is being run through a square root, any time that you see two overlapping factorial segments in your product, eliminate one and pull the other outside of the root.

Your sum denominator also risks overflowing. Instead of dividing by a huge product, consider successive division by each term.

I'm not going to bother puzzling out how to express this in constexpr C++, but a reasonable algorithm sketched out in Python could look like

from collections import defaultdict
from functools import reduce
from math import sqrt, isclose
from operator import mul

def product(sequence):
return reduce(mul, sequence, 1)

def f(k):
return product(range(1, int(k)+1))

def clebsch_gordan(j1, m1, j2, m2, J, M):
# The numerator arguments to factorial
num_args = [int(k) for k in (
J + j1 - j2,
J - j1 + j2,
j1 + j2 - J,
J + M,
J - M,
j1 - m1,
j1 + m1,
j2 - m2,
j2 + m2,
)]

# All factors, starting by adding from numerator
facs = defaultdict(int)
for k in num_args:
for i in range(2, k + 1):
facs[i] += 1

# Subtract factors from denominator
den_arg = int(j1 + j2 + J + 1)
for i in range(2, den_arg + 1):
facs[i] -= 1

# Subtract two from any orders within the square root that need it
outside_root = 1
for x, n in facs.items():
if n >= 2:
quotient, facs[x] = divmod(n, 2)
outside_root *= x**quotient

# Factorial-equivalent expression replacing the old numerator, denominator
# and applying the square root
radicand = (2*J + 1) * product(k**n for k, n in facs.items())

jmin = int(max(0, j2 - J - m1, j1 + m2 - J))
jmax = int(min(j2 + m2, j1 - m1, j1 + j2 - J))
jsum = 1 if jmin > jmax else 0
sign = -1 if jmin % 2 else 1

for k in range(jmin, jmax + 1):
jsum += (
sign
/ f(k)
/ f(j1 + j2 - J - k)
/ f(j1 - m1 - k)
/ f(j2 + m2 - k)
/ f(J - j2 + m1 + k)
/ f(J - j1 - m2 + k)
)
sign = -sign

fac_expr = outside_root * sqrt(radicand)

# For comparison; delete me
old_numerator = product(k for n in num_args for k in range(2, n+1))
old_denominator = product(range(2, den_arg+1))
old_fac_expr = sqrt(old_numerator * (2*J + 1) / old_denominator)
print(f'Old: {old_fac_expr:.3f} = sqrt({2*J + 1} * {old_numerator} / {old_denominator})')
new_numerator = product(x**n for x, n in facs.items() if n > 0)
new_denominator = product(x**-n for x, n in facs.items() if n < 0)
print(f'New: {fac_expr:.3f} = {outside_root} * sqrt({2*J + 1} * {new_numerator} / {new_denominator})\n')

return fac_expr * jsum

assert isclose(clebsch_gordan(0, 0, 0, 0, 0, 0), 1)
assert isclose(clebsch_gordan(2, 0, 3, 0, 5, 0), sqrt(10 / 21))
assert isclose(clebsch_gordan(2, 1, 3, 0, 5, 1), 2 * sqrt(2 / 21))
assert isclose(clebsch_gordan(2, 0, 3, 1, 5, 1), sqrt(3 / 7))
assert isclose(clebsch_gordan(2, 1, 3, 1, 5, 2), sqrt(0.5))
assert isclose(clebsch_gordan(0.5, 0.5, 0.5, -0.5, 1, 0), sqrt(0.5))
assert isclose(clebsch_gordan(0.5, 0.5, 0.5, -0.5, 0, 0), sqrt(0.5))
assert isclose(clebsch_gordan(2.5, 0.5, 1.5, 0.5, 2, 1), -5 / (2 * sqrt(21)))


Output:

Old: 1.000 = sqrt(1 * 1 / 1)
New: 1.000 = 1 * sqrt(1 * 1 / 1)

Old: 99.369 = sqrt(11 * 35831808000 / 39916800)
New: 99.369 = 1440 * sqrt(11 * 24 / 55440)

Old: 133.318 = sqrt(11 * 64497254400 / 39916800)
New: 133.318 = 864 * sqrt(11 * 120 / 55440)

Old: 125.694 = sqrt(11 * 57330892800 / 39916800)
New: 125.694 = 1152 * sqrt(11 * 60 / 55440)

Old: 203.647 = sqrt(11 * 150493593600 / 39916800)
New: 203.647 = 288 * sqrt(11 * 360 / 7920)

Old: 0.707 = sqrt(3 * 1 / 6)
New: 0.707 = 1 * sqrt(3 * 1 / 6)

Old: 0.707 = sqrt(1 * 1 / 2)
New: 0.707 = 1 * sqrt(1 * 1 / 2)

Old: 1.309 = sqrt(5 * 1728 / 5040)
New: 1.309 = 12 * sqrt(5 * 2 / 840)

• So if I understood correctly possible improvements would be: Use double instead of int for factorials. Also possibly storing the factorials in a vector, pre-computation to try to do some prior manipulation like cancelation. Jan 31, 2022 at 14:25
• Using double replaces one problem with another: instead of overflow risk, you'll have accuracy loss risk. Ideally, keep your ints, but do cancellation - not necessarily of the kind you've described with vectors, though I guess that's possible; instead do analysis of the ranges of your terms, identify the largest factorial numerator term, and replace that term and your denominator with a pi-expression that uses the remaining range as in my example Jan 31, 2022 at 14:28
• How would I do such an analysis though, if not with vectors? Jan 31, 2022 at 14:30
• You have a fixed number of numerators - find the one whose factorial argument is highest. This can just be a series of nested max calls. Jan 31, 2022 at 14:34
• Seems I was a bit slow on the uptake. Sense, that makes perfect sense! Jan 31, 2022 at 14:36

It may be worth introducing a new type to represent spin values; this can store twice the spin as an integer, and convert to double as necessary by dividing by two. That simplifies the is_integral() and is_half_integral() functions (which can be members of the Spin class). A similar class for angular momentum would also simplify the code.

A minimal version looks something like

// binary fixed-point quantity for multiples of ½
class semi_integer
{
int x2;                     // twice the actual value

public:
constexpr explicit semi_integer(double x)
: x2{static_cast<int>(2*x)}
{
if (std::fmod(x, 0.5)) {
throw std::domain_error("Spin value must be a multiple of 0.5");
}
}

constexpr semi_integer(int x)
: x2{2*x}
{}

constexpr operator double() const {
return x2 / 2.0;
}

constexpr bool is_integer() const {
return x2 % 2 == 0;
}
constexpr bool is_half_integer() const {
return !is_integer();
}

constexpr semi_integer& operator+=(semi_integer other) {
x2 += other.x2;
return *this;
}
};

class spin : public semi_integer
{
public:
constexpr spin(double d)
: semi_integer{d}
{}
constexpr spin operator+(spin other) {
return static_cast<spin&>(other += *this);
}
};

class angular_momentum : public semi_integer
{
public:
constexpr angular_momentum(double d)
: semi_integer{d}
{}
};


We'll use it like this:

constexpr double clebsch_gordan(spin j1, angular_momentum m1,
spin j2, angular_momentum m2,
spin J,  angular_momentum M)
{
if (j1.is_integer() != m1.is_integer()) {
throw std::domain_error("clebsch_gordan: j1 and m1 must both be integral or both half integral.");
}
if (j2.is_integer() != m2.is_integer()) {
throw std::domain_error("clebsch_gordan: j2 and m2 must both be integral or both half integral.");
}
if ((j1 + j2).is_integer() != J.is_integer()) {
throw std::domain_error("clebsch_gordan: J is no valid value for j1 and j2 combination.");
}
if (j1 < 0 || j2 < 0) {
throw std::domain_error("clebsch_gordan: j1 and j2 must be non-negative.");
}
if (J > j1 + j2 || J < std::abs(j1 - j2)) {
throw std::domain_error("clebsch_gordan: J must be |j1-j2| <= J <= j1+j2.");
}
if (std::abs(m1) > j1 || std::abs(m2) > j2 || std::abs(M) > J) {
throw std::domain_error("clebsch_gordan: m1 and m2 must be |m1| <= j1 and |m2| <= j2 and |m1+m2| <= J");
}
if (m1 + m2 != M) {
throw std::domain_error("clebsch_gordan: m1 + m2 != M.");
}


We no longer need to check that the values are multiples of ½, because the spin and angular_momentum constructors enforce that.

Other things - we seem to be missing some necessary headers - it helps reviewers if the code is complete.

Also, we seem to be missing using std::abs; somewhere (or we're including the deprecated C header <math.h>) - prefer using std::abs normally.

constexpr_sqrt() could have its helper function as a local. It should also use the standard implementation at run-time.

constexpr double constexpr_sqrt(double x)
{
if (std::is_constant_evaluated() && x >= 0 && std::isfinite(x)) {
const std::function<double(double,double,double)>& f =
[&f](double x, double curr, double prev){
return curr == prev
? curr
: f(x, 0.5 * (curr + x / curr), curr);
};
return f(x, x, 0);
} else {
return std::sqrt(x);
}
}


It's not clear why we have to pass M as a parameter, given that we throw an exception if it's not equal to the sum of m1 and m2. Why not just compute it from those two values?

• I followed the implementation of Mathematica, which also has the M parameter present. It might serve as a fail safe, if someone tries to accidentally compute a wrong combination. I actually have a Angular_Momentum class, so it's a nice idea to implement this function with that class. Do you have any comment about the actual implementation itself? Jan 31, 2022 at 13:06