/* specfunc/hyperg.c

*

* Copyright (C) 1996, 1997, 1998, 1999, 2000 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 3 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 */

/* Miscellaneous implementations of use

* for evaluation of hypergeometric functions.

*/

#include <config.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_sf_exp.h>

#include <gsl/gsl_sf_gamma.h>

#include "error.h"

#include "hyperg.h"

#define SUM_LARGE (1.0e-5*GSL_DBL_MAX)

int

gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x,

gsl_sf_result * result

)

{

double an = a;

double bn = b;

double n = 1.0;

double del = 1.0;

double abs_del = 1.0;

double max_abs_del = 1.0;

double sum_val = 1.0;

double sum_err = 0.0;

while(abs_del/fabs(sum_val) > 0.25*GSL_DBL_EPSILON) {

double u, abs_u;

if(bn == 0.0) {

DOMAIN_ERROR(result);

}

if(an == 0.0) {

result->val = sum_val;

result->err = sum_err;

result->err = 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);

return GSL_SUCCESS;

}

if (n > 10000.0) {

result->val = sum_val;

result->err = sum_err;

GSL_ERROR ("hypergeometric series failed to converge", GSL_EFAILED);

}

u = x * (an/(bn*n));

abs_u = fabs(u);

if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) {

result->val = sum_val;

result->err = fabs(sum_val);

GSL_ERROR ("overflow", GSL_EOVRFLW);

}

del *= u;

sum_val = del;

if(fabs(sum_val) > SUM_LARGE) {

result->val = sum_val;

result->err = fabs(sum_val);

GSL_ERROR ("overflow", GSL_EOVRFLW);

}

abs_del = fabs(del);

max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del);

sum_err = 2.0*GSL_DBL_EPSILON*abs_del;

an = 1.0;

bn = 1.0;

n = 1.0;

}

result->val = sum_val;

result->err = sum_err;

result->err = abs_del;

result->err = 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);

return GSL_SUCCESS;

}

int

gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result)

{

if(fabs(x/b) < 1.0) {

const double u = x/b;

const double v = 1.0/(1.0-u);

const double pre = pow(v,a);

const double uv = u*v;

const double uv2 = uv*uv;

const double t1 = a*(a 1.0)/(2.0*b)*uv2;

const double t2a = a*(a 1.0)/(24.0*b*b)*uv2;

const double t2b = 12.0 16.0*(a 2.0)*uv 3.0*(a 2.0)*(a 3.0)*uv2;

const double t2 = t2a*t2b;

result->val = pre * (1.0 - t1 t2);

result->err = pre * GSL_DBL_EPSILON * (1.0 fabs(t1) fabs(t2));

result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

return GSL_SUCCESS;

}

else {

DOMAIN_ERROR(result);

}

}

int

gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x,

gsl_sf_result * result,

double * ln_multiplier

)

{

double N = floor(b); /* b = N eps */

double eps = b - N;

if(fabs(eps) < GSL_SQRT_DBL_EPSILON) {

double lnpre_val;

double lnpre_err;

gsl_sf_result M;

if(b > 1.0) {

double tmp = (1.0-b)*log(x);

gsl_sf_result lg_bm1;

gsl_sf_result lg_a;

gsl_sf_lngamma_e(b-1.0, &lg_bm1);

gsl_sf_lngamma_e(a, &lg_a);

lnpre_val = tmp x lg_bm1.val - lg_a.val;

lnpre_err = lg_bm1.err lg_a.err GSL_DBL_EPSILON * (fabs(x) fabs(tmp));

gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, -x, &M);

}

else {

gsl_sf_result lg_1mb;

gsl_sf_result lg_1pamb;

gsl_sf_lngamma_e(1.0-b, &lg_1mb);

gsl_sf_lngamma_e(1.0 a-b, &lg_1pamb);

lnpre_val = lg_1mb.val - lg_1pamb.val;

lnpre_err = lg_1mb.err lg_1pamb.err;

gsl_sf_hyperg_1F1_large_b_e(a, b, x, &M);

}

if(lnpre_val > GSL_LOG_DBL_MAX-10.0) {

result->val = M.val;

result->err = M.err;

*ln_multiplier = lnpre_val;

GSL_ERROR ("overflow", GSL_EOVRFLW);

}

else {

gsl_sf_result epre;

int stat_e = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &epre);

result->val = epre.val * M.val;

result->err = epre.val * M.err epre.err * fabs(M.val);

result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

*ln_multiplier = 0.0;

return stat_e;

}

}

else {

double omb_lnx = (1.0-b)*log(x);

gsl_sf_result lg_1mb; double sgn_1mb;

gsl_sf_result lg_1pamb; double sgn_1pamb;

gsl_sf_result lg_bm1; double sgn_bm1;

gsl_sf_result lg_a; double sgn_a;

gsl_sf_result M1, M2;

double lnpre1_val, lnpre2_val;

double lnpre1_err, lnpre2_err;

double sgpre1, sgpre2;

gsl_sf_hyperg_1F1_large_b_e( a, b, x, &M1);

gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, x, &M2);

gsl_sf_lngamma_sgn_e(1.0-b, &lg_1mb, &sgn_1mb);

gsl_sf_lngamma_sgn_e(1.0 a-b, &lg_1pamb, &sgn_1pamb);

gsl_sf_lngamma_sgn_e(b-1.0, &lg_bm1, &sgn_bm1);

gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);

lnpre1_val = lg_1mb.val - lg_1pamb.val;

lnpre1_err = lg_1mb.err lg_1pamb.err;

lnpre2_val = lg_bm1.val - lg_a.val - omb_lnx - x;

lnpre2_err = lg_bm1.err lg_a.err GSL_DBL_EPSILON * (fabs(omb_lnx) fabs(x));

sgpre1 = sgn_1mb * sgn_1pamb;

sgpre2 = sgn_bm1 * sgn_a;

if(lnpre1_val > GSL_LOG_DBL_MAX-10.0 || lnpre2_val > GSL_LOG_DBL_MAX-10.0) {

double max_lnpre_val = GSL_MAX(lnpre1_val,lnpre2_val);

double max_lnpre_err = GSL_MAX(lnpre1_err,lnpre2_err);

double lp1 = lnpre1_val - max_lnpre_val;

double lp2 = lnpre2_val - max_lnpre_val;

double t1 = sgpre1*exp(lp1);

double t2 = sgpre2*exp(lp2);

result->val = t1*M1.val t2*M2.val;

result->err = fabs(t1)*M1.err fabs(t2)*M2.err;

result->err = GSL_DBL_EPSILON * exp(max_lnpre_err) * (fabs(t1*M1.val) fabs(t2*M2.val));

result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

*ln_multiplier = max_lnpre_val;

GSL_ERROR ("overflow", GSL_EOVRFLW);

}

else {

double t1 = sgpre1*exp(lnpre1_val);

double t2 = sgpre2*exp(lnpre2_val);

result->val = t1*M1.val t2*M2.val;

result->err = fabs(t1) * M1.err fabs(t2)*M2.err;

result->err = GSL_DBL_EPSILON * (exp(lnpre1_err)*fabs(t1*M1.val) exp(lnpre2_err)*fabs(t2*M2.val));

result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);

*ln_multiplier = 0.0;

return GSL_SUCCESS;

}

}

}

/* [Carlson, p.109] says the error in truncating this asymptotic series

* is less than the absolute value of the first neglected term.

*

* A termination argument is provided, so that the series will

* be summed at most up to n=n_trunc. If n_trunc is set negative,

* then the series is summed until it appears to start diverging.

*/

int

gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x,

int n_trunc,

gsl_sf_result * result

)

{

const int maxiter = 2000;

double an = a;

double bn = b;

double n = 1.0;

double sum = 1.0;

double del = 1.0;

double abs_del = 1.0;

double max_abs_del = 1.0;

double last_abs_del = 1.0;

while(abs_del/fabs(sum) > GSL_DBL_EPSILON && n < maxiter) {

double u = an * (bn/n * x);

double abs_u = fabs(u);

if(abs_u > 1.0 && (max_abs_del > GSL_DBL_MAX/abs_u)) {

result->val = sum;

result->err = fabs(sum);

GSL_ERROR ("overflow", GSL_EOVRFLW);

}

del *= u;

sum = del;

abs_del = fabs(del);

if(abs_del > last_abs_del) break; /* series is probably starting to grow */

last_abs_del = abs_del;

max_abs_del = GSL_MAX(abs_del, max_abs_del);

an = 1.0;

bn = 1.0;

n = 1.0;

if(an == 0.0 || bn == 0.0) break; /* series terminated */

if(n_trunc >= 0 && n >= n_trunc) break; /* reached requested timeout */

}

result->val = sum;

result->err = GSL_DBL_EPSILON * n abs_del;

if(n >= maxiter)

GSL_ERROR ("error", GSL_EMAXITER);

else

return GSL_SUCCESS;

}