From 53ca34dcb47966dcdcad8584967dfcc4cb9457d3 Mon Sep 17 00:00:00 2001 From: Alexey Alexandrovich Illarionov Date: Wed, 6 Nov 2013 03:40:01 -0500 Subject: [PATCH] coupling3j rewritten for better precision --- specfunc/Makefile.am | 2 +- specfunc/coupling.c | 296 ++++++++++++++++++++++++++++++++++---------- specfunc/coupling3j_regge.c | 76 ++++++++++++ 3 files changed, 311 insertions(+), 63 deletions(-) create mode 100644 specfunc/coupling3j_regge.c diff --git a/specfunc/Makefile.am b/specfunc/Makefile.am index 11b1fbc..07bcbc3 100644 --- a/specfunc/Makefile.am +++ b/specfunc/Makefile.am @@ -6,7 +6,7 @@ noinst_HEADERS = bessel_amp_phase.h bessel_olver.h bessel_temme.h bessel.h hyper INCLUDES = -I$(top_srcdir) -libgslspecfunc_la_SOURCES = airy.c airy_der.c airy_zero.c atanint.c bessel.c bessel.h bessel_I0.c bessel_I1.c bessel_In.c bessel_Inu.c bessel_J0.c bessel_J1.c bessel_Jn.c bessel_Jnu.c bessel_K0.c bessel_K1.c bessel_Kn.c bessel_Knu.c bessel_Y0.c bessel_Y1.c bessel_Yn.c bessel_Ynu.c bessel_amp_phase.c bessel_amp_phase.h bessel_i.c bessel_j.c bessel_k.c bessel_olver.c bessel_temme.c bessel_y.c bessel_zero.c bessel_sequence.c beta.c beta_inc.c clausen.c coulomb.c coupling.c coulomb_bound.c dawson.c debye.c dilog.c elementary.c ellint.c elljac.c erfc.c exp.c expint.c expint3.c fermi_dirac.c gegenbauer.c gamma.c gamma_inc.c hyperg_0F1.c hyperg_2F0.c hyperg_1F1.c hyperg_2F1.c hyperg_U.c hyperg.c laguerre.c lambert.c legendre_H3d.c legendre_Qn.c legendre_con.c legendre_poly.c log.c mathieu_angfunc.c mathieu_charv.c mathieu_coeff.c mathieu_radfunc.c mathieu_workspace.c poch.c pow_int.c psi.c recurse.h result.c shint.c sinint.c synchrotron.c transport.c trig.c zeta.c +libgslspecfunc_la_SOURCES = airy.c airy_der.c airy_zero.c atanint.c bessel.c bessel.h bessel_I0.c bessel_I1.c bessel_In.c bessel_Inu.c bessel_J0.c bessel_J1.c bessel_Jn.c bessel_Jnu.c bessel_K0.c bessel_K1.c bessel_Kn.c bessel_Knu.c bessel_Y0.c bessel_Y1.c bessel_Yn.c bessel_Ynu.c bessel_amp_phase.c bessel_amp_phase.h bessel_i.c bessel_j.c bessel_k.c bessel_olver.c bessel_temme.c bessel_y.c bessel_zero.c bessel_sequence.c beta.c beta_inc.c clausen.c coulomb.c coupling.c coupling3j_regge.c coulomb_bound.c dawson.c debye.c dilog.c elementary.c ellint.c elljac.c erfc.c exp.c expint.c expint3.c fermi_dirac.c gegenbauer.c gamma.c gamma_inc.c hyperg_0F1.c hyperg_2F0.c hyperg_1F1.c hyperg_2F1.c hyperg_U.c hyperg.c laguerre.c lambert.c legendre_H3d.c legendre_Qn.c legendre_con.c legendre_poly.c log.c mathieu_angfunc.c mathieu_charv.c mathieu_coeff.c mathieu_radfunc.c mathieu_workspace.c poch.c pow_int.c psi.c recurse.h result.c shint.c sinint.c synchrotron.c transport.c trig.c zeta.c TESTS = $(check_PROGRAMS) diff --git a/specfunc/coupling.c b/specfunc/coupling.c index 6d68391..7c91342 100644 --- a/specfunc/coupling.c +++ b/specfunc/coupling.c @@ -17,7 +17,8 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ -/* Author: G. Jungman */ +/* Author: G. Jungman + * A. A. Illarionov */ #include #include @@ -105,9 +106,177 @@ m_selection_fails(int two_ja, int two_jb, int two_jc, ); } +/* Calculation of 3j-symbol is based on edge-recursion algorithm + * described by Tuzun, Comp Phys Comm 112 (1998) */ -/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ +extern int regge3j_index_permutations_table[720]; + +static inline +int +regge3j_find_ordered_params( + int twoj1, int twoj2, int twoj3, int twom1, int twom2, int twom3, + int regge[9], int * sign_param) +{ + int tmp[9] = { + (-twoj1 + twoj2 + twoj3)/2, (twoj1 - twoj2 + twoj3)/ 2, (twoj1 + twoj2 - twoj3)/2, + (twoj1 - twom1)/2, (twoj2 - twom2)/2, (twoj3 - twom3)/2, + (twoj1 + twom1)/2, (twoj2 + twom2)/2 , (twoj3 + twom3)/2 + }; + int j; + int * i = regge3j_index_permutations_table, + * ie = regge3j_index_permutations_table + 720; + + for (; i < ie ; i += 10) { + if (tmp[i[0]] > tmp[i[4]]) continue; + if (tmp[i[4]] > tmp[i[8]]) continue; + if (tmp[i[8]] > tmp[i[3]]) continue; + if (tmp[i[3]] > tmp[i[1]]) continue; + if (tmp[i[4]] > tmp[i[7]]) continue; + if ((tmp[i[4]] == tmp[i[7]]) && (tmp[i[5]] > tmp[i[8]])) continue; + + for(j = 0; j < 9 ; ++j) + regge[j] = tmp[i[j]]; + *sign_param = i[9]; + return 0; + } + return 1; +} + +static inline +int +m1_pow(int n) +{ + return 1 - 2*(n&1); +} + +static inline +void +coupling_3j_sum_e(int r[9], gsl_sf_result * result) +{ + gsl_sf_result pp, p; + double tmp; + int k; + + gsl_sf_choose_e(r[2], r[7], &pp); + pp.val *= m1_pow(r[7]); + + if (r[0] == 0) { + *result = pp; + return ; + } + + tmp = ((double)(r[7] * r[5])) / (r[2] - r[7] + 1) ; + p.val = pp.val * (r[8] - tmp); + p.err = pp.err * fabs(r[8] - tmp) + fabs(pp.val) * + (abs(r[8]) + fabs(tmp)) * 2.0 * GSL_DBL_EPSILON; + if (r[0] == 1) { + *result = p; + return; + } + + for (k = r[0] - 2; k >= 0; --k) { + tmp = (r[0] - k)*(r[4] - k); + double tmp1 = (-r[2]*(r[5] + 1) + (r[4] - k - 1)*(r[1] + k + 2) + + (r[3] + k + 1)*(r[0] - k)) / tmp; + double tmp2 = (r[3] + k + 2)*(r[1] + k + 2) / tmp; + + result->val = p.val * tmp1 - pp.val * tmp2; + result->err = p.err * tmp1 + pp.err * tmp2 + + 2.0 * GSL_DBL_EPSILON * (fabs(p.val * tmp1) + fabs(pp.val*tmp2)); + + pp = p; + p = *result; + } + + return; +} + +static inline +int +plus_log_pochhammer_ratio_e(gsl_sf_result * result, const int a, const int b, int n) +{ + double r = 1.; + int status; + gsl_sf_result f; + + int i; + for(i = n ; i > 0 ; --i) + r *= ((double)(a + i - 1))/((double)(b + i - 1)); + + status = gsl_sf_log_e(r, &f); + f.err += (n * 2) * GSL_DBL_EPSILON; + + result->val += f.val; + result->err += f.err; // here we use the fact that all values have the same sign + + return status; +} + + +static inline +int +coupling_3j_log_factor_e(int r[9], gsl_sf_result * result) +{ + int status; + + status = gsl_sf_log_e(r[0] + r[1] + r[2] + 1., result); + result->val *= -1; + + status += plus_log_pochhammer_ratio_e(result, r[8] + 1, r[6] + r[8] + 1, r[2] - r[4]); + status += plus_log_pochhammer_ratio_e(result, r[5] + 1, r[7] + 1, r[6] - r[5]); + status += plus_log_pochhammer_ratio_e(result, r[0] + 1, r[8] + r[5] + 1, r[6] - r[5]); + status += plus_log_pochhammer_ratio_e(result, 1, r[1] + 1, r[0]); + status += plus_log_pochhammer_ratio_e(result, 1, r[3] + r[6] + 1, r[0]); + + return status; +} + +static inline +void +coupling_3j_0_0_0_abs_e(int a, int b, int c, gsl_sf_result * result) +/* + * a = (twoj1 + twoj3 - twoj2) / 2 and even + * b = (twoj1 + twoj2 - twoj3) / 2 and even + * c = (twoj2 + twoj3 - twoj1) / 2 and even + * a >= b >= c + */ +{ + int a2 = a / 2, b2 = b / 2, c2 = c / 2; + int i; + double val = 1.; + + for (i = 1 ; i <= b2 ; ++i) { + val *= (b2 + i) * (a2 + i) * (a2 + i); + val /= (a + i) * (a + b2 + i) * i; + } + for (i = 1 ; i <= c2 ; ++i) { + val *= (c2 + i) * (a2 + b2 + i) * (a2 + b2 + i); + val /= (a + b + i) * (a + b + c2 + i) * i; + } + + result->val = sqrt(val / (a + b + c + 1.)); + result->err = (b2 + c2 + 1) * GSL_DBL_EPSILON * fabs(result->val); + + return; +} + +static inline +void +sort_3(int * a0, int * a1 , int * a2) +{ +#ifndef locswap +#define locswap(a, b) {int t = (b); (b) = (a); (a) = t;} + if ( *a1 > *a2 ) locswap(*a1, *a2); + if ( *a0 > *a1 ) locswap(*a0, *a1); + if ( *a1 > *a2 ) locswap(*a1, *a2); +#undef locswap +#endif +} + + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ int gsl_sf_coupling_3j_e (int two_ja, int two_jb, int two_jc, @@ -126,72 +295,75 @@ gsl_sf_coupling_3j_e (int two_ja, int two_jb, int two_jc, result->err = 0.0; return GSL_SUCCESS; } - else if ( two_ma == 0 && two_mb == 0 && two_mc == 0 - && ((two_ja + two_jb + two_jc) % 4 == 2) ) { - /* Special case for (ja jb jc; 0 0 0) = 0 when ja+jb+jc=odd */ - result->val = 0.0; - result->err = 0.0; - return GSL_SUCCESS; - } else { - int jca = (-two_ja + two_jb + two_jc) / 2, - jcb = ( two_ja - two_jb + two_jc) / 2, - jcc = ( two_ja + two_jb - two_jc) / 2, - jmma = ( two_ja - two_ma) / 2, - jmmb = ( two_jb - two_mb) / 2, - jmmc = ( two_jc - two_mc) / 2, - jpma = ( two_ja + two_ma) / 2, - jpmb = ( two_jb + two_mb) / 2, - jpmc = ( two_jc + two_mc) / 2, - jsum = ( two_ja + two_jb + two_jc) / 2, - kmin = locMax3 (0, jpmb - jmmc, jmma - jpmc), - kmax = locMin3 (jcc, jmma, jpmb), - k, sign = GSL_IS_ODD (kmin - jpma + jmmb) ? -1 : 1, - status = 0; - double sum_pos = 0.0, sum_neg = 0.0, sum_err = 0.0; - gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4, term, lnorm; - - status += gsl_sf_lnchoose_e (two_ja, jcc , &bcn1); - status += gsl_sf_lnchoose_e (two_jb, jcc , &bcn2); - status += gsl_sf_lnchoose_e (jsum+1, jcc , &bcd1); - status += gsl_sf_lnchoose_e (two_ja, jmma, &bcd2); - status += gsl_sf_lnchoose_e (two_jb, jmmb, &bcd3); - status += gsl_sf_lnchoose_e (two_jc, jpmc, &bcd4); - - lnorm.val = 0.5 * (bcn1.val + bcn2.val - bcd1.val - bcd2.val - bcd3.val - - bcd4.val - log(two_jc + 1.0)); - lnorm.err = 0.5 * (bcn1.err + bcn2.err + bcd1.err + bcd2.err + bcd3.err - + bcd4.err + GSL_DBL_EPSILON * log(two_jc + 1.0)); - - for (k = kmin; k <= kmax; k++) { - status += gsl_sf_lnchoose_e (jcc, k, &bc1); - status += gsl_sf_lnchoose_e (jcb, jmma - k, &bc2); - status += gsl_sf_lnchoose_e (jca, jpmb - k, &bc3); - status += gsl_sf_exp_err_e(bc1.val + bc2.val + bc3.val + lnorm.val, - bc1.err + bc2.err + bc3.err + lnorm.err, - &term); - - if (status != 0) { - OVERFLOW_ERROR (result); - } - - if (sign < 0) { - sum_neg += term.val; + int j_sum = (two_ja + two_jb + two_jc) / 2; + if ( two_ma == 0 && two_mb == 0 && two_mc == 0 ) { + if (j_sum % 2 != 0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; } else { - sum_pos += term.val; + int a = (two_ja + two_jb - two_jc) / 2, + b = (two_jb + two_jc - two_ja) / 2, + c = (two_jc + two_ja - two_jb) / 2; + sort_3(&c, &b, &a); + + coupling_3j_0_0_0_abs_e(a, b, c, result); + result->val *= m1_pow(j_sum/2); + return GSL_SUCCESS; } + } else { + int regge[9] = {0}, sign_param; + + regge3j_find_ordered_params(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc, + regge, &sign_param); + + if (regge[3] == regge[6] && regge[4] == regge[7] && regge[5] == regge[8]) { + if (j_sum % 2 != 0) { + result->val = 0.; + result->err = 0.; + return GSL_SUCCESS; + } else { + coupling_3j_0_0_0_abs_e(regge[1], regge[2], regge[0], result); + result->val *= m1_pow(j_sum/2 + sign_param * + (regge[0] + regge[1] + regge[2])); + return GSL_SUCCESS; + } + } else { + int status; + gsl_sf_result factor; - sum_err += term.err; + /* we have to compute factor anyway for proper error estimation */ + status = coupling_3j_log_factor_e(regge, &factor); + factor.val *= 0.5; + factor.err *= 0.5; - sign = -sign; - } - - result->val = sum_pos - sum_neg; - result->err = sum_err; - result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg); - result->err += 2.0 * GSL_DBL_EPSILON * (kmax - kmin) * fabs(result->val); + coupling_3j_sum_e(regge, result); - return GSL_SUCCESS; + if (result->val == 0.) { + result->err *= exp(factor.val); + return status; + } else { + double rel_err = result->err / fabs(result->val); + + sign_param = GSL_SIGN(result->val) * m1_pow((regge[1] - regge[8]) + + (regge[0] + regge[1] + regge[2]) * sign_param); + /* now contains the sign */ + + status += gsl_sf_log_abs_e(result->val, result); + result->err += rel_err + factor.err + 2.0 * GSL_DBL_EPSILON * + (fabs(result->val) + fabs(factor.val)); + result->val += factor.val; + + result->val = exp(result->val); + result->err *= result->val; + + result->val *= sign_param; + + return status; + } + } + } } } diff --git a/specfunc/coupling3j_regge.c b/specfunc/coupling3j_regge.c new file mode 100644 index 0000000..439bf9d --- /dev/null +++ b/specfunc/coupling3j_regge.c @@ -0,0 +1,76 @@ +/* automatically generated */ + +extern int regge3j_index_permutations_table [720] = { +0, 1, 2, 3, 4, 5, 6, 7, 8, 0, +0, 3, 6, 1, 4, 7, 2, 5, 8, 0, +0, 1, 2, 6, 7, 8, 3, 4, 5, 1, +0, 3, 6, 2, 5, 8, 1, 4, 7, 1, +3, 4, 5, 0, 1, 2, 6, 7, 8, 1, +1, 4, 7, 0, 3, 6, 2, 5, 8, 1, +3, 4, 5, 6, 7, 8, 0, 1, 2, 0, +1, 4, 7, 2, 5, 8, 0, 3, 6, 0, +6, 7, 8, 0, 1, 2, 3, 4, 5, 0, +2, 5, 8, 0, 3, 6, 1, 4, 7, 0, +6, 7, 8, 3, 4, 5, 0, 1, 2, 1, +2, 5, 8, 1, 4, 7, 0, 3, 6, 1, +0, 2, 1, 3, 5, 4, 6, 8, 7, 1, +0, 6, 3, 1, 7, 4, 2, 8, 5, 1, +0, 2, 1, 6, 8, 7, 3, 5, 4, 0, +0, 6, 3, 2, 8, 5, 1, 7, 4, 0, +3, 5, 4, 0, 2, 1, 6, 8, 7, 0, +1, 7, 4, 0, 6, 3, 2, 8, 5, 0, +3, 5, 4, 6, 8, 7, 0, 2, 1, 1, +1, 7, 4, 2, 8, 5, 0, 6, 3, 1, +6, 8, 7, 0, 2, 1, 3, 5, 4, 1, +2, 8, 5, 0, 6, 3, 1, 7, 4, 1, +6, 8, 7, 3, 5, 4, 0, 2, 1, 0, +2, 8, 5, 1, 7, 4, 0, 6, 3, 0, +1, 0, 2, 4, 3, 5, 7, 6, 8, 1, +3, 0, 6, 4, 1, 7, 5, 2, 8, 1, +1, 0, 2, 7, 6, 8, 4, 3, 5, 0, +3, 0, 6, 5, 2, 8, 4, 1, 7, 0, +4, 3, 5, 1, 0, 2, 7, 6, 8, 0, +4, 1, 7, 3, 0, 6, 5, 2, 8, 0, +4, 3, 5, 7, 6, 8, 1, 0, 2, 1, +4, 1, 7, 5, 2, 8, 3, 0, 6, 1, +7, 6, 8, 1, 0, 2, 4, 3, 5, 1, +5, 2, 8, 3, 0, 6, 4, 1, 7, 1, +7, 6, 8, 4, 3, 5, 1, 0, 2, 0, +5, 2, 8, 4, 1, 7, 3, 0, 6, 0, +1, 2, 0, 4, 5, 3, 7, 8, 6, 0, +3, 6, 0, 4, 7, 1, 5, 8, 2, 0, +1, 2, 0, 7, 8, 6, 4, 5, 3, 1, +3, 6, 0, 5, 8, 2, 4, 7, 1, 1, +4, 5, 3, 1, 2, 0, 7, 8, 6, 1, +4, 7, 1, 3, 6, 0, 5, 8, 2, 1, +4, 5, 3, 7, 8, 6, 1, 2, 0, 0, +4, 7, 1, 5, 8, 2, 3, 6, 0, 0, +7, 8, 6, 1, 2, 0, 4, 5, 3, 0, +5, 8, 2, 3, 6, 0, 4, 7, 1, 0, +7, 8, 6, 4, 5, 3, 1, 2, 0, 1, +5, 8, 2, 4, 7, 1, 3, 6, 0, 1, +2, 0, 1, 5, 3, 4, 8, 6, 7, 0, +6, 0, 3, 7, 1, 4, 8, 2, 5, 0, +2, 0, 1, 8, 6, 7, 5, 3, 4, 1, +6, 0, 3, 8, 2, 5, 7, 1, 4, 1, +5, 3, 4, 2, 0, 1, 8, 6, 7, 1, +7, 1, 4, 6, 0, 3, 8, 2, 5, 1, +5, 3, 4, 8, 6, 7, 2, 0, 1, 0, +7, 1, 4, 8, 2, 5, 6, 0, 3, 0, +8, 6, 7, 2, 0, 1, 5, 3, 4, 0, +8, 2, 5, 6, 0, 3, 7, 1, 4, 0, +8, 6, 7, 5, 3, 4, 2, 0, 1, 1, +8, 2, 5, 7, 1, 4, 6, 0, 3, 1, +2, 1, 0, 5, 4, 3, 8, 7, 6, 1, +6, 3, 0, 7, 4, 1, 8, 5, 2, 1, +2, 1, 0, 8, 7, 6, 5, 4, 3, 0, +6, 3, 0, 8, 5, 2, 7, 4, 1, 0, +5, 4, 3, 2, 1, 0, 8, 7, 6, 0, +7, 4, 1, 6, 3, 0, 8, 5, 2, 0, +5, 4, 3, 8, 7, 6, 2, 1, 0, 1, +7, 4, 1, 8, 5, 2, 6, 3, 0, 1, +8, 7, 6, 2, 1, 0, 5, 4, 3, 1, +8, 5, 2, 6, 3, 0, 7, 4, 1, 1, +8, 7, 6, 5, 4, 3, 2, 1, 0, 0, +8, 5, 2, 7, 4, 1, 6, 3, 0, 0 +}; -- 1.8.1.5