| /* specfunc/coupling.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or (at |
| * your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, but |
| * WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| /* Author: G. Jungman */ |
| |
| #include <config.h> |
| #include <stdlib.h> |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_errno.h> |
| #include <gsl/gsl_sf_gamma.h> |
| #include <gsl/gsl_sf_coupling.h> |
| |
| #include "error.h" |
| |
| inline |
| static |
| int locMax3(const int a, const int b, const int c) |
| { |
| int d = GSL_MAX(a, b); |
| return GSL_MAX(d, c); |
| } |
| |
| inline |
| static |
| int locMin3(const int a, const int b, const int c) |
| { |
| int d = GSL_MIN(a, b); |
| return GSL_MIN(d, c); |
| } |
| |
| inline |
| static |
| int locMin5(const int a, const int b, const int c, const int d, const int e) |
| { |
| int f = GSL_MIN(a, b); |
| int g = GSL_MIN(c, d); |
| int h = GSL_MIN(f, g); |
| return GSL_MIN(e, h); |
| } |
| |
| |
| /* See: [Thompson, Atlas for Computing Mathematical Functions] */ |
| |
| static |
| int |
| delta(int ta, int tb, int tc, gsl_sf_result * d) |
| { |
| gsl_sf_result f1, f2, f3, f4; |
| int status = 0; |
| status += gsl_sf_fact_e((ta + tb - tc)/2, &f1); |
| status += gsl_sf_fact_e((ta + tc - tb)/2, &f2); |
| status += gsl_sf_fact_e((tb + tc - ta)/2, &f3); |
| status += gsl_sf_fact_e((ta + tb + tc)/2 + 1, &f4); |
| if(status != 0) { |
| OVERFLOW_ERROR(d); |
| } |
| d->val = f1.val * f2.val * f3.val / f4.val; |
| d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val); |
| return GSL_SUCCESS; |
| } |
| |
| |
| static |
| int |
| triangle_selection_fails(int two_ja, int two_jb, int two_jc) |
| { |
| return ((two_jb < abs(two_ja - two_jc)) || (two_jb > two_ja + two_jc)); |
| } |
| |
| |
| static |
| int |
| m_selection_fails(int two_ja, int two_jb, int two_jc, |
| int two_ma, int two_mb, int two_mc) |
| { |
| return ( |
| abs(two_ma) > two_ja |
| || abs(two_mb) > two_jb |
| || abs(two_mc) > two_jc |
| || GSL_IS_ODD(two_ja + two_ma) |
| || GSL_IS_ODD(two_jb + two_mb) |
| || GSL_IS_ODD(two_jc + two_mc) |
| || (two_ma + two_mb + two_mc) != 0 |
| ); |
| } |
| |
| |
| /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
| |
| |
| int |
| gsl_sf_coupling_3j_e (int two_ja, int two_jb, int two_jc, |
| int two_ma, int two_mb, int two_mc, |
| gsl_sf_result * result) |
| { |
| /* CHECK_POINTER(result) */ |
| |
| if(two_ja < 0 || two_jb < 0 || two_jc < 0) { |
| DOMAIN_ERROR(result); |
| } |
| else if ( triangle_selection_fails(two_ja, two_jb, two_jc) |
| || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc) |
| ) { |
| 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, norm, term; |
| gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4; |
| |
| status += gsl_sf_choose_e (two_ja, jcc , &bcn1); |
| status += gsl_sf_choose_e (two_jb, jcc , &bcn2); |
| status += gsl_sf_choose_e (jsum+1, jcc , &bcd1); |
| status += gsl_sf_choose_e (two_ja, jmma, &bcd2); |
| status += gsl_sf_choose_e (two_jb, jmmb, &bcd3); |
| status += gsl_sf_choose_e (two_jc, jpmc, &bcd4); |
| |
| if (status != 0) { |
| OVERFLOW_ERROR (result); |
| } |
| |
| norm = sqrt (bcn1.val * bcn2.val) |
| / sqrt (bcd1.val * bcd2.val * bcd3.val * bcd4.val * ((double) two_jc + 1.0)); |
| |
| for (k = kmin; k <= kmax; k++) { |
| status += gsl_sf_choose_e (jcc, k, &bc1); |
| status += gsl_sf_choose_e (jcb, jmma - k, &bc2); |
| status += gsl_sf_choose_e (jca, jpmb - k, &bc3); |
| |
| if (status != 0) { |
| OVERFLOW_ERROR (result); |
| } |
| |
| term = bc1.val * bc2.val * bc3.val; |
| |
| if (sign < 0) { |
| sum_neg += norm * term; |
| } else { |
| sum_pos += norm * term; |
| } |
| |
| sign = -sign; |
| } |
| |
| result->val = sum_pos - sum_neg; |
| result->err = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg); |
| result->err += 2.0 * GSL_DBL_EPSILON * (kmax - kmin) * fabs(result->val); |
| |
| return GSL_SUCCESS; |
| } |
| } |
| |
| int |
| gsl_sf_coupling_6j_INCORRECT_e(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf, |
| gsl_sf_result * result) |
| { |
| return gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result); |
| } |
| |
| |
| int |
| gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf, |
| gsl_sf_result * result) |
| { |
| /* CHECK_POINTER(result) */ |
| |
| if( two_ja < 0 || two_jb < 0 || two_jc < 0 |
| || two_jd < 0 || two_je < 0 || two_jf < 0 |
| ) { |
| DOMAIN_ERROR(result); |
| } |
| else if( triangle_selection_fails(two_ja, two_jb, two_jc) |
| || 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) |
| ) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else { |
| gsl_sf_result n1; |
| gsl_sf_result d1, d2, d3, d4, d5, d6; |
| double norm; |
| int tk, tkmin, tkmax; |
| double phase; |
| double sum_pos = 0.0; |
| double sum_neg = 0.0; |
| double sumsq_err = 0.0; |
| int status = 0; |
| status += delta(two_ja, two_jb, two_jc, &d1); |
| status += delta(two_ja, two_je, two_jf, &d2); |
| status += delta(two_jb, two_jd, two_jf, &d3); |
| status += delta(two_je, two_jd, two_jc, &d4); |
| if(status != GSL_SUCCESS) { |
| OVERFLOW_ERROR(result); |
| } |
| norm = sqrt(d1.val) * sqrt(d2.val) * sqrt(d3.val) * sqrt(d4.val); |
| |
| tkmin = locMax3(0, |
| two_ja + two_jd - two_jc - two_jf, |
| two_jb + two_je - two_jc - two_jf); |
| |
| tkmax = locMin5(two_ja + two_jb + two_je + two_jd + 2, |
| two_ja + two_jb - two_jc, |
| two_je + two_jd - two_jc, |
| two_ja + two_je - two_jf, |
| two_jb + two_jd - two_jf); |
| |
| phase = GSL_IS_ODD((two_ja + two_jb + two_je + two_jd + tkmin)/2) |
| ? -1.0 |
| : 1.0; |
| |
| for(tk=tkmin; tk<=tkmax; tk += 2) { |
| double term; |
| double term_err; |
| gsl_sf_result den_1, den_2; |
| gsl_sf_result d1_a, d1_b; |
| status = 0; |
| |
| status += gsl_sf_fact_e((two_ja + two_jb + two_je + two_jd - tk)/2 + 1, &n1); |
| status += gsl_sf_fact_e(tk/2, &d1_a); |
| status += gsl_sf_fact_e((two_jc + two_jf - two_ja - two_jd + tk)/2, &d1_b); |
| status += gsl_sf_fact_e((two_jc + two_jf - two_jb - two_je + tk)/2, &d2); |
| status += gsl_sf_fact_e((two_ja + two_jb - two_jc - tk)/2, &d3); |
| status += gsl_sf_fact_e((two_je + two_jd - two_jc - tk)/2, &d4); |
| status += gsl_sf_fact_e((two_ja + two_je - two_jf - tk)/2, &d5); |
| status += gsl_sf_fact_e((two_jb + two_jd - two_jf - tk)/2, &d6); |
| |
| if(status != GSL_SUCCESS) { |
| OVERFLOW_ERROR(result); |
| } |
| |
| d1.val = d1_a.val * d1_b.val; |
| d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err; |
| |
| den_1.val = d1.val*d2.val*d3.val; |
| den_1.err = d1.err * fabs(d2.val*d3.val); |
| den_1.err += d2.err * fabs(d1.val*d3.val); |
| den_1.err += d3.err * fabs(d1.val*d2.val); |
| |
| den_2.val = d4.val*d5.val*d6.val; |
| den_2.err = d4.err * fabs(d5.val*d6.val); |
| den_2.err += d5.err * fabs(d4.val*d6.val); |
| den_2.err += d6.err * fabs(d4.val*d5.val); |
| |
| term = phase * n1.val / den_1.val / den_2.val; |
| phase = -phase; |
| term_err = n1.err / fabs(den_1.val) / fabs(den_2.val); |
| term_err += fabs(term / den_1.val) * den_1.err; |
| term_err += fabs(term / den_2.val) * den_2.err; |
| |
| if(term >= 0.0) { |
| sum_pos += norm*term; |
| } |
| else { |
| sum_neg -= norm*term; |
| } |
| |
| sumsq_err += norm*norm * term_err*term_err; |
| } |
| |
| result->val = sum_pos - sum_neg; |
| result->err = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg); |
| result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0)); |
| result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val); |
| |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| int |
| gsl_sf_coupling_RacahW_e(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf, |
| gsl_sf_result * result) |
| { |
| int status = gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result); |
| int phase_sum = (two_ja + two_jb + two_jc + two_jd)/2; |
| result->val *= ( GSL_IS_ODD(phase_sum) ? -1.0 : 1.0 ); |
| return status; |
| } |
| |
| |
| int |
| gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf, |
| int two_jg, int two_jh, int two_ji, |
| gsl_sf_result * result) |
| { |
| /* CHECK_POINTER(result) */ |
| |
| if( two_ja < 0 || two_jb < 0 || two_jc < 0 |
| || two_jd < 0 || two_je < 0 || two_jf < 0 |
| || two_jg < 0 || two_jh < 0 || two_ji < 0 |
| ) { |
| DOMAIN_ERROR(result); |
| } |
| else if( triangle_selection_fails(two_ja, two_jb, two_jc) |
| || triangle_selection_fails(two_jd, two_je, two_jf) |
| || triangle_selection_fails(two_jg, two_jh, two_ji) |
| || 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) |
| ) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else { |
| int tk; |
| int tkmin = locMax3(abs(two_ja-two_ji), abs(two_jh-two_jd), abs(two_jb-two_jf)); |
| int tkmax = locMin3(two_ja + two_ji, two_jh + two_jd, two_jb + two_jf); |
| double sum_pos = 0.0; |
| double sum_neg = 0.0; |
| double sumsq_err = 0.0; |
| double phase; |
| for(tk=tkmin; tk<=tkmax; tk += 2) { |
| gsl_sf_result s1, s2, s3; |
| double term; |
| double term_err; |
| int status = 0; |
| |
| status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk, two_jh, two_jd, two_jg, &s1); |
| status += gsl_sf_coupling_6j_e(two_jb, two_jf, tk, two_jd, two_jh, two_je, &s2); |
| status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk, two_jf, two_jb, two_jc, &s3); |
| |
| if(status != GSL_SUCCESS) { |
| OVERFLOW_ERROR(result); |
| } |
| term = s1.val * s2.val * s3.val; |
| term_err = s1.err * fabs(s2.val*s3.val); |
| term_err += s2.err * fabs(s1.val*s3.val); |
| term_err += s3.err * fabs(s1.val*s2.val); |
| |
| if(term >= 0.0) { |
| sum_pos += (tk + 1) * term; |
| } |
| else { |
| sum_neg -= (tk + 1) * term; |
| } |
| |
| sumsq_err += ((tk+1) * term_err) * ((tk+1) * term_err); |
| } |
| |
| phase = GSL_IS_ODD(tkmin) ? -1.0 : 1.0; |
| |
| result->val = phase * (sum_pos - sum_neg); |
| result->err = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg); |
| result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0)); |
| result->err += 2.0 * GSL_DBL_EPSILON * (tkmax-tkmin + 2.0) * fabs(result->val); |
| |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
| |
| #include "eval.h" |
| |
| double gsl_sf_coupling_3j(int two_ja, int two_jb, int two_jc, |
| int two_ma, int two_mb, int two_mc) |
| { |
| EVAL_RESULT(gsl_sf_coupling_3j_e(two_ja, two_jb, two_jc, |
| two_ma, two_mb, two_mc, |
| &result)); |
| } |
| |
| |
| double gsl_sf_coupling_6j_INCORRECT(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf) |
| { |
| EVAL_RESULT(gsl_sf_coupling_6j_INCORRECT_e(two_ja, two_jb, two_jc, |
| two_jd, two_je, two_jf, |
| &result)); |
| } |
| |
| |
| double gsl_sf_coupling_6j(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf) |
| { |
| EVAL_RESULT(gsl_sf_coupling_6j_e(two_ja, two_jb, two_jc, |
| two_jd, two_je, two_jf, |
| &result)); |
| } |
| |
| |
| double gsl_sf_coupling_RacahW(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf) |
| { |
| EVAL_RESULT(gsl_sf_coupling_RacahW_e(two_ja, two_jb, two_jc, |
| two_jd, two_je, two_jf, |
| &result)); |
| } |
| |
| |
| double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc, |
| int two_jd, int two_je, int two_jf, |
| int two_jg, int two_jh, int two_ji) |
| { |
| EVAL_RESULT(gsl_sf_coupling_9j_e(two_ja, two_jb, two_jc, |
| two_jd, two_je, two_jf, |
| two_jg, two_jh, two_ji, |
| &result)); |
| } |
