|
Packit |
67cb25 |
/* specfunc/hyperg_1F1.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
* Copyright (C) 2010 Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Author: G. Jungman */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_elementary.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_exp.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_bessel.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_gamma.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_laguerre.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_sf_hyperg.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "error.h"
|
|
Packit |
67cb25 |
#include "hyperg.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Asymptotic result for 1F1(a, b, x) x -> -Infinity.
|
|
Packit |
67cb25 |
* Assumes b-a != neg integer and b != neg integer.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_asymp_negx(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result lg_b;
|
|
Packit |
67cb25 |
gsl_sf_result lg_bma;
|
|
Packit |
67cb25 |
double sgn_b;
|
|
Packit |
67cb25 |
double sgn_bma;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
|
|
Packit |
67cb25 |
int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) {
|
|
Packit |
67cb25 |
gsl_sf_result F;
|
|
Packit |
67cb25 |
int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F);
|
|
Packit |
67cb25 |
if(F.val != 0) {
|
|
Packit |
67cb25 |
double ln_term_val = a*log(-x);
|
|
Packit |
67cb25 |
double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val));
|
|
Packit |
67cb25 |
double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val;
|
|
Packit |
67cb25 |
double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
|
|
Packit |
67cb25 |
sgn_bma*sgn_b*F.val, F.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_F);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return stat_F;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Asymptotic result for 1F1(a, b, x) x -> +Infinity
|
|
Packit |
67cb25 |
* Assumes b != neg integer and a != neg integer
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_asymp_posx(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result lg_b;
|
|
Packit |
67cb25 |
gsl_sf_result lg_a;
|
|
Packit |
67cb25 |
double sgn_b;
|
|
Packit |
67cb25 |
double sgn_a;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
|
|
Packit |
67cb25 |
int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) {
|
|
Packit |
67cb25 |
gsl_sf_result F;
|
|
Packit |
67cb25 |
int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F);
|
|
Packit |
67cb25 |
if(stat_F == GSL_SUCCESS && F.val != 0) {
|
|
Packit |
67cb25 |
double lnx = log(x);
|
|
Packit |
67cb25 |
double ln_term_val = (a-b)*lnx;
|
|
Packit |
67cb25 |
double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx)
|
|
Packit |
67cb25 |
+ 2.0 * GSL_DBL_EPSILON * fabs(a-b);
|
|
Packit |
67cb25 |
double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x;
|
|
Packit |
67cb25 |
double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
|
|
Packit |
67cb25 |
sgn_a*sgn_b*F.val, F.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_F);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return stat_F;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Asymptotic result from Slater 4.3.7
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* To get the general series, write M(a,b,x) as
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* M(a,b,x)=sum ((a)_n/(b)_n) (x^n / n!)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* and expand (b)_n in inverse powers of b as follows
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* -log(1/(b)_n) = sum_(k=0)^(n-1) log(b+k)
|
|
Packit |
67cb25 |
* = n log(b) + sum_(k=0)^(n-1) log(1+k/b)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Do a taylor expansion of the log in 1/b and sum the resulting terms
|
|
Packit |
67cb25 |
* using the standard algebraic formulas for finite sums of powers of
|
|
Packit |
67cb25 |
* k. This should then give
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* M(a,b,x) = sum_(n=0)^(inf) (a_n/n!) (x/b)^n * (1 - n(n-1)/(2b)
|
|
Packit |
67cb25 |
* + (n-1)n(n+1)(3n-2)/(24b^2) + ...
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* which can be summed explicitly. The trick for summing it is to take
|
|
Packit |
67cb25 |
* derivatives of sum_(i=0)^(inf) a_n*y^n/n! = (1-y)^(-a);
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [BJG 16/01/2007]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_largebx(const double a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double y = x/b;
|
|
Packit |
67cb25 |
double f = exp(-a*log1p(-y));
|
|
Packit |
67cb25 |
double t1 = -((a*(a+1.0))/(2*b))*pow((y/(1.0-y)),2.0);
|
|
Packit |
67cb25 |
double t2 = (1/(24*b*b))*((a*(a+1)*y*y)/pow(1-y,4))*(12+8*(2*a+1)*y+(3*a*a-a-2)*y*y);
|
|
Packit |
67cb25 |
double t3 = (-1/(48*b*b*b*pow(1-y,6)))*a*((a + 1)*((y*((a + 1)*(a*(y*(y*((y*(a - 2) + 16)*(a - 1)) + 72)) + 96)) + 24)*pow(y, 2)));
|
|
Packit |
67cb25 |
result->val = f * (1 + t1 + t2 + t3);
|
|
Packit |
67cb25 |
result->err = 2*fabs(f*t3) + 2*GSL_DBL_EPSILON*fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Asymptotic result for x < 2b-4a, 2b-4a large.
|
|
Packit |
67cb25 |
* [Abramowitz+Stegun, 13.5.21]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* assumes 0 <= x/(2b-4a) <= 1
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double eta = 2.0*b - 4.0*a;
|
|
Packit |
67cb25 |
double cos2th = x/eta;
|
|
Packit |
67cb25 |
double sin2th = 1.0 - cos2th;
|
|
Packit |
67cb25 |
double th = acos(sqrt(cos2th));
|
|
Packit |
67cb25 |
double pre_h = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
|
|
Packit |
67cb25 |
gsl_sf_result lg_b;
|
|
Packit |
67cb25 |
int stat_lg = gsl_sf_lngamma_e(b, &lg_b);
|
|
Packit |
67cb25 |
double t1 = 0.5*(1.0-b)*log(0.25*x*eta);
|
|
Packit |
67cb25 |
double t2 = 0.25*log(pre_h);
|
|
Packit |
67cb25 |
double lnpre_val = lg_b.val + 0.5*x + t1 - t2;
|
|
Packit |
67cb25 |
double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2));
|
|
Packit |
67cb25 |
#if SMALL_ANGLE
|
|
Packit |
67cb25 |
const double eps = asin(sqrt(cos2th)); /* theta = pi/2 - eps */
|
|
Packit |
67cb25 |
double s1 = (fmod(a, 1.0) == 0.0) ? 0.0 : sin(a*M_PI);
|
|
Packit |
67cb25 |
double eta_reduc = (fmod(eta + 1, 4.0) == 0.0) ? 0.0 : fmod(eta + 1, 8.0);
|
|
Packit |
67cb25 |
double phi1 = 0.25*eta_reduc*M_PI;
|
|
Packit |
67cb25 |
double phi2 = 0.25*eta*(2*eps + sin(2.0*eps));
|
|
Packit |
67cb25 |
double s2 = sin(phi1 - phi2);
|
|
Packit |
67cb25 |
#else
|
|
Packit |
67cb25 |
double s1 = sin(a*M_PI);
|
|
Packit |
67cb25 |
double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
double ser_val = s1 + s2;
|
|
Packit |
67cb25 |
double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2));
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
|
|
Packit |
67cb25 |
ser_val, ser_err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_lg);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Luke's rational approximation.
|
|
Packit |
67cb25 |
* See [Luke, Algorithms for the Computation of Mathematical Functions, p.182]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Like the case of the 2F1 rational approximations, these are
|
|
Packit |
67cb25 |
* probably guaranteed to converge for x < 0, barring gross
|
|
Packit |
67cb25 |
* numerical instability in the pre-asymptotic regime.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_luke(const double a, const double c, const double xin,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double RECUR_BIG = 1.0e+50;
|
|
Packit |
67cb25 |
const int nmax = 5000;
|
|
Packit |
67cb25 |
int n = 3;
|
|
Packit |
67cb25 |
const double x = -xin;
|
|
Packit |
67cb25 |
const double x3 = x*x*x;
|
|
Packit |
67cb25 |
const double t0 = a/c;
|
|
Packit |
67cb25 |
const double t1 = (a+1.0)/(2.0*c);
|
|
Packit |
67cb25 |
const double t2 = (a+2.0)/(2.0*(c+1.0));
|
|
Packit |
67cb25 |
double F = 1.0;
|
|
Packit |
67cb25 |
double prec;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double Bnm3 = 1.0; /* B0 */
|
|
Packit |
67cb25 |
double Bnm2 = 1.0 + t1 * x; /* B1 */
|
|
Packit |
67cb25 |
double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x); /* B2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double Anm3 = 1.0; /* A0 */
|
|
Packit |
67cb25 |
double Anm2 = Bnm2 - t0 * x; /* A1 */
|
|
Packit |
67cb25 |
double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x; /* A2 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(1) {
|
|
Packit |
67cb25 |
double npam1 = n + a - 1;
|
|
Packit |
67cb25 |
double npcm1 = n + c - 1;
|
|
Packit |
67cb25 |
double npam2 = n + a - 2;
|
|
Packit |
67cb25 |
double npcm2 = n + c - 2;
|
|
Packit |
67cb25 |
double tnm1 = 2*n - 1;
|
|
Packit |
67cb25 |
double tnm3 = 2*n - 3;
|
|
Packit |
67cb25 |
double tnm5 = 2*n - 5;
|
|
Packit |
67cb25 |
double F1 = (n-a-2) / (2*tnm3*npcm1);
|
|
Packit |
67cb25 |
double F2 = (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1);
|
|
Packit |
67cb25 |
double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
|
|
Packit |
67cb25 |
double E = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
|
|
Packit |
67cb25 |
double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
|
|
Packit |
67cb25 |
double r = An/Bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
prec = fabs((F - r)/F);
|
|
Packit |
67cb25 |
F = r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(prec < GSL_DBL_EPSILON || n > nmax) break;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
|
|
Packit |
67cb25 |
An /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm3 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm3 /= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
|
|
Packit |
67cb25 |
An *= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn *= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 *= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 *= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 *= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 *= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm3 *= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm3 *= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Bnm3 = Bnm2;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
Anm3 = Anm2;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = F;
|
|
Packit |
67cb25 |
result->err = 2.0 * fabs(F * prec);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Series for 1F1(1,b,x)
|
|
Packit |
67cb25 |
* b > 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double sum_val = 1.0;
|
|
Packit |
67cb25 |
double sum_err = 0.0;
|
|
Packit |
67cb25 |
double term = 1.0;
|
|
Packit |
67cb25 |
double n = 1.0;
|
|
Packit |
67cb25 |
while(fabs(term/sum_val) > 0.25*GSL_DBL_EPSILON) {
|
|
Packit |
67cb25 |
term *= x/(b+n-1);
|
|
Packit |
67cb25 |
sum_val += term;
|
|
Packit |
67cb25 |
sum_err += 8.0*GSL_DBL_EPSILON*fabs(term) + GSL_DBL_EPSILON*fabs(sum_val);
|
|
Packit |
67cb25 |
n += 1.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = sum_val;
|
|
Packit |
67cb25 |
result->err = sum_err;
|
|
Packit |
67cb25 |
result->err += 2.0 * fabs(term);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(1,b,x)
|
|
Packit |
67cb25 |
* b >= 1, b integer
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(b < 1) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 1) {
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 2) {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 3) {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_2_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_n_e(b-1, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(1,b,x)
|
|
Packit |
67cb25 |
* b >=1, b real
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* checked OK: [GJ] Thu Oct 1 16:46:35 MDT 1998
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_1(const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ax = fabs(x);
|
|
Packit |
67cb25 |
double ib = floor(b + 0.1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(b < 1.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 1.0) {
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b >= 1.4*ax) {
|
|
Packit |
67cb25 |
return hyperg_1F1_1_series(b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) {
|
|
Packit |
67cb25 |
return hyperg_1F1_1_int((int)ib, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.0) {
|
|
Packit |
67cb25 |
if(x > 100.0 && b < 0.75*x) {
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_posx(1.0, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b < 1.0e+05) {
|
|
Packit |
67cb25 |
/* Recurse backward on b, from a
|
|
Packit |
67cb25 |
* chosen offset point. For x > 0,
|
|
Packit |
67cb25 |
* which holds here, this should
|
|
Packit |
67cb25 |
* be a stable direction.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double off = ceil(1.4*x-b) + 1.0;
|
|
Packit |
67cb25 |
double bp = b + off;
|
|
Packit |
67cb25 |
gsl_sf_result M;
|
|
Packit |
67cb25 |
int stat_s = hyperg_1F1_1_series(bp, x, &M);
|
|
Packit |
67cb25 |
const double err_rat = M.err / fabs(M.val);
|
|
Packit |
67cb25 |
while(bp > b+0.1) {
|
|
Packit |
67cb25 |
/* M(1,b-1) = x/(b-1) M(1,b) + 1 */
|
|
Packit |
67cb25 |
bp -= 1.0;
|
|
Packit |
67cb25 |
M.val = 1.0 + x/bp * M.val;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = M.val;
|
|
Packit |
67cb25 |
result->err = err_rat * fabs(M.val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val);
|
|
Packit |
67cb25 |
return stat_s;
|
|
Packit |
67cb25 |
} else if (fabs(x) < fabs(b) && fabs(x) < sqrt(fabs(b)) * fabs(b-x)) {
|
|
Packit |
67cb25 |
return hyperg_1F1_largebx(1.0, b, x, result);
|
|
Packit |
67cb25 |
} else if (fabs(x) > fabs(b)) {
|
|
Packit |
67cb25 |
return hyperg_1F1_1_series(b, x, result);
|
|
Packit |
67cb25 |
} else {
|
|
Packit |
67cb25 |
return hyperg_1F1_large2bm4a(1.0, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x <= 0 and b not large compared to |x|
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(ax < 10.0 && b < 10.0) {
|
|
Packit |
67cb25 |
return hyperg_1F1_1_series(b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) {
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_negx(1.0, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return hyperg_1F1_luke(1.0, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(a,b,x)/Gamma(b) for b->0
|
|
Packit |
67cb25 |
* [limit of Abramowitz+Stegun 13.3.7]
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double eta = a*x;
|
|
Packit |
67cb25 |
if(eta > 0.0) {
|
|
Packit |
67cb25 |
double root_eta = sqrt(eta);
|
|
Packit |
67cb25 |
gsl_sf_result I1_scaled;
|
|
Packit |
67cb25 |
int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled);
|
|
Packit |
67cb25 |
if(I1_scaled.val <= 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Note that 13.3.7 contains higher terms which are zeroth order
|
|
Packit |
67cb25 |
in b. These make a non-negligible contribution to the sum.
|
|
Packit |
67cb25 |
With the first correction term, the I1 above is replaced by
|
|
Packit |
67cb25 |
I1 + (2/3)*a*(x/(4a))**(3/2)*I2(2*root_eta). We will add
|
|
Packit |
67cb25 |
this as part of the result and error estimate. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double corr1 =(2.0/3.0)*a*pow(x/(4.0*a),1.5)*gsl_sf_bessel_In_scaled(2, 2.0*root_eta)
|
|
Packit |
67cb25 |
;
|
|
Packit |
67cb25 |
const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(2.0*root_eta) + log(I1_scaled.val+corr1);
|
|
Packit |
67cb25 |
const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs((I1_scaled.err+corr1)/I1_scaled.val);
|
|
Packit |
67cb25 |
return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(eta == 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* eta < 0 */
|
|
Packit |
67cb25 |
double root_eta = sqrt(-eta);
|
|
Packit |
67cb25 |
gsl_sf_result J1;
|
|
Packit |
67cb25 |
int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1;;
|
|
Packit |
67cb25 |
if(J1.val <= 0.0) {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double t1 = 0.5*x;
|
|
Packit |
67cb25 |
const double t2 = 0.5*log(-eta);
|
|
Packit |
67cb25 |
const double t3 = fabs(x);
|
|
Packit |
67cb25 |
const double t4 = log(J1.val);
|
|
Packit |
67cb25 |
const double lnr_val = t1 + t2 + t3 + t4;
|
|
Packit |
67cb25 |
const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val);
|
|
Packit |
67cb25 |
gsl_sf_result ex;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex);
|
|
Packit |
67cb25 |
result->val = -ex.val;
|
|
Packit |
67cb25 |
result->err = ex.err;
|
|
Packit |
67cb25 |
return stat_e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1'(a,b,x)/1F1(a,b,x)
|
|
Packit |
67cb25 |
* Uses Gautschi's version of the CF.
|
|
Packit |
67cb25 |
* [Gautschi, Math. Comp. 31, 994 (1977)]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Supposedly this suffers from the "anomalous convergence"
|
|
Packit |
67cb25 |
* problem when b < x. I have seen anomalous convergence
|
|
Packit |
67cb25 |
* in several of the continued fractions associated with
|
|
Packit |
67cb25 |
* 1F1(a,b,x). This particular CF formulation seems stable
|
|
Packit |
67cb25 |
* for b > x. However, it does display a painful artifact
|
|
Packit |
67cb25 |
* of the anomalous convergence; the convergence plateaus
|
|
Packit |
67cb25 |
* unless b >>> x. For example, even for b=1000, x=1, this
|
|
Packit |
67cb25 |
* method locks onto a ratio which is only good to about
|
|
Packit |
67cb25 |
* 4 digits. Apparently the rest of the digits are hiding
|
|
Packit |
67cb25 |
* way out on the plateau, but finite-precision lossage
|
|
Packit |
67cb25 |
* means you will never get them.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double RECUR_BIG = GSL_SQRT_DBL_MAX;
|
|
Packit |
67cb25 |
const int maxiter = 5000;
|
|
Packit |
67cb25 |
int n = 1;
|
|
Packit |
67cb25 |
double Anm2 = 1.0;
|
|
Packit |
67cb25 |
double Bnm2 = 0.0;
|
|
Packit |
67cb25 |
double Anm1 = 0.0;
|
|
Packit |
67cb25 |
double Bnm1 = 1.0;
|
|
Packit |
67cb25 |
double a1 = 1.0;
|
|
Packit |
67cb25 |
double b1 = 1.0;
|
|
Packit |
67cb25 |
double An = b1*Anm1 + a1*Anm2;
|
|
Packit |
67cb25 |
double Bn = b1*Bnm1 + a1*Bnm2;
|
|
Packit |
67cb25 |
double an, bn;
|
|
Packit |
67cb25 |
double fn = An/Bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(n < maxiter) {
|
|
Packit |
67cb25 |
double old_fn;
|
|
Packit |
67cb25 |
double del;
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
an = (a+n)*x/((b-x+n-1)*(b-x+n));
|
|
Packit |
67cb25 |
bn = 1.0;
|
|
Packit |
67cb25 |
An = bn*Anm1 + an*Anm2;
|
|
Packit |
67cb25 |
Bn = bn*Bnm1 + an*Bnm2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
|
|
Packit |
67cb25 |
An /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
old_fn = fn;
|
|
Packit |
67cb25 |
fn = An/Bn;
|
|
Packit |
67cb25 |
del = old_fn/fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = a/(b-x) * fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(n == maxiter)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif /* 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1'(a,b,x)/1F1(a,b,x)
|
|
Packit |
67cb25 |
* Uses Gautschi's series transformation of the
|
|
Packit |
67cb25 |
* continued fraction. This is apparently the best
|
|
Packit |
67cb25 |
* method for getting this ratio in the stable region.
|
|
Packit |
67cb25 |
* The convergence is monotone and supergeometric
|
|
Packit |
67cb25 |
* when b > x.
|
|
Packit |
67cb25 |
* Assumes a >= -1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a == 0.0) {
|
|
Packit |
67cb25 |
*result = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const int maxiter = 5000;
|
|
Packit |
67cb25 |
double sum = 1.0;
|
|
Packit |
67cb25 |
double pk = 1.0;
|
|
Packit |
67cb25 |
double rhok = 0.0;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=1; k
|
|
Packit |
67cb25 |
double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
|
|
Packit |
67cb25 |
rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
|
|
Packit |
67cb25 |
pk *= rhok;
|
|
Packit |
67cb25 |
sum += pk;
|
|
Packit |
67cb25 |
if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
*result = a/(b-x) * sum;
|
|
Packit |
67cb25 |
if(k == maxiter)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(a+1,b,x)/1F1(a,b,x)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* I think this suffers from typical "anomalous convergence".
|
|
Packit |
67cb25 |
* I could not find a region where it was truly useful.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_CF1(const double a, const double b, const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double RECUR_BIG = GSL_SQRT_DBL_MAX;
|
|
Packit |
67cb25 |
const int maxiter = 5000;
|
|
Packit |
67cb25 |
int n = 1;
|
|
Packit |
67cb25 |
double Anm2 = 1.0;
|
|
Packit |
67cb25 |
double Bnm2 = 0.0;
|
|
Packit |
67cb25 |
double Anm1 = 0.0;
|
|
Packit |
67cb25 |
double Bnm1 = 1.0;
|
|
Packit |
67cb25 |
double a1 = b - a - 1.0;
|
|
Packit |
67cb25 |
double b1 = b - x - 2.0*(a+1.0);
|
|
Packit |
67cb25 |
double An = b1*Anm1 + a1*Anm2;
|
|
Packit |
67cb25 |
double Bn = b1*Bnm1 + a1*Bnm2;
|
|
Packit |
67cb25 |
double an, bn;
|
|
Packit |
67cb25 |
double fn = An/Bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(n < maxiter) {
|
|
Packit |
67cb25 |
double old_fn;
|
|
Packit |
67cb25 |
double del;
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
an = (a + n - 1.0) * (b - a - n);
|
|
Packit |
67cb25 |
bn = b - x - 2.0*(a+n);
|
|
Packit |
67cb25 |
An = bn*Anm1 + an*Anm2;
|
|
Packit |
67cb25 |
Bn = bn*Bnm1 + an*Bnm2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
|
|
Packit |
67cb25 |
An /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
old_fn = fn;
|
|
Packit |
67cb25 |
fn = An/Bn;
|
|
Packit |
67cb25 |
del = old_fn/fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = fn;
|
|
Packit |
67cb25 |
if(n == maxiter)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif /* 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(a,b+1,x)/1F1(a,b,x)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This seemed to suffer from "anomalous convergence".
|
|
Packit |
67cb25 |
* However, I have no theory for this recurrence.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double RECUR_BIG = GSL_SQRT_DBL_MAX;
|
|
Packit |
67cb25 |
const int maxiter = 5000;
|
|
Packit |
67cb25 |
int n = 1;
|
|
Packit |
67cb25 |
double Anm2 = 1.0;
|
|
Packit |
67cb25 |
double Bnm2 = 0.0;
|
|
Packit |
67cb25 |
double Anm1 = 0.0;
|
|
Packit |
67cb25 |
double Bnm1 = 1.0;
|
|
Packit |
67cb25 |
double a1 = b + 1.0;
|
|
Packit |
67cb25 |
double b1 = (b + 1.0) * (b - x);
|
|
Packit |
67cb25 |
double An = b1*Anm1 + a1*Anm2;
|
|
Packit |
67cb25 |
double Bn = b1*Bnm1 + a1*Bnm2;
|
|
Packit |
67cb25 |
double an, bn;
|
|
Packit |
67cb25 |
double fn = An/Bn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
while(n < maxiter) {
|
|
Packit |
67cb25 |
double old_fn;
|
|
Packit |
67cb25 |
double del;
|
|
Packit |
67cb25 |
n++;
|
|
Packit |
67cb25 |
Anm2 = Anm1;
|
|
Packit |
67cb25 |
Bnm2 = Bnm1;
|
|
Packit |
67cb25 |
Anm1 = An;
|
|
Packit |
67cb25 |
Bnm1 = Bn;
|
|
Packit |
67cb25 |
an = (b + n) * (b + n - 1.0 - a) * x;
|
|
Packit |
67cb25 |
bn = (b + n) * (b + n - 1.0 - x);
|
|
Packit |
67cb25 |
An = bn*Anm1 + an*Anm2;
|
|
Packit |
67cb25 |
Bn = bn*Bnm1 + an*Bnm2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
|
|
Packit |
67cb25 |
An /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bn /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm1 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Anm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
Bnm2 /= RECUR_BIG;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
old_fn = fn;
|
|
Packit |
67cb25 |
fn = An/Bn;
|
|
Packit |
67cb25 |
del = old_fn/fn;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*result = fn;
|
|
Packit |
67cb25 |
if(n == maxiter)
|
|
Packit |
67cb25 |
GSL_ERROR ("error", GSL_EMAXITER);
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif /* 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(a,b,x)
|
|
Packit |
67cb25 |
* |a| <= 1, b > 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double bma = b-a;
|
|
Packit |
67cb25 |
const double oma = 1.0-a;
|
|
Packit |
67cb25 |
const double ap1mb = 1.0+a-b;
|
|
Packit |
67cb25 |
const double abs_bma = fabs(bma);
|
|
Packit |
67cb25 |
const double abs_oma = fabs(oma);
|
|
Packit |
67cb25 |
const double abs_ap1mb = fabs(ap1mb);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double ax = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(a == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 1.0 && b >= 1.0) {
|
|
Packit |
67cb25 |
return hyperg_1F1_1(b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == -1.0) {
|
|
Packit |
67cb25 |
result->val = 1.0 + a/b * x;
|
|
Packit |
67cb25 |
result->err = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b >= 1.4*ax) {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.0) {
|
|
Packit |
67cb25 |
if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_posx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b < 5.0e+06) {
|
|
Packit |
67cb25 |
/* Recurse backward on b from
|
|
Packit |
67cb25 |
* a suitably high point.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double b_del = ceil(1.4*x-b) + 1.0;
|
|
Packit |
67cb25 |
double bp = b + b_del;
|
|
Packit |
67cb25 |
gsl_sf_result r_Mbp1;
|
|
Packit |
67cb25 |
gsl_sf_result r_Mb;
|
|
Packit |
67cb25 |
double Mbp1;
|
|
Packit |
67cb25 |
double Mb;
|
|
Packit |
67cb25 |
double Mbm1;
|
|
Packit |
67cb25 |
int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
|
|
Packit |
67cb25 |
int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp, x, &r_Mb);
|
|
Packit |
67cb25 |
const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
|
|
Packit |
67cb25 |
Mbp1 = r_Mbp1.val;
|
|
Packit |
67cb25 |
Mb = r_Mb.val;
|
|
Packit |
67cb25 |
while(bp > b+0.1) {
|
|
Packit |
67cb25 |
/* Do backward recursion. */
|
|
Packit |
67cb25 |
Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
|
|
Packit |
67cb25 |
bp -= 1.0;
|
|
Packit |
67cb25 |
Mbp1 = Mb;
|
|
Packit |
67cb25 |
Mb = Mbm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = Mb;
|
|
Packit |
67cb25 |
result->err = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_0, stat_1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (fabs(x) < fabs(b) && fabs(a*x) < sqrt(fabs(b)) * fabs(b-x)) {
|
|
Packit |
67cb25 |
return hyperg_1F1_largebx(a, b, x, result);
|
|
Packit |
67cb25 |
} else {
|
|
Packit |
67cb25 |
return hyperg_1F1_large2bm4a(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x < 0 and b not large compared to |x|
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(ax < 10.0 && b < 10.0) {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_negx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return hyperg_1F1_luke(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(b+eps,b,x)
|
|
Packit |
67cb25 |
* |eps|<=1, b > 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
|
|
Packit |
67cb25 |
/* If b-a is very small and x/b is not too large we can
|
|
Packit |
67cb25 |
* use this explicit approximation.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* v2 = a/(2b^2(b+1))
|
|
Packit |
67cb25 |
* v3 = a(b-2a)/(3b^3(b+1)(b+2))
|
|
Packit |
67cb25 |
* ...
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* See [Luke, Mathematical Functions and Their Approximations, p.292]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This cannot be used for b near a negative integer or zero.
|
|
Packit |
67cb25 |
* Also, if x/b is large the deviation from exp(x) behaviour grows.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double a = b + eps;
|
|
Packit |
67cb25 |
gsl_sf_result exab;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_e(a*x/b, &exab);
|
|
Packit |
67cb25 |
double v2 = a/(2.0*b*b*(b+1.0));
|
|
Packit |
67cb25 |
double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
|
|
Packit |
67cb25 |
double v = v2 + v3 * x;
|
|
Packit |
67cb25 |
double f = (1.0 - eps*x*x*v);
|
|
Packit |
67cb25 |
result->val = exab.val * f;
|
|
Packit |
67cb25 |
result->err = exab.err * fabs(f);
|
|
Packit |
67cb25 |
result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
|
|
Packit |
67cb25 |
result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Otherwise use a Kummer transformation to reduce
|
|
Packit |
67cb25 |
* it to the small a case.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Kummer_1F1;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
|
|
Packit |
67cb25 |
if(Kummer_1F1.val != 0.0) {
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
|
|
Packit |
67cb25 |
Kummer_1F1.val, Kummer_1F1.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
result->val = 0.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return stat_K;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* E(x) = exp(x) x > 0
|
|
Packit |
67cb25 |
* = 1 x < 0
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* a >= 1/2
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result I;
|
|
Packit |
67cb25 |
int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
|
|
Packit |
67cb25 |
gsl_sf_result lg;
|
|
Packit |
67cb25 |
int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
|
|
Packit |
67cb25 |
double ln_term = (0.5-a)*log(0.25*fabs(x));
|
|
Packit |
67cb25 |
double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
|
|
Packit |
67cb25 |
double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
|
|
Packit |
67cb25 |
I.val, I.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Determine middle parts of diagonal recursion along b=2a
|
|
Packit |
67cb25 |
* from two endpoints, i.e.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* given: M(a,b) and M(a+1,b+2)
|
|
Packit |
67cb25 |
* get: M(a+1,b+1) and M(a,b+1)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
inline
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_diag_step(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
const double Mab, const double Map1bp2,
|
|
Packit |
67cb25 |
double * Map1bp1, double * Mabp1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a == b) {
|
|
Packit |
67cb25 |
*Map1bp1 = Mab;
|
|
Packit |
67cb25 |
*Mabp1 = Mab - x/(b+1.0) * Map1bp2;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
*Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
|
|
Packit |
67cb25 |
*Mabp1 = (a * *Map1bp1 - b * Mab)/(a-b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif /* 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Determine endpoint of diagonal recursion.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* given: M(a,b) and M(a+1,b+2)
|
|
Packit |
67cb25 |
* get: M(a+1,b) and M(a+1,b+1)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
inline
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_diag_end_step(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
const double Mab, const double Map1bp2,
|
|
Packit |
67cb25 |
double * Map1b, double * Map1bp1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
*Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
|
|
Packit |
67cb25 |
*Map1b = Mab + x/b * *Map1bp1;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif /* 0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle the case of a and b both positive integers.
|
|
Packit |
67cb25 |
* Assumes a > 0 and b > 0.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ax = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(a == b) {
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result); /* 1F1(a,a,x) */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 1) {
|
|
Packit |
67cb25 |
return gsl_sf_exprel_n_e(b-1, x, result); /* 1F1(1,b,x) */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == a + 1) {
|
|
Packit |
67cb25 |
gsl_sf_result K;
|
|
Packit |
67cb25 |
int stat_K = gsl_sf_exprel_n_e(a, -x, &K); /* 1F1(1,1+a,-x) */
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
K.val, K.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == b + 1) {
|
|
Packit |
67cb25 |
gsl_sf_result ex;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_e(x, &ex);
|
|
Packit |
67cb25 |
result->val = ex.val * (1.0 + x/b);
|
|
Packit |
67cb25 |
result->err = ex.err * (1.0 + x/b);
|
|
Packit |
67cb25 |
result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == b + 2) {
|
|
Packit |
67cb25 |
gsl_sf_result ex;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_e(x, &ex);
|
|
Packit |
67cb25 |
double poly = (1.0 + x/b*(2.0 + x/(b+1.0)));
|
|
Packit |
67cb25 |
result->val = ex.val * poly;
|
|
Packit |
67cb25 |
result->err = ex.err * fabs(poly);
|
|
Packit |
67cb25 |
result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
|
|
Packit |
67cb25 |
return stat_e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 2*a) {
|
|
Packit |
67cb25 |
return hyperg_1F1_beq2a_pos(a, x, result); /* 1F1(a,2a,x) */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( ( b < 10 && a < 10 && ax < 5.0 )
|
|
Packit |
67cb25 |
|| ( b > a*ax )
|
|
Packit |
67cb25 |
|| ( b > a && ax < 5.0 )
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b > a && b >= 2*a + x) {
|
|
Packit |
67cb25 |
/* Use the Gautschi CF series, then
|
|
Packit |
67cb25 |
* recurse backward to a=0 for normalization.
|
|
Packit |
67cb25 |
* This will work for either sign of x.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double rap;
|
|
Packit |
67cb25 |
int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap;;
|
|
Packit |
67cb25 |
double ra = 1.0 + x/a * rap;
|
|
Packit |
67cb25 |
double Ma = GSL_SQRT_DBL_MIN;
|
|
Packit |
67cb25 |
double Map1 = ra * Ma;
|
|
Packit |
67cb25 |
double Mnp1 = Map1;
|
|
Packit |
67cb25 |
double Mn = Ma;
|
|
Packit |
67cb25 |
double Mnm1;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
for(n=a; n>0; n--) {
|
|
Packit |
67cb25 |
Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
|
|
Packit |
67cb25 |
Mnp1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = Ma/Mn;
|
|
Packit |
67cb25 |
result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
|
|
Packit |
67cb25 |
return stat_CF1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b > a && b < 2*a + x && b > x) {
|
|
Packit |
67cb25 |
/* Use the Gautschi series representation of
|
|
Packit |
67cb25 |
* the continued fraction. Then recurse forward
|
|
Packit |
67cb25 |
* to the a=b line for normalization. This will
|
|
Packit |
67cb25 |
* work for either sign of x, although we do need
|
|
Packit |
67cb25 |
* to check for b > x, for when x is positive.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double rap;
|
|
Packit |
67cb25 |
int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap;;
|
|
Packit |
67cb25 |
double ra = 1.0 + x/a * rap;
|
|
Packit |
67cb25 |
gsl_sf_result ex;
|
|
Packit |
67cb25 |
int stat_ex;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double Ma = GSL_SQRT_DBL_MIN;
|
|
Packit |
67cb25 |
double Map1 = ra * Ma;
|
|
Packit |
67cb25 |
double Mnm1 = Ma;
|
|
Packit |
67cb25 |
double Mn = Map1;
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
for(n=a+1; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
stat_ex = gsl_sf_exp_e(x, &ex); /* 1F1(b,b,x) */
|
|
Packit |
67cb25 |
result->val = ex.val * Ma/Mn;
|
|
Packit |
67cb25 |
result->err = ex.err * fabs(Ma/Mn);
|
|
Packit |
67cb25 |
result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x >= 0.0) {
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(b < a) {
|
|
Packit |
67cb25 |
/* The point b,b is below the b=2a+x line.
|
|
Packit |
67cb25 |
* Forward recursion on a from b,b+1 is possible.
|
|
Packit |
67cb25 |
* Note that a > b + 1 as well, since we already tried a = b + 1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
|
|
Packit |
67cb25 |
double ex = exp(x);
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
double Mnm1 = ex; /* 1F1(b,b,x) */
|
|
Packit |
67cb25 |
double Mn = ex * (1.0 + x/b); /* 1F1(b+1,b,x) */
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
for(n=b+1; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = Mn;
|
|
Packit |
67cb25 |
result->err = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
|
|
Packit |
67cb25 |
result->err *= fabs(a-b)+1.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* b > a
|
|
Packit |
67cb25 |
* b < 2a + x
|
|
Packit |
67cb25 |
* b <= x (otherwise we would have finished above)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Gautschi anomalous convergence region. However, we can
|
|
Packit |
67cb25 |
* recurse forward all the way from a=0,1 because we are
|
|
Packit |
67cb25 |
* always underneath the b=2a+x line.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result r_Mn;
|
|
Packit |
67cb25 |
double Mnm1 = 1.0; /* 1F1(0,b,x) */
|
|
Packit |
67cb25 |
double Mn; /* 1F1(1,b,x) */
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
gsl_sf_exprel_n_e(b-1, x, &r_Mn);
|
|
Packit |
67cb25 |
Mn = r_Mn.val;
|
|
Packit |
67cb25 |
for(n=1; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = Mn;
|
|
Packit |
67cb25 |
result->err = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x < 0
|
|
Packit |
67cb25 |
* b < a (otherwise we would have tripped one of the above)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(a <= 0.5*(b-x) || a >= -x) {
|
|
Packit |
67cb25 |
/* Gautschi continued fraction is in the anomalous region,
|
|
Packit |
67cb25 |
* so we must find another way. We recurse down in b,
|
|
Packit |
67cb25 |
* from the a=b line.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double ex = exp(x);
|
|
Packit |
67cb25 |
double Manp1 = ex;
|
|
Packit |
67cb25 |
double Man = ex * (1.0 + x/(a-1.0));
|
|
Packit |
67cb25 |
double Manm1;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
for(n=a-1; n>b; n--) {
|
|
Packit |
67cb25 |
Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
|
|
Packit |
67cb25 |
Manp1 = Man;
|
|
Packit |
67cb25 |
Man = Manm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = Man;
|
|
Packit |
67cb25 |
result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
|
|
Packit |
67cb25 |
result->err *= fabs(b-a)+1.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Pick a0 such that b ~= 2a0 + x, then
|
|
Packit |
67cb25 |
* recurse down in b from a0,a0 to determine
|
|
Packit |
67cb25 |
* the values near the line b=2a+x. Then recurse
|
|
Packit |
67cb25 |
* forward on a from a0.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
int a0 = (int) ceil(0.5*(b-x));
|
|
Packit |
67cb25 |
double Ma0b; /* M(a0,b) */
|
|
Packit |
67cb25 |
double Ma0bp1; /* M(a0,b+1) */
|
|
Packit |
67cb25 |
double Ma0p1b; /* M(a0+1,b) */
|
|
Packit |
67cb25 |
double Mnm1;
|
|
Packit |
67cb25 |
double Mn;
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
int n;
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double ex = exp(x);
|
|
Packit |
67cb25 |
double Ma0np1 = ex;
|
|
Packit |
67cb25 |
double Ma0n = ex * (1.0 + x/(a0-1.0));
|
|
Packit |
67cb25 |
double Ma0nm1;
|
|
Packit |
67cb25 |
for(n=a0-1; n>b; n--) {
|
|
Packit |
67cb25 |
Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
|
|
Packit |
67cb25 |
Ma0np1 = Ma0n;
|
|
Packit |
67cb25 |
Ma0n = Ma0nm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
Ma0bp1 = Ma0np1;
|
|
Packit |
67cb25 |
Ma0b = Ma0n;
|
|
Packit |
67cb25 |
Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialise the recurrence correctly BJG */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (a0 >= a)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mn = Ma0b;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (a0 + 1>= a)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mn = Ma0p1b;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mnm1 = Ma0b;
|
|
Packit |
67cb25 |
Mn = Ma0p1b;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(n=a0+1; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = Mn;
|
|
Packit |
67cb25 |
result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
|
|
Packit |
67cb25 |
result->err *= fabs(b-a)+1.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
|
|
Packit |
67cb25 |
* When the terms are all positive, this
|
|
Packit |
67cb25 |
* must work. We will assume this here.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(a == 0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 1.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
int N = -a;
|
|
Packit |
67cb25 |
double poly = 1.0;
|
|
Packit |
67cb25 |
int k;
|
|
Packit |
67cb25 |
for(k=N-1; k>=0; k--) {
|
|
Packit |
67cb25 |
double t = (a+k)/(b+k) * (x/(k+1));
|
|
Packit |
67cb25 |
double r = t + 1.0/poly;
|
|
Packit |
67cb25 |
if(r > 0.9*GSL_DBL_MAX/poly) {
|
|
Packit |
67cb25 |
OVERFLOW_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
poly *= r; /* P_n = 1 + t_n P_{n-1} */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
result->val = poly;
|
|
Packit |
67cb25 |
result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Evaluate negative integer a case by relation
|
|
Packit |
67cb25 |
* to Laguerre polynomials. This is more general than
|
|
Packit |
67cb25 |
* the direct polynomial evaluation, but is safe
|
|
Packit |
67cb25 |
* for all values of x.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
|
|
Packit |
67cb25 |
* = n B(b,n) Laguerre[n,b-1,x]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* assumes b is not a negative integer
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int n = -a;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result lag;
|
|
Packit |
67cb25 |
const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
|
|
Packit |
67cb25 |
if(b < 0.0) {
|
|
Packit |
67cb25 |
gsl_sf_result lnfact;
|
|
Packit |
67cb25 |
gsl_sf_result lng1;
|
|
Packit |
67cb25 |
gsl_sf_result lng2;
|
|
Packit |
67cb25 |
double s1, s2;
|
|
Packit |
67cb25 |
const int stat_f = gsl_sf_lnfact_e(n, &lnfact);
|
|
Packit |
67cb25 |
const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1;;
|
|
Packit |
67cb25 |
const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2;;
|
|
Packit |
67cb25 |
const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
|
|
Packit |
67cb25 |
const double lnpre_err = lnfact.err + lng1.err + lng2.err
|
|
Packit |
67cb25 |
+ 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
|
|
Packit |
67cb25 |
const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
|
|
Packit |
67cb25 |
s1*s2*lag.val, lag.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
gsl_sf_result lnbeta;
|
|
Packit |
67cb25 |
gsl_sf_lnbeta_e(b, n, &lnbeta);
|
|
Packit |
67cb25 |
if(fabs(lnbeta.val) < 0.1) {
|
|
Packit |
67cb25 |
/* As we have noted, when B(x,y) is near 1,
|
|
Packit |
67cb25 |
* evaluating log(B(x,y)) is not accurate.
|
|
Packit |
67cb25 |
* Instead we evaluate B(x,y) directly.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double ln_term_val = log(1.25*n);
|
|
Packit |
67cb25 |
const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
|
|
Packit |
67cb25 |
gsl_sf_result beta;
|
|
Packit |
67cb25 |
int stat_b = gsl_sf_beta_e(b, n, &beta);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
|
|
Packit |
67cb25 |
lag.val, lag.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
result->val *= beta.val/1.25;
|
|
Packit |
67cb25 |
result->err *= beta.val/1.25;
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* B(x,y) was not near 1, so it is safe to use
|
|
Packit |
67cb25 |
* the logarithmic values.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double ln_n = log(n);
|
|
Packit |
67cb25 |
const double ln_term_val = lnbeta.val + ln_n;
|
|
Packit |
67cb25 |
const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
|
|
Packit |
67cb25 |
lag.val, lag.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_l);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle negative integer a case for x > 0 and
|
|
Packit |
67cb25 |
* generic b.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
|
|
Packit |
67cb25 |
* M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const int n = -a;
|
|
Packit |
67cb25 |
const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
|
|
Packit |
67cb25 |
double sgpoch;
|
|
Packit |
67cb25 |
gsl_sf_result lnpoch;
|
|
Packit |
67cb25 |
gsl_sf_result U;
|
|
Packit |
67cb25 |
const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
|
|
Packit |
67cb25 |
const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
|
|
Packit |
67cb25 |
const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
|
|
Packit |
67cb25 |
sgn * sgpoch * U.val, U.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Assumes a <= -1, b <= -1, and b <= a.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.0) {
|
|
Packit |
67cb25 |
return hyperg_1F1_a_negint_poly(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Apply a Kummer transformation to make x > 0 so
|
|
Packit |
67cb25 |
* we can evaluate the polynomial safely. Of course,
|
|
Packit |
67cb25 |
* this assumes b <= a, which must be true for
|
|
Packit |
67cb25 |
* a<0 and b<0, since otherwise the thing is undefined.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result K;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
K.val, K.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* [Abramowitz+Stegun, 13.1.3]
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
|
|
Packit |
67cb25 |
* { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* b not an integer >= 2
|
|
Packit |
67cb25 |
* a-b not a negative integer
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double bp = 2.0 - b;
|
|
Packit |
67cb25 |
const double ap = a - b + 1.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result lg_ap, lg_bp;
|
|
Packit |
67cb25 |
double sg_ap;
|
|
Packit |
67cb25 |
int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
|
|
Packit |
67cb25 |
int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
|
|
Packit |
67cb25 |
int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
|
|
Packit |
67cb25 |
double t1 = (bp-1.0) * log(x);
|
|
Packit |
67cb25 |
double lnpre_val = lg_ap.val - lg_bp.val + t1;
|
|
Packit |
67cb25 |
double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result lg_2mbp, lg_1papmbp;
|
|
Packit |
67cb25 |
double sg_2mbp, sg_1papmbp;
|
|
Packit |
67cb25 |
int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp, &lg_2mbp, &sg_2mbp);
|
|
Packit |
67cb25 |
int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
|
|
Packit |
67cb25 |
int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
|
|
Packit |
67cb25 |
double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
|
|
Packit |
67cb25 |
double lnc1_err = lg_2mbp.err + lg_1papmbp.err
|
|
Packit |
67cb25 |
+ GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result M;
|
|
Packit |
67cb25 |
gsl_sf_result_e10 U;
|
|
Packit |
67cb25 |
int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
|
|
Packit |
67cb25 |
int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
|
|
Packit |
67cb25 |
int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result_e10 term_M;
|
|
Packit |
67cb25 |
int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
|
|
Packit |
67cb25 |
sg_2mbp*sg_1papmbp*M.val, M.err,
|
|
Packit |
67cb25 |
&term_M);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const double ombp = 1.0 - bp;
|
|
Packit |
67cb25 |
const double Uee_val = U.e10*M_LN10;
|
|
Packit |
67cb25 |
const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
|
|
Packit |
67cb25 |
const double Mee_val = term_M.e10*M_LN10;
|
|
Packit |
67cb25 |
const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
|
|
Packit |
67cb25 |
int stat_e1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Do a little dance with the exponential prefactors
|
|
Packit |
67cb25 |
* to avoid overflows in intermediate results.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(Uee_val > Mee_val) {
|
|
Packit |
67cb25 |
const double factorM_val = exp(Mee_val-Uee_val);
|
|
Packit |
67cb25 |
const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
|
|
Packit |
67cb25 |
const double inner_val = term_M.val*factorM_val - ombp*U.val;
|
|
Packit |
67cb25 |
const double inner_err =
|
|
Packit |
67cb25 |
term_M.err*factorM_val + fabs(ombp) * U.err
|
|
Packit |
67cb25 |
+ fabs(term_M.val) * factorM_err
|
|
Packit |
67cb25 |
+ GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
|
|
Packit |
67cb25 |
stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
|
|
Packit |
67cb25 |
sg_ap*inner_val, inner_err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
const double factorU_val = exp(Uee_val - Mee_val);
|
|
Packit |
67cb25 |
const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
|
|
Packit |
67cb25 |
const double inner_val = term_M.val - ombp*factorU_val*U.val;
|
|
Packit |
67cb25 |
const double inner_err =
|
|
Packit |
67cb25 |
term_M.err + fabs(ombp*factorU_val*U.err)
|
|
Packit |
67cb25 |
+ fabs(ombp*factorU_err*U.val)
|
|
Packit |
67cb25 |
+ GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
|
|
Packit |
67cb25 |
stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
|
|
Packit |
67cb25 |
sg_ap*inner_val, inner_err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Handle case of generic positive a, b.
|
|
Packit |
67cb25 |
* Assumes b-a is not a negative integer.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_ab_pos(const double a, const double b,
|
|
Packit |
67cb25 |
const double x,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double ax = fabs(x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if( ( b < 10.0 && a < 10.0 && ax < 5.0 )
|
|
Packit |
67cb25 |
|| ( b > a*ax )
|
|
Packit |
67cb25 |
|| ( b > a && ax < 5.0 )
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( x < -100.0
|
|
Packit |
67cb25 |
&& GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Large negative x asymptotic.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_negx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( x > 100.0
|
|
Packit |
67cb25 |
&& GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Large positive x asymptotic.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_posx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(fabs(b-a) <= 1.0) {
|
|
Packit |
67cb25 |
/* Directly handle b near a.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_beps_bgt0(a-b, b, x, result); /* a = b + eps */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
else if(b > a && b >= 2*a + x) {
|
|
Packit |
67cb25 |
/* Use the Gautschi CF series, then
|
|
Packit |
67cb25 |
* recurse backward to a near 0 for normalization.
|
|
Packit |
67cb25 |
* This will work for either sign of x.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double rap;
|
|
Packit |
67cb25 |
int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap;;
|
|
Packit |
67cb25 |
double ra = 1.0 + x/a * rap;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double Ma = GSL_SQRT_DBL_MIN;
|
|
Packit |
67cb25 |
double Map1 = ra * Ma;
|
|
Packit |
67cb25 |
double Mnp1 = Map1;
|
|
Packit |
67cb25 |
double Mn = Ma;
|
|
Packit |
67cb25 |
double Mnm1;
|
|
Packit |
67cb25 |
gsl_sf_result Mn_true;
|
|
Packit |
67cb25 |
int stat_Mt;
|
|
Packit |
67cb25 |
double n;
|
|
Packit |
67cb25 |
for(n=a; n>0.5; n -= 1.0) {
|
|
Packit |
67cb25 |
Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
|
|
Packit |
67cb25 |
Mnp1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = (Ma/Mn) * Mn_true.val;
|
|
Packit |
67cb25 |
result->err = fabs(Ma/Mn) * Mn_true.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b > a && b < 2*a + x && b > x) {
|
|
Packit |
67cb25 |
/* Use the Gautschi series representation of
|
|
Packit |
67cb25 |
* the continued fraction. Then recurse forward
|
|
Packit |
67cb25 |
* to near the a=b line for normalization. This will
|
|
Packit |
67cb25 |
* work for either sign of x, although we do need
|
|
Packit |
67cb25 |
* to check for b > x, which is relevant when x is positive.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Mn_true;
|
|
Packit |
67cb25 |
int stat_Mt;
|
|
Packit |
67cb25 |
double rap;
|
|
Packit |
67cb25 |
int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap;;
|
|
Packit |
67cb25 |
double ra = 1.0 + x/a * rap;
|
|
Packit |
67cb25 |
double Ma = GSL_SQRT_DBL_MIN;
|
|
Packit |
67cb25 |
double Mnm1 = Ma;
|
|
Packit |
67cb25 |
double Mn = ra * Mnm1;
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
double n;
|
|
Packit |
67cb25 |
for(n=a+1.0; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
|
|
Packit |
67cb25 |
result->val = Ma/Mn * Mn_true.val;
|
|
Packit |
67cb25 |
result->err = fabs(Ma/Mn) * Mn_true.err;
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x >= 0.0) {
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(b < a) {
|
|
Packit |
67cb25 |
/* Forward recursion on a from a=b+eps-1,b+eps.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double N = floor(a-b);
|
|
Packit |
67cb25 |
double eps = a - b - N;
|
|
Packit |
67cb25 |
gsl_sf_result r_M0;
|
|
Packit |
67cb25 |
gsl_sf_result r_M1;
|
|
Packit |
67cb25 |
int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
|
|
Packit |
67cb25 |
int stat_1 = hyperg_1F1_beps_bgt0(eps, b, x, &r_M1);
|
|
Packit |
67cb25 |
double M0 = r_M0.val;
|
|
Packit |
67cb25 |
double M1 = r_M1.val;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double Mam1 = M0;
|
|
Packit |
67cb25 |
double Ma = M1;
|
|
Packit |
67cb25 |
double Map1;
|
|
Packit |
67cb25 |
double ap;
|
|
Packit |
67cb25 |
double start_pair = fabs(M0) + fabs(M1);
|
|
Packit |
67cb25 |
double minim_pair = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
double pair_ratio;
|
|
Packit |
67cb25 |
double rat_0 = fabs(r_M0.err/r_M0.val);
|
|
Packit |
67cb25 |
double rat_1 = fabs(r_M1.err/r_M1.val);
|
|
Packit |
67cb25 |
for(ap=b+eps; ap
|
|
Packit |
67cb25 |
Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
|
|
Packit |
67cb25 |
Mam1 = Ma;
|
|
Packit |
67cb25 |
Ma = Map1;
|
|
Packit |
67cb25 |
minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
pair_ratio = start_pair/minim_pair;
|
|
Packit |
67cb25 |
result->val = Ma;
|
|
Packit |
67cb25 |
result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
|
|
Packit |
67cb25 |
result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_0, stat_1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* b > a
|
|
Packit |
67cb25 |
* b < 2a + x
|
|
Packit |
67cb25 |
* b <= x
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Recurse forward on a from a=eps,eps+1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double eps = a - floor(a);
|
|
Packit |
67cb25 |
gsl_sf_result r_Mnm1;
|
|
Packit |
67cb25 |
gsl_sf_result r_Mn;
|
|
Packit |
67cb25 |
int stat_0 = hyperg_1F1_small_a_bgt0(eps, b, x, &r_Mnm1);
|
|
Packit |
67cb25 |
int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
|
|
Packit |
67cb25 |
double Mnm1 = r_Mnm1.val;
|
|
Packit |
67cb25 |
double Mn = r_Mn.val;
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double n;
|
|
Packit |
67cb25 |
double start_pair = fabs(Mn) + fabs(Mnm1);
|
|
Packit |
67cb25 |
double minim_pair = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
double pair_ratio;
|
|
Packit |
67cb25 |
double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
|
|
Packit |
67cb25 |
double rat_1 = fabs(r_Mn.err/r_Mn.val);
|
|
Packit |
67cb25 |
for(n=eps+1.0; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
pair_ratio = start_pair/minim_pair;
|
|
Packit |
67cb25 |
result->val = Mn;
|
|
Packit |
67cb25 |
result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
|
|
Packit |
67cb25 |
result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
|
|
Packit |
67cb25 |
result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_0, stat_1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* x < 0
|
|
Packit |
67cb25 |
* b < a
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(a <= 0.5*(b-x) || a >= -x) {
|
|
Packit |
67cb25 |
/* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double N = floor(a - b);
|
|
Packit |
67cb25 |
double eps = 1.0 + N - a + b;
|
|
Packit |
67cb25 |
gsl_sf_result r_Manp1;
|
|
Packit |
67cb25 |
gsl_sf_result r_Man;
|
|
Packit |
67cb25 |
int stat_0 = hyperg_1F1_beps_bgt0(-eps, a+eps, x, &r_Manp1);
|
|
Packit |
67cb25 |
int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
|
|
Packit |
67cb25 |
double Manp1 = r_Manp1.val;
|
|
Packit |
67cb25 |
double Man = r_Man.val;
|
|
Packit |
67cb25 |
double Manm1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double n;
|
|
Packit |
67cb25 |
double start_pair = fabs(Manp1) + fabs(Man);
|
|
Packit |
67cb25 |
double minim_pair = GSL_DBL_MAX;
|
|
Packit |
67cb25 |
double pair_ratio;
|
|
Packit |
67cb25 |
double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
|
|
Packit |
67cb25 |
double rat_1 = fabs(r_Man.err/r_Man.val);
|
|
Packit |
67cb25 |
for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
|
|
Packit |
67cb25 |
Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
|
|
Packit |
67cb25 |
Manp1 = Man;
|
|
Packit |
67cb25 |
Man = Manm1;
|
|
Packit |
67cb25 |
minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* FIXME: this is a nasty little hack; there is some
|
|
Packit |
67cb25 |
(transient?) instability in this recurrence for some
|
|
Packit |
67cb25 |
values. I can tell when it happens, which is when
|
|
Packit |
67cb25 |
this pair_ratio is large. But I do not know how to
|
|
Packit |
67cb25 |
measure the error in terms of it. I guessed quadratic
|
|
Packit |
67cb25 |
below, but it is probably worse than that.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
pair_ratio = start_pair/minim_pair;
|
|
Packit |
67cb25 |
result->val = Man;
|
|
Packit |
67cb25 |
result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
|
|
Packit |
67cb25 |
result->err *= pair_ratio*pair_ratio + 1.0;
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_0, stat_1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Pick a0 such that b ~= 2a0 + x, then
|
|
Packit |
67cb25 |
* recurse down in b from a0,a0 to determine
|
|
Packit |
67cb25 |
* the values near the line b=2a+x. Then recurse
|
|
Packit |
67cb25 |
* forward on a from a0.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
double epsa = a - floor(a);
|
|
Packit |
67cb25 |
double a0 = floor(0.5*(b-x)) + epsa;
|
|
Packit |
67cb25 |
double N = floor(a0 - b);
|
|
Packit |
67cb25 |
double epsb = 1.0 + N - a0 + b;
|
|
Packit |
67cb25 |
double Ma0b;
|
|
Packit |
67cb25 |
double Ma0bp1;
|
|
Packit |
67cb25 |
double Ma0p1b;
|
|
Packit |
67cb25 |
int stat_a0;
|
|
Packit |
67cb25 |
double Mnm1;
|
|
Packit |
67cb25 |
double Mn;
|
|
Packit |
67cb25 |
double Mnp1;
|
|
Packit |
67cb25 |
double n;
|
|
Packit |
67cb25 |
double err_rat;
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_sf_result r_Ma0np1;
|
|
Packit |
67cb25 |
gsl_sf_result r_Ma0n;
|
|
Packit |
67cb25 |
int stat_0 = hyperg_1F1_beps_bgt0(-epsb, a0+epsb, x, &r_Ma0np1);
|
|
Packit |
67cb25 |
int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
|
|
Packit |
67cb25 |
double Ma0np1 = r_Ma0np1.val;
|
|
Packit |
67cb25 |
double Ma0n = r_Ma0n.val;
|
|
Packit |
67cb25 |
double Ma0nm1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
|
|
Packit |
67cb25 |
Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
|
|
Packit |
67cb25 |
Ma0np1 = Ma0n;
|
|
Packit |
67cb25 |
Ma0n = Ma0nm1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
Ma0bp1 = Ma0np1;
|
|
Packit |
67cb25 |
Ma0b = Ma0n;
|
|
Packit |
67cb25 |
Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */
|
|
Packit |
67cb25 |
stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialise the recurrence correctly BJG */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (a0 >= a - 0.1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mn = Ma0b;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (a0 + 1>= a - 0.1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mn = Ma0p1b;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
Mnm1 = Ma0b;
|
|
Packit |
67cb25 |
Mn = Ma0p1b;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(n=a0+1.0; n
|
|
Packit |
67cb25 |
Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
|
|
Packit |
67cb25 |
Mnm1 = Mn;
|
|
Packit |
67cb25 |
Mn = Mnp1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
result->val = Mn;
|
|
Packit |
67cb25 |
result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
|
|
Packit |
67cb25 |
return stat_a0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Assumes b != integer
|
|
Packit |
67cb25 |
* Assumes a != integer when x > 0
|
|
Packit |
67cb25 |
* Assumes b-a != neg integer when x < 0
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
static
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
hyperg_1F1_ab_neg(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double bma = b - a;
|
|
Packit |
67cb25 |
const double abs_x = fabs(x);
|
|
Packit |
67cb25 |
const double abs_a = fabs(a);
|
|
Packit |
67cb25 |
const double abs_b = fabs(b);
|
|
Packit |
67cb25 |
const double size_a = GSL_MAX(abs_a, 1.0);
|
|
Packit |
67cb25 |
const double size_b = GSL_MAX(abs_b, 1.0);
|
|
Packit |
67cb25 |
const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if( (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
|
|
Packit |
67cb25 |
|| (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( x > 0.0
|
|
Packit |
67cb25 |
&& size_b > size_a
|
|
Packit |
67cb25 |
&& size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Series terms are positive definite up until
|
|
Packit |
67cb25 |
* there is a sign change. But by then the
|
|
Packit |
67cb25 |
* terms are small due to the last condition.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
|
|
Packit |
67cb25 |
|| (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Use Kummer transformation to render series safe.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Kummer_1F1;
|
|
Packit |
67cb25 |
int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
Kummer_1F1.val, Kummer_1F1.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( x < -30.0
|
|
Packit |
67cb25 |
&& GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Large negative x asymptotic.
|
|
Packit |
67cb25 |
* Note that we do not check if b-a is a negative integer.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_negx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if( x > 100.0
|
|
Packit |
67cb25 |
&& GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
|
|
Packit |
67cb25 |
) {
|
|
Packit |
67cb25 |
/* Large positive x asymptotic.
|
|
Packit |
67cb25 |
* Note that we do not check if a is a negative integer.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_posx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
|
|
Packit |
67cb25 |
return hyperg_1F1_U(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* FIXME: if all else fails, try the series... BJG */
|
|
Packit |
67cb25 |
if (x < 0.0) {
|
|
Packit |
67cb25 |
/* Apply Kummer Transformation */
|
|
Packit |
67cb25 |
int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
|
|
Packit |
67cb25 |
double K_factor = exp(x);
|
|
Packit |
67cb25 |
result->val *= K_factor;
|
|
Packit |
67cb25 |
result->err *= K_factor;
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
} else {
|
|
Packit |
67cb25 |
int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
return status;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Sadness... */
|
|
Packit |
67cb25 |
/* result->val = 0.0; */
|
|
Packit |
67cb25 |
/* result->err = 0.0; */
|
|
Packit |
67cb25 |
/* GSL_ERROR ("error", GSL_EUNIMPL); */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == b) {
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b < 0 && (a < b || a > 0)) {
|
|
Packit |
67cb25 |
/* Standard domain error due to singularity. */
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x > 100.0 && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
|
|
Packit |
67cb25 |
/* x -> +Inf asymptotic */
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_posx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) {
|
|
Packit |
67cb25 |
/* x -> -Inf asymptotic */
|
|
Packit |
67cb25 |
return hyperg_1F1_asymp_negx(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a < 0 && b < 0) {
|
|
Packit |
67cb25 |
return hyperg_1F1_ab_negint(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a < 0 && b > 0) {
|
|
Packit |
67cb25 |
/* Use Kummer to reduce it to the positive integer case.
|
|
Packit |
67cb25 |
* Note that b > a, strictly, since we already trapped b = a.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Kummer_1F1;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
Kummer_1F1.val, Kummer_1F1.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* a > 0 and b > 0 */
|
|
Packit |
67cb25 |
return hyperg_1F1_ab_posint(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
|
|
Packit |
67cb25 |
gsl_sf_result * result
|
|
Packit |
67cb25 |
)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
const double bma = b - a;
|
|
Packit |
67cb25 |
const double rinta = floor(a + 0.5);
|
|
Packit |
67cb25 |
const double rintb = floor(b + 0.5);
|
|
Packit |
67cb25 |
const double rintbma = floor(bma + 0.5);
|
|
Packit |
67cb25 |
const int a_integer = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
|
|
Packit |
67cb25 |
const int b_integer = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
|
|
Packit |
67cb25 |
const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
|
|
Packit |
67cb25 |
const int b_neg_integer = ( b < -0.1 && b_integer );
|
|
Packit |
67cb25 |
const int a_neg_integer = ( a < -0.1 && a_integer );
|
|
Packit |
67cb25 |
const int bma_neg_integer = ( bma < -0.1 && bma_integer );
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* CHECK_POINTER(result) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(x == 0.0) {
|
|
Packit |
67cb25 |
/* Testing for this before testing a and b
|
|
Packit |
67cb25 |
* is somewhat arbitrary. The result is that
|
|
Packit |
67cb25 |
* we have 1F1(a,0,0) = 1.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b == 0.0) {
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == 0.0) {
|
|
Packit |
67cb25 |
result->val = 1.0;
|
|
Packit |
67cb25 |
result->err = 0.0;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a == b) {
|
|
Packit |
67cb25 |
/* case: a==b; exp(x)
|
|
Packit |
67cb25 |
* It's good to test exact equality now.
|
|
Packit |
67cb25 |
* We also test approximate equality later.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gsl_sf_exp_e(x, result);
|
|
Packit |
67cb25 |
} else if(fabs(b) < _1F1_INT_THRESHOLD && fabs(a) < _1F1_INT_THRESHOLD) {
|
|
Packit |
67cb25 |
/* a and b near zero: 1 + a/b (exp(x)-1)
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Note that neither a nor b is zero, since
|
|
Packit |
67cb25 |
* we eliminated that with the above tests.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_sf_result exm1;
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_expm1_e(x, &exm1);
|
|
Packit |
67cb25 |
double sa = ( a > 0.0 ? 1.0 : -1.0 );
|
|
Packit |
67cb25 |
double sb = ( b > 0.0 ? 1.0 : -1.0 );
|
|
Packit |
67cb25 |
double lnab = log(fabs(a/b)); /* safe */
|
|
Packit |
67cb25 |
gsl_sf_result hx;
|
|
Packit |
67cb25 |
int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
|
|
Packit |
67cb25 |
sa * sb * exm1.val, exm1.err,
|
|
Packit |
67cb25 |
&hx;;
|
|
Packit |
67cb25 |
result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val); /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
|
|
Packit |
67cb25 |
result->err = hx.err;
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_hx, stat_e);
|
|
Packit |
67cb25 |
} else if (fabs(b) < _1F1_INT_THRESHOLD && fabs(x*a) < 1) {
|
|
Packit |
67cb25 |
/* b near zero and a not near zero
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
const double m_arg = 1.0/(0.5*b);
|
|
Packit |
67cb25 |
gsl_sf_result F_renorm;
|
|
Packit |
67cb25 |
int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
|
|
Packit |
67cb25 |
int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
|
|
Packit |
67cb25 |
0.5*F_renorm.val, 0.5*F_renorm.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_m, stat_F);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a_integer && b_integer) {
|
|
Packit |
67cb25 |
/* Check for reduction to the integer case.
|
|
Packit |
67cb25 |
* Relies on the arbitrary "near an integer" test.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b_neg_integer && !(a_neg_integer && a > b)) {
|
|
Packit |
67cb25 |
/* Standard domain error due to
|
|
Packit |
67cb25 |
* uncancelled singularity.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
DOMAIN_ERROR(result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a_neg_integer) {
|
|
Packit |
67cb25 |
return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(b > 0.0) {
|
|
Packit |
67cb25 |
if(-1.0 <= a && a <= 1.0) {
|
|
Packit |
67cb25 |
/* Handle small a explicitly.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
return hyperg_1F1_small_a_bgt0(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(bma_neg_integer) {
|
|
Packit |
67cb25 |
/* Catch this now, to avoid problems in the
|
|
Packit |
67cb25 |
* generic evaluation code.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Kummer_1F1;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
Kummer_1F1.val, Kummer_1F1.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a < 0.0 && fabs(x) < 2*GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
/* Use Kummer to reduce it to the generic positive case.
|
|
Packit |
67cb25 |
* Note that b > a, strictly, since we already trapped b = a.
|
|
Packit |
67cb25 |
* Also b-(b-a)=a, and a is not a negative integer here,
|
|
Packit |
67cb25 |
* so the generic evaluation is safe.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result Kummer_1F1;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
Kummer_1F1.val, Kummer_1F1.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if (a > 0) {
|
|
Packit |
67cb25 |
/* a > 0.0 */
|
|
Packit |
67cb25 |
return hyperg_1F1_ab_pos(a, b, x, result);
|
|
Packit |
67cb25 |
} else {
|
|
Packit |
67cb25 |
return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* b < 0.0 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if(bma_neg_integer && x < 0.0) {
|
|
Packit |
67cb25 |
/* Handle this now to prevent problems
|
|
Packit |
67cb25 |
* in the generic evaluation.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result K;
|
|
Packit |
67cb25 |
int stat_K;
|
|
Packit |
67cb25 |
int stat_e;
|
|
Packit |
67cb25 |
if(a < 0.0) {
|
|
Packit |
67cb25 |
/* Kummer transformed version of safe polynomial.
|
|
Packit |
67cb25 |
* The condition a < 0 is equivalent to b < b-a,
|
|
Packit |
67cb25 |
* which is the condition required for the series
|
|
Packit |
67cb25 |
* to be positive definite here.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
/* Generic eval for negative integer a. */
|
|
Packit |
67cb25 |
stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
K.val, K.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(a > 0.0) {
|
|
Packit |
67cb25 |
/* Use Kummer to reduce it to the generic negative case.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
gsl_sf_result K;
|
|
Packit |
67cb25 |
int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
|
|
Packit |
67cb25 |
K.val, K.err,
|
|
Packit |
67cb25 |
result);
|
|
Packit |
67cb25 |
return GSL_ERROR_SELECT_2(stat_e, stat_K);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
return hyperg_1F1_ab_neg(a, b, x, result);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#if 0
|
|
Packit |
67cb25 |
/* Luke in the canonical case.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
|
|
Packit |
67cb25 |
double prec;
|
|
Packit |
67cb25 |
return hyperg_1F1_luke(a, b, x, result, &prec);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Luke with Kummer transformation.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
|
|
Packit |
67cb25 |
double prec;
|
|
Packit |
67cb25 |
double Kummer_1F1;
|
|
Packit |
67cb25 |
double ex;
|
|
Packit |
67cb25 |
int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
|
|
Packit |
67cb25 |
int stat_e = gsl_sf_exp_e(x, &ex);
|
|
Packit |
67cb25 |
if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
|
|
Packit |
67cb25 |
double lnr = log(fabs(Kummer_1F1)) + x;
|
|
Packit |
67cb25 |
if(lnr < GSL_LOG_DBL_MAX) {
|
|
Packit |
67cb25 |
*result = ex * Kummer_1F1;
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
*result = GSL_POSINF;
|
|
Packit |
67cb25 |
GSL_ERROR ("overflow", GSL_EOVRFLW);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else if(stat_F != GSL_SUCCESS) {
|
|
Packit |
67cb25 |
*result = 0.0;
|
|
Packit |
67cb25 |
return stat_F;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else {
|
|
Packit |
67cb25 |
*result = 0.0;
|
|
Packit |
67cb25 |
return stat_e;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "eval.h"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double gsl_sf_hyperg_1F1(double a, double b, double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
|
|
Packit |
67cb25 |
}
|