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

   1 /* specfunc/hyperg.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 /* Miscellaneous implementations of use
  23  * for evaluation of hypergeometric functions.
  24  */
  25 #include <config.h>
  26 #include <gsl/gsl_math.h>
  27 #include <gsl/gsl_errno.h>
  28 #include <gsl/gsl_sf_exp.h>
  29 #include <gsl/gsl_sf_gamma.h>
  30
  31 #include "error.h"
  32 #include "hyperg.h"
  33
  34 #define SUM_LARGE  (1.0e-5*GSL_DBL_MAX)
  35
  36
  37 int
  38 gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x,
  39                               gsl_sf_result * result
  40                               )
  41 {
  42   double an  = a;
  43   double bn  = b;
  44   double n   = 1.0;
  45   double del = 1.0;
  46   double abs_del = 1.0;
  47   double max_abs_del = 1.0;
  48   double sum_val = 1.0;
  49   double sum_err = 0.0;
  50
  51   while(abs_del/fabs(sum_val) > 0.25*GSL_DBL_EPSILON) {
  52     double u, abs_u;
  53
  54     if(bn == 0.0) {
  55       DOMAIN_ERROR(result);
  56     }
  57
  58     if(an == 0.0) {
  59       result->val  = sum_val;
  60       result->err  = sum_err;
  61       result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);
  62       return GSL_SUCCESS;
  63     }
  64
  65     if (n > 10000.0) {
  66       result->val  = sum_val;
  67       result->err  = sum_err;
  68       GSL_ERROR ("hypergeometric series failed to converge", GSL_EFAILED);
  69     }
  70
  71     u = x * (an/(bn*n));
  72     abs_u = fabs(u);
  73     if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) {
  74       result->val = sum_val;
  75       result->err = fabs(sum_val);
  76       GSL_ERROR ("overflow", GSL_EOVRFLW);
  77     }
  78     del *= u;
  79     sum_val += del;
  80     if(fabs(sum_val) > SUM_LARGE) {
  81       result->val = sum_val;
  82       result->err = fabs(sum_val);
  83       GSL_ERROR ("overflow", GSL_EOVRFLW);
  84     }
  85
  86     abs_del = fabs(del);
  87     max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del);
  88     sum_err += 2.0*GSL_DBL_EPSILON*abs_del;
  89
  90     an += 1.0;
  91     bn += 1.0;
  92     n  += 1.0;
  93   }
  94
  95   result->val  = sum_val;
  96   result->err  = sum_err;
  97   result->err += abs_del;
  98   result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);
  99
 100   return GSL_SUCCESS;
 101 }
 102
 103
 104 int
 105 gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result)
 106 {
 107   if(fabs(x/b) < 1.0) {
 108     const double u = x/b;
 109     const double v = 1.0/(1.0-u);
 110     const double pre = pow(v,a);
 111     const double uv  = u*v;
 112     const double uv2 = uv*uv;
 113     const double t1  = a*(a+1.0)/(2.0*b)*uv2;
 114     const double t2a = a*(a+1.0)/(24.0*b*b)*uv2;
 115     const double t2b = 12.0 + 16.0*(a+2.0)*uv + 3.0*(a+2.0)*(a+3.0)*uv2;
 116     const double t2  = t2a*t2b;
 117     result->val  = pre * (1.0 - t1 + t2);
 118     result->err  = pre * GSL_DBL_EPSILON * (1.0 + fabs(t1) + fabs(t2));
 119     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 120     return GSL_SUCCESS;
 121   }
 122   else {
 123     DOMAIN_ERROR(result);
 124   }
 125 }
 126
 127
 128 int
 129 gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x,
 130                              gsl_sf_result * result,
 131                              double * ln_multiplier
 132                              )
 133 {
 134   double N   = floor(b);  /* b = N + eps */
 135   double eps = b - N;
 136
 137   if(fabs(eps) < GSL_SQRT_DBL_EPSILON) {
 138     double lnpre_val;
 139     double lnpre_err;
 140     gsl_sf_result M;
 141     if(b > 1.0) {
 142       double tmp = (1.0-b)*log(x);
 143       gsl_sf_result lg_bm1;
 144       gsl_sf_result lg_a;
 145       gsl_sf_lngamma_e(b-1.0, &lg_bm1);
 146       gsl_sf_lngamma_e(a, &lg_a);
 147       lnpre_val = tmp + x + lg_bm1.val - lg_a.val;
 148       lnpre_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(x) + fabs(tmp));
 149       gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, -x, &M);
 150     }
 151     else {
 152       gsl_sf_result lg_1mb;
 153       gsl_sf_result lg_1pamb;
 154       gsl_sf_lngamma_e(1.0-b, &lg_1mb);
 155       gsl_sf_lngamma_e(1.0+a-b, &lg_1pamb);
 156       lnpre_val = lg_1mb.val - lg_1pamb.val;
 157       lnpre_err = lg_1mb.err + lg_1pamb.err;
 158       gsl_sf_hyperg_1F1_large_b_e(a, b, x, &M);
 159     }
 160
 161     if(lnpre_val > GSL_LOG_DBL_MAX-10.0) {
 162       result->val  = M.val;
 163       result->err  = M.err;
 164       *ln_multiplier = lnpre_val;
 165       GSL_ERROR ("overflow", GSL_EOVRFLW);
 166     }
 167     else {
 168       gsl_sf_result epre;
 169       int stat_e = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &epre);
 170       result->val  = epre.val * M.val;
 171       result->err  = epre.val * M.err + epre.err * fabs(M.val);
 172       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 173       *ln_multiplier = 0.0;
 174       return stat_e;
 175     }
 176   }
 177   else {
 178     double omb_lnx = (1.0-b)*log(x);
 179     gsl_sf_result lg_1mb;    double sgn_1mb;
 180     gsl_sf_result lg_1pamb;  double sgn_1pamb;
 181     gsl_sf_result lg_bm1;    double sgn_bm1;
 182     gsl_sf_result lg_a;      double sgn_a;
 183     gsl_sf_result M1, M2;
 184     double lnpre1_val, lnpre2_val;
 185     double lnpre1_err, lnpre2_err;
 186     double sgpre1, sgpre2;
 187     gsl_sf_hyperg_1F1_large_b_e(    a,     b, x, &M1);
 188     gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, x, &M2);
 189
 190     gsl_sf_lngamma_sgn_e(1.0-b,   &lg_1mb,   &sgn_1mb);
 191     gsl_sf_lngamma_sgn_e(1.0+a-b, &lg_1pamb, &sgn_1pamb);
 192
 193     gsl_sf_lngamma_sgn_e(b-1.0, &lg_bm1, &sgn_bm1);
 194     gsl_sf_lngamma_sgn_e(a,     &lg_a,   &sgn_a);
 195
 196     lnpre1_val = lg_1mb.val - lg_1pamb.val;
 197     lnpre1_err = lg_1mb.err + lg_1pamb.err;
 198     lnpre2_val = lg_bm1.val - lg_a.val - omb_lnx - x;
 199     lnpre2_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(omb_lnx)+fabs(x));
 200     sgpre1 = sgn_1mb * sgn_1pamb;
 201     sgpre2 = sgn_bm1 * sgn_a;
 202
 203     if(lnpre1_val > GSL_LOG_DBL_MAX-10.0 || lnpre2_val > GSL_LOG_DBL_MAX-10.0) {
 204       double max_lnpre_val = GSL_MAX(lnpre1_val,lnpre2_val);
 205       double max_lnpre_err = GSL_MAX(lnpre1_err,lnpre2_err);
 206       double lp1 = lnpre1_val - max_lnpre_val;
 207       double lp2 = lnpre2_val - max_lnpre_val;
 208       double t1  = sgpre1*exp(lp1);
 209       double t2  = sgpre2*exp(lp2);
 210       result->val  = t1*M1.val + t2*M2.val;
 211       result->err  = fabs(t1)*M1.err + fabs(t2)*M2.err;
 212       result->err += GSL_DBL_EPSILON * exp(max_lnpre_err) * (fabs(t1*M1.val) + fabs(t2*M2.val));
 213       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 214       *ln_multiplier = max_lnpre_val;
 215       GSL_ERROR ("overflow", GSL_EOVRFLW);
 216     }
 217     else {
 218       double t1 = sgpre1*exp(lnpre1_val);
 219       double t2 = sgpre2*exp(lnpre2_val);
 220       result->val  = t1*M1.val + t2*M2.val;
 221       result->err  = fabs(t1) * M1.err + fabs(t2)*M2.err;
 222       result->err += GSL_DBL_EPSILON * (exp(lnpre1_err)*fabs(t1*M1.val) + exp(lnpre2_err)*fabs(t2*M2.val));
 223       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 224       *ln_multiplier = 0.0;
 225       return GSL_SUCCESS;
 226     }
 227   }
 228 }
 229
 230
 231
 232 /* [Carlson, p.109] says the error in truncating this asymptotic series
 233  * is less than the absolute value of the first neglected term.
 234  *
 235  * A termination argument is provided, so that the series will
 236  * be summed at most up to n=n_trunc. If n_trunc is set negative,
 237  * then the series is summed until it appears to start diverging.
 238  */
 239 int
 240 gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x,
 241                               int n_trunc,
 242                               gsl_sf_result * result
 243                               )
 244 {
 245   const int maxiter = 2000;
 246   double an = a;
 247   double bn = b;
 248   double n   = 1.0;
 249   double sum = 1.0;
 250   double del = 1.0;
 251   double abs_del = 1.0;
 252   double max_abs_del = 1.0;
 253   double last_abs_del = 1.0;
 254
 255   while(abs_del/fabs(sum) > GSL_DBL_EPSILON && n < maxiter) {
 256
 257     double u = an * (bn/n * x);
 258     double abs_u = fabs(u);
 259
 260     if(abs_u > 1.0 && (max_abs_del > GSL_DBL_MAX/abs_u)) {
 261       result->val = sum;
 262       result->err = fabs(sum);
 263       GSL_ERROR ("overflow", GSL_EOVRFLW);
 264     }
 265
 266     del *= u;
 267     sum += del;
 268
 269     abs_del = fabs(del);
 270
 271     if(abs_del > last_abs_del) break; /* series is probably starting to grow */
 272
 273     last_abs_del = abs_del;
 274     max_abs_del  = GSL_MAX(abs_del, max_abs_del);
 275
 276     an += 1.0;
 277     bn += 1.0;
 278     n  += 1.0;
 279
 280     if(an == 0.0 || bn == 0.0) break;        /* series terminated */
 281
 282     if(n_trunc >= 0 && n >= n_trunc) break;  /* reached requested timeout */
 283   }
 284
 285   result->val = sum;
 286   result->err = GSL_DBL_EPSILON * n + abs_del;
 287   if(n >= maxiter)
 288     GSL_ERROR ("error", GSL_EMAXITER);
 289   else
 290     return GSL_SUCCESS;
 291 }