hyperg_u.c

开放gsl矩阵运算

/* specfunc/hyperg_U.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_exp.h"
#include "gsl_sf_gamma.h"
#include "gsl_sf_bessel.h"
#include "gsl_sf_laguerre.h"
#include "gsl_sf_pow_int.h"
#include "gsl_sf_hyperg.h"

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

#define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON)

#define SERIES_EVAL_OK(a,b,x) ((fabs(a) < 5 && b < 5 && x < 2.0) || (fabs(a) <  10 && b < 10 && x < 1.0))

#define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x))


/* Log[U(a,2a,x)]
 * [Abramowitz+stegun, 13.6.21]
 * Assumes x > 0, a > 1/2.
 */
static
int
hyperg_lnU_beq2a(const double a, const double x, gsl_sf_result * result)
{
  const double lx = log(x);
  const double nu = a - 0.5;
  const double lnpre = 0.5*(x - M_LNPI) - nu*lx;
  gsl_sf_result lnK;
  gsl_sf_bessel_lnKnu_e(nu, 0.5*x, &lnK);
  result->val  = lnpre + lnK.val;
  result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + 0.5*M_LNPI + fabs(nu*lx));
  result->err += lnK.err;
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  return GSL_SUCCESS;
}


/* Evaluate u_{N+1}/u_N by Steed's continued fraction method.
 *
 * u_N := Gamma[a+N]/Gamma[a] U(a + N, b, x)
 *
 * u_{N+1}/u_N = (a+N) U(a+N+1,b,x)/U(a+N,b,x)
 */
static
int
hyperg_U_CF1(const double a, const double b, const int N, const double x,
             double * result, int * count)
{
  const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  const int maxiter = 20000;
  int n = 1;
  double Anm2 = 1.0;
  double Bnm2 = 0.0;
  double Anm1 = 0.0;
  double Bnm1 = 1.0;
  double a1 = -(a + N);
  double b1 =  (b - 2.0*a - x - 2.0*(N+1));
  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 + n - b)*(a + N + n - 1.0);
    bn =  (b - 2.0*a - x - 2.0*(N+n));
    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 = fn;
  *count  = n;

  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}


/* Large x asymptotic for  x^a U(a,b,x)
 * Based on SLATEC D9CHU() [W. Fullerton]
 *
 * Uses a rational approximation due to Luke.
 * See [Luke, Algorithms for the Computation of Special Functions, p. 252]
 *     [Luke, Utilitas Math. (1977)]
 *
 * z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
 *
 * This assumes that a is not a negative integer and
 * that 1+a-b is not a negative integer. If one of them
 * is, then the 2F0 actually terminates, the above
 * relation is an equality, and the sum should be
 * evaluated directly [see below].
 */
static
int
d9chu(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double EPS   = 8.0 * GSL_DBL_EPSILON;  /* EPS = 4.0D0*D1MACH(4)   */
  const int maxiter = 500;
  double aa[4], bb[4];
  int i;

  double bp = 1.0 + a - b;
  double ab = a*bp;
  double ct2 = 2.0 * (x - ab);
  double sab = a + bp;
  
  double ct3 = sab + 1.0 + ab;
  double anbn = ct3 + sab + 3.0;
  double ct1 = 1.0 + 2.0*x/anbn;

  bb[0] = 1.0;
  aa[0] = 1.0;

  bb[1] = 1.0 + 2.0*x/ct3;
  aa[1] = 1.0 + ct2/ct3;
  
  bb[2] = 1.0 + 6.0*ct1*x/ct3;
  aa[2] = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3;

  for(i=4; i<maxiter; i++) {
    int j;
    double c2;
    double d1z;
    double g1, g2, g3;
    double x2i1 = 2*i - 3;
    ct1   = x2i1/(x2i1-2.0);
    anbn += x2i1 + sab;
    ct2   = (x2i1 - 1.0)/anbn;
    c2    = x2i1*ct2 - 1.0;
    d1z   = 2.0*x2i1*x/anbn;
    
    ct3 = sab*ct2;
    g1  = d1z + ct1*(c2+ct3);
    g2  = d1z - c2;
    g3  = ct1*(1.0 - ct3 - 2.0*ct2);
    
    bb[3] = g1*bb[2] + g2*bb[1] + g3*bb[0];
    aa[3] = g1*aa[2] + g2*aa[1] + g3*aa[0];
    
    if(fabs(aa[3]*bb[0]-aa[0]*bb[3]) < EPS*fabs(bb[3]*bb[0])) break;
    
    for(j=0; j<3; j++) {
      aa[j] = aa[j+1];
      bb[j] = bb[j+1];
    }
  }
  
  result->val = aa[3]/bb[3];
  result->err = 8.0 * GSL_DBL_EPSILON * fabs(result->val);
  
  if(i == maxiter) {
    GSL_ERROR ("error", GSL_EMAXITER);
  }
  else {
    return GSL_SUCCESS;
  }
}


/* Evaluate asymptotic for z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
 * We check for termination of the 2F0 as a special case.
 * Assumes x > 0.
 * Also assumes a,b are not too large compared to x.
 */
static
int
hyperg_zaU_asymp(const double a, const double b, const double x, gsl_sf_result *result)
{
  const double ap = a;
  const double bp = 1.0 + a - b;
  const double rintap = floor(ap + 0.5);
  const double rintbp = floor(bp + 0.5);
  const int ap_neg_int = ( ap < 0.0 && fabs(ap - rintap) < INT_THRESHOLD );
  const int bp_neg_int = ( bp < 0.0 && fabs(bp - rintbp) < INT_THRESHOLD );

  if(ap_neg_int || bp_neg_int) {
    /* Evaluate 2F0 polynomial.
     */
    double mxi = -1.0/x;
    double nmax = -(int)(GSL_MIN(ap,bp) - 0.1);
    double tn  = 1.0;
    double sum = 1.0;
    double n   = 1.0;
    double sum_err = 0.0;
    while(n <= nmax) {
      double apn = (ap+n-1.0);
      double bpn = (bp+n-1.0);
      tn  *= ((apn/n)*mxi)*bpn;
      sum += tn;
      sum_err += 2.0 * GSL_DBL_EPSILON * fabs(tn);
      n += 1.0;
    }
    result->val  = sum;
    result->err  = sum_err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(nmax)+1.0) * fabs(sum);
    return GSL_SUCCESS;
  }
  else {
    return d9chu(a,b,x,result);
  }
}


/* Evaluate finite sum which appears below.
 */
static
int
hyperg_U_finite_sum(int N, double a, double b, double x, double xeps,
                    gsl_sf_result * result)
{
  int i;
  double sum_val;
  double sum_err;

  if(N <= 0) {
    double t_val = 1.0;
    double t_err = 0.0;
    gsl_sf_result poch;
    int stat_poch;

    sum_val = 1.0;
    sum_err = 0.0;
    for(i=1; i<= -N; i++) {
      const double xi1  = i - 1;
      const double mult = (a+xi1)*x/((b+xi1)*(xi1+1.0));
      t_val *= mult;
      t_err += fabs(mult) * t_err + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
      sum_val += t_val;
      sum_err += t_err;
    }

    stat_poch = gsl_sf_poch_e(1.0+a-b, -a, &poch);

    result->val  = sum_val * poch.val;
    result->err  = fabs(sum_val) * poch.err + sum_err * fabs(poch.val);
    result->err += fabs(poch.val) * (fabs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
    return stat_poch;
  }
  else {
    const int M = N - 2;
    if(M < 0) {
      result->val = 0.0;
      result->err = 0.0;
      return GSL_SUCCESS;
    }
    else {
      gsl_sf_result gbm1;
      gsl_sf_result gamr;
      int stat_gbm1;
      int stat_gamr;
      double t_val = 1.0;
      double t_err = 0.0;

      sum_val = 1.0;
      sum_err = 0.0;
      for(i=1; i<=M; i++) {
        const double mult = (a-b+i)*x/((1.0-b+i)*i);
        t_val *= mult;
	t_err += t_err * fabs(mult) + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
        sum_val += t_val;
	sum_err += t_err;
      }

      stat_gbm1 = gsl_sf_gamma_e(b-1.0, &gbm1);
      stat_gamr = gsl_sf_gammainv_e(a,  &gamr);

      if(stat_gbm1 == GSL_SUCCESS) {
        gsl_sf_result powx1N;
	int stat_p = gsl_sf_pow_int_e(x, 1-N, &powx1N);
        double pe_val = powx1N.val * xeps;
	double pe_err = powx1N.err * fabs(xeps) + 2.0 * GSL_DBL_EPSILON * fabs(pe_val);
        double coeff_val = gbm1.val * gamr.val * pe_val;
	double coeff_err = gbm1.err * fabs(gamr.val * pe_val)
                         + gamr.err * fabs(gbm1.val * pe_val)
                         + fabs(gbm1.val * gamr.val) * pe_err
                         + 2.0 * GSL_DBL_EPSILON * fabs(coeff_val);

	result->val  = sum_val * coeff_val;
	result->err  = fabs(sum_val) * coeff_err + sum_err * fabs(coeff_val);
	result->err += 2.0 * GSL_DBL_EPSILON * (M+2.0) * fabs(result->val);
	result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
	return stat_p;
      }
      else {
        result->val = 0.0;
	result->err = 0.0;
        return stat_gbm1;
      }
    }
  }
}


/* Based on SLATEC DCHU() [W. Fullerton]
 * Assumes x > 0.
 * This is just a series summation method, and
 * it is not good for large a.
 *
 * I patched up the window for 1+a-b near zero. [GJ]
 */
static
int
hyperg_U_series(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */
  const double SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON;

  if(fabs(1.0 + a - b) < SQRT_EPS) {
    /* Original Comment: ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X
     */
    /* We can however do the following:
     * U(a,b,x) = U(a,a+1,x) when 1+a-b=0
     * and U(a,a+1,x) = x^(-a).
     */
    double lnr = -a * log(x);
    int stat_e =  gsl_sf_exp_e(lnr, result);
    result->err += 2.0 * SQRT_EPS * fabs(result->val);
    return stat_e;
  }
  else {
    double aintb = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) );
    double beps  = b - aintb;
    int N = aintb;
    
    double lnx  = log(x);
    double xeps = exp(-beps*lnx);

    /* Evaluate finite sum.
     */
    gsl_sf_result sum;
    int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum);


    /* Evaluate infinite sum. */

    int istrt = ( N < 1 ? 1-N : 0 );
    double xi = istrt;

    gsl_sf_result gamr;
    gsl_sf_result powx;
    int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr);
    int stat_powx = gsl_sf_pow_int_e(x, istrt, &powx);
    double sarg   = beps*M_PI;
    double sfact  = ( sarg != 0.0 ? sarg/sin(sarg) : 1.0 );
    double factor_val = sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * gamr.val * powx.val;
    double factor_err = fabs(gamr.val) * powx.err + fabs(powx.val) * gamr.err
                      + 2.0 * GSL_DBL_EPSILON * fabs(factor_val);

    gsl_sf_result pochai;
    gsl_sf_result gamri1;
    gsl_sf_result gamrni;
    int stat_pochai = gsl_sf_poch_e(a, xi, &pochai);
    int stat_gamri1 = gsl_sf_gammainv_e(xi + 1.0, &gamri1);
    int stat_gamrni = gsl_sf_gammainv_e(aintb + xi, &gamrni);
    int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);
    int stat_gamall = GSL_ERROR_SELECT_4(stat_sum, stat_gam123, stat_pochai, stat_powx);

    gsl_sf_result pochaxibeps;
    gsl_sf_result gamrxi1beps;
    int stat_pochaxibeps = gsl_sf_poch_e(a, xi-beps, &pochaxibeps);
    int stat_gamrxi1beps = gsl_sf_gammainv_e(xi + 1.0 - beps, &gamrxi1beps);

    int stat_all = GSL_ERROR_SELECT_3(stat_gamall, stat_pochaxibeps, stat_gamrxi1beps);

    double b0_val = factor_val * pochaxibeps.val * gamrni.val * gamrxi1beps.val;
    double b0_err =  fabs(factor_val * pochaxibeps.val * gamrni.val) * gamrxi1beps.err
                   + fabs(factor_val * pochaxibeps.val * gamrxi1beps.val) * gamrni.err
		   + fabs(factor_val * gamrni.val * gamrxi1beps.val) * pochaxibeps.err
		   + fabs(pochaxibeps.val * gamrni.val * gamrxi1beps.val) * factor_err
                   + 2.0 * GSL_DBL_EPSILON * fabs(b0_val);

    if(fabs(xeps-1.0) < 0.5) {
      /*
       C  X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE
       C  CAREFUL IN EVALUATING THE DIFFERENCES.
       */
      int i;
      gsl_sf_result pch1ai;
      gsl_sf_result pch1i;
      gsl_sf_result poch1bxibeps;
      int stat_pch1ai = gsl_sf_pochrel_e(a + xi, -beps, &pch1ai);
      int stat_pch1i  = gsl_sf_pochrel_e(xi + 1.0 - beps, beps, &pch1i);
      int stat_poch1bxibeps = gsl_sf_pochrel_e(b+xi, -beps, &poch1bxibeps);
      double c0_t1_val = beps*pch1ai.val*pch1i.val;
      double c0_t1_err = fabs(beps) * fabs(pch1ai.val) * pch1i.err
                       + fabs(beps) * fabs(pch1i.val)  * pch1ai.err
                       + 2.0 * GSL_DBL_EPSILON * fabs(c0_t1_val);
      double c0_t2_val = -poch1bxibeps.val + pch1ai.val - pch1i.val + c0_t1_val;
      double c0_t2_err =  poch1bxibeps.err + pch1ai.err + pch1i.err + c0_t1_err
                       + 2.0 * GSL_DBL_EPSILON * fabs(c0_t2_val);
      double c0_val = factor_val * pochai.val * gamrni.val * gamri1.val * c0_t2_val;
      double c0_err =  fabs(factor_val * pochai.val * gamrni.val * gamri1.val) * c0_t2_err
                     + fabs(factor_val * pochai.val * gamrni.val * c0_t2_val) * gamri1.err
		     + fabs(factor_val * pochai.val * gamri1.val * c0_t2_val) * gamrni.err
		     + fabs(factor_val * gamrni.val * gamri1.val * c0_t2_val) * pochai.err
		     + fabs(pochai.val * gamrni.val * gamri1.val * c0_t2_val) * factor_err
                     + 2.0 * GSL_DBL_EPSILON * fabs(c0_val);
      /*
       C  XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS)
       */
      gsl_sf_result dexprl;
      int stat_dexprl = gsl_sf_exprel_e(-beps*lnx, &dexprl);
      double xeps1_val = lnx * dexprl.val;
      double xeps1_err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(beps*lnx)) * fabs(dexprl.val)
                       + fabs(lnx) * dexprl.err
                       + 2.0 * GSL_DBL_EPSILON * fabs(xeps1_val);
      double dchu_val = sum.val + c0_val + xeps1_val*b0_val;
      double dchu_err = sum.err + c0_err
                      + fabs(xeps1_val)*b0_err + xeps1_err * fabs(b0_val)
                      + fabs(b0_val*lnx)*dexprl.err
                      + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(c0_val) + fabs(xeps1_val*b0_val));
      double xn = N;
      double t_val;
      double t_err;

      stat_all = GSL_ERROR_SELECT_5(stat_all, stat_dexprl, stat_poch1bxibeps, stat_pch1i, stat_pch1ai);

      for(i=1; i<2000; i++) {