Wigcpp: A Wigner-3j 6j and 9j Symbol Calculation Library Written in C++ 17

Wigcpp is a high precision and high performance C++ library for computing Wigner 3j, 6j and 9j symbols using prime factorization and multi word integer arithmetic based on the precalculated prime factorization table.

The computational methods implemented in this project are derived from WIGXJPF.

Features

In Numerical Analysis

It could be proved that wigner $xj$ symbols can be presented as:

$$ W_{x} = \frac{n\sqrt{s}}{q} $$

where $n$, $s$ and $q$ are integers. Algorithm in wigcpp leverages this principle to mitigate precision loss in floating-point arithmetic: it decomposes the mathematical expression for $xj$ symbol into integer components $n$, $s$ and $q$. These integers are then converted to floating-point numbers and used in only three floating-point operations: multiply, square root and division. In overall, six floating-point operations are used in the whole procedure of calculating $xj$ symbol, which keeps the relative error never exceeded $6\varepsilon$. $\varepsilon$ is the machine epsilon. For 80-bit long double of the x87 floating-point unit, $\varepsilon$ is $2^{-64}$. For more deatils, please look through the source code, and Citation.

In Language

1. RAII in C++ automates the acquisition and release of resources, including thread-local and global ones, so users don't need to manage them manually.

2. Exception-free. This library is compatible with -fno-exceptions option.

Build

Wigcpp uses CMake as its build system. To build it, users need to change directory into the project root, using

cmake -S. -B build

to generate compile configurations, and

cmake --build build

to compile the target, then use

cmake --install build --prefix <Install_prefix>

to install the product of the compilation, <Install_prefix> must be substituted as an actual path.

The CMakeLists.txt of this project provides three options to control the products of building: BUILD_SHARED_LIBS, BUILD_FORTRAN_INTERFACE and BUILD_TEST.

These three options were set defaultly as:

options	status
BUILD_FORTRAN_INTERFACE	ON
BUILD_SHARED_LIBS	OFF
BUILD_TEST	OFF

Which means that static library is built defaultly.

API

Wigcpp provides C, C++ and Fortran interface. In C interface, it provides these several functions:

void wigcpp_global_init(int max_two_j, int wigner_type);

void wigcpp_reset_tls();

double clebsch_gordan(int two_j1, int two_j2, int two_m1, int two_m2, int two_J, int two_M);

double wigner3j(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3);

double wigner6j(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6);

double wigner9j(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6, int two_j7, int two_j8, int two_j9);

All of the angular-momentum quantum numbers and magnetic quantum numbers must be passed to the functions in their doubled form, that means if you have a physical value $j$, you must pass $2j$ to these functions.

For the same reason, wigcpp_global_init function accepts twice the maximum physical angular momentum value as its first parameter. And wigcpp_global_init accepts the maximum wigner symbol type that will be used in the whole calculation process as its second parameter, which is must be 3, 6 or 9.

Calling of wigcpp_global_init must be done in the Main Thread, if you call this function in other threads, the behavior of it is undefined.

wigcpp_reset_tls's duty is to reset the Thread Local Storage which is used by wigner3j, wigner6j and wigner9j, then provide a clean thread local state when a thread begins to execute a new task. An example using OpenMP with wigcpp will be provided in the Multi-threaded calling of functuions, which contains the calling of wigcpp_reset_tls.

A simple example in C is as followed.

/* test.c */

#include "stdio.h"
#include "float.h"
#include "wigcpp/wigcpp.h"

int main(void){
	wigcpp_global_init(2 * 100, 9);
	double result = wigner3j(2 * 1, 2 * 1, 2 * 2, 0, 0, 0);
	printf("result is %.*g\n", (int)DBL_DECIMAL_DIG, result);
}

For gcc, the compilation command will be:

gcc -I<Install_prefix>/include \
    -L<Install_prefix>/lib\
    test.c -lwigcpp -lstdc++ -lm

Run the executable, you will see result is 0.36514837167011072 on 8 bytes double type platform.

In C++ interface, we use namespace to encapsulate these functions. Declarations of these functions are:

namespace wigcpp{

	void global_init(int max_two_j, int wigner_type);

	void reset_tls();

	double cg(int two_j1, int two_j2, int two_m1, int two_m2, int two_J, int two_M);

	double three_j(int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3);

	double six_j(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6);

	double nine_j(int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6, int two_j7, int two_j8, int two_j9);

}

A simple example is as followed:

/* test.cpp */

#include "wigcpp/wigcpp.hpp"
#include <cstdio>
#include <cfloat>

auto main(void) -> int{
	wigcpp::global_init(2 * 100, 9);
	double result = wigcpp::three_j(2 * 1, 2 * 1, 2 * 2, 0, 0, 0);
	std::printf("result is %.*g\n", (int)DBL_DECIMAL_DIG, result);
}

Compile command is:

g++ -I<Install_prefix>/include \
    -L<Install_prefix>/lib \
    test.cpp -lwigcpp

For Fortran interface, it maintains the same name as C interface:

subroutine wigcpp_global_init(max_two_j, wigner_type)
	integer :: max_two_j, wigner_type
end subroutine

subroutine wigcpp_reset_tls()
end subroutine

function clebsch_gordan(two_j1, two_j2, two_m1, two_m2, two_J, two_M)
	integer :: two_j1, two_j2, two_m1, two_m2, two_J, two_M
	real :: clebsch_gordan
end function

function wigner3j(two_j1, two_j2, two_j3, two_m1, two_m2, two_m3)
	integer :: two_j1, two_j2, two_j3, two_m1, two_m2, two_m3
	real :: wigner3j
end function

function wigner6j(two_j1, two_j2, two_j3, two_j4, two_j5, two_j6)
	integer :: two_j1, two_j2, two_j3, two_j4, two_j5, two_j6
	real :: wigner6j
end function

function wigner9j(two_j1, two_j2, two_j3, two_j4, two_j5, two_j6, two_j7, two_j8, two_j9)
	integer :: two_j1, two_j2, two_j3, two_j4, two_j5, two_j6, two_j7, two_j8, two_j9
	real :: wigner9j
end function

A simple example is:

!test.f90

program main
	use wigcpp
	implicit none
	real :: result

	call wigcpp_global_init(2 * 100, 9)

	result = clebsch_gordan(35, 37, 3, 5, 66, 8)
	print *, "result of cg is", result

	result = wigner3j(2 * 1, 2 * 1, 2 * 2, 0, 0, 0)
	print *, "result of 3j is", result

	result = wigner6j(2*2, 2*2, 2*2, 2*2, 2*2, 2*2)
	print *, "result of 6j is", result

	result = wigner9j(2*20, 2*20, 2*20, 2*20, 2*20, 2*20, 2*20, 2*20, 2*20)
	print *, "result of 9j is", result

end program main

Compile command will be:

gfortran -I<Install_prefix>/include \
         -L<Install_prefix>/lib \
         test.f90 -lwigcpp -lstdc++

Then run the executable, you will see

 result of cg is  0.109003529
 result of 3j is  0.365148365
 result of 6j is  -4.28571440E-02
 result of 9j is   5.73250327E-05

in the stdout.

Multi-threaded calling of functions

The use of thread_local in C++11 eliminates the need for users to explicitly initialize thread-local resources at thread startup. So when calling these functions in pthread-like(1:1 threading model) multi-threads, you can proceed just as you would in a single-threaded environment.

#include "wigcpp/wigcpp.hpp"
#include <iostream>
#include <thread>
#include <vector>

struct ThreadData{
	int two_j1, two_j2, two_j3, two_m1, two_m2, two_m3;
	double result;
};

auto main(void) -> int{
	constexpr int kThreads = 4;

	std::vector<std::thread> threads;
	std::vector<ThreadData> thread_data(kThreads);

	thread_data[0] = {2 * 400, 2 * 80, 2 * 480, 2 * 1, -1 * 2, 0, 0.0};
	thread_data[1] = {3, 7, 10, 1,-1, 0, 0.0};
	thread_data[2] = {2, 4, 6, 0, 0, 0, 0.0};
	thread_data[3] = {2 * 300, 2 * 400, 2 * 700, 2 *50, -25 * 2, -25 * 2, 0.0};

	wigcpp::global_init(2 * 1000, 3);

	for(int i = 0; i < kThreads; ++i){
		threads.emplace_back([&thread_data, i]{
			thread_data[i].result = wigcpp::three_j(
				thread_data[i].two_j1, thread_data[i].two_j2, thread_data[i].two_j3, 
				thread_data[i].two_m1, thread_data[i].two_m2, thread_data[i].two_m3);
		});
	}

	for(auto &t : threads){
		t.join();
	}

  for(int i = 0; i < kThreads; ++i){
    std::cout << "Result in Thread" << i << ': ' << thread_data[i].result << std::endl;
  }
}

Output in stdout will be:

Result in Thread0: 0.00840976
Result in Thread1: 0.194625
Result in Thread2: -0.29277
Result in Thread3: -6.07534e-05

But when using thread pool, users must make sure that TLS has a certain and correct initial state when a worker thread begins processing a new task. wigcpp_reset_tls is designed to handle this problem. An example using OpenMP in Fortran is

!$OMP PARALLEL PRIVATE(iJtot, imj, imjp, iLr, iLp) NUM_THREADS(24)

call wigcpp_reset_tls

!$OMP DO COLLAPSE(5) SCHEDULE(DYNAMIC)
  do iJtot = 0, nJtot
      do imj = -qn_j, qn_j
          do imjp = -qn_jp, qn_jp
              do iLr = 0, nLr
                  do iLp = 0, nLp
                      iml = imj - imjp
                      if (abs(iml) > iLp) cycle
                      if ((iJtot > qn_j+iLr .or. iJtot < abs(qn_j-iLr)) .or. (iJtot > qn_jp+iLp .or. iJtot < abs(qn_jp-iLp))) cycle
                      CG1(iJtot,imj,iLr) = clebsch_gordan(qn_j*2,iLr*2,imj*2,0,iJtot*2,imj*2)
                      CG2(iJtot,imjp,iLp,iml) = clebsch_gordan(qn_jp*2,iLp*2,imjp*2,iml*2,iJtot*2,(imjp + iml)*2)
                  end do
              end do
          end do
      end do
  end do
!$OMP END DO
!$OMP END PARALLEL

Also, when you using wigcpp in any M:N threading model, make sure that using wigcpp_tls_reset to reset the state of Thread Local Storage at the begining of the task process in every threads.

ToDo

This project is in progress. In further, more benchmarks will be implemented, and performance optimization will be conducted through techniques such as:

1. LTO/PGO builds (LTO using GNU has added).
2. More tests for wigner6j and wigner9j.

License

This project is licensed under the GPL-3.0 license, see the LICENSE for details.

Contact

If you have any question, feel free to open an issue, or Email to the maintainer: liuhaotian0406@163.com

Citation

This project uses prime factorization and multi word integer arithmetic to calculate wigner-3j, 6j and 9j symbols, which is derived from the article:

H. T. Johansson and C. Forssén, "Fast and Accurate Evaluation of Wigner 3j, 6j, and 9j Symbols Using Prime Factorization and Multiword Integer Arithmetic" , SIAM J. Sci. Comput., 38(1) (2016), A376-A384.

Acknowledgments

Algorithm in this project was inspired by the WIGXJPF library, developed by the main developer Håkan T. Johansson, for high-precision computation of Wigner symbols. WIGXJPF is licensed under the GNU Lesser General Public License v3.0 (LGPL-3.0), and a copy of the license is provided in COPYING.LESSER. We thank the authors for their open contribution to computational physics.