The coupling of three spins requires the total sum to be integer. This is
not checked on gsl_sf_coupling_6j_e or gsl_sf_coupling_9j_e.
We discovered it when trying to evaluate the 9j symbol:
(1/2 1/2 1/2)
(1/2 1/2 1/2)
( 0 0 0 )
and got a non-zero result. We believe that the recursive calculations of
6j should have given the zero, but found also that lacking similar checks.
The attached patch also include some test cases. The 6js have been
cross-checked with Mathematica.
Best regards,
Håkan Johansson
diff -ru gsl-1.15-orig/specfunc//coupling.c gsl-1.15/specfunc//coupling.c
--- gsl-1.15-orig/specfunc//coupling.c 2010-12-26 18:57:08.000000000 +0100
+++ gsl-1.15/specfunc//coupling.c 2013-07-11 15:49:06.000000000 +0200
@@ -206,6 +206,10 @@
|| triangle_selection_fails(two_ja, two_je, two_jf)
|| triangle_selection_fails(two_jb, two_jd, two_jf)
|| triangle_selection_fails(two_je, two_jd, two_jc)
+ || GSL_IS_ODD(two_ja + two_jb + two_jc)
+ || GSL_IS_ODD(two_ja + two_je + two_jf)
+ || GSL_IS_ODD(two_jb + two_jd + two_jf)
+ || GSL_IS_ODD(two_je + two_jd + two_jc)
) {
result->val = 0.0;
result->err = 0.0;
@@ -335,6 +339,12 @@
|| triangle_selection_fails(two_ja, two_jd, two_jg)
|| triangle_selection_fails(two_jb, two_je, two_jh)
|| triangle_selection_fails(two_jc, two_jf, two_ji)
+ || GSL_IS_ODD(two_ja + two_jb + two_jc)
+ || GSL_IS_ODD(two_jd + two_je + two_jf)
+ || GSL_IS_ODD(two_jg + two_jh + two_ji)
+ || GSL_IS_ODD(two_ja + two_jd + two_jg)
+ || GSL_IS_ODD(two_jb + two_je + two_jh)
+ /* two_jc + two_jf + two_ji even as a consequence of above tests */
) {
result->val = 0.0;
result->err = 0.0;
Endast i gsl-1.15/specfunc/: coupling.c~
Endast i gsl-1.15/specfunc/: .deps
Endast i gsl-1.15/specfunc/: Makefile
diff -ru gsl-1.15-orig/specfunc//test_sf.c gsl-1.15/specfunc//test_sf.c
--- gsl-1.15-orig/specfunc//test_sf.c 2010-12-26 21:37:41.000000000 +0100
+++ gsl-1.15/specfunc//test_sf.c 2013-07-11 15:48:14.000000000 +0200
@@ -414,6 +414,8 @@
TEST_SF(s, gsl_sf_coupling_6j_e, (4, 4, 4, 2, 2, 2, &r), sqrt(7.0/3.0)/10.0, TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_coupling_6j_e, (6, 6, 6, 4, 4, 4, &r), -sqrt(3.0/5.0)/14.0, TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_coupling_6j_e, (6, 6, 6, 4, 4, 2, &r), -sqrt(3.0/5.0)/7.0, TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_coupling_6j_e, (1, 0, 1, 0, 1, 0, &r), -sqrt(1.0/2.0), TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_coupling_6j_e, (1, 0, 1, 1, 0, 1, &r), -1.0/2.0, TEST_TOL0, GSL_SUCCESS);
/* Test 6j error checking */
@@ -433,6 +435,12 @@
TEST_SF(s, gsl_sf_coupling_6j_e, (2, 7, 4, 2, 2, 2, &r), 0, 0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_coupling_6j_e, (7, 2, 4, 2, 2, 2, &r), 0, 0, GSL_SUCCESS);
+ /* Test 6j half-integer/integer coupling conditions */
+
+ TEST_SF(s, gsl_sf_coupling_6j_e, (1, 1, 1, 0, 1, 1, &r), 0, 0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_coupling_6j_e, (1, 1, 1, 1, 0, 1, &r), 0, 0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_coupling_6j_e, (1, 1, 1, 1, 1, 0, &r), 0, 0, GSL_SUCCESS);
+
/* Test 9j */
TEST_SF(s, gsl_sf_coupling_9j_e, (4, 2, 4, 3, 3, 2, 1, 1, 2, &r), -sqrt(1.0/6.0)/10.0, TEST_TOL2, GSL_SUCCESS);
@@ -459,6 +467,11 @@
TEST_SF(s, gsl_sf_coupling_9j_e, (4, 2, 4, 3, 3, 2, 10, 1, 2, &r), 0, 0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_coupling_9j_e, (4, 2, 4, 3, 3, 2, 1, 10, 2, &r), 0, 0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_coupling_9j_e, (4, 2, 4, 3, 3, 2, 1, 1, 10, &r), 0, 0, GSL_SUCCESS);
+
+ /* Test 9j half-integer/integer coupling conditions */
+
+ TEST_SF(s, gsl_sf_coupling_9j_e, (1, 1, 1, 1, 1, 1, 0, 0, 0, &r), 0, 0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_coupling_9j_e, (1, 1, 0, 1, 1, 0, 1, 1, 0, &r), 0, 0, GSL_SUCCESS);
return s;
}
Endast i gsl-1.15/specfunc/: test_sf.c~
