psi.c

math library from gnu

      result->err  = GSL_DBL_EPSILON * (fabs(t1) + fabs(x/(t2*t2)));
      result->err += result_c.err;
      result->err += GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else if(x < 1.0) { /* x = v */
      const double t1 = 1.0/x;
      cheb_eval_e(&psi_cs, 2.0*x-1.0, &result_c);
      
      result->val  = -t1 + result_c.val;
      result->err  = GSL_DBL_EPSILON * t1;
      result->err += result_c.err;
      result->err += GSL_DBL_EPSILON * fabs(result->val);
      return GSL_SUCCESS;
    }
    else { /* x = 1 + v */
      const double v = x - 1.0;
      return cheb_eval_e(&psi_cs, 2.0*v-1.0, result);
    }
  }
}


/* psi(z) for large |z| in the right half-plane; [Abramowitz + Stegun, 6.3.18] */
static
gsl_complex
psi_complex_asymp(gsl_complex z)
{
  /* coefficients in the asymptotic expansion for large z;
   * let w = z^(-2) and write the expression in the form
   *
   *   ln(z) - 1/(2z) - 1/12 w (1 + c1 w + c2 w + c3 w + ... )
   */
  static const double c1 = -0.1;
  static const double c2 =  1.0/21.0;
  static const double c3 = -0.05;

  gsl_complex zi = gsl_complex_inverse(z);
  gsl_complex w  = gsl_complex_mul(zi, zi);
  gsl_complex cs;

  /* Horner method evaluation of term in parentheses */
  gsl_complex sum;
  sum = gsl_complex_mul_real(w, c3/c2);
  sum = gsl_complex_add_real(sum, 1.0);
  sum = gsl_complex_mul_real(sum, c2/c1);
  sum = gsl_complex_mul(sum, w);
  sum = gsl_complex_add_real(sum, 1.0);
  sum = gsl_complex_mul_real(sum, c1);
  sum = gsl_complex_mul(sum, w);
  sum = gsl_complex_add_real(sum, 1.0);

  /* correction added to log(z) */
  cs = gsl_complex_mul(sum, w);
  cs = gsl_complex_mul_real(cs, -1.0/12.0);
  cs = gsl_complex_add(cs, gsl_complex_mul_real(zi, -0.5));

  return gsl_complex_add(gsl_complex_log(z), cs);
}



/* psi(z) for complex z in the right half-plane */
static int
psi_complex_rhp(
  gsl_complex z,
  gsl_sf_result * result_re,
  gsl_sf_result * result_im
  )
{
  int n_recurse = 0;
  int i;
  gsl_complex a;

  if(GSL_REAL(z) == 0.0 && GSL_IMAG(z) == 0.0)
  {
    result_re->val = 0.0;
    result_im->val = 0.0;
    result_re->err = 0.0;
    result_im->err = 0.0;
    return GSL_EDOM;
  }

  /* compute the number of recurrences to apply */
  if(GSL_REAL(z) < 20.0 && fabs(GSL_IMAG(z)) < 20.0)
  {
    const double sp = sqrt(20.0 + GSL_IMAG(z));
    const double sn = sqrt(20.0 - GSL_IMAG(z));
    const double rhs = sp*sn - GSL_REAL(z);
    if(rhs > 0.0) n_recurse = ceil(rhs);
  }

  /* compute asymptotic at the large value z + n_recurse */
  a = psi_complex_asymp(gsl_complex_add_real(z, n_recurse));

  result_re->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(a));
  result_im->err = 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(a));

  /* descend recursively, if necessary */
  for(i = n_recurse; i >= 1; --i)
  {
    gsl_complex zn = gsl_complex_add_real(z, i - 1.0);
    gsl_complex zn_inverse = gsl_complex_inverse(zn);
    a = gsl_complex_sub(a, zn_inverse);

    /* accumulate the error, to catch cancellations */
    result_re->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_REAL(zn_inverse));
    result_im->err += 2.0 * GSL_DBL_EPSILON * fabs(GSL_IMAG(zn_inverse));
  }

  result_re->val = GSL_REAL(a);
  result_im->val = GSL_IMAG(a);

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

  return GSL_SUCCESS;
}



/* generic polygamma; assumes n >= 0 and x > 0
 */
static int
psi_n_xg0(const int n, const double x, gsl_sf_result * result)
{
  if(n == 0) {
    return gsl_sf_psi_e(x, result);
  }
  else {
    /* Abramowitz + Stegun 6.4.10 */
    gsl_sf_result ln_nf;
    gsl_sf_result hzeta;
    int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
    int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
    int stat_e  = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
                                           hzeta.val, hzeta.err,
                                           result);
    if(GSL_IS_EVEN(n)) result->val = -result->val;
    return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
  }
}



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

int gsl_sf_psi_int_e(const int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n <= 0) {
    DOMAIN_ERROR(result);
  }
  else if(n <= PSI_TABLE_NMAX) {
    result->val = psi_table[n];
    result->err = GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    /* Abramowitz+Stegun 6.3.18 */
    const double c2 = -1.0/12.0;
    const double c3 =  1.0/120.0;
    const double c4 = -1.0/252.0;
    const double c5 =  1.0/240.0;
    const double ni2 = (1.0/n)*(1.0/n);
    const double ser = ni2 * (c2 + ni2 * (c3 + ni2 * (c4 + ni2*c5)));
    result->val  = log(n) - 0.5/n + ser;
    result->err  = GSL_DBL_EPSILON * (fabs(log(n)) + fabs(0.5/n) + fabs(ser));
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}


int gsl_sf_psi_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */
  return psi_x(x, result);
}


int
gsl_sf_psi_1piy_e(const double y, gsl_sf_result * result)
{
  const double ay = fabs(y);

  /* CHECK_POINTER(result) */

  if(ay > 1000.0) {
    /* [Abramowitz+Stegun, 6.3.19] */
    const double yi2 = 1.0/(ay*ay);
    const double lny = log(ay);
    const double sum = yi2 * (1.0/12.0 + 1.0/120.0 * yi2 + 1.0/252.0 * yi2*yi2);
    result->val = lny + sum;
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
    return GSL_SUCCESS;
  }
  else if(ay > 10.0) {
    /* [Abramowitz+Stegun, 6.3.19] */
    const double yi2 = 1.0/(ay*ay);
    const double lny = log(ay);
    const double sum = yi2 * (1.0/12.0 +
                         yi2 * (1.0/120.0 +
                           yi2 * (1.0/252.0 +
                             yi2 * (1.0/240.0 +
                               yi2 * (1.0/132.0 + 691.0/32760.0 * yi2)))));
    result->val = lny + sum;
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lny) + fabs(sum));
    return GSL_SUCCESS;
  }
  else if(ay > 1.0){
    const double y2 = ay*ay;
    const double x  = (2.0*ay - 11.0)/9.0;
    const double v  = y2*(1.0/(1.0+y2) + 0.5/(4.0+y2));
    gsl_sf_result result_c;
    cheb_eval_e(&r1py_cs, x, &result_c);
    result->val  = result_c.val - M_EULER + v;
    result->err  = result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(v) + M_EULER + fabs(result_c.val));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    result->err *= 5.0; /* FIXME: losing a digit somewhere... maybe at x=... ? */
    return GSL_SUCCESS;
  }
  else {
    /* [Abramowitz+Stegun, 6.3.17]
     *
     * -M_EULER + y^2 Sum[1/n 1/(n^2 + y^2), {n,1,M}]
     *   +     Sum[1/n^3, {n,M+1,Infinity}]
     *   - y^2 Sum[1/n^5, {n,M+1,Infinity}]
     *   + y^4 Sum[1/n^7, {n,M+1,Infinity}]
     *   - y^6 Sum[1/n^9, {n,M+1,Infinity}]
     *   + O(y^8)
     *
     * We take M=50 for at least 15 digit precision.
     */
    const int M = 50;
    const double y2 = y*y;
    const double c0 = 0.00019603999466879846570;
    const double c2 = 3.8426659205114376860e-08;
    const double c4 = 1.0041592839497643554e-11;
    const double c6 = 2.9516743763500191289e-15;
    const double p  = c0 + y2 *(-c2 + y2*(c4 - y2*c6));
    double sum = 0.0;
    double v;
    
    int n;
    for(n=1; n<=M; n++) {
      sum += 1.0/(n * (n*n + y*y));
    }

    v = y2 * (sum + p);
    result->val  = -M_EULER + v;
    result->err  = GSL_DBL_EPSILON * (M_EULER + fabs(v));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}


int gsl_sf_psi_1_int_e(const int n, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */
  if(n <= 0) {
    DOMAIN_ERROR(result);
  }
  else if(n <= PSI_1_TABLE_NMAX) {
    result->val = psi_1_table[n];
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    /* Abramowitz+Stegun 6.4.12
     * double-precision for n > 100
     */
    const double c0 = -1.0/30.0;
    const double c1 =  1.0/42.0;
    const double c2 = -1.0/30.0;
    const double ni2 = (1.0/n)*(1.0/n);
    const double ser =  ni2*ni2 * (c0 + ni2*(c1 + c2*ni2));
    result->val = (1.0 + 0.5/n + 1.0/(6.0*n*n) + ser) / n;
    result->err = GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
}


int gsl_sf_psi_1_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x == 0.0 || x == -1.0 || x == -2.0) {
    DOMAIN_ERROR(result);
  }
  else if(x > 0.0)
  {
    return psi_n_xg0(1, x, result);
  }
  else if(x > -5.0)
  {
    /* Abramowitz + Stegun 6.4.6 */
    int M = -floor(x);
    double fx = x + M;
    double sum = 0.0;
    int m;

    if(fx == 0.0)
      DOMAIN_ERROR(result);

    for(m = 0; m < M; ++m)
      sum += 1.0/((x+m)*(x+m));

    {
      int stat_psi = psi_n_xg0(1, fx, result);
      result->val += sum;
      result->err += M * GSL_DBL_EPSILON * sum;
      return stat_psi;
    }
  }
  else
  {
    /* Abramowitz + Stegun 6.4.7 */
    const double sin_px = sin(M_PI * x);
    const double d = M_PI*M_PI/(sin_px*sin_px);
    gsl_sf_result r;
    int stat_psi = psi_n_xg0(1, 1.0-x, &r);
    result->val = d - r.val;
    result->err = r.err + 2.0*GSL_DBL_EPSILON*d;
    return stat_psi;
  }
}


int gsl_sf_psi_n_e(const int n, const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(n == 0)
  {
    return gsl_sf_psi_e(x, result);
  }
  else if(n == 1)
  {
    return gsl_sf_psi_1_e(x, result);
  }
  else if(n < 0 || x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else {
    gsl_sf_result ln_nf;
    gsl_sf_result hzeta;
    int stat_hz = gsl_sf_hzeta_e(n+1.0, x, &hzeta);
    int stat_nf = gsl_sf_lnfact_e((unsigned int) n, &ln_nf);
    int stat_e  = gsl_sf_exp_mult_err_e(ln_nf.val, ln_nf.err,
                                           hzeta.val, hzeta.err,
                                           result);
    if(GSL_IS_EVEN(n)) result->val = -result->val;
    return GSL_ERROR_SELECT_3(stat_e, stat_nf, stat_hz);
  }
}


int
gsl_sf_complex_psi_e(
  const double x,
  const double y,
  gsl_sf_result * result_re,
  gsl_sf_result * result_im
  )
{
  if(x >= 0.0)
  {
    gsl_complex z = gsl_complex_rect(x, y);
    return psi_complex_rhp(z, result_re, result_im);
  }
  else
  {
    /* reflection formula [Abramowitz+Stegun, 6.3.7] */
    gsl_complex z = gsl_complex_rect(x, y);
    gsl_complex omz = gsl_complex_rect(1.0 - x, -y);
    gsl_complex zpi = gsl_complex_mul_real(z, M_PI);
    gsl_complex cotzpi = gsl_complex_cot(zpi);
    int ret_val = psi_complex_rhp(omz, result_re, result_im);

    if(GSL_IS_REAL(GSL_REAL(cotzpi)) && GSL_IS_REAL(GSL_IMAG(cotzpi)))
    {
      result_re->val -= M_PI * GSL_REAL(cotzpi);
      result_im->val -= M_PI * GSL_IMAG(cotzpi);
      return ret_val;
    }
    else
    {
      GSL_ERROR("singularity", GSL_EDOM);
    }
  }
}



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

#include "eval.h"

double gsl_sf_psi_int(const int n)
{
  EVAL_RESULT(gsl_sf_psi_int_e(n, &result));
}

double gsl_sf_psi(const double x)
{
  EVAL_RESULT(gsl_sf_psi_e(x, &result));
}

double gsl_sf_psi_1piy(const double x)
{
  EVAL_RESULT(gsl_sf_psi_1piy_e(x, &result));
}

double gsl_sf_psi_1_int(const int n)
{
  EVAL_RESULT(gsl_sf_psi_1_int_e(n, &result));
}

double gsl_sf_psi_1(const double x)
{
  EVAL_RESULT(gsl_sf_psi_1_e(x, &result));
}

double gsl_sf_psi_n(const int n, const double x)
{
  EVAL_RESULT(gsl_sf_psi_n_e(n, x, &result));
}