/* 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., 675 Mass Ave, Cambridge, MA 02139, 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"inlinestaticint locMax3(const int a, const int b, const int c){  int d = GSL_MAX(a, b);  return GSL_MAX(d, c);}inlinestaticint locMin3(const int a, const int b, const int c){  int d = GSL_MIN(a, b);  return GSL_MIN(d, c);}inlinestaticint 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] */staticintdelta(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;}staticinttriangle_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));}staticintm_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 *-*-*-*-*-*-*-*-*-*-*-*/intgsl_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;  }}intgsl_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);}intgsl_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;  }}intgsl_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;}intgsl_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));}