/* specfunc/hyperg_U.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman

 * Copyright (C) 2009, 2010 Brian Gough

 * 

 * 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_exp.h>

#include <gsl/gsl_sf_gamma.h>

#include <gsl/gsl_sf_bessel.h>

#include <gsl/gsl_sf_laguerre.h>

#include <gsl/gsl_sf_pow_int.h>

#include <gsl/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

    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;

      }

    }

  }

}

/* Evaluate infinite sum which appears below.

 */

static

int

hyperg_U_infinite_sum_stable(int N, double a, double bint, double b, double beps, double x, double xeps, gsl_sf_result sum,

                             gsl_sf_result * result)

{

  const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */

    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(bint + xi, &gamrni);

    int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);

    int stat_gamall = GSL_ERROR_SELECT_3(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);

      /*

       C  X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE

       C  STRAIGHTFORWARD FORMULATION IS STABLE.

       */

      int i;

      double dchu_val;

      double dchu_err;

      double t_val;

      double t_err;

      gsl_sf_result dgamrbxi;

      int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);

      double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps;

      double a0_err =  fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err

                     + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err

                     + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err

                     + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err

                     + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);

      stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi);

      b0_val = xeps * b0_val / beps;

      b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);

      dchu_val = sum.val + a0_val - b0_val;

      dchu_err = sum.err + a0_err + b0_err

        + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));

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

        double xi = istrt + i;

        double xi1 = istrt + i - 1;

        double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);

        double b0_multiplier = (a+xi1-beps)*x/((bint+xi1)*(xi-beps));

        a0_val *= a0_multiplier;

        a0_err += fabs(a0_multiplier) * a0_err;

        b0_val *= b0_multiplier;

        b0_err += fabs(b0_multiplier) * b0_err;

        t_val = a0_val - b0_val;

        t_err = a0_err + b0_err;

        dchu_val += t_val;

        dchu_err += t_err;

        if(fabs(t_val) < EPS*fabs(dchu_val)) break;

      }

      result->val  = dchu_val;

      result->err  = 2.0 * dchu_err;

      result->err += 2.0 * fabs(t_val);

      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);

      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {

        GSL_ERROR ("error", GSL_EMAXITER);

      }

      else {

        return stat_all;

      }

}

static

int

hyperg_U_infinite_sum_simple(int N, double a, double bint, double b, double beps, double x, double xeps, gsl_sf_result sum,

                             gsl_sf_result * result)

{

  const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */

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

    double xi = istrt;

    gsl_sf_result powx;

    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 ) * powx.val;

    double factor_err = fabs(powx.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(bint + xi, &gamrni);

    int stat_gam123 = GSL_ERROR_SELECT_2(stat_gamri1, stat_gamrni);

    int stat_gamall = GSL_ERROR_SELECT_3(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 X =  sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * powx.val * gsl_sf_poch(1 + a - b, xi - 1 + b - beps) * gsl_sf_gammainv(a);

    double b0_val = X * 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);

      /*

       C  X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE

       C  STRAIGHTFORWARD FORMULATION IS STABLE.

       */

      int i;

      double dchu_val;

      double dchu_err;

      double t_val;

      double t_err;

      gsl_sf_result gamr;

      gsl_sf_result dgamrbxi;

      int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr;;

      int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);

      double a0_val = factor_val * gamr.val * pochai.val * dgamrbxi.val * gamri1.val / beps;

      double a0_err =  fabs(factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps) * gamr.err

                     + fabs(factor_val * gamr.val * dgamrbxi.val * gamri1.val / beps) * pochai.err

                     + fabs(factor_val * gamr.val * pochai.val * gamri1.val / beps) * dgamrbxi.err

                     + fabs(factor_val * gamr.val * pochai.val * dgamrbxi.val / beps) * gamri1.err

                     + fabs(pochai.val * gamr.val * dgamrbxi.val * gamri1.val / beps) * factor_err

                     + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);

      stat_all = GSL_ERROR_SELECT_3(stat_all, stat_gamr, stat_dgamrbxi);

      b0_val = xeps * b0_val / beps;

      b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);

      dchu_val = sum.val + a0_val - b0_val;

      dchu_err = sum.err + a0_err + b0_err

        + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));

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

        double xi = istrt + i;

        double xi1 = istrt + i - 1;

        double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);

        double b0_multiplier = (a+xi1-beps)*x/((bint+xi1)*(xi-beps));

        a0_val *= a0_multiplier;

        a0_err += fabs(a0_multiplier) * a0_err;

        b0_val *= b0_multiplier;

        b0_err += fabs(b0_multiplier) * b0_err;

        t_val = a0_val - b0_val;

        t_err = a0_err + b0_err;

        dchu_val += t_val;

        dchu_err += t_err;

        if(!gsl_finite(t_val) || fabs(t_val) < EPS*fabs(dchu_val)) break;

      }

      result->val  = dchu_val;

      result->err  = 2.0 * dchu_err;

      result->err += 2.0 * fabs(t_val);

      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);

      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {

        GSL_ERROR ("error", GSL_EMAXITER);

      }

      else {

        return stat_all;

      }

}

static

int

hyperg_U_infinite_sum_improved(int N, double a, double bint, double b, double beps, double x, double xeps, gsl_sf_result sum,

                               gsl_sf_result * result)

{

  const double EPS      = 2.0 * GSL_DBL_EPSILON;  /* EPS = D1MACH(3) */

  const double lnx = log(x);

    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(bint + xi, &gamrni);

    int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);

    int stat_gamall = GSL_ERROR_SELECT_3(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);

      /*

       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++) {

        const double xi  = istrt + i;

        const double xi1 = istrt + i - 1;

        const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps);

        const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps));

        const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi);

        const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps));

        b0_val *= b0_multiplier;

        b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;

        c0_val  = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val;

        c0_err  =  fabs(c0_multiplier_1) * c0_err

                 + fabs(c0_multiplier_2) * b0_err

                 + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON

                 + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON;

        t_val  = c0_val + xeps1_val*b0_val;

        t_err  = c0_err + fabs(xeps1_val)*b0_err;

        t_err += fabs(b0_val*lnx) * dexprl.err;

        t_err += fabs(b0_val)*xeps1_err;

        dchu_val += t_val;

        dchu_err += t_err;

        if(fabs(t_val) < EPS*fabs(dchu_val)) break;

      }

      result->val  = dchu_val;

      result->err  = 2.0 * dchu_err;

      result->err += 2.0 * fabs(t_val);

      result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);

      result->err *= 2.0; /* FIXME: fudge factor */

      if(i >= 2000) {

        GSL_ERROR ("error", GSL_EMAXITER);

      }

      else {

        return stat_all;

      }

}

/* 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 SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON;

double bint = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) );

  double beps  = b - bint;

  double a_beps = a - beps;

  double r_a_beps = floor(a_beps + 0.5);

  double a_beps_int = ( fabs(a_beps - r_a_beps) < INT_THRESHOLD );

/*  double a_b_1 = a-b+1;

  double r_a_b_1 = floor(a_b_1+0.5);

  double r_a_b_1_int = (fabs(a_b_1-r_a_b_1)< INT_THRESHOLD);

  Check for (a-beps) being a member of -N; N being 0,1,... */

  if (a_beps_int && a_beps <= 0)

  { 	beps=beps - 1 + floor(a_beps);bint=bint + 1 - floor(a_beps);

  }

  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 {

     int N = (int) bint;

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

    int stat_inf;

    /* Evaluate infinite sum. */

    if(fabs(xeps-1.0) > 0.5 ) {

      stat_inf = hyperg_U_infinite_sum_stable(N, a, bint, b, beps, x, xeps, sum, result);

    } else if (1+a-b < 0 && 1+a-b==floor(1+a-b) && beps != 0) {

       stat_inf = hyperg_U_infinite_sum_simple(N, a, bint, b, beps, x, xeps, sum, result);

    } else {

      stat_inf = hyperg_U_infinite_sum_improved(N, a, bint, b, beps, x, xeps, sum, result);

    }

    return GSL_ERROR_SELECT_2(stat_sum, stat_inf);

  }

}

/* Assumes b > 0 and x > 0.

 */

static

int

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

{

  if(a == -1.0) {

    /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x

     */

    result->val  = -b + x;

    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));

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

    return GSL_SUCCESS;

  }

  else if(a == 0.0) {

    /* U(0,c+1,x) = Laguerre[c,0,x] = 1

     */

    result->val = 1.0;

    result->err = 0.0;

    return GSL_SUCCESS;

  }

  else if(ASYMP_EVAL_OK(a,b,x)) {

    double p = pow(x, -a);

    gsl_sf_result asymp;

    int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);

    result->val  = asymp.val * p;

    result->err  = asymp.err * p;

    result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p;

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

    return stat_asymp;

  }

  else {

    return hyperg_U_series(a, b, x, result);

  }

}

/* Assumes b > 0 and x > 0.

 */

static

int

hyperg_U_small_a_bgt0(const double a, const double b, const double x,

                      gsl_sf_result * result,

                      double * ln_multiplier

                      )

{

  if(a == 0.0) {

    result->val = 1.0;

    result->err = 0.0;

    *ln_multiplier = 0.0;

    return GSL_SUCCESS;

  }

  else if(   (b > 5000.0 && x < 0.90 * fabs(b))

          || (b >  500.0 && x < 0.50 * fabs(b))

    ) {

    int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier);

    if(stat == GSL_EOVRFLW)

      return GSL_SUCCESS;

    else

      return stat;

  }

  else if(b > 15.0) {

    /* Recurse up from b near 1.

     */

    double eps = b - floor(b);

    double b0  = 1.0 + eps;

    gsl_sf_result r_Ubm1;

    gsl_sf_result r_Ub;

    int stat_0 = hyperg_U_small_ab(a, b0,     x, &r_Ubm1);

    int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub);

    double Ubm1 = r_Ubm1.val;

    double Ub   = r_Ub.val;

    double Ubp1;

    double bp;

    for(bp = b0+1.0; bp

      Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;

      Ubm1 = Ub;

      Ub   = Ubp1;

    }

    result->val  = Ub;

    result->err  = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub);

    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub);

    *ln_multiplier = 0.0;

    return GSL_ERROR_SELECT_2(stat_0, stat_1);

  }

  else {

    *ln_multiplier = 0.0;

    return hyperg_U_small_ab(a, b, x, result);

  }

}

/* We use this to keep track of large

 * dynamic ranges in the recursions.

 * This can be important because sometimes

 * we want to calculate a very large and

 * a very small number and the answer is

 * the product, of order 1. This happens,

 * for instance, when we apply a Kummer

 * transform to make b positive and

 * both x and b are large.

 */

#define RESCALE_2(u0,u1,factor,count)      \

do {                                       \

  double au0 = fabs(u0);                   \

  if(au0 > factor) {                       \

    u0 /= factor;                          \

    u1 /= factor;                          \

    count++;                               \

  }                                        \

  else if(au0 < 1.0/factor) {              \

    u0 *= factor;                          \

    u1 *= factor;                          \

    count--;                               \

  }                                        \

} while (0)

/* Specialization to b >= 1, for integer parameters.

 * Assumes x > 0.

 */

static

int

hyperg_U_int_bge1(const int a, const int b, const double x,

                  gsl_sf_result_e10 * result)

{

  if(a == 0) {

    result->val = 1.0;

    result->err = 0.0;

    result->e10 = 0;

    return GSL_SUCCESS;

  }

  else if(a == -1) {

    result->val  = -b + x;

    result->err  = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));

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

    result->e10  = 0;

    return GSL_SUCCESS;

  }

  else if(b == a + 1) {

    /* U(a,a+1,x) = x^(-a)

     */

    return gsl_sf_exp_e10_e(-a*log(x), result);

  }