hyperg_1f1.c

math library from gnu

/* specfunc/hyperg_1F1.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_elementary.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_laguerre.h>
#include <gsl/gsl_sf_hyperg.h>

#include "error.h"
#include "hyperg.h"

#define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON)


/* Asymptotic result for 1F1(a, b, x)  x -> -Infinity.
 * Assumes b-a != neg integer and b != neg integer.
 */
static
int
hyperg_1F1_asymp_negx(const double a, const double b, const double x,
                     gsl_sf_result * result)
{
  gsl_sf_result lg_b;
  gsl_sf_result lg_bma;
  double sgn_b;
  double sgn_bma;

  int stat_b   = gsl_sf_lngamma_sgn_e(b,   &lg_b,   &sgn_b);
  int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma);

  if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) {
    gsl_sf_result F;
    int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F);
    if(F.val != 0) {
      double ln_term_val = a*log(-x);
      double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val));
      double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val;
      double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err;
      int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
                                            sgn_bma*sgn_b*F.val, F.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_F);
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
      return stat_F;
    }
  }
  else {
    DOMAIN_ERROR(result);
  }
}


/* Asymptotic result for 1F1(a, b, x)  x -> +Infinity
 * Assumes b != neg integer and a != neg integer
 */
static
int
hyperg_1F1_asymp_posx(const double a, const double b, const double x,
                      gsl_sf_result * result)
{
  gsl_sf_result lg_b;
  gsl_sf_result lg_a;
  double sgn_b;
  double sgn_a;

  int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
  int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);

  if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) {
    gsl_sf_result F;
    int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F);
    if(stat_F == GSL_SUCCESS && F.val != 0) {
      double lnx = log(x);
      double ln_term_val = (a-b)*lnx;
      double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx)
                         + 2.0 * GSL_DBL_EPSILON * fabs(a-b);
      double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x;
      double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x);
      int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
                                            sgn_a*sgn_b*F.val, F.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_F);
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
      return stat_F;
    }
  }
  else {
    DOMAIN_ERROR(result);
  }
}

/* Asymptotic result from Slater 4.3.7 
 * 
 * To get the general series, write M(a,b,x) as
 *
 *  M(a,b,x)=sum ((a)_n/(b)_n) (x^n / n!)
 *
 * and expand (b)_n in inverse powers of b as follows
 *
 * -log(1/(b)_n) = sum_(k=0)^(n-1) log(b+k)
 *             = n log(b) + sum_(k=0)^(n-1) log(1+k/b)
 *
 * Do a taylor expansion of the log in 1/b and sum the resulting terms
 * using the standard algebraic formulas for finite sums of powers of
 * k.  This should then give
 *
 * M(a,b,x) = sum_(n=0)^(inf) (a_n/n!) (x/b)^n * (1 - n(n-1)/(2b) 
 *                          + (n-1)n(n+1)(3n-2)/(24b^2) + ...
 *
 * which can be summed explicitly. The trick for summing it is to take
 * derivatives of sum_(i=0)^(inf) a_n*y^n/n! = (1-y)^(-a);
 *
 * [BJG 16/01/2007]
 */

static 
int
hyperg_1F1_largebx(const double a, const double b, const double x, gsl_sf_result * result)
{
  double y = x/b;
  double f = exp(-a*log1p(-y));
  double t1 = -((a*(a+1.0))/(2*b))*pow((y/(1.0-y)),2.0);
  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);
  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)));
  result->val = f * (1 + t1 + t2 + t3);
  result->err = 2*fabs(f*t3) + 2*GSL_DBL_EPSILON*fabs(result->val);
  return GSL_SUCCESS;
}
 
/* Asymptotic result for x < 2b-4a, 2b-4a large.
 * [Abramowitz+Stegun, 13.5.21]
 *
 * assumes 0 <= x/(2b-4a) <= 1
 */
static
int
hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result)
{
  double eta    = 2.0*b - 4.0*a;
  double cos2th = x/eta;
  double sin2th = 1.0 - cos2th;
  double th = acos(sqrt(cos2th));
  double pre_h  = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
  gsl_sf_result lg_b;
  int stat_lg = gsl_sf_lngamma_e(b, &lg_b);
  double t1 = 0.5*(1.0-b)*log(0.25*x*eta);
  double t2 = 0.25*log(pre_h);
  double lnpre_val = lg_b.val + 0.5*x + t1 - t2;
  double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2));
#if SMALL_ANGLE
  const double eps = asin(sqrt(cos2th));  /* theta = pi/2 - eps */
  double s1 = (fmod(a, 1.0) == 0.0) ? 0.0 : sin(a*M_PI);
  double eta_reduc = (fmod(eta + 1, 4.0) == 0.0) ? 0.0 : fmod(eta + 1, 8.0);
  double phi1 = 0.25*eta_reduc*M_PI;
  double phi2 = 0.25*eta*(2*eps + sin(2.0*eps));
  double s2 = sin(phi1 - phi2);
#else
  double s1 = sin(a*M_PI);
  double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
#endif
  double ser_val = s1 + s2;
  double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2));
  int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
                                        ser_val, ser_err,
                                        result);
  return GSL_ERROR_SELECT_2(stat_e, stat_lg);
}


/* Luke's rational approximation.
 * See [Luke, Algorithms for the Computation of Mathematical Functions, p.182]
 *
 * Like the case of the 2F1 rational approximations, these are
 * probably guaranteed to converge for x < 0, barring gross
 * numerical instability in the pre-asymptotic regime.
 */
static
int
hyperg_1F1_luke(const double a, const double c, const double xin,
                gsl_sf_result * result)
{
  const double RECUR_BIG = 1.0e+50;
  const int nmax = 5000;
  int n = 3;
  const double x  = -xin;
  const double x3 = x*x*x;
  const double t0 = a/c;
  const double t1 = (a+1.0)/(2.0*c);
  const double t2 = (a+2.0)/(2.0*(c+1.0));
  double F = 1.0;
  double prec;

  double Bnm3 = 1.0;                                  /* B0 */
  double Bnm2 = 1.0 + t1 * x;                         /* B1 */
  double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */
 
  double Anm3 = 1.0;                                                      /* A0 */
  double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
  double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */

  while(1) {
    double npam1 = n + a - 1;
    double npcm1 = n + c - 1;
    double npam2 = n + a - 2;
    double npcm2 = n + c - 2;
    double tnm1  = 2*n - 1;
    double tnm3  = 2*n - 3;
    double tnm5  = 2*n - 5;
    double F1 =  (n-a-2) / (2*tnm3*npcm1);
    double F2 =  (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1);
    double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
    double E  = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1);

    double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
    double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
    double r = An/Bn;

    prec = fabs((F - r)/F);
    F = r;

    if(prec < GSL_DBL_EPSILON || n > nmax) break;

    if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
      An   /= RECUR_BIG;
      Bn   /= RECUR_BIG;
      Anm1 /= RECUR_BIG;
      Bnm1 /= RECUR_BIG;
      Anm2 /= RECUR_BIG;
      Bnm2 /= RECUR_BIG;
      Anm3 /= RECUR_BIG;
      Bnm3 /= RECUR_BIG;
    }
    else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
      An   *= RECUR_BIG;
      Bn   *= RECUR_BIG;
      Anm1 *= RECUR_BIG;
      Bnm1 *= RECUR_BIG;
      Anm2 *= RECUR_BIG;
      Bnm2 *= RECUR_BIG;
      Anm3 *= RECUR_BIG;
      Bnm3 *= RECUR_BIG;
    }

    n++;
    Bnm3 = Bnm2;
    Bnm2 = Bnm1;
    Bnm1 = Bn;
    Anm3 = Anm2;
    Anm2 = Anm1;
    Anm1 = An;
  }

  result->val  = F;
  result->err  = 2.0 * fabs(F * prec);
  result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F);

  return GSL_SUCCESS;
}

/* Series for 1F1(1,b,x)
 * b > 0
 */
static
int
hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result)
{
  double sum_val = 1.0;
  double sum_err = 0.0;
  double term = 1.0;
  double n    = 1.0;
  while(fabs(term/sum_val) > 0.25*GSL_DBL_EPSILON) {
    term *= x/(b+n-1);
    sum_val += term;
    sum_err += 8.0*GSL_DBL_EPSILON*fabs(term) + GSL_DBL_EPSILON*fabs(sum_val);
    n += 1.0;
  }
  result->val  = sum_val;
  result->err  = sum_err;
  result->err += 2.0 *  fabs(term);
  return GSL_SUCCESS;
}


/* 1F1(1,b,x)
 * b >= 1, b integer
 */
static
int
hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result)
{
  if(b < 1) {
    DOMAIN_ERROR(result);
  }
  else if(b == 1) {
    return gsl_sf_exp_e(x, result);
  }
  else if(b == 2) {
    return gsl_sf_exprel_e(x, result);
  }
  else if(b == 3) {
    return gsl_sf_exprel_2_e(x, result);
  }
  else {
    return gsl_sf_exprel_n_e(b-1, x, result);
  }
}


/* 1F1(1,b,x)
 * b >=1, b real
 *
 * checked OK: [GJ] Thu Oct  1 16:46:35 MDT 1998
 */
static
int
hyperg_1F1_1(const double b, const double x, gsl_sf_result * result)
{
  double ax = fabs(x);
  double ib = floor(b + 0.1);

  if(b < 1.0) {
    DOMAIN_ERROR(result);
  }
  else if(b == 1.0) {
    return gsl_sf_exp_e(x, result);
  }
  else if(b >= 1.4*ax) {
    return hyperg_1F1_1_series(b, x, result);
  }
  else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) {
    return hyperg_1F1_1_int((int)ib, x, result);
  }
  else if(x > 0.0) {
    if(x > 100.0 && b < 0.75*x) {
      return hyperg_1F1_asymp_posx(1.0, b, x, result);
    }
    else if(b < 1.0e+05) {
      /* Recurse backward on b, from a
       * chosen offset point. For x > 0,
       * which holds here, this should
       * be a stable direction.
       */
      const double off = ceil(1.4*x-b) + 1.0;
      double bp = b + off;
      gsl_sf_result M;
      int stat_s = hyperg_1F1_1_series(bp, x, &M);
      const double err_rat = M.err / fabs(M.val);
      while(bp > b+0.1) {
        /* M(1,b-1) = x/(b-1) M(1,b) + 1 */
        bp -= 1.0;
        M.val  = 1.0 + x/bp * M.val;
      }
      result->val  = M.val;
      result->err  = err_rat * fabs(M.val);
      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val);
      return stat_s;
    } else if (fabs(x) < fabs(b) && fabs(x) < sqrt(fabs(b)) * fabs(b-x)) {
      return hyperg_1F1_largebx(1.0, b, x, result);
    } else if (fabs(x) > fabs(b)) {
      return hyperg_1F1_1_series(b, x, result);
    } else {
      return hyperg_1F1_large2bm4a(1.0, b, x, result);
    }
  }
  else {
    /* x <= 0 and b not large compared to |x|
     */
    if(ax < 10.0 && b < 10.0) {
      return hyperg_1F1_1_series(b, x, result);
    }
    else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) {
      return hyperg_1F1_asymp_negx(1.0, b, x, result);
    }
    else {
      return hyperg_1F1_luke(1.0, b, x, result);
    }
  }
}


/* 1F1(a,b,x)/Gamma(b) for b->0
 * [limit of Abramowitz+Stegun 13.3.7]
 */
static