htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/coupling.c

   1 /* specfunc/coupling.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <stdlib.h>
  24 #include <gsl/gsl_math.h>
  25 #include <gsl/gsl_errno.h>
  26 #include <gsl/gsl_sf_gamma.h>
  27 #include <gsl/gsl_sf_coupling.h>
  28
  29 #include "error.h"
  30
  31 inline
  32 static
  33 int locMax3(const int a, const int b, const int c)
  34 {
  35   int d = GSL_MAX(a, b);
  36   return GSL_MAX(d, c);
  37 }
  38
  39 inline
  40 static
  41 int locMin3(const int a, const int b, const int c)
  42 {
  43   int d = GSL_MIN(a, b);
  44   return GSL_MIN(d, c);
  45 }
  46
  47 inline
  48 static
  49 int locMin5(const int a, const int b, const int c, const int d, const int e)
  50 {
  51   int f = GSL_MIN(a, b);
  52   int g = GSL_MIN(c, d);
  53   int h = GSL_MIN(f, g);
  54   return GSL_MIN(e, h);
  55 }
  56
  57
  58 /* See: [Thompson, Atlas for Computing Mathematical Functions] */
  59
  60 static
  61 int
  62 delta(int ta, int tb, int tc, gsl_sf_result * d)
  63 {
  64   gsl_sf_result f1, f2, f3, f4;
  65   int status = 0;
  66   status += gsl_sf_fact_e((ta + tb - tc)/2, &f1);
  67   status += gsl_sf_fact_e((ta + tc - tb)/2, &f2);
  68   status += gsl_sf_fact_e((tb + tc - ta)/2, &f3);
  69   status += gsl_sf_fact_e((ta + tb + tc)/2 + 1, &f4);
  70   if(status != 0) {
  71     OVERFLOW_ERROR(d);
  72   }
  73   d->val = f1.val * f2.val * f3.val / f4.val;
  74   d->err = 4.0 * GSL_DBL_EPSILON * fabs(d->val);
  75   return GSL_SUCCESS;
  76 }
  77
  78
  79 static
  80 int
  81 triangle_selection_fails(int two_ja, int two_jb, int two_jc)
  82 {
  83   return ((two_jb < abs(two_ja - two_jc)) || (two_jb > two_ja + two_jc));
  84 }
  85
  86
  87 static
  88 int
  89 m_selection_fails(int two_ja, int two_jb, int two_jc,
  90                   int two_ma, int two_mb, int two_mc)
  91 {
  92   return (
  93          abs(two_ma) > two_ja
  94       || abs(two_mb) > two_jb
  95       || abs(two_mc) > two_jc
  96       || GSL_IS_ODD(two_ja + two_ma)
  97       || GSL_IS_ODD(two_jb + two_mb)
  98       || GSL_IS_ODD(two_jc + two_mc)
  99       || (two_ma + two_mb + two_mc) != 0
 100           );
 101 }
 102
 103
 104 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 105
 106
 107 int
 108 gsl_sf_coupling_3j_e (int two_ja, int two_jb, int two_jc,
 109                       int two_ma, int two_mb, int two_mc,
 110                       gsl_sf_result * result)
 111 {
 112   /* CHECK_POINTER(result) */
 113
 114   if(two_ja < 0 || two_jb < 0 || two_jc < 0) {
 115     DOMAIN_ERROR(result);
 116   }
 117   else if (   triangle_selection_fails(two_ja, two_jb, two_jc)
 118            || m_selection_fails(two_ja, two_jb, two_jc, two_ma, two_mb, two_mc)
 119      ) {
 120     result->val = 0.0;
 121     result->err = 0.0;
 122     return GSL_SUCCESS;
 123   }
 124   else {
 125     int jca  = (-two_ja + two_jb + two_jc) / 2,
 126         jcb  = ( two_ja - two_jb + two_jc) / 2,
 127         jcc  = ( two_ja + two_jb - two_jc) / 2,
 128         jmma = ( two_ja - two_ma) / 2,
 129         jmmb = ( two_jb - two_mb) / 2,
 130         jmmc = ( two_jc - two_mc) / 2,
 131         jpma = ( two_ja + two_ma) / 2,
 132         jpmb = ( two_jb + two_mb) / 2,
 133         jpmc = ( two_jc + two_mc) / 2,
 134         jsum = ( two_ja + two_jb + two_jc) / 2,
 135         kmin = locMax3 (0, jpmb - jmmc, jmma - jpmc),
 136         kmax = locMin3 (jcc, jmma, jpmb),
 137         k, sign = GSL_IS_ODD (kmin - jpma + jmmb) ? -1 : 1,
 138         status = 0;
 139     double sum_pos = 0.0, sum_neg = 0.0, norm, term;
 140     gsl_sf_result bc1, bc2, bc3, bcn1, bcn2, bcd1, bcd2, bcd3, bcd4;
 141
 142     status += gsl_sf_choose_e (two_ja, jcc , &bcn1);
 143     status += gsl_sf_choose_e (two_jb, jcc , &bcn2);
 144     status += gsl_sf_choose_e (jsum+1, jcc , &bcd1);
 145     status += gsl_sf_choose_e (two_ja, jmma, &bcd2);
 146     status += gsl_sf_choose_e (two_jb, jmmb, &bcd3);
 147     status += gsl_sf_choose_e (two_jc, jpmc, &bcd4);
 148
 149     if (status != 0) {
 150       OVERFLOW_ERROR (result);
 151     }
 152
 153     norm = sqrt (bcn1.val * bcn2.val)
 154            / sqrt (bcd1.val * bcd2.val * bcd3.val * bcd4.val * ((double) two_jc + 1.0));
 155
 156     for (k = kmin; k <= kmax; k++) {
 157       status += gsl_sf_choose_e (jcc, k, &bc1);
 158       status += gsl_sf_choose_e (jcb, jmma - k, &bc2);
 159       status += gsl_sf_choose_e (jca, jpmb - k, &bc3);
 160
 161       if (status != 0) {
 162         OVERFLOW_ERROR (result);
 163       }
 164
 165       term = bc1.val * bc2.val * bc3.val;
 166
 167       if (sign < 0) {
 168         sum_neg += norm * term;
 169       } else {
 170         sum_pos += norm * term;
 171       }
 172
 173       sign = -sign;
 174     }
 175
 176     result->val  = sum_pos - sum_neg;
 177     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 178     result->err += 2.0 * GSL_DBL_EPSILON * (kmax - kmin) * fabs(result->val);
 179
 180     return GSL_SUCCESS;
 181   }
 182 }
 183
 184 int
 185 gsl_sf_coupling_6j_INCORRECT_e(int two_ja, int two_jb, int two_jc,
 186                                int two_jd, int two_je, int two_jf,
 187                                gsl_sf_result * result)
 188 {
 189   return gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result);
 190 }
 191
 192
 193 int
 194 gsl_sf_coupling_6j_e(int two_ja, int two_jb, int two_jc,
 195                      int two_jd, int two_je, int two_jf,
 196                      gsl_sf_result * result)
 197 {
 198   /* CHECK_POINTER(result) */
 199
 200   if(   two_ja < 0 || two_jb < 0 || two_jc < 0
 201      || two_jd < 0 || two_je < 0 || two_jf < 0
 202      ) {
 203     DOMAIN_ERROR(result);
 204   }
 205   else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
 206           || triangle_selection_fails(two_ja, two_je, two_jf)
 207           || triangle_selection_fails(two_jb, two_jd, two_jf)
 208           || triangle_selection_fails(two_je, two_jd, two_jc)
 209      ) {
 210     result->val = 0.0;
 211     result->err = 0.0;
 212     return GSL_SUCCESS;
 213   }
 214   else {
 215     gsl_sf_result n1;
 216     gsl_sf_result d1, d2, d3, d4, d5, d6;
 217     double norm;
 218     int tk, tkmin, tkmax;
 219     double phase;
 220     double sum_pos = 0.0;
 221     double sum_neg = 0.0;
 222     double sumsq_err = 0.0;
 223     int status = 0;
 224     status += delta(two_ja, two_jb, two_jc, &d1);
 225     status += delta(two_ja, two_je, two_jf, &d2);
 226     status += delta(two_jb, two_jd, two_jf, &d3);
 227     status += delta(two_je, two_jd, two_jc, &d4);
 228     if(status != GSL_SUCCESS) {
 229       OVERFLOW_ERROR(result);
 230     }
 231     norm = sqrt(d1.val) * sqrt(d2.val) * sqrt(d3.val) * sqrt(d4.val);
 232
 233     tkmin = locMax3(0,
 234                    two_ja + two_jd - two_jc - two_jf,
 235                    two_jb + two_je - two_jc - two_jf);
 236
 237     tkmax = locMin5(two_ja + two_jb + two_je + two_jd + 2,
 238                     two_ja + two_jb - two_jc,
 239                     two_je + two_jd - two_jc,
 240                     two_ja + two_je - two_jf,
 241                     two_jb + two_jd - two_jf);
 242
 243     phase = GSL_IS_ODD((two_ja + two_jb + two_je + two_jd + tkmin)/2)
 244             ? -1.0
 245             :  1.0;
 246
 247     for(tk=tkmin; tk<=tkmax; tk += 2) {
 248       double term;
 249       double term_err;
 250       gsl_sf_result den_1, den_2;
 251       gsl_sf_result d1_a, d1_b;
 252       status = 0;
 253
 254       status += gsl_sf_fact_e((two_ja + two_jb + two_je + two_jd - tk)/2 + 1, &n1);
 255       status += gsl_sf_fact_e(tk/2, &d1_a);
 256       status += gsl_sf_fact_e((two_jc + two_jf - two_ja - two_jd + tk)/2, &d1_b);
 257       status += gsl_sf_fact_e((two_jc + two_jf - two_jb - two_je + tk)/2, &d2);
 258       status += gsl_sf_fact_e((two_ja + two_jb - two_jc - tk)/2, &d3);
 259       status += gsl_sf_fact_e((two_je + two_jd - two_jc - tk)/2, &d4);
 260       status += gsl_sf_fact_e((two_ja + two_je - two_jf - tk)/2, &d5);
 261       status += gsl_sf_fact_e((two_jb + two_jd - two_jf - tk)/2, &d6);
 262
 263       if(status != GSL_SUCCESS) {
 264         OVERFLOW_ERROR(result);
 265       }
 266
 267       d1.val = d1_a.val * d1_b.val;
 268       d1.err = d1_a.err * fabs(d1_b.val) + fabs(d1_a.val) * d1_b.err;
 269
 270       den_1.val  = d1.val*d2.val*d3.val;
 271       den_1.err  = d1.err * fabs(d2.val*d3.val);
 272       den_1.err += d2.err * fabs(d1.val*d3.val);
 273       den_1.err += d3.err * fabs(d1.val*d2.val);
 274
 275       den_2.val  = d4.val*d5.val*d6.val;
 276       den_2.err  = d4.err * fabs(d5.val*d6.val);
 277       den_2.err += d5.err * fabs(d4.val*d6.val);
 278       den_2.err += d6.err * fabs(d4.val*d5.val);
 279
 280       term  = phase * n1.val / den_1.val / den_2.val;
 281       phase = -phase;
 282       term_err  = n1.err / fabs(den_1.val) / fabs(den_2.val);
 283       term_err += fabs(term / den_1.val) * den_1.err;
 284       term_err += fabs(term / den_2.val) * den_2.err;
 285
 286       if(term >= 0.0) {
 287         sum_pos += norm*term;
 288       }
 289       else {
 290         sum_neg -= norm*term;
 291       }
 292
 293       sumsq_err += norm*norm * term_err*term_err;
 294     }
 295
 296     result->val  = sum_pos - sum_neg;
 297     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 298     result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
 299     result->err += 2.0 * GSL_DBL_EPSILON * (tkmax - tkmin + 2.0) * fabs(result->val);
 300
 301     return GSL_SUCCESS;
 302   }
 303 }
 304
 305
 306 int
 307 gsl_sf_coupling_RacahW_e(int two_ja, int two_jb, int two_jc,
 308                          int two_jd, int two_je, int two_jf,
 309                          gsl_sf_result * result)
 310 {
 311   int status = gsl_sf_coupling_6j_e(two_ja, two_jb, two_je, two_jd, two_jc, two_jf, result);
 312   int phase_sum = (two_ja + two_jb + two_jc + two_jd)/2;
 313   result->val *= ( GSL_IS_ODD(phase_sum) ? -1.0 : 1.0 );
 314   return status;
 315 }
 316
 317
 318 int
 319 gsl_sf_coupling_9j_e(int two_ja, int two_jb, int two_jc,
 320                      int two_jd, int two_je, int two_jf,
 321                      int two_jg, int two_jh, int two_ji,
 322                      gsl_sf_result * result)
 323 {
 324   /* CHECK_POINTER(result) */
 325
 326   if(   two_ja < 0 || two_jb < 0 || two_jc < 0
 327      || two_jd < 0 || two_je < 0 || two_jf < 0
 328      || two_jg < 0 || two_jh < 0 || two_ji < 0
 329      ) {
 330     DOMAIN_ERROR(result);
 331   }
 332   else if(   triangle_selection_fails(two_ja, two_jb, two_jc)
 333           || triangle_selection_fails(two_jd, two_je, two_jf)
 334           || triangle_selection_fails(two_jg, two_jh, two_ji)
 335           || triangle_selection_fails(two_ja, two_jd, two_jg)
 336           || triangle_selection_fails(two_jb, two_je, two_jh)
 337           || triangle_selection_fails(two_jc, two_jf, two_ji)
 338      ) {
 339     result->val = 0.0;
 340     result->err = 0.0;
 341     return GSL_SUCCESS;
 342   }
 343   else {
 344     int tk;
 345     int tkmin = locMax3(abs(two_ja-two_ji), abs(two_jh-two_jd), abs(two_jb-two_jf));
 346     int tkmax = locMin3(two_ja + two_ji, two_jh + two_jd, two_jb + two_jf);
 347     double sum_pos = 0.0;
 348     double sum_neg = 0.0;
 349     double sumsq_err = 0.0;
 350     double phase;
 351     for(tk=tkmin; tk<=tkmax; tk += 2) {
 352       gsl_sf_result s1, s2, s3;
 353       double term;
 354       double term_err;
 355       int status = 0;
 356
 357       status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk,  two_jh, two_jd, two_jg,  &s1);
 358       status += gsl_sf_coupling_6j_e(two_jb, two_jf, tk,  two_jd, two_jh, two_je,  &s2);
 359       status += gsl_sf_coupling_6j_e(two_ja, two_ji, tk,  two_jf, two_jb, two_jc,  &s3);
 360
 361       if(status != GSL_SUCCESS) {
 362         OVERFLOW_ERROR(result);
 363       }
 364       term = s1.val * s2.val * s3.val;
 365       term_err  = s1.err * fabs(s2.val*s3.val);
 366       term_err += s2.err * fabs(s1.val*s3.val);
 367       term_err += s3.err * fabs(s1.val*s2.val);
 368
 369       if(term >= 0.0) {
 370         sum_pos += (tk + 1) * term;
 371       }
 372       else {
 373         sum_neg -= (tk + 1) * term;
 374       }
 375
 376       sumsq_err += ((tk+1) * term_err) * ((tk+1) * term_err);
 377     }
 378
 379     phase = GSL_IS_ODD(tkmin) ? -1.0 : 1.0;
 380
 381     result->val  = phase * (sum_pos - sum_neg);
 382     result->err  = 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
 383     result->err += sqrt(sumsq_err / (0.5*(tkmax-tkmin)+1.0));
 384     result->err += 2.0 * GSL_DBL_EPSILON * (tkmax-tkmin + 2.0) * fabs(result->val);
 385
 386     return GSL_SUCCESS;
 387   }
 388 }
 389
 390
 391 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 392
 393 #include "eval.h"
 394
 395 double gsl_sf_coupling_3j(int two_ja, int two_jb, int two_jc,
 396                           int two_ma, int two_mb, int two_mc)
 397 {
 398   EVAL_RESULT(gsl_sf_coupling_3j_e(two_ja, two_jb, two_jc,
 399                                    two_ma, two_mb, two_mc,
 400                                    &result));
 401 }
 402
 403
 404 double gsl_sf_coupling_6j_INCORRECT(int two_ja, int two_jb, int two_jc,
 405                                     int two_jd, int two_je, int two_jf)
 406 {
 407   EVAL_RESULT(gsl_sf_coupling_6j_INCORRECT_e(two_ja, two_jb, two_jc,
 408                                              two_jd, two_je, two_jf,
 409                                              &result));
 410 }
 411
 412
 413 double gsl_sf_coupling_6j(int two_ja, int two_jb, int two_jc,
 414                           int two_jd, int two_je, int two_jf)
 415 {
 416   EVAL_RESULT(gsl_sf_coupling_6j_e(two_ja, two_jb, two_jc,
 417                                    two_jd, two_je, two_jf,
 418                                    &result));
 419 }
 420
 421
 422 double gsl_sf_coupling_RacahW(int two_ja, int two_jb, int two_jc,
 423                           int two_jd, int two_je, int two_jf)
 424 {
 425   EVAL_RESULT(gsl_sf_coupling_RacahW_e(two_ja, two_jb, two_jc,
 426                                       two_jd, two_je, two_jf,
 427                                       &result));
 428 }
 429
 430
 431 double gsl_sf_coupling_9j(int two_ja, int two_jb, int two_jc,
 432                           int two_jd, int two_je, int two_jf,
 433                           int two_jg, int two_jh, int two_ji)
 434 {
 435   EVAL_RESULT(gsl_sf_coupling_9j_e(two_ja, two_jb, two_jc,
 436                                    two_jd, two_je, two_jf,
 437                                    two_jg, two_jh, two_ji,
 438                                    &result));
 439 }