libwignernj

Exact evaluation of Wigner 3j, 6j, and 9j symbols, Clebsch-Gordan coefficients, Racah W-coefficients, Fano X-coefficients, and Gaunt coefficients in C99.

All intermediate arithmetic is exact integer arithmetic using a prime-factorization representation; floating-point conversion happens only at the final step. Results are accurate to the last bit of the chosen output precision.

Algorithm: prime factorization of factorials (Dodds & Wiechers, Comput. Phys. Commun. 4, 268, 1972; doi:10.1016/0010-4655(72)90019-7; refined in subsequent work) combined with the multiword-integer Racah sum of Johansson & Forssén, SIAM J. Sci. Comput. 38(1), A376–A384, 2016 (doi:10.1137/15M1021908).

Language interfaces: C (primary), C++ (header-only wrapper, links against libwignernj), Python (CPython extension), Fortran 90 (iso_c_binding).

Argument convention

Every coupling-coefficient routine (3j, 6j, 9j, Clebsch-Gordan, Racah W, Fano X, Gaunt, real-Gaunt) takes its angular-momentum arguments as twice their value so that half-integers are represented exactly as odd integers:

j = 3/2  →  tj = 3
m = -1/2 →  tm = -1

This applies to the C, C++, and Fortran interfaces. The Python interface also accepts plain floats (e.g. 0.5) and fractions.Fraction objects.

The single exception is wignernj_real_ylm_in_complex_ylm, whose argument is always an integer orbital angular momentum and is therefore taken as plain l rather than 2*l.

Building

cmake -B build && cmake --build build
ctest --test-dir build

CMake options:

Option	Default	Description
BUILD_SHARED_LIBS	ON	Shared library
WIGNERNJ_BUILD_FORTRAN	ON	Fortran interface
WIGNERNJ_BUILD_TESTS	ON	C/Fortran test suite
WIGNERNJ_BUILD_CXX_TESTS	ON	C++ header tests
WIGNERNJ_BUILD_EXAMPLES	ON	Build and run the language-binding example programs as ctest tests
WIGNERNJ_BUILD_LTO	ON	Link-time optimisation; auto-disabled if the toolchain does not support it (and on MSVC, where it conflicts with WINDOWS_EXPORT_ALL_SYMBOLS)
WIGNERNJ_BUILD_PYTHON	OFF	Python extension (dynamically linked against libwignernj)
WIGNERNJ_BUILD_QUADMATH	OFF	libquadmath / IEEE 754 binary128 (__float128) interface
WIGNERNJ_BUILD_MPFR	OFF	MPFR arbitrary-precision interface
WIGNERNJ_BUILD_FLINT	OFF	Use FLINT/GMP/MPFR for the bigint backend (instead of the in-tree schoolbook with Karatsuba) — sub-quadratic multiplication at large j via Toom-Cook / Schönhage--Strassen
WIGNERNJ_BUILD_COVERAGE	OFF	Build with --coverage -O0 for lcov / Codecov (gcc, clang)

A separate preprocessor switch -DBIGINT_FORCE_PORTABLE (passed via CMAKE_C_FLAGS) forces the multiword-integer back-end onto its pure-C99 fallback path even on compilers that support __uint128_t; this is exercised in the CI matrix to verify that the native and fallback code paths produce bit-identical output.

C API

#include "wignernj.h"

/* Wigner 3j:  ( j1  j2  j3 ) */
/*             ( m1  m2  m3 ) */
double      wigner3j  (int tj1, int tj2, int tj3, int tm1, int tm2, int tm3);
float       wigner3j_f(int tj1, int tj2, int tj3, int tm1, int tm2, int tm3);
long double wigner3j_l(int tj1, int tj2, int tj3, int tm1, int tm2, int tm3);

/* Wigner 6j:  { j1 j2 j3 } */
/*             { j4 j5 j6 } */
double      wigner6j  (int tj1, int tj2, int tj3, int tj4, int tj5, int tj6);
float       wigner6j_f(int tj1, int tj2, int tj3, int tj4, int tj5, int tj6);
long double wigner6j_l(int tj1, int tj2, int tj3, int tj4, int tj5, int tj6);

/* Wigner 9j (row-major order) */
double      wigner9j  (int tj11, int tj12, int tj13,
                       int tj21, int tj22, int tj23,
                       int tj31, int tj32, int tj33);

/* Clebsch-Gordan:  <j1 m1; j2 m2 | J M> */
double      clebsch_gordan  (int tj1, int tm1, int tj2, int tm2, int tJ, int tM);

/* Racah W-coefficient:  W(j1 j2 J j3; j12 j23) */
double      racah_w  (int tj1, int tj2, int tJ, int tj3, int tj12, int tj23);

/* Fano X-coefficient:  X(j1 j2 j12; j3 j4 j34; j13 j24 J)
 *   = sqrt[(2j12+1)(2j34+1)(2j13+1)(2j24+1)] * {9j} */
double      fano_x   (int tj1, int tj2, int tj12,
                      int tj3, int tj4, int tj34,
                      int tj13, int tj24, int tJ);

/* Gaunt coefficient:  integral Y_{l1}^{m1} Y_{l2}^{m2} Y_{l3}^{m3} dΩ */
double      gaunt  (int tl1, int tm1, int tl2, int tm2, int tl3, int tm3);

/* Real-spherical-harmonic Gaunt coefficient (Wikipedia/Condon-Shortley) */
double      gaunt_real(int tl1, int tm1, int tl2, int tm2, int tl3, int tm3);

/* Real Y_lm in the complex Y_lm basis:
 *   S_{l,m_r} = sum_{m_c} C[m_r, m_c] Y_l^{m_c}
 * fills a column-major (2l+1)x(2l+1) complex matrix.  The element
 * type `wignernj_cdouble_t` is a typedef shim in wignernj.h that
 * maps to the native complex-double type on every supported
 * toolchain: `double _Complex` on gcc/clang/Intel, `_Dcomplex` on
 * MSVC, and a layout-compat `struct {double _pair[2];}` in C++.
 * Callers holding `double _Complex *` or `std::complex<double> *`
 * (via the C++ wrapper) pass them directly without a cast. */
void        wignernj_real_ylm_in_complex_ylm(int l, wignernj_cdouble_t *C_out);

Each function is available in three precisions: double (no suffix), float (_f), and long double (_l). When the library is built with -DWIGNERNJ_BUILD_QUADMATH=ON, an additional _q variant returning __float128 is exposed in wignernj_quadmath.h (see below). Functions return 0 for symbols that vanish by selection rules; selection-rule violations are not errors.

Linking: pkg-config --libs libwignernj or -lwignernj -lm.

libquadmath API

Build with -DWIGNERNJ_BUILD_QUADMATH=ON (requires a compiler with __float128 support: GCC, Clang, or Intel ICC/ICX on Linux/macOS; not Apple Clang or MSVC). Include wignernj_quadmath.h in addition to wignernj.h. Each coupling-coefficient routine gains a _q variant returning __float128 (IEEE 754 binary128, 113-bit mantissa); wignernj_real_ylm_in_complex_ylm gains a _q variant that fills a wignernj_cfloat128_t matrix.

#include "wignernj.h"
#include "wignernj_quadmath.h"

__float128 v = wigner6j_q(4, 4, 4, 4, 4, 4);

The Fortran module also exposes the corresponding wigner3j_q, wigner6j_q, ..., gaunt_real_q bind(c) interfaces and the real-valued convenience wrappers w3jq, w6jq, w9jq, wcgq, wracahwq, wgauntq, wgaunt_realq that take real(real128) arguments, plus wreal_ylm_in_complex_ylmq(l, c_out) filling a complex(c_float128_complex) matrix.

C++ callers get the same wignernj::symbol3j<T>() / cg<T>() / ... template surface at __float128 by including wignernj_quadmath.hpp in addition to (or instead of) wignernj.hpp:

#include "wignernj_quadmath.hpp"
__float128 v = wignernj::symbol3j<__float128>(2, 2, 0, 0, 0, 0);

MPFR API

Build with -DWIGNERNJ_BUILD_MPFR=ON (requires libmpfr). Include wignernj_mpfr.h in addition to wignernj.h. Set the output precision on rop via mpfr_init2 before calling; the rounding mode is the last argument and may be any of the standard MPFR modes (MPFR_RNDN, MPFR_RNDZ, MPFR_RNDD, MPFR_RNDU, MPFR_RNDA).

#include "wignernj.h"
#include "wignernj_mpfr.h"

mpfr_t v;
mpfr_init2(v, 256);   /* 256-bit precision */

wigner3j_mpfr(v, 4, 4, 0,  2, -2, 0,  MPFR_RNDN);
wigner6j_mpfr(v, 2, 2, 0,  2,  2, 0,  MPFR_RNDN);
wigner9j_mpfr(v, 2, 2, 2,  2, 2, 2,  2, 2, 2,  MPFR_RNDN);
clebsch_gordan_mpfr(v, 2, 2,  2, -2,  4, 0,  MPFR_RNDN);
racah_w_mpfr(v, 2, 2, 4,  2,  4, 4,  MPFR_RNDN);
gaunt_mpfr(v, 4, 2,  4, -2,  4, 0,  MPFR_RNDN);
gaunt_real_mpfr(v, 4, 2,  4, -2,  0, 0,  MPFR_RNDN);

mpfr_clear(v);

The real ↔ complex Y_lm basis-overlap matrix uses a different signature because its output is a (2l+1) × (2l+1) complex matrix rather than a scalar: wignernj_real_ylm_in_complex_ylm_mpfr(l, C_re, C_im, rnd) fills two parallel mpfr_t arrays of length (2l+1)² (real and imaginary parts) that the caller mpfr_init2's and mpfr_clear's.

Link with -lwignernj -lmpfr -lm.

C++ API

The header-only wrapper wignernj.hpp provides a template interface that accepts either 2*j integers or real-valued doubles (half-integers). The wrapper has no separate translation unit, but every function forwards to a C symbol in libwignernj, so you still need to link the C library (-lwignernj -lm):

#include "wignernj.hpp"

// Integer 2*j form
double v = wignernj::symbol3j<double>(2, 2, 0,  0, 0, 0);
float  f = wignernj::symbol6j<float> (2, 2, 2,  2, 2, 2);

// Real-valued form (throws std::invalid_argument if not a half-integer)
double v = wignernj::symbol3j(1.0, 1.0, 0.0,  0.0, 0.0, 0.0);
double c = wignernj::cg(0.5, 0.5, 0.5, -0.5, 1.0, 0.0);

// Available functions: symbol3j, symbol6j, symbol9j, cg, racahw, fanox, gaunt,
// gauntreal, real_ylm_in_complex_ylm

Link with -lwignernj -lm (and -lmpfr if WIGNERNJ_BUILD_MPFR=ON).

Python API

pip install wignernj                  # from PyPI (pre-built wheels for
                                      # Linux x86_64/aarch64, macOS arm64,
                                      # and Windows x86_64; sdist fallback
                                      # for any other target)
pip install -e .                      # or build from a checkout

import wignernj

wignernj.wigner3j(1, 1, 0, 0, 0, 0)             # integer 2*j form
wignernj.wigner3j(0.5, 0.5, 1, 0.5, -0.5, 0)   # real half-integer form
wignernj.wigner6j(1, 1, 2, 1, 1, 2)
wignernj.wigner9j(1, 1, 2, 1, 1, 2, 2, 2, 4)
wignernj.clebsch_gordan(1, 1, 1, -1, 2, 0)
wignernj.racah_w(1, 1, 2, 1, 2, 2)
wignernj.fano_x(1, 1, 2,  1, 1, 2,  2, 2, 4)
wignernj.gaunt(2, 1, 2, -1, 2, 0, precision='longdouble')
wignernj.gaunt_real(2, 1, 2, -1, 0, 0)
wignernj.real_ylm_in_complex_ylm(1)                  # 3x3 unitary real->complex Y matrix

The optional precision= keyword selects 'float', 'double' (default), or 'longdouble'. Arguments may be integers, floats, or fractions.Fraction.

Fortran API

The wignernj module provides real-valued wrappers w3j, w6j, w9j, wcg, wracahw, wfanox, wgaunt, and wgaunt_real that accept double-precision real arguments, plus a typed real-to-complex-Y subroutine wreal_ylm_in_complex_ylm(l, c_out) that fills a complex(c_double_complex) (2l+1) x (2l+1) matrix:

use wignernj
use iso_c_binding, only: c_double_complex
real(8) :: v
complex(c_double_complex), allocatable :: C(:,:)
v = w3j(1.0d0, 1.0d0, 0.0d0,  0.0d0, 0.0d0, 0.0d0)
v = w6j(1.0d0, 1.0d0, 2.0d0,  1.0d0, 1.0d0, 2.0d0)
v = wcg(0.5d0, 0.5d0, 0.5d0, -0.5d0, 1.0d0, 0.0d0)
v = wfanox(1.0d0, 1.0d0, 2.0d0,  1.0d0, 1.0d0, 2.0d0,  2.0d0, 2.0d0, 4.0d0)
v = wgaunt(2.0d0, 1.0d0, 2.0d0, -1.0d0, 2.0d0, 0.0d0)
v = wgaunt_real(2.0d0, 1.0d0, 2.0d0, -1.0d0, 0.0d0, 0.0d0)
allocate(C(3, 3))
call wreal_ylm_in_complex_ylm(1, C)

Raw bind(c) interfaces using 2*j integers are also available for all three precisions (the long double variants require gfortran or Cray Fortran), plus the _q (real(real128)) variants when the library is built with WIGNERNJ_BUILD_QUADMATH=ON.

Link with -lwignernj_f03 -lwignernj -lm.

Limits

The prime list and its inverse-lookup index are hard-coded into the compiled library (≈49 kB of read-only data), which is what lets the library work with no caller-side initialization step. The default sieve limit was chosen on a rule of thumb that ~50 kB is a reasonable upper bound for compile-time constant tables; the resulting prime table covers factorials up to 20020!, which translates to:

• 3j / 6j / CG / Racah W / complex Gaunt / real Gaunt: j1+j2+j3 ≤ 20019 (equal-j: j ≤ 6673)
• 9j / Fano X: equal-j j ≤ 5004 (k-dependent triangle denominators reach (4j+1)!; Fano X delegates to the 9j pipeline)

Exceeding these limits prints a diagnostic to stderr and aborts. The ceiling is not architectural: it is set by the size of the compile-time prime table (PRIME_SIEVE_LIMIT in the auto-generated src/prime_table_macros.h) and can be raised by regenerating the table with tools/gen_prime_table.py and rebuilding, at the cost of a proportionally larger compiled-in table. MAX_FACTORIAL_ARG in the same header is derived from the sieve limit, so contributors do not need to keep the two values in sync by hand.

The 9j is also O(j⁴) in computation time; evaluations with j > a few hundred can be slow. See docs/reference.md for details.

Documentation

Full API reference with mathematical definitions, selection rules, and per-language examples: docs/reference.md.

Citing libwignernj

If libwignernj contributes to published work, please cite:

S. Lehtola, libwignernj: a reusable C/C++/Fortran/Python library for exact Wigner symbols and related coefficients, arXiv:2605.06634 (2026). doi:10.48550/arXiv.2605.06634.

@article{Lehtola2026libwignernj,
  author  = {Lehtola, Susi},
  title   = {{libwignernj}: a reusable {C}/{C}++/{F}ortran/{P}ython
             library for exact {W}igner symbols and related coefficients},
  journal = {arXiv preprint},
  year    = {2026},
  eprint  = {2605.06634},
  archivePrefix = {arXiv},
  doi     = {10.48550/arXiv.2605.06634}
}

A machine-readable CITATION.cff is provided at the repository root; GitHub's "Cite this repository" button reads it.

Repository layout

• include/, src/ — C99 core library and language wrappers
• tests/ — C/C++/Fortran unit tests, OOM-injection harness, symmetry oracles; Python tests under tests/python/
• tests/gen_refs.py — regenerates the sympy-based reference tables consumed by tests/test_3j.c, test_6j.c, etc.
• tests/cmake_downstream/ — minimal out-of-tree project demonstrating consumption via find_package(wignernj REQUIRED COMPONENTS Fortran)
• benchmarks/ — library-versioning microbenches and profile drivers (bench_term_cache.c, bench_sweep.c, bench_mul.c, bench_div128.c, profile_3j_4000.c, …) used to validate libwignernj changes against prior versions; see benchmarks/README.md. Comparative benchmarks against external libraries (WIGXJPF, GSL) are published with the paper as supplementary material rather than living in this repo.
• examples/ — single-file demonstrations of every public symbol in C, C++, Fortran, and Python (all_symbols.*) plus focused per-topic demos that mirror across all four bindings (e.g. real_basis_lz.* building the orbital l_z matrix in the real-Y basis); all built and run as ctest tests when WIGNERNJ_BUILD_EXAMPLES=ON (the default).
• tools/ — prime-table and source-list generators run at build time
• docs/ — extended reference and the descriptor paper

License

BSD 3-Clause — see LICENSE.