/* specfunc/poch.c

 * 

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

 * Copyright (C) 2009 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_log.h>

#include <gsl/gsl_sf_pow_int.h>

#include <gsl/gsl_sf_psi.h>

#include <gsl/gsl_sf_gamma.h>

#include "error.h"

static const double bern[21] = {

   0.0   /* no element 0 */,  

  +0.833333333333333333333333333333333e-01,

  -0.138888888888888888888888888888888e-02,

  +0.330687830687830687830687830687830e-04,

  -0.826719576719576719576719576719576e-06,

  +0.208767569878680989792100903212014e-07,

  -0.528419013868749318484768220217955e-09,

  +0.133825365306846788328269809751291e-10,

  -0.338968029632258286683019539124944e-12,

  +0.858606205627784456413590545042562e-14,

  -0.217486869855806187304151642386591e-15,

  +0.550900282836022951520265260890225e-17,

  -0.139544646858125233407076862640635e-18,

  +0.353470703962946747169322997780379e-20,

  -0.895351742703754685040261131811274e-22,

  +0.226795245233768306031095073886816e-23,

  -0.574472439520264523834847971943400e-24,

  +0.145517247561486490186626486727132e-26,

  -0.368599494066531017818178247990866e-28,

  +0.933673425709504467203255515278562e-30,

  -0.236502241570062993455963519636983e-31

};

/* ((a)_x - 1)/x in the "small x" region where

 * cancellation must be controlled.

 *

 * Based on SLATEC DPOCH1().

 */

/*

C When ABS(X) is so small that substantial cancellation will occur if

C the straightforward formula is used, we use an expansion due

C to Fields and discussed by Y. L. Luke, The Special Functions and Their

C Approximations, Vol. 1, Academic Press, 1969, page 34.

C

C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as

C        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .

C In order to maintain significance in POCH1, we write for positive a

C        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)

C                       = 1.0 + Q*EXPREL(Q) .

C Likewise the polynomial is written

C        POLY = 1.0 + X*POLY1(A,X) .

C Thus,

C        POCH1(A,X) = (POCH(A,X) - 1) / X

C                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)

C

*/

static

int

pochrel_smallx(const double a, const double x, gsl_sf_result * result)

{

  /*

   SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))

   ALNEPS = LOG(D1MACH(3))

   */

  const double SQTBIG = 1.0/(2.0*M_SQRT2*M_SQRT3*GSL_SQRT_DBL_MIN);

  const double ALNEPS = GSL_LOG_DBL_EPSILON - M_LN2;

  if(x == 0.0) {

    return gsl_sf_psi_e(a, result);

  }

  else {

    const double bp   = (  (a < -0.5) ? 1.0-a-x : a );

    const int    incr = ( (bp < 10.0) ? 11.0-bp : 0 );

    const double b    = bp + incr;

    double dpoch1;

    gsl_sf_result dexprl;

    int stat_dexprl;

    int i;

    double var    = b + 0.5*(x-1.0);

    double alnvar = log(var);

    double q = x*alnvar;

    double poly1 = 0.0;

    if(var < SQTBIG) {

      const int nterms = (int)(-0.5*ALNEPS/alnvar + 1.0);

      const double var2 = (1.0/var)/var;

      const double rho  = 0.5 * (x + 1.0);

      double term = var2;

      double gbern[24];

      int k, j;

      gbern[1] = 1.0;

      gbern[2] = -rho/12.0;

      poly1 = gbern[2] * term;

      if(nterms > 20) {

        /* NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD */

        /* nterms = 20; */

        result->val = 0.0;

        result->err = 0.0;

        GSL_ERROR ("error", GSL_ESANITY);

      }

      for(k=2; k<=nterms; k++) {

        double gbk = 0.0;

        for(j=1; j<=k; j++) {

          gbk += bern[k-j+1]*gbern[j];

        }

        gbern[k+1] = -rho*gbk/k;

        term  *= (2*k-2-x)*(2*k-1-x)*var2;

        poly1 += gbern[k+1]*term;

      }

    }

    stat_dexprl = gsl_sf_expm1_e(q, &dexprl);

    if(stat_dexprl != GSL_SUCCESS) {

      result->val = 0.0;

      result->err = 0.0;

      return stat_dexprl;

    }

    dexprl.val = dexprl.val/q;

    poly1 *= (x - 1.0);

    dpoch1 = dexprl.val * (alnvar + q * poly1) + poly1;

    for(i=incr-1; i >= 0; i--) {

      /*

       C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION

       C TO OBTAIN DPOCH1(BP,X).

       */

      double binv = 1.0/(bp+i);

      dpoch1 = (dpoch1 - binv) / (1.0 + x*binv);

    }

    if(bp == a) {

      result->val = dpoch1;

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

      return GSL_SUCCESS;

    }

    else {

      /*

       C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A

       C REFLECTION FORMULA TO OBTAIN DPOCH1(A,X).

       */

      double sinpxx = sin(M_PI*x)/x;

      double sinpx2 = sin(0.5*M_PI*x);

      double t1 = sinpxx/tan(M_PI*b);

      double t2 = 2.0*sinpx2*(sinpx2/x);

      double trig  = t1 - t2;

      result->val  = dpoch1 * (1.0 + x*trig) + trig;

      result->err  = (fabs(dpoch1*x) + 1.0) * GSL_DBL_EPSILON * (fabs(t1) + fabs(t2));

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

      return GSL_SUCCESS;

    }    

  }

}

/* Assumes a>0 and a+x>0.

 */

static

int

lnpoch_pos(const double a, const double x, gsl_sf_result * result)

{

  double absx = fabs(x);

  if(absx > 0.1*a || absx*log(GSL_MAX_DBL(a,2.0)) > 0.1) {

    if(a < GSL_SF_GAMMA_XMAX && a+x < GSL_SF_GAMMA_XMAX) {

      /* If we can do it by calculating the gamma functions

       * directly, then that will be more accurate than

       * doing the subtraction of the logs.

       */

      gsl_sf_result g1;

      gsl_sf_result g2;

      gsl_sf_gammainv_e(a,   &g1;;

      gsl_sf_gammainv_e(a+x, &g2;;

      result->val  = -log(g2.val/g1.val);

      result->err  = g1.err/fabs(g1.val) + g2.err/fabs(g2.val);

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

      return GSL_SUCCESS;

    }

    else {

      /* Otherwise we must do the subtraction.

       */

      gsl_sf_result lg1;

      gsl_sf_result lg2;

      int stat_1 = gsl_sf_lngamma_e(a,   &lg1;;

      int stat_2 = gsl_sf_lngamma_e(a+x, &lg2;;

      result->val  = lg2.val - lg1.val;

      result->err  = lg2.err + lg1.err;

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

      return GSL_ERROR_SELECT_2(stat_1, stat_2);

    }

  }

  else if(absx < 0.1*a && a > 15.0) {

    /* Be careful about the implied subtraction.

     * Note that both a+x and and a must be

     * large here since a is not small

     * and x is not relatively large.

     * So we calculate using Stirling for Log[Gamma(z)].

     *

     *   Log[Gamma(a+x)/Gamma(a)] = x(Log[a]-1) + (x+a-1/2)Log[1+x/a]

     *                              + (1/(1+eps)   - 1) / (12 a)

     *                              - (1/(1+eps)^3 - 1) / (360 a^3)

     *                              + (1/(1+eps)^5 - 1) / (1260 a^5)

     *                              - (1/(1+eps)^7 - 1) / (1680 a^7)

     *                              + ...

     */

    const double eps = x/a;

    const double den = 1.0 + eps;

    const double d3 = den*den*den;

    const double d5 = d3*den*den;

    const double d7 = d5*den*den;

    const double c1 = -eps/den;

    const double c3 = -eps*(3.0+eps*(3.0+eps))/d3;

    const double c5 = -eps*(5.0+eps*(10.0+eps*(10.0+eps*(5.0+eps))))/d5;

    const double c7 = -eps*(7.0+eps*(21.0+eps*(35.0+eps*(35.0+eps*(21.0+eps*(7.0+eps))))))/d7;

    const double p8 = gsl_sf_pow_int(1.0+eps,8);

    const double c8 = 1.0/p8             - 1.0;  /* these need not   */

    const double c9 = 1.0/(p8*(1.0+eps)) - 1.0;  /* be very accurate */

    const double a4 = a*a*a*a;

    const double a6 = a4*a*a;

    const double ser_1 = c1 + c3/(30.0*a*a) + c5/(105.0*a4) + c7/(140.0*a6);

    const double ser_2 = c8/(99.0*a6*a*a) - 691.0/360360.0 * c9/(a6*a4);

    const double ser = (ser_1 + ser_2)/ (12.0*a);

    double term1 = x * log(a/M_E);

    double term2;

    gsl_sf_result ln_1peps;

    gsl_sf_log_1plusx_e(eps, &ln_1peps);  /* log(1 + x/a) */

    term2 = (x + a - 0.5) * ln_1peps.val;

    result->val  = term1 + term2 + ser;

    result->err  = GSL_DBL_EPSILON*fabs(term1);

    result->err += fabs((x + a - 0.5)*ln_1peps.err);

    result->err += fabs(ln_1peps.val) * GSL_DBL_EPSILON * (fabs(x) + fabs(a) + 0.5);

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

    return GSL_SUCCESS;

  }

  else {

    gsl_sf_result poch_rel;

    int stat_p = pochrel_smallx(a, x, &poch_rel);

    double eps = x*poch_rel.val;

    int stat_e = gsl_sf_log_1plusx_e(eps, result);

    result->err  = 2.0 * fabs(x * poch_rel.err / (1.0 + eps));

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

    return GSL_ERROR_SELECT_2(stat_e, stat_p);

  }

}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int

gsl_sf_lnpoch_e(const double a, const double x, gsl_sf_result * result)

{

  /* CHECK_POINTER(result) */

  if(a <= 0.0 || a+x <= 0.0) {

    DOMAIN_ERROR(result);

  }

  else if(x == 0.0) {

    result->val = 0.0;

    result->err = 0.0;

    return GSL_SUCCESS;

  }

  else {

    return lnpoch_pos(a, x, result);

  }

}

int

gsl_sf_lnpoch_sgn_e(const double a, const double x,

                       gsl_sf_result * result, double * sgn)

{

  if(x == 0.0) {

    *sgn = 1.0;

    result->val = 0.0;

    result->err = 0.0;

    return GSL_SUCCESS;

  }

  else if(a > 0.0 && a+x > 0.0) {

    *sgn = 1.0;

    return lnpoch_pos(a, x, result);

  }

  else if (a <= 0 && a == floor(a)) {

    /* Special cases for infinite denominator Gamma(a) */

    if (a + x < 0  && x == floor(x)) {

      /* Handle the case where both a and a+x are negative integers. */

      gsl_sf_result result_pos;

      /* Use the reflection formula AMS6.1.17

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

      int stat = lnpoch_pos (-a, -x, &result_pos);

      double f = log (a / (a + x));

      double s = (fmod(x, 2) == 0) ? 1 : -1;

      result->val = f - result_pos.val;

      result->err = result_pos.err + 2.0 * GSL_DBL_EPSILON * f;

      *sgn = s;

      return stat;

    } else if (a + x == 0) {

      /* Handle a+x = 0 i.e. Gamma(0)/Gamma(a) */

      /* poch (-a,a) == (-1)^a Gamma(a+1) */

      int stat = gsl_sf_lngamma_sgn_e (-a + 1, result, sgn);

      double s = (fmod(-a, 2) == 0) ? 1 : -1;

      *sgn *= s;

      return stat;

    } else {

      /* Handle finite numerator, Gamma(a+x) for a+x != 0 or neg int */

      result->val = GSL_NEGINF;

      result->err = 0.0;

      *sgn = 1;

      return GSL_SUCCESS;

    }

  } else if(a < 0.0 && a+x < 0.0) {

    /* Reduce to positive case using reflection.

     */

    double sin_1 = sin(M_PI * (1.0 - a));

    double sin_2 = sin(M_PI * (1.0 - a - x));

    if(sin_1 == 0.0 || sin_2 == 0.0) {

      *sgn = 0.0;

      DOMAIN_ERROR(result);

    }

    else {

      gsl_sf_result lnp_pos;

      int stat_pp   = lnpoch_pos(1.0-a, -x, &lnp_pos);

      double lnterm = log(fabs(sin_1/sin_2));

      result->val  = lnterm - lnp_pos.val;

      result->err  = lnp_pos.err;

      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(1.0-a) + fabs(1.0-a-x)) * fabs(lnterm);

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

      *sgn = GSL_SIGN(sin_1*sin_2);

      return stat_pp;

    }

  }

  else {

    /* Evaluate gamma ratio directly.

     */

    gsl_sf_result lg_apn;

    gsl_sf_result lg_a;

    double s_apn, s_a;

    int stat_apn = gsl_sf_lngamma_sgn_e(a+x, &lg_apn, &s_apn);

    int stat_a   = gsl_sf_lngamma_sgn_e(a,   &lg_a,   &s_a);

    if(stat_apn == GSL_SUCCESS && stat_a == GSL_SUCCESS) {

      result->val  = lg_apn.val - lg_a.val;

      result->err  = lg_apn.err + lg_a.err;

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

      *sgn = s_a * s_apn;

      return GSL_SUCCESS;

    }

    else if(stat_apn == GSL_EDOM || stat_a == GSL_EDOM){

      *sgn = 0.0;

      DOMAIN_ERROR(result);

    }

    else {

      result->val = 0.0;

      result->err = 0.0;

      *sgn = 0.0;

      return GSL_FAILURE;

    }

  }

}

int

gsl_sf_poch_e(const double a, const double x, gsl_sf_result * result)

{

  /* CHECK_POINTER(result) */

  if(x == 0.0) {

    result->val = 1.0;

    result->err = 0.0;

    return GSL_SUCCESS;

  } else {

    gsl_sf_result lnpoch;

    double sgn;

    int stat_lnpoch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);

    if (lnpoch.val == GSL_NEGINF) {

      result->val = 0;

      result->err = 0;

      return stat_lnpoch;

    } else {

      int stat_exp    = gsl_sf_exp_err_e(lnpoch.val, lnpoch.err, result);

      result->val *= sgn;

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

      return GSL_ERROR_SELECT_2(stat_exp, stat_lnpoch);

    }

  }

}

int

gsl_sf_pochrel_e(const double a, const double x, gsl_sf_result * result)

{

  const double absx = fabs(x);

  const double absa = fabs(a);

  /* CHECK_POINTER(result) */

  if(absx > 0.1*absa || absx*log(GSL_MAX(absa,2.0)) > 0.1) {

    gsl_sf_result lnpoch;

    double sgn;

    int stat_poch = gsl_sf_lnpoch_sgn_e(a, x, &lnpoch, &sgn);

    if(lnpoch.val > GSL_LOG_DBL_MAX) {

      OVERFLOW_ERROR(result);

    }

    else {

      const double el = exp(lnpoch.val);

      result->val  = (sgn*el - 1.0)/x;

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

      result->err += 2.0 * GSL_DBL_EPSILON * (fabs(sgn*el) + 1.0) / fabs(x);

      return stat_poch;

    }

  }

  else {

    return pochrel_smallx(a, x, result);

  }

}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_lnpoch(const double a, const double x)

{

  EVAL_RESULT(gsl_sf_lnpoch_e(a, x, &result));

}

double gsl_sf_poch(const double a, const double x)

{

  EVAL_RESULT(gsl_sf_poch_e(a, x, &result));

}

double gsl_sf_pochrel(const double a, const double x)

{

  EVAL_RESULT(gsl_sf_pochrel_e(a, x, &result));

}