hyperg_1f1.c

开放gsl矩阵运算

/* 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 2 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_elementary.h"
#include "gsl_sf_exp.h"
#include "gsl_sf_bessel.h"
#include "gsl_sf_gamma.h"
#include "gsl_sf_laguerre.h"
#include "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 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));
  double s1 = sin(a*M_PI);
  double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
  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) > 2.0*GSL_DBL_EPSILON) {
    term *= x/(b+n-1);
    sum_val += term;
    sum_err += 2.0 * 4.0 * GSL_DBL_EPSILON * fabs(term);
    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 {
      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
int
hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result)
{
  double eta = a*x;
  if(eta > 0.0) {
    double root_eta = sqrt(eta);
    gsl_sf_result I1_scaled;
    int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled);
    if(I1_scaled.val <= 0.0) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM);
    }
    else {
      const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(x) + log(I1_scaled.val);
      const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(I1_scaled.err/I1_scaled.val);
      return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
    }
  }
  else if(eta == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* eta < 0 */
    double root_eta = sqrt(-eta);
    gsl_sf_result J1;
    int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1);
    if(J1.val <= 0.0) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM);
    }
    else {
      const double t1 = 0.5*x;
      const double t2 = 0.5*log(-eta);
      const double t3 = fabs(x);
      const double t4 = log(J1.val);
      const double lnr_val = t1 + t2 + t3 + t4;
      const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val);
      gsl_sf_result ex;
      int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex);
      result->val = -ex.val;
      result->err =  ex.err;
      return stat_e;
    }
  }
  
}


/* 1F1'(a,b,x)/1F1(a,b,x)
 * Uses Gautschi's version of the CF.
 * [Gautschi, Math. Comp. 31, 994 (1977)]
 *
 * Supposedly this suffers from the "anomalous convergence"
 * problem when b < x. I have seen anomalous convergence
 * in several of the continued fractions associated with
 * 1F1(a,b,x). This particular CF formulation seems stable
 * for b > x. However, it does display a painful artifact
 * of the anomalous convergence; the convergence plateaus
 * unless b >>> x. For example, even for b=1000, x=1, this
 * method locks onto a ratio which is only good to about
 * 4 digits. Apparently the rest of the digits are hiding
 * way out on the plateau, but finite-precision lossage
 * means you will never get them.
 */
#if 0
static
int
hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 5000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = 1.0;
  double b1 = 1.0;
  double An = b1*Anm1 + a1*Anm2;
  double Bn = b1*Bnm1 + a1*Bnm2;
  double an, bn;
  double fn = An/Bn;

  while(n < maxiter) {
    double old_fn;
    double del;
    n++;
    Anm2 = Anm1;
    Bnm2 = Bnm1;
    Anm1 = An;
    Bnm1 = Bn;
    an = (a+n)*x/((b-x+n-1)*(b-x+n));
    bn = 1.0;
    An = bn*Anm1 + an*Anm2;
    Bn = bn*Bnm1 + an*Bnm2;

    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;
    }

    old_fn = fn;
    fn = An/Bn;
    del = old_fn/fn;
    
    if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
  }

  *result = a/(b-x) * fn;

  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */


/* 1F1'(a,b,x)/1F1(a,b,x)
 * Uses Gautschi's series transformation of the
 * continued fraction. This is apparently the best
 * method for getting this ratio in the stable region.
 * The convergence is monotone and supergeometric
 * when b > x.
 * Assumes a >= -1.