diff --git a/src/soap/functions.cpp b/src/soap/functions.cpp
index d80c846e1fbe8770dfe62d7a0b23e49f973bc539..59f8ec5290ab4105d9ac474ddc1cec09c8a2d0fd 100644
--- a/src/soap/functions.cpp
+++ b/src/soap/functions.cpp
@@ -247,14 +247,15 @@ std::complex<double> pow_nnan(std::complex<double> z, double a) {
}
}
-const int FACTORIAL_CACHE_SIZE = 16;
+const int FACTORIAL_CACHE_SIZE = 19;
const long int FACTORIAL_CACHE[] = {
1,
1, 2, 6,
24, 120, 720,
5040, 40320, 362880,
3628800, 39916800, 479001600,
- 6227020800, 87178291200, 1307674368000 };
+ 6227020800, 87178291200, 1307674368000,
+ 20922789888000, 355687428096000, 6402373705728000};
long int factorial(int x) {
if (x < FACTORIAL_CACHE_SIZE) {
@@ -270,14 +271,15 @@ long int factorial(int x) {
}
}
-const int FACTORIAL2_CACHE_SIZE = 16;
+const int FACTORIAL2_CACHE_SIZE = 19;
const long int FACTORIAL2_CACHE[] = {
1,
1, 2, 3,
8, 15, 48,
105, 384, 945,
3840, 10395, 46080,
- 135135, 645120, 2027025 };
+ 135135, 645120, 2027025,
+ 10321920, 34459425, 185794560};
long int factorial2(int x) {
if (x < 0) {
diff --git a/src/soap/linalg/wigner.cpp b/src/soap/linalg/wigner.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..efe2bbcd2458ba6d326e2b860c017c1e748610ab
--- /dev/null
+++ b/src/soap/linalg/wigner.cpp
@@ -0,0 +1,505 @@
+/*******************************************************-/
+ * This source code is subject to the terms of the GNU -/
+ * Lesser Public License. If a copy of the LGPL was not -/
+ * distributed with this file, you can obtain one at -/
+ * https://www.gnu.org/licenses/lgpl.html. -/
+ ********************************************************/
+
+#include "soap/linalg/wigner.hpp"
+
+namespace WignerSymbols {
+std::vector<double> wigner3j(double l2, double l3,
+ double m1, double m2, double m3)
+{
+ // We compute the numeric limits of double precision.
+ double huge = sqrt(std::numeric_limits<double>::max()/20.0);
+ double srhuge = sqrt(huge);
+ double tiny = std::numeric_limits<double>::min();
+ double srtiny = sqrt(tiny);
+ double eps = std::numeric_limits<double>::epsilon();
+
+ // We enforce the selection rules.
+ bool select(true);
+ select = (
+ std::fabs(m1+m2+m3)<eps
+ && std::fabs(m2) <= l2+eps
+ && std::fabs(m3) <= l3+eps
+ );
+
+ if (!select) return std::vector<double>(1,0.0);
+
+ // We compute the limits of l1.
+ double l1min = std::max(std::fabs(l2-l3),std::fabs(m1));
+ double l1max = l2+l3;
+
+ // We compute the size of the resulting array.
+ int size = (int)std::floor(l1max-l1min+1.0+eps);
+ std::vector<double> thrcof(size,0.0);
+
+ // If l1min=l1max, we have an analytical formula.
+ if (size==1)
+ {
+ thrcof[0] = pow(-1.0,std::floor(std::fabs(l2+m2-l3+m3)))/sqrt(l1min+l2+l3+1.0);
+ }
+
+ // Another special case where the recursion relation fails.
+ else
+ {
+ // We start with an arbitrary value.
+ thrcof[0] = srtiny;
+
+ // From now on, we check the variation of |alpha(l1)|.
+ double alphaNew, l1(l1min);
+ if (l1min==0.0)
+ alphaNew = -(m3-m2+2.0*wigner3j_auxB(l1,l2,l3,m1,m2,m3))/wigner3j_auxA(1.0,l2,l3,m1,m2,m3);
+ else
+ alphaNew = -wigner3j_auxB(l1min,l2,l3,m1,m2,m3)
+ /(l1min*wigner3j_auxA(l1min+1.0,l2,l3,m1,m2,m3));
+
+ // We compute the two-term recursion.
+ thrcof[1] = alphaNew*thrcof[0];
+
+ // If size > 2, we start the forward recursion.
+ if (size>2)
+ {
+ // We start with an arbitrary value.
+ thrcof[0] = srtiny;
+
+ // From now on, we check the variation of |alpha(l1)|.
+ double alphaOld, alphaNew, beta, l1(l1min);
+ if (l1min==0.0)
+ alphaNew = -(m3-m2+2.0*wigner3j_auxB(l1,l2,l3,m1,m2,m3))/wigner3j_auxA(1.0,l2,l3,m1,m2,m3);
+ else
+ alphaNew = -wigner3j_auxB(l1min,l2,l3,m1,m2,m3)
+ /(l1min*wigner3j_auxA(l1min+1.0,l2,l3,m1,m2,m3));
+
+ // We compute the two-term recursion.
+ thrcof[1] = alphaNew*thrcof[0];
+
+ // We compute the rest of the recursion.
+ int i = 1;
+ bool alphaVar = false;
+ do
+ {
+ // Bookkeeping:
+ i++; // Next term in recursion
+ alphaOld = alphaNew; // Monitoring of |alpha(l1)|.
+ l1 += 1.0; // l1 = l1+1
+
+ // New coefficients in recursion.
+ alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+ /(l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3));
+
+ beta = -(l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3)
+ /(l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3));
+
+ // Application of the recursion.
+ thrcof[i] = alphaNew*thrcof[i-1]+beta*thrcof[i-2];
+
+ // We check if we are overflowing.
+ if (std::fabs(thrcof[i])>srhuge)
+ {
+ std::cout << "We renormalized the forward recursion." << std::endl;
+ for (std::vector<double>::iterator it = thrcof.begin(); it != thrcof.begin()+i; ++it)
+ {
+ //if (std::fabs(*it) < srtiny) *it = 0;
+ //else
+ *it /= srhuge;
+ }
+ }
+
+ // This piece of code checks whether we have reached
+ // the classical region. If we have, the second if
+ // sets alphaVar to true and we break this loop at the
+ // next iteration because we need thrcof(l1mid+1) to
+ // compute the scalar lambda.
+ if (alphaVar) break;
+
+ if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0)
+ alphaVar=true;
+
+ } while (i<(size-1)); // Loop stops when we have computed all values.
+
+ // If this is the case, we have stumbled upon a classical region.
+ // We start the backwards recursion.
+ if (i!=size-1)
+ {
+ // We keep the two terms around l1mid to compute the factor later.
+ double l1midm1(thrcof[i-2]),l1mid(thrcof[i-1]),l1midp1(thrcof[i]);
+
+ // We compute the backward recursion by providing an arbitrary
+ // startint value.
+ thrcof[size-1] = srtiny;
+
+ // We compute the two-term recursion.
+ l1 = l1max;
+ alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+ /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+ thrcof[size-2] = alphaNew*thrcof[size-1];
+
+ // We compute the rest of the backward recursion.
+ int j = size-2;
+ do
+ {
+ // Bookkeeping
+ j--; // Previous term in recursion.
+ l1 -= 1.0; // l1 = l1-1
+
+ // New coefficients in recursion.
+ alphaNew = -wigner3j_auxB(l1,l2,l3,m1,m2,m3)
+ /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+ beta = -l1*wigner3j_auxA(l1+1.0,l2,l3,m1,m2,m3)
+ /((l1+1.0)*wigner3j_auxA(l1,l2,l3,m1,m2,m3));
+
+ // Application of the recursion.
+ thrcof[j] = alphaNew*thrcof[j+1]+beta*thrcof[j+2];
+
+ // We check if we are overflowing.
+ if (std::fabs(thrcof[j]>srhuge))
+ {
+ std::cout << "We renormalized the backward recursion." << std::endl;
+ for (std::vector<double>::iterator it = thrcof.begin()+j; it != thrcof.end(); ++it)
+ {
+ //if (std::fabs(*it) < srtiny) *it = 0;
+ //else
+ *it /= srhuge;
+ }
+ }
+
+ } while (j>(i-2)); // Loop stops when we are at l1=l1mid-1.
+
+ // We now compute the scaling factor for the forward recursion.
+ double lambda = (l1midp1*thrcof[j+2]+l1mid*thrcof[j+1]+l1midm1*thrcof[j])
+ /(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1);
+
+ // We scale the forward recursion.
+ for (std::vector<double>::iterator it = thrcof.begin(); it != thrcof.begin()+j; ++it)
+ {
+ *it *= lambda;
+ }
+ }
+ }
+ }
+
+ // We compute the overall factor.
+ double sum = 0.0;
+ for (int k=0;k<size;k++)
+ {
+ sum += (2.0*(l1min+k)+1.0)*thrcof[k]*thrcof[k];
+ }
+ //std::cout << sum << std::endl;
+
+ //std::cout << "(-1)^(l2-l3-m1): " << pow(-1.0,l2-l3-m1) << " sgn:" << sgn(thrcof[size-1]) << std::endl;
+ double c1 = pow(-1.0,l2-l3-m1)*sgn(thrcof[size-1]);
+ //std::cout << "c1: " << c1 << std::endl;
+ for (std::vector<double>::iterator it = thrcof.begin(); it != thrcof.end(); ++it)
+ {
+ //std::cout << *it << ", " << c1 << ", ";
+ *it *= c1/sqrt(sum);
+ //std::cout << *it << std::endl;
+ }
+ return thrcof;
+}
+
+double wigner3j(double l1, double l2, double l3,
+ double m1, double m2, double m3)
+{
+ // We enforce the selection rules.
+ bool select(true);
+ select = (
+ std::fabs(m1+m2+m3)<1.0e-10
+ && std::floor(l1+l2+l3)==(l1+l2+l3)
+ && l3 >= std::fabs(l1-l2)
+ && l3 <= l1+l2
+ && std::fabs(m1) <= l1
+ && std::fabs(m2) <= l2
+ && std::fabs(m3) <= l3
+ );
+
+ if (!select) return 0.0;
+
+ // We compute l1min and the position of the array we will want.
+ double l1min = std::max(std::fabs(l2-l3),std::fabs(m1));
+
+ // We fetch the proper value in the array.
+ int index = (int)(l1-l1min);
+
+ return wigner3j(l2,l3,m1,m2,m3)[index];
+}
+
+std::vector<double> wigner6j(double l2, double l3,
+ double l4, double l5, double l6)
+{
+ // We compute the numeric limits of double precision.
+ double huge = std::numeric_limits<double>::max();
+ double srhuge = sqrt(huge);
+ double tiny = std::numeric_limits<double>::min();
+ double srtiny = sqrt(tiny);
+ double eps = std::numeric_limits<double>::epsilon();
+
+ // We enforce the selection rules.
+ bool select(true);
+
+ // Triangle relations for the four tryads
+ select = (
+ std::fabs(l4-l2) <= l6 && l6 <= l4+l2
+ && std::fabs(l4-l5) <= l3 && l3 <= l4+l5
+ );
+
+ // Sum rule of the tryads
+ select = (
+ std::floor(l4+l2+l6)==(l4+l2+l6)
+ && std::floor(l4+l5+l3)==(l4+l5+l3)
+ );
+
+ if (!select) return std::vector<double>(1,0.0);
+
+ // We compute the limits of l1.
+ double l1min = std::max(std::fabs(l2-l3),std::fabs(l5-l6));
+ double l1max = std::min(l2+l3,l5+l6);
+
+ // We compute the size of the resulting array.
+ unsigned int size = (int)std::floor(l1max-l1min+1.0+eps);
+ std::vector<double> sixcof(size,0.0);
+
+ // If l1min=l1max, we have an analytical formula.
+ if (size==1)
+ {
+ sixcof[0] = 1.0/sqrt((l1min+l1min+1.0)*(l4+l4+1.0));
+ sixcof[0] *= ((int)std::floor(l2+l3+l5+l6+eps) & 1 ? -1.0 : 1.0);
+ }
+
+ // Otherwise, we start the forward recursion.
+ else
+ {
+ // We start with an arbitrary value.
+ sixcof[0] = srtiny;
+
+ // From now on, we check the variation of |alpha(l1)|.
+ double alphaNew, l1(l1min);
+
+ if (l1min==0)
+ alphaNew = -(l2*(l2+1.0)+l3*(l3+1.0)+l5*(l5+1.0)+l6*(l6+1.0)-2.0*l4*(l4+1.0))/wigner6j_auxA(1.0,l2,l3,l4,l5,l6);
+
+ else
+ alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+ /(l1min*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+
+ // We compute the two-term recursion.
+ sixcof[1] = alphaNew*sixcof[0];
+
+ if (size>2)
+ {
+ // We start with an arbitrary value.
+ sixcof[0] = srtiny;
+
+ // From now on, we check the variation of |alpha(l1)|.
+ double alphaOld, alphaNew, beta, l1(l1min);
+ if (l1min==0)
+ alphaNew = -(l2*(l2+1.0)+l3*(l3+1.0)+l5*(l5+1.0)+l6*(l6+1.0)-2.0*l4*(l4+1.0))/wigner6j_auxA(1.0,l2,l3,l4,l5,l6);
+
+ else
+ alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+ /(l1min*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+
+ // We compute the two-term recursion.
+ sixcof[1] = alphaNew*sixcof[0];
+
+ // We compute the rest of the recursion.
+ unsigned int i = 1;
+ bool alphaVar = false;
+ do
+ {
+ // Bookkeeping:
+ i++; // Next term in recursion
+ alphaOld = alphaNew; // Monitoring of |alpha(l1)|.
+ l1 += 1.0; // l1 = l1+1
+
+ // New coefficients in recursion.
+ alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+ /(l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+
+ beta = -(l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6)
+ /(l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6));
+
+ // Application of the recursion.
+ sixcof[i] = alphaNew*sixcof[i-1]+beta*sixcof[i-2];
+
+ // We check if we are overflowing.
+ if (std::fabs(sixcof[i]>srhuge))
+ {
+ std::cout << "We renormalized the forward recursion." << std::endl;
+ for (std::vector<double>::iterator it = sixcof.begin(); it != sixcof.begin()+i; ++it)
+ {
+ *it /= srhuge;
+ }
+ }
+
+ // This piece of code checks whether we have reached
+ // the classical region. If we have, the second if
+ // sets alphaVar to true and we break this loop at the
+ // next iteration because we need sixcof(l1mid+1) to
+ // compute the scalar.
+ if (alphaVar) break;
+
+ if (std::fabs(alphaNew)-std::fabs(alphaOld)>0.0)
+ alphaVar=true;
+
+ } while (i<(size-1)); // Loop stops when we have computed all values.
+
+ // If this is the case, we have stumbled upon a classical region.
+ // We start the backwards recursion.
+ if (i!=size-1)
+ {
+ // We keep the two terms around l1mid to compute the factor later.
+ double l1midm1(sixcof[i-2]),l1mid(sixcof[i-1]),l1midp1(sixcof[i]);
+
+ // We compute the backward recursion by providing an arbitrary
+ // startint value.
+ sixcof[size-1] = srtiny;
+
+ // We compute the two-term recursion.
+ l1 = l1max;
+ alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+ /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+ sixcof[size-2] = alphaNew*sixcof[size-1];
+
+ // We compute the rest of the backward recursion.
+ unsigned int j = size-2;
+ do
+ {
+ // Bookkeeping
+ j--; // Previous term in recursion.
+ l1 -= 1.0; // l1 = l1-1
+
+ // New coefficients in recursion.
+ alphaNew = -wigner6j_auxB(l1,l2,l3,l4,l5,l6)
+ /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+ beta = -l1*wigner6j_auxA(l1+1.0,l2,l3,l4,l5,l6)
+ /((l1+1.0)*wigner6j_auxA(l1,l2,l3,l4,l5,l6));
+
+ // Application of the recursion.
+ sixcof[j] = alphaNew*sixcof[j+1]+beta*sixcof[j+2];
+
+ // We check if we are overflowing.
+ if (std::fabs(sixcof[j]>srhuge))
+ {
+ std::cout << "We renormalized the backward recursion." << std::endl;
+ for (std::vector<double>::iterator it = sixcof.begin()+j; it != sixcof.end(); ++it)
+ {
+ *it /= srhuge;
+ }
+ }
+
+ } while (j>(i-2)); // Loop stops when we are at l1=l1mid-1.
+
+ // We now compute the scaling factor for the forward recursion.
+ double lambda = (l1midp1*sixcof[j+2]+l1mid*sixcof[j+1]+l1midm1*sixcof[j])
+ /(l1midp1*l1midp1+l1mid*l1mid+l1midm1*l1midm1);
+
+ // We scale the forward recursion.
+ for (std::vector<double>::iterator it = sixcof.begin(); it != sixcof.begin()+j; ++it)
+ {
+ *it *= lambda;
+ }
+ }
+ }
+ }
+
+ // We compute the overall factor.
+ double sum = 0.0;
+ for (unsigned int k=0;k<size;k++)
+ {
+ sum += (2.0*(l1min+k)+1.0)*(2.0*l4+1.0)*sixcof[k]*sixcof[k];
+ }
+ double c1 = pow(-1.0,std::floor(l2+l3+l5+l6+eps))*sgn(sixcof[size-1])/sqrt(sum);
+
+ for (std::vector<double>::iterator it = sixcof.begin(); it != sixcof.end(); ++it)
+ {
+ *it *= c1;
+ }
+ return sixcof;
+}
+
+double wigner6j(double l1, double l2, double l3,
+ double l4, double l5, double l6)
+{
+ // We enforce the selection rules.
+ bool select(true);
+
+ // Triangle relations for the four tryads
+ select = (
+ std::fabs(l1-l2) <= l3 && l3 <= l1+l2
+ && std::fabs(l1-l5) <= l6 && l6 <= l1+l5
+ && std::fabs(l4-l2) <= l6 && l6 <= l4+l2
+ && std::fabs(l4-l5) <= l3 && l3 <= l4+l5
+ );
+
+ // Sum rule of the tryads
+ select = (
+ std::floor(l1+l2+l3)==(l1+l2+l3)
+ && std::floor(l1+l5+l6)==(l1+l5+l6)
+ && std::floor(l4+l2+l6)==(l4+l2+l6)
+ && std::floor(l4+l5+l3)==(l4+l5+l3)
+ );
+
+ if (!select) return 0.0;
+
+ // We compute l1min and the position of the array we will want.
+ double l1min = std::max(std::fabs(l2-l3),std::fabs(l5-l6));
+ int index = (int)(l1-l1min);
+
+ return wigner6j(l2,l3,l4,l5,l6)[index];
+}
+
+double wigner3j_auxA(double l1, double l2, double l3,
+ double m1, double m2, double m3)
+{
+ double T1 = l1*l1-pow(l2-l3,2.0);
+ double T2 = pow(l2+l3+1.0,2.0)-l1*l1;
+ double T3 = l1*l1-m1*m1;
+
+ return sqrt(T1*T2*T3);
+}
+
+double wigner3j_auxB(double l1, double l2, double l3,
+ double m1, double m2, double m3)
+{
+ double T1 = -(2.0*l1+1.0);
+ double T2 = l2*(l2+1.0)*m1;
+ double T3 = l3*(l3+1.0)*m1;
+ double T4 = l1*(l1+1.0)*(m3-m2);
+
+ return T1*(T2-T3-T4);
+}
+
+double wigner6j_auxA(double l1, double l2, double l3,
+ double l4, double l5, double l6)
+{
+ double T1 = l1*l1-pow(l2-l3,2.0);
+ double T2 = pow(l2+l3+1.0,2.0)-l1*l1;
+ double T3 = l1*l1-pow(l5-l6,2.0);
+ double T4 = pow(l5+l6+1.0,2.0)-l1*l1;
+
+ return sqrt(T1*T2*T3*T4);
+}
+
+double wigner6j_auxB(double l1, double l2, double l3,
+ double l4, double l5, double l6)
+{
+ double T0 = 2.0*l1+1.0;
+
+ double T1 = l1*(l1+1.0);
+ double T2 = -l1*(l1+1.0)+l2*(l2+1.0)+l3*(l3+1.0);
+
+ double T3 = l5*(l5+1.0);
+ double T4 = l1*(l1+1.0)+l2*(l2+1.0)-l3*(l3+1.0);
+
+ double T5 = l6*(l6+1.0);
+ double T6 = l1*(l1+1.0)-l2*(l2+1.0)+l3*(l3+1.0);
+
+ double T7 = 2.0*l1*(l1+1.0)*l4*(l4+1.0);
+
+ return (T0*(T1*T2+T3*T4+T5*T6-T7));
+}
+}
diff --git a/src/soap/linalg/wigner.hpp b/src/soap/linalg/wigner.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cb90d534379b0261a7a2f33811e71e54c06ae3a6
--- /dev/null
+++ b/src/soap/linalg/wigner.hpp
@@ -0,0 +1,87 @@
+/*******************************************************-/
+ * This source code is subject to the terms of the GNU -/
+ * Lesser Public License. If a copy of the LGPL was not -/
+ * distributed with this file, you can obtain one at -/
+ * https://www.gnu.org/licenses/lgpl.html. -/
+ ********************************************************/
+
+#ifndef WIGNER_SYMBOLS_CPP_H
+#define WIGNER_SYMBOLS_CPP_H
+
+/** \file wignerSymbols-cpp.h
+ *
+ * \author Joey Dumont <joey.dumont@gmail.com>
+ *
+ * \since 2013-08-16
+ *
+ * \brief Defines utility functions for the evaluation of Wigner-3j and -6j symbols.
+ *
+ * We compute the Wigner-3j and -6j symbols
+ * f(L1) = ( L1 L2 L3)
+ * (-M2-M3 M2 M3)
+ * for all allowed values of L1, the other parameters
+ * being held fixed. The algorithm is based on the work
+ * by Schulten and Gordon.
+ * K. Schulten, "Exact recursive evaluation of 3j- and 6j-coefficients for quantum-mechanical coupling of angular momenta,"
+ * J. Math. Phys. 16, 1961 (1975).
+ * K. Schulten and R. G. Gordon, "Recursive evaluation of 3j and 6j coefficients,"
+ * Comput. Phys. Commun. 11, 269–278 (1976).
+ */
+
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <vector>
+#include <iostream>
+
+namespace WignerSymbols {
+
+template <typename T>
+double sgn(T val)
+{
+ int sgn = (T(0) < val) - (val < T(0));
+ if (sgn == 0)
+ return 1.0;
+ else
+ return (double)sgn;
+}
+
+/*! @name Evaluation of Wigner-3j and -6j symbols.
+ * We implement Schulten's algorithm in C++.
+ */
+///@{
+std::vector<double> wigner3j(double l2, double l3,
+ double m1, double m2, double m3);
+
+double wigner3j(double l1, double l2, double l3,
+ double m1, double m2, double m3);
+
+double wigner3j_auxA(double l1, double l2, double l3,
+ double m1, double m2, double m3);
+double wigner3j_auxB(double l1, double l2, double l3,
+ double m1, double m2, double m3);
+
+std::vector<double> wigner6j(double l2, double l3,
+ double l4, double l5, double l6);
+
+double wigner6j(double l1, double l2, double l3,
+ double l4, double l5, double l6);
+
+double wigner6j_auxA(double l1, double l2, double l3,
+ double l4, double l5, double l6);
+
+double wigner6j_auxB(double l1, double l2, double l3,
+ double l4, double l5, double l6);
+
+/*! Computes the Clebsch-Gordan coefficient by relating it to the
+ * Wigner 3j symbol. It sometimes eases the notation to use the
+ * Clebsch-Gordan coefficients directly. */
+inline double clebschGordan(double l1, double l2, double l3,
+ double m1, double m2, double m3)
+{
+ // We simply compute it via the 3j symbol.
+ return (pow(-1.0,l1-l2+m3)*sqrt(2.0*l3+1.0)*wigner3j(l1,l2,l3,m1,m2,-m3));
+}
+}
+
+#endif // WIGNER_SYMBOLS_CPP_H