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

   1 /* specfunc/hyperg_0F1.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 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_gamma.h>
  27 #include <gsl/gsl_sf_bessel.h>
  28 #include <gsl/gsl_sf_hyperg.h>
  29
  30 #include "error.h"
  31
  32 #define locEPS  (1000.0*GSL_DBL_EPSILON)
  33
  34
  35 /* Evaluate bessel_I(nu, x), allowing nu < 0.
  36  * This is fine here because we do not not allow
  37  * nu to be a negative integer.
  38  * x > 0.
  39  */
  40 static
  41 int
  42 hyperg_0F1_bessel_I(const double nu, const double x, gsl_sf_result * result)
  43 {
  44   if(x > GSL_LOG_DBL_MAX) {
  45     OVERFLOW_ERROR(result);
  46   }
  47
  48   if(nu < 0.0) {
  49     const double anu = -nu;
  50     const double s   = 2.0/M_PI * sin(anu*M_PI);
  51     const double ex  = exp(x);
  52     gsl_sf_result I;
  53     gsl_sf_result K;
  54     int stat_I = gsl_sf_bessel_Inu_scaled_e(anu, x, &I);
  55     int stat_K = gsl_sf_bessel_Knu_scaled_e(anu, x, &K);
  56     result->val  = ex * I.val + s * (K.val / ex);
  57     result->err  = ex * I.err + fabs(s * K.err/ex);
  58     result->err += fabs(s * (K.val/ex)) * GSL_DBL_EPSILON * anu * M_PI;
  59     return GSL_ERROR_SELECT_2(stat_K, stat_I);
  60   }
  61   else {
  62     const double ex  = exp(x);
  63     gsl_sf_result I;
  64     int stat_I = gsl_sf_bessel_Inu_scaled_e(nu, x, &I);
  65     result->val = ex * I.val;
  66     result->err = ex * I.err + GSL_DBL_EPSILON * fabs(result->val);
  67     return stat_I;
  68   }
  69 }
  70
  71
  72 /* Evaluate bessel_J(nu, x), allowing nu < 0.
  73  * This is fine here because we do not not allow
  74  * nu to be a negative integer.
  75  * x > 0.
  76  */
  77 static
  78 int
  79 hyperg_0F1_bessel_J(const double nu, const double x, gsl_sf_result * result)
  80 {
  81   if(nu < 0.0) {
  82     const double anu = -nu;
  83     const double s   = sin(anu*M_PI);
  84     const double c   = cos(anu*M_PI);
  85     gsl_sf_result J;
  86     gsl_sf_result Y;
  87     int stat_J = gsl_sf_bessel_Jnu_e(anu, x, &J);
  88     int stat_Y = gsl_sf_bessel_Ynu_e(anu, x, &Y);
  89     result->val  = c * J.val - s * Y.val;
  90     result->err  = fabs(c * J.err) + fabs(s * Y.err);
  91     result->err += fabs(anu * M_PI) * GSL_DBL_EPSILON * fabs(J.val + Y.val);
  92     return GSL_ERROR_SELECT_2(stat_Y, stat_J);
  93   }
  94   else {
  95     return gsl_sf_bessel_Jnu_e(nu, x, result);
  96   }
  97 }
  98
  99
 100 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 101
 102 int
 103 gsl_sf_hyperg_0F1_e(double c, double x, gsl_sf_result * result)
 104 {
 105   const double rintc = floor(c + 0.5);
 106   const int c_neg_integer = (c < 0.0 && fabs(c - rintc) < locEPS);
 107
 108   /* CHECK_POINTER(result) */
 109
 110   if(c == 0.0 || c_neg_integer) {
 111     DOMAIN_ERROR(result);
 112   }
 113   else if(x < 0.0) {
 114     gsl_sf_result Jcm1;
 115     gsl_sf_result lg_c;
 116     double sgn;
 117     int stat_g = gsl_sf_lngamma_sgn_e(c, &lg_c, &sgn);
 118     int stat_J = hyperg_0F1_bessel_J(c-1.0, 2.0*sqrt(-x), &Jcm1);
 119     if(stat_g != GSL_SUCCESS) {
 120       result->val = 0.0;
 121       result->err = 0.0;
 122       return stat_g;
 123     }
 124     else if(Jcm1.val == 0.0) {
 125       result->val = 0.0;
 126       result->err = 0.0;
 127       return stat_J;
 128     }
 129     else {
 130       const double tl = log(-x)*0.5*(1.0-c);
 131       double ln_pre_val = lg_c.val + tl;
 132       double ln_pre_err = lg_c.err + 2.0 * GSL_DBL_EPSILON * fabs(tl);
 133       return gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
 134                                       sgn*Jcm1.val, Jcm1.err,
 135                                       result);
 136     }
 137   }
 138   else if(x == 0.0) {
 139     result->val = 1.0;
 140     result->err = 1.0;
 141     return GSL_SUCCESS;
 142   }
 143   else {
 144     gsl_sf_result Icm1;
 145     gsl_sf_result lg_c;
 146     double sgn;
 147     int stat_g = gsl_sf_lngamma_sgn_e(c, &lg_c, &sgn);
 148     int stat_I = hyperg_0F1_bessel_I(c-1.0, 2.0*sqrt(x), &Icm1);
 149     if(stat_g != GSL_SUCCESS) {
 150       result->val = 0.0;
 151       result->err = 0.0;
 152       return stat_g;
 153     }
 154     else if(Icm1.val == 0.0) {
 155       result->val = 0.0;
 156       result->err = 0.0;
 157       return stat_I;
 158     }
 159     else {
 160       const double tl = log(x)*0.5*(1.0-c);
 161       const double ln_pre_val = lg_c.val + tl;
 162       const double ln_pre_err = lg_c.err + 2.0 * GSL_DBL_EPSILON * fabs(tl);
 163       return gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
 164                                       sgn*Icm1.val, Icm1.err,
 165                                       result);
 166     }
 167   }
 168 }
 169
 170
 171 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 172
 173 #include "eval.h"
 174
 175 double gsl_sf_hyperg_0F1(const double c, const double x)
 176 {
 177   EVAL_RESULT(gsl_sf_hyperg_0F1_e(c, x, &result));
 178 }