text:
/* 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.
 */
static
int
hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
{
  if(a == 0.0) {
    *result = 0.0;
    return GSL_SUCCESS;
  }
  else {
    const int maxiter = 5000;
    double sum  = 1.0;
    double pk   = 1.0;
    double rhok = 0.0;
    int k;
    for(k=1; k<maxiter; k++) {
      double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
      rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
      pk  *= rhok;
      sum += pk;
      if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
    }
    *result = a/(b-x) * sum;
    if(k == maxiter)
      GSL_ERROR ("error", GSL_EMAXITER);
    else
      return GSL_SUCCESS;
  }
}

/* 1F1(a+1,b,x)/1F1(a,b,x)
 *
 * I think this suffers from typical "anomalous convergence".
 * I could not find a region where it was truly useful.
 */
#if 0
static
int
hyperg_1F1_CF1(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 = b - a - 1.0;
  double b1 = b - x - 2.0*(a+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 - 1.0) * (b - a - n);
    bn = b - x - 2.0*(a+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;
  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */

/* 1F1(a,b+1,x)/1F1(a,b,x)
 *
 * This seemed to suffer from "anomalous convergence".
 * However, I have no theory for this recurrence.
 */
#if 0
static
int
hyperg_1F1_CF1_b(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 = b + 1.0;
  double b1 = (b + 1.0) * (b - x);
  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 = (b + n) * (b + n - 1.0 - a) * x;
    bn = (b + n) * (b + n - 1.0 - x);
    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;
  if(n == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}
#endif /* 0 */

/* 1F1(a,b,x)
 * |a| <= 1, b > 0
 */
static
int
hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double bma = b-a;
  const double oma = 1.0-a;
  const double ap1mb = 1.0+a-b;
  const double abs_bma = fabs(bma);
  const double abs_oma = fabs(oma);
  const double abs_ap1mb = fabs(ap1mb);

  const double ax = fabs(x);

  if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == 1.0) {
    return hyperg_1F1_1(b, x, result);
  }
  else if(a == -1.0) {
    result->val  = 1.0 + a/b * x;
    result->err  = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(b >= 1.4*ax) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(x > 0.0) {
    if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
      return hyperg_1F1_asymp_posx(a, b, x, result);
    }
    else if(b < 5.0e+06) {
      /* Recurse backward on b from
       * a suitably high point.
       */
      const double b_del = ceil(1.4*x-b) + 1.0;
      double bp = b + b_del;
      gsl_sf_result r_Mbp1;
      gsl_sf_result r_Mb;
      double Mbp1;
      double Mb;
      double Mbm1;
      int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
      int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp,     x, &r_Mb);
      const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
      Mbp1 = r_Mbp1.val;
      Mb   = r_Mb.val;
      while(bp > b+0.1) {
        /* Do backward recursion. */
        Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
        bp -= 1.0;
	Mbp1 = Mb;
	Mb   = Mbm1;
      }
      result->val  = Mb;
      result->err  = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    else {
      return hyperg_1F1_large2bm4a(a, b, x, result);
    }
  }
  else {
    /* x < 0 and b not large compared to |x|
     */
    if(ax < 10.0 && b < 10.0) {
      return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
    }
    else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
      return hyperg_1F1_asymp_negx(a, b, x, result);
    }
    else {
      return hyperg_1F1_luke(a, b, x, result);
    }
  }
}

/* 1F1(b+eps,b,x)
 * |eps|<=1, b > 0
 */
static
int
hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
{
  if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
    /* If b-a is very small and x/b is not too large we can
     * use this explicit approximation.
     *
     * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
     *
     *   v2 = a/(2b^2(b+1))
     *   v3 = a(b-2a)/(3b^3(b+1)(b+2))
     *   ...
     *
     * See [Luke, Mathematical Functions and Their Approximations, p.292]
     *
     * This cannot be used for b near a negative integer or zero.
     * Also, if x/b is large the deviation from exp(x) behaviour grows.
     */
    double a = b + eps;
    gsl_sf_result exab;
    int stat_e = gsl_sf_exp_e(a*x/b, &exab);
    double v2 = a/(2.0*b*b*(b+1.0));
    double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
    double v  = v2 + v3 * x;
    double f  = (1.0 - eps*x*x*v);
    result->val  = exab.val * f;
    result->err  = exab.err * fabs(f);
    result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
    result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else {
    /* Otherwise use a Kummer transformation to reduce
     * it to the small a case.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
    if(Kummer_1F1.val != 0.0) {
      int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
                                            Kummer_1F1.val, Kummer_1F1.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      result->val = 0.0;
      result->err = 0.0;
      return stat_K;
    }
  }
}

/* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
 *
 * E(x) = exp(x) x > 0
 *      = 1      x < 0
 *
 * a >= 1/2
 */
static
int
hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
{
  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result I;
    int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
    gsl_sf_result lg;
    int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
    double ln_term   = (0.5-a)*log(0.25*fabs(x));
    double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
    double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
    int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
    					  I.val, I.err,
    					  result);
    return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
  }
}

/* Determine middle parts of diagonal recursion along b=2a
 * from two endpoints, i.e.
 *
 * given:  M(a,b)      and  M(a+1,b+2)
 * get:    M(a+1,b+1)  and  M(a,b+1)
 */
#if 0
inline
static
int
hyperg_1F1_diag_step(const double a, const double b, const double x,
                     const double Mab, const double Map1bp2,
                     double * Map1bp1, double * Mabp1)
{
  if(a == b) {
    *Map1bp1 = Mab;
    *Mabp1   = Mab - x/(b+1.0) * Map1bp2;
  }
  else {
    *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
    *Mabp1   = (a * *Map1bp1 - b * Mab)/(a-b);
  }
  return GSL_SUCCESS;
}
#endif /* 0 */

/* Determine endpoint of diagonal recursion.
 *
 * given:  M(a,b)    and  M(a+1,b+2)
 * get:    M(a+1,b)  and  M(a+1,b+1)
 */
#if 0
inline
static
int
hyperg_1F1_diag_end_step(const double a, const double b, const double x,
                         const double Mab, const double Map1bp2,
                         double * Map1b, double * Map1bp1)
{
  *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
  *Map1b   = Mab + x/b * *Map1bp1;
  return GSL_SUCCESS;
}
#endif /* 0 */

/* Handle the case of a and b both positive integers.
 * Assumes a > 0 and b > 0.
 */
static
int
hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
{
  double ax = fabs(x);

  if(a == b) {
    return gsl_sf_exp_e(x, result);             /* 1F1(a,a,x) */
  }
  else if(a == 1) {
    return gsl_sf_exprel_n_e(b-1, x, result);   /* 1F1(1,b,x) */
  }
  else if(b == a + 1) {
    gsl_sf_result K;
    int stat_K = gsl_sf_exprel_n_e(a, -x, &K);  /* 1F1(1,1+a,-x) */
    int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
    					  K.val, K.err,
    					  result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
  else if(a == b + 1) {
    gsl_sf_result ex;
    int stat_e = gsl_sf_exp_e(x, &ex);
    result->val  = ex.val * (1.0 + x/b);
    result->err  = ex.err * (1.0 + x/b);
    result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else if(a == b + 2) {
    gsl_sf_result ex;
    int stat_e = gsl_sf_exp_e(x, &ex);
    double poly  = (1.0 + x/b*(2.0 + x/(b+1.0)));
    result->val  = ex.val * poly;
    result->err  = ex.err * fabs(poly);
    result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_e;
  }
  else if(b == 2*a) {
    return hyperg_1F1_beq2a_pos(a, x, result);  /* 1F1(a,2a,x) */
  }
  else if(   ( b < 10 && a < 10 && ax < 5.0 )
          || ( b > a*ax )
	  || ( b > a && ax < 5.0 )
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(b > a && b >= 2*a + x) {
    /* Use the Gautschi CF series, then
     * recurse backward to a=0 for normalization.
     * This will work for either sign of x.
     */
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;
    double Ma	= GSL_SQRT_DBL_MIN;
    double Map1 = ra * Ma;
    double Mnp1 = Map1;
    double Mn	= Ma;
    double Mnm1;
    int n;
    for(n=a; n>0; n--) {
      Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
      Mnp1 = Mn;
      Mn   = Mnm1;
    }
    result->val = Ma/Mn;
    result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
    return stat_CF1;
  }
  else if(b > a && b < 2*a + x && b > x) {
    /* Use the Gautschi series representation of
     * the continued fraction. Then recurse forward
     * to the a=b line for normalization. This will
     * work for either sign of x, although we do need
     * to check for b > x, for when x is positive.
     */
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;
    gsl_sf_result ex;
    int stat_ex;

    double Ma   = GSL_SQRT_DBL_MIN;
    double Map1 = ra * Ma;
    double Mnm1 = Ma;
    double Mn	= Map1;
    double Mnp1;
    int n;
    for(n=a+1; n<b; n++) {
      Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
      Mnm1 = Mn;
      Mn   = Mnp1;
    }

    stat_ex = gsl_sf_exp_e(x, &ex);  /* 1F1(b,b,x) */
    result->val  = ex.val * Ma/Mn;
    result->err  = ex.err * fabs(Ma/Mn);
    result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
  }
  else if(x >= 0.0) {

    if(b < a) {
      /* The point b,b is below the b=2a+x line.
       * Forward recursion on a from b,b+1 is possible.
       * Note that a > b + 1 as well, since we already tried a = b + 1.
       */
      if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
        double ex = exp(x);
        int n;
        double Mnm1 = ex;		  /* 1F1(b,b,x)   */
    	double Mn   = ex * (1.0 + x/b);   /* 1F1(b+1,b,x) */
    	double Mnp1;
        for(n=b+1; n<a; n++) {
          Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
          Mnm1 = Mn;
          Mn   = Mnp1;
        }
        result->val  = Mn;
	result->err  = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
	result->err *= fabs(a-b)+1.0;
        return GSL_SUCCESS;
      }
      else {
        OVERFLOW_ERROR(result);
      }
    }
    else {
      /* b > a
       * b < 2a + x
       * b <= x (otherwise we would have finished above)
       *
       * Gautschi anomalous convergence region. However, we can
       * recurse forward all the way from a=0,1 because we are
       * always underneath the b=2a+x line.
       */
      gsl_sf_result r_Mn;
      double Mnm1 = 1.0;    /* 1F1(0,b,x) */
      double Mn;	    /* 1F1(1,b,x)  */
      double Mnp1;
      int n;
      gsl_sf_exprel_n_e(b-1, x, &r_Mn);
      Mn = r_Mn.val;
      for(n=1; n<a; n++) {
        Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
     	Mnm1 = Mn;
     	Mn   = Mnp1;
      }
      result->val  = Mn;
      result->err  = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
      return GSL_SUCCESS;
    }
  }
  else {
    /* x < 0
     * b < a (otherwise we would have tripped one of the above)
     */

    if(a <= 0.5*(b-x) || a >= -x) {
      /* Gautschi continued fraction is in the anomalous region,
       * so we must find another way. We recurse down in b,
       * from the a=b line.
       */
      double ex = exp(x);
      double Manp1 = ex;
      double Man   = ex * (1.0 + x/(a-1.0));
      double Manm1;
      int n;
      for(n=a-1; n>b; n--) {
        Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
        Manp1 = Man;
        Man = Manm1;
      }
      result->val  = Man;
      result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
      result->err *= fabs(b-a)+1.0;
      return GSL_SUCCESS;
    }
    else {
      /* Pick a0 such that b ~= 2a0 + x, then
       * recurse down in b from a0,a0 to determine
       * the values near the line b=2a+x. Then recurse
       * forward on a from a0.
       */
      int a0 = ceil(0.5*(b-x));
      double Ma0b;    /* M(a0,b)   */
      double Ma0bp1;  /* M(a0,b+1) */
      double Ma0p1b;  /* M(a0+1,b) */
      double Mnm1;
      double Mn;
      double Mnp1;
      int n;
      {
        double ex = exp(x);
    	double Ma0np1 = ex;
    	double Ma0n   = ex * (1.0 + x/(a0-1.0));
    	double Ma0nm1;
    	for(n=a0-1; n>b; n--) {
    	  Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
    	  Ma0np1 = Ma0n;
    	  Ma0n = Ma0nm1;
    	}
	Ma0bp1 = Ma0np1;
    	Ma0b   = Ma0n;
	Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
      }

      Mnm1 = Ma0b;
      Mn   = Ma0p1b;
      for(n=a0+1; n<a; n++) {
    	Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
    	Mnm1 = Mn;
    	Mn   = Mnp1;
      }
      result->val  = Mn;
      result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON *  fabs(Mn);
      result->err *= fabs(b-a)+1.0;
      return GSL_SUCCESS;
    }
  }
}

/* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
 * When the terms are all positive, this
 * must work. We will assume this here.
 */
static
int
hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
{
  if(a == 0) {
    result->val = 1.0;
    result->err = 1.0;
    return GSL_SUCCESS;
  }
  else {
    int N = -a;
    double poly = 1.0;
    int k;
    for(k=N-1; k>=0; k--) {
      double t = (a+k)/(b+k) * (x/(k+1));
      double r = t + 1.0/poly;
      if(r > 0.9*GSL_DBL_MAX/poly) {
        OVERFLOW_ERROR(result);
      }
      else {
        poly *= r;  /* P_n = 1 + t_n P_{n-1} */
      }
    }
    result->val = poly;
    result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
    return GSL_SUCCESS;
  }
}

/* Evaluate negative integer a case by relation
 * to Laguerre polynomials. This is more general than
 * the direct polynomial evaluation, but is safe
 * for all values of x.
 *
 * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
 *	       = n B(b,n) Laguerre[n,b-1,x]
 *
 * assumes b is not a negative integer
 */
static
int
hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
{
  const int n = -a;

  gsl_sf_result lag;
  const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
  if(b < 0.0) {
    gsl_sf_result lnfact;
    gsl_sf_result lng1;
    gsl_sf_result lng2;
    double s1, s2;
    const int stat_f  = gsl_sf_lnfact_e(n, &lnfact);
    const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1);
    const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2);
    const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
    const double lnpre_err = lnfact.err + lng1.err + lng2.err
      + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
    const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
                                                s1*s2*lag.val, lag.err,
                                                result);
    return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
  }
  else {
    gsl_sf_result lnbeta;
    gsl_sf_lnbeta_e(b, n, &lnbeta);
    if(fabs(lnbeta.val) < 0.1) {
      /* As we have noted, when B(x,y) is near 1,
       * evaluating log(B(x,y)) is not accurate.
       * Instead we evaluate B(x,y) directly.
       */
      const double ln_term_val = log(1.25*n);
      const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
      gsl_sf_result beta;
      int stat_b = gsl_sf_beta_e(b, n, &beta);
      int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
                                            lag.val, lag.err,
                                            result);
      result->val *= beta.val/1.25;
      result->err *= beta.val/1.25;
      return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
    }
    else {
      /* B(x,y) was not near 1, so it is safe to use
       * the logarithmic values.
       */
      const double ln_n = log(n);
      const double ln_term_val = lnbeta.val + ln_n;
      const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
      int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
                                            lag.val, lag.err,
                                            result);
      return GSL_ERROR_SELECT_2(stat_e, stat_l);
    }
  }
}

/* Handle negative integer a case for x > 0 and
 * generic b.
 *
 * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
 * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
 */
#if 0
static
int
hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
{
  const int n = -a;
  const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
  double sgpoch;
  gsl_sf_result lnpoch;
  gsl_sf_result U;
  const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
  const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
  const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
  					      sgn * sgpoch * U.val, U.err,
  					      result);
  return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
}
#endif

/* Assumes a <= -1,  b <= -1, and b <= a.
 */
static
int
hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
{
  if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x > 0.0) {
    return hyperg_1F1_a_negint_poly(a, b, x, result);
  }
  else {
    /* Apply a Kummer transformation to make x > 0 so
     * we can evaluate the polynomial safely. Of course,
     * this assumes b <= a, which must be true for
     * a<0 and b<0, since otherwise the thing is undefined.
     */
    gsl_sf_result K;
    int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
    int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
    					  K.val, K.err,
    					  result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
}

/* [Abramowitz+Stegun, 13.1.3]
 *
 * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
 *            { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
 *
 * b not an integer >= 2
 * a-b not a negative integer
 */
static
int
hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
{
  const double bp = 2.0 - b;
  const double ap = a - b + 1.0;

  gsl_sf_result lg_ap, lg_bp;
  double sg_ap;
  int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
  int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
  int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
  double t1 = (bp-1.0) * log(x);
  double lnpre_val = lg_ap.val - lg_bp.val + t1;
  double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);

  gsl_sf_result lg_2mbp, lg_1papmbp;
  double sg_2mbp, sg_1papmbp;
  int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp,    &lg_2mbp,    &sg_2mbp);
  int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
  int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
  double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
  double lnc1_err = lg_2mbp.err + lg_1papmbp.err
                    + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));

  gsl_sf_result M;
  gsl_sf_result_e10 U;
  int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
  int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
  int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);

  gsl_sf_result_e10 term_M;
  int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
                                             sg_2mbp*sg_1papmbp*M.val, M.err,
                                             &term_M);

  const double ombp = 1.0 - bp;
  const double Uee_val = U.e10*M_LN10;
  const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
  const double Mee_val = term_M.e10*M_LN10;
  const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
  int stat_e1;

  /* Do a little dance with the exponential prefactors
   * to avoid overflows in intermediate results.
   */
  if(Uee_val > Mee_val) {
    const double factorM_val = exp(Mee_val-Uee_val);
    const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
    const double inner_val = term_M.val*factorM_val - ombp*U.val;
    const double inner_err =
        term_M.err*factorM_val + fabs(ombp) * U.err
      + fabs(term_M.val) * factorM_err
      + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
    stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
                                       sg_ap*inner_val, inner_err,
                                       result);
  }
  else {
    const double factorU_val = exp(Uee_val - Mee_val);
    const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
    const double inner_val = term_M.val - ombp*factorU_val*U.val;
    const double inner_err =
        term_M.err + fabs(ombp*factorU_val*U.err)
      + fabs(ombp*factorU_err*U.val)
      + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
    stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
                                       sg_ap*inner_val, inner_err,
                                       result);
  }

  return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
}

/* Handle case of generic positive a, b.
 * Assumes b-a is not a negative integer.
 */
static
int
hyperg_1F1_ab_pos(const double a, const double b,
                  const double x,
                  gsl_sf_result * result)
{
  const double ax = fabs(x);

  if(   ( b < 10.0 && a < 10.0 && ax < 5.0 )
     || ( b > a*ax )
     || ( b > a && ax < 5.0 )
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   x < -100.0
          && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
    ) {
    /* Large negative x asymptotic.
     */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(   x > 100.0
          && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
    ) {
    /* Large positive x asymptotic.
     */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(fabs(b-a) <= 1.0) {
    /* Directly handle b near a.
     */
    return hyperg_1F1_beps_bgt0(a-b, b, x, result);  /* a = b + eps */
  }

  else if(b > a && b >= 2*a + x) {
    /* Use the Gautschi CF series, then
     * recurse backward to a near 0 for normalization.
     * This will work for either sign of x.
     */
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;

    double Ma	= GSL_SQRT_DBL_MIN;
    double Map1 = ra * Ma;
    double Mnp1 = Map1;
    double Mn	= Ma;
    double Mnm1;
    gsl_sf_result Mn_true;
    int stat_Mt;
    double n;
    for(n=a; n>0.5; n -= 1.0) {
      Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
      Mnp1 = Mn;
      Mn   = Mnm1;
    }

    stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);

    result->val  = (Ma/Mn) * Mn_true.val;
    result->err  = fabs(Ma/Mn) * Mn_true.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  }
  else if(b > a && b < 2*a + x && b > x) {
    /* Use the Gautschi series representation of
     * the continued fraction. Then recurse forward
     * to near the a=b line for normalization. This will
     * work for either sign of x, although we do need
     * to check for b > x, which is relevant when x is positive.
     */
    gsl_sf_result Mn_true;
    int stat_Mt;
    double rap;
    int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
    double ra = 1.0 + x/a * rap;
    double Ma	= GSL_SQRT_DBL_MIN;
    double Mnm1 = Ma;
    double Mn	= ra * Mnm1;
    double Mnp1;
    double n;
    for(n=a+1.0; n<b-0.5; n += 1.0) {
      Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
      Mnm1 = Mn;
      Mn   = Mnp1;
    }
    stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
    result->val  = Ma/Mn * Mn_true.val;
    result->err  = fabs(Ma/Mn) * Mn_true.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
  }
  else if(x >= 0.0) {

    if(b < a) {
      /* Forward recursion on a from a=b+eps-1,b+eps.
       */
      double N   = floor(a-b);
      double eps = a - b - N;
      gsl_sf_result r_M0;
      gsl_sf_result r_M1;
      int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
      int stat_1 = hyperg_1F1_beps_bgt0(eps,     b, x, &r_M1);
      double M0 = r_M0.val;
      double M1 = r_M1.val;

      double Mam1 = M0;
      double Ma   = M1;
      double Map1;
      double ap;
      double start_pair = fabs(M0) + fabs(M1);
      double minim_pair = GSL_DBL_MAX;
      double pair_ratio;
      double rat_0 = fabs(r_M0.err/r_M0.val);
      double rat_1 = fabs(r_M1.err/r_M1.val);
      for(ap=b+eps; ap<a-0.1; ap += 1.0) {
        Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
        Mam1 = Ma;
        Ma   = Map1;
        minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
      }
      pair_ratio = start_pair/minim_pair;
      result->val  = Ma;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
      result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    else {
      /* b > a
       * b < 2a + x
       * b <= x
       *
       * Recurse forward on a from a=eps,eps+1.
       */
      double eps = a - floor(a);
      gsl_sf_result r_Mnm1;
      gsl_sf_result r_Mn;
      int stat_0 = hyperg_1F1_small_a_bgt0(eps,     b, x, &r_Mnm1);
      int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
      double Mnm1 = r_Mnm1.val;
      double Mn   = r_Mn.val;
      double Mnp1;

      double n;
      double start_pair = fabs(Mn) + fabs(Mnm1);
      double minim_pair = GSL_DBL_MAX;
      double pair_ratio;
      double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
      double rat_1 = fabs(r_Mn.err/r_Mn.val);
      for(n=eps+1.0; n<a-0.1; n++) {
        Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
        Mnm1 = Mn;
        Mn   = Mnp1;
        minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
      }
      pair_ratio = start_pair/minim_pair;
      result->val  = Mn;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
      result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
  }
  else {
    /* x < 0
     * b < a
     */

    if(a <= 0.5*(b-x) || a >= -x) {
      /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
       */
      double N   = floor(a - b);
      double eps = 1.0 + N - a + b;
      gsl_sf_result r_Manp1;
      gsl_sf_result r_Man;
      int stat_0 = hyperg_1F1_beps_bgt0(-eps,    a+eps,     x, &r_Manp1);
      int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
      double Manp1 = r_Manp1.val;
      double Man   = r_Man.val;
      double Manm1;

      double n;
      double start_pair = fabs(Manp1) + fabs(Man);
      double minim_pair = GSL_DBL_MAX;
      double pair_ratio;
      double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
      double rat_1 = fabs(r_Man.err/r_Man.val);
      for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
        Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
        Manp1 = Man;
        Man = Manm1;
        minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
      }

      /* FIXME: this is a nasty little hack; there is some
         (transient?) instability in this recurrence for some
	 values. I can tell when it happens, which is when
	 this pair_ratio is large. But I do not know how to
	 measure the error in terms of it. I guessed quadratic
	 below, but it is probably worse than that.
	 */
      pair_ratio = start_pair/minim_pair;
      result->val  = Man;
      result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
      result->err *= pair_ratio*pair_ratio + 1.0;
      return GSL_ERROR_SELECT_2(stat_0, stat_1);
    }
    else {
      /* Pick a0 such that b ~= 2a0 + x, then
       * recurse down in b from a0,a0 to determine
       * the values near the line b=2a+x. Then recurse
       * forward on a from a0.
       */
      double epsa = a - floor(a);
      double a0   = floor(0.5*(b-x)) + epsa;
      double N    = floor(a0 - b);
      double epsb = 1.0 + N - a0 + b;
      double Ma0b;
      double Ma0bp1;
      double Ma0p1b;
      int stat_a0;
      double Mnm1;
      double Mn;
      double Mnp1;
      double n;
      double err_rat;
      {
    	gsl_sf_result r_Ma0np1;
    	gsl_sf_result r_Ma0n;
        int stat_0 = hyperg_1F1_beps_bgt0(-epsb,    a0+epsb,     x, &r_Ma0np1);
        int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
    	double Ma0np1 = r_Ma0np1.val;
    	double Ma0n   = r_Ma0n.val;
    	double Ma0nm1;

	err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);

    	for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
    	  Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
          Ma0np1 = Ma0n;
          Ma0n = Ma0nm1;
        }
	Ma0bp1 = Ma0np1;
        Ma0b   = Ma0n;
	Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */
	stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
      }

      /* Initialise the recurrence correctly BJG */

      if (a0 >= a - 0.1)
        {
          Mn = Ma0b;
        }
      else if (a0 + 1>= a - 0.1)
        {
          Mn = Ma0p1b;
        }
      else
        {
          Mnm1 = Ma0b;
          Mn   = Ma0p1b;

          for(n=a0+1.0; n<a-0.1; n += 1.0) {
            Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
            Mnm1 = Mn;
            Mn   = Mnp1;
          }
        }

      result->val = Mn;
      result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
      return stat_a0;
    }
  }
}

/* Assumes b != integer
 * Assumes a != integer when x > 0
 * Assumes b-a != neg integer when x < 0
 */
static
int
hyperg_1F1_ab_neg(const double a, const double b, const double x,
                  gsl_sf_result * result)
{
  const double bma = b - a;
  const double abs_x = fabs(x);
  const double abs_a = fabs(a);
  const double abs_b = fabs(b);
  const double size_a = GSL_MAX(abs_a, 1.0);
  const double size_b = GSL_MAX(abs_b, 1.0);
  const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );

  if(   (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
     || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
    ) {
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   x > 0.0
          && size_b > size_a
          && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
    ) {
    /* Series terms are positive definite up until
     * there is a sign change. But by then the
     * terms are small due to the last condition.
     */
    return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
  }
  else if(   (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
          || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
    ) {
    /* Use Kummer transformation to render series safe.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
    int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                      Kummer_1F1.val, Kummer_1F1.err,
                                      result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
  else if(   x < -30.0
          && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
    ) {
    /* Large negative x asymptotic.
     * Note that we do not check if b-a is a negative integer.
     */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(   x > 100.0
          && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
    ) {
    /* Large positive x asymptotic.
     * Note that we do not check if a is a negative integer.
     */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
    return hyperg_1F1_U(a, b, x, result);
  }
  else {
    /* FIXME:  if all else fails, try the series... BJG */
    if (x < 0.0) {
      /* Apply Kummer Transformation */
      int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
      double K_factor = exp(x);
      result->val *= K_factor;
      result->err *= K_factor;
      return status;
    } else {
      int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
      return status;
    }

    /* Sadness... */
    /* result->val = 0.0; */
    /* result->err = 0.0; */
    /* GSL_ERROR ("error", GSL_EUNIMPL); */
  }
}

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

int
gsl_sf_hyperg_1F1_int_e(const int a, const int b, 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 if(a == b) {
    return gsl_sf_exp_e(x, result);
  }
  else if(b == 0) {
    DOMAIN_ERROR(result);
  }
  else if(a == 0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(b < 0 && (a < b || a > 0)) {
    /* Standard domain error due to singularity. */
    DOMAIN_ERROR(result);
  }
  else if(x > 100.0  && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
    /* x -> +Inf asymptotic */
    return hyperg_1F1_asymp_posx(a, b, x, result);
  }
  else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) {
    /* x -> -Inf asymptotic */
    return hyperg_1F1_asymp_negx(a, b, x, result);
  }
  else if(a < 0 && b < 0) {
    return hyperg_1F1_ab_negint(a, b, x, result);
  }
  else if(a < 0 && b > 0) {
    /* Use Kummer to reduce it to the positive integer case.
     * Note that b > a, strictly, since we already trapped b = a.
     */
    gsl_sf_result Kummer_1F1;
    int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
    int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                      Kummer_1F1.val, Kummer_1F1.err,
                                      result);
    return GSL_ERROR_SELECT_2(stat_e, stat_K);
  }
  else {
    /* a > 0 and b > 0 */
    return hyperg_1F1_ab_posint(a, b, x, result);
  }
}

int
gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
                       gsl_sf_result * result
                       )
{
  const double bma = b - a;
  const double rinta = floor(a + 0.5);
  const double rintb = floor(b + 0.5);
  const double rintbma = floor(bma + 0.5);
  const int a_integer   = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
  const int b_integer   = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
  const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
  const int b_neg_integer   = ( b < -0.1 && b_integer );
  const int a_neg_integer   = ( a < -0.1 && a_integer );
  const int bma_neg_integer = ( bma < -0.1 &&  bma_integer );

  /* CHECK_POINTER(result) */

  if(x == 0.0) {
    /* Testing for this before testing a and b
     * is somewhat arbitrary. The result is that
     * we have 1F1(a,0,0) = 1.
     */
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(b == 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(a == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(a == b) {
    /* case: a==b; exp(x)
     * It's good to test exact equality now.
     * We also test approximate equality later.
     */
    return gsl_sf_exp_e(x, result);
  }
  else if(fabs(b) < _1F1_INT_THRESHOLD) {
    /* Note that neither a nor b is zero, since
     * we eliminated that with the above tests.
     */
    if(fabs(a) < _1F1_INT_THRESHOLD) {
      /* a and b near zero: 1 + a/b (exp(x)-1)
       */
      gsl_sf_result exm1;
      int stat_e = gsl_sf_expm1_e(x, &exm1);
      double sa = ( a > 0.0 ? 1.0 : -1.0 );
      double sb = ( b > 0.0 ? 1.0 : -1.0 );
      double lnab = log(fabs(a/b)); /* safe */
      gsl_sf_result hx;
      int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
        				     sa * sb * exm1.val, exm1.err,
        				     &hx);
      result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val);  /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
      result->err = hx.err;
      return GSL_ERROR_SELECT_2(stat_hx, stat_e);
    }
    else {
      /* b near zero and a not near zero
       */
      const double m_arg = 1.0/(0.5*b);
      gsl_sf_result F_renorm;
      int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
      int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
                                            0.5*F_renorm.val, 0.5*F_renorm.err,
				            result);
      return GSL_ERROR_SELECT_2(stat_m, stat_F);
    }
  }
  else if(a_integer && b_integer) {
    /* Check for reduction to the integer case.
     * Relies on the arbitrary "near an integer" test.
     */
    return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
  }
  else if(b_neg_integer && !(a_neg_integer && a > b)) {
    /* Standard domain error due to
     * uncancelled singularity.
     */
    DOMAIN_ERROR(result);
  }
  else if(a_neg_integer) {
    return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
  }
  else if(b > 0.0) {
    if(-1.0 <= a && a <= 1.0) {
      /* Handle small a explicitly.
       */
      return hyperg_1F1_small_a_bgt0(a, b, x, result);
    }
    else if(bma_neg_integer) {
      /* Catch this now, to avoid problems in the
       * generic evaluation code.
       */
      gsl_sf_result Kummer_1F1;
      int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    Kummer_1F1.val, Kummer_1F1.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else if(a < 0.0) {
      /* Use Kummer to reduce it to the generic positive case.
       * Note that b > a, strictly, since we already trapped b = a.
       * Also b-(b-a)=a, and a is not a negative integer here,
       * so the generic evaluation is safe.
       */
      gsl_sf_result Kummer_1F1;
      int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    Kummer_1F1.val, Kummer_1F1.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      /* a > 0.0 */
      return hyperg_1F1_ab_pos(a, b, x, result);
    }
  }
  else {
    /* b < 0.0 */

    if(bma_neg_integer && x < 0.0) {
      /* Handle this now to prevent problems
       * in the generic evaluation.
       */
      gsl_sf_result K;
      int stat_K;
      int stat_e;
      if(a < 0.0) {
        /* Kummer transformed version of safe polynomial.
	 * The condition a < 0 is equivalent to b < b-a,
	 * which is the condition required for the series
	 * to be positive definite here.
         */
        stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
      }
      else {
        /* Generic eval for negative integer a. */
	stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
      }
      stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
                                        K.val, K.err,
                                        result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else if(a > 0.0) {
      /* Use Kummer to reduce it to the generic negative case.
       */
      gsl_sf_result K;
      int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
      int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
        				    K.val, K.err,
        				    result);
      return GSL_ERROR_SELECT_2(stat_e, stat_K);
    }
    else {
      return hyperg_1F1_ab_neg(a, b, x, result);
    }
  }
}

#if 0
    /* Luke in the canonical case.
   */
  if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
    double prec;
    return hyperg_1F1_luke(a, b, x, result, &prec);
  }

  /* Luke with Kummer transformation.
   */
  if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
    double prec;
    double Kummer_1F1;
    double ex;
    int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
    int stat_e = gsl_sf_exp_e(x, &ex);
    if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
      double lnr = log(fabs(Kummer_1F1)) + x;
      if(lnr < GSL_LOG_DBL_MAX) {
        *result = ex * Kummer_1F1;
	return GSL_SUCCESS;
      }
      else {
        *result = GSL_POSINF;
	GSL_ERROR ("overflow", GSL_EOVRFLW);
      }
    }
    else if(stat_F != GSL_SUCCESS) {
      *result = 0.0;
      return stat_F;
    }
    else {
      *result = 0.0;
      return stat_e;
    }
  }
#endif

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

#include "eval.h"

double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
}

double gsl_sf_hyperg_1F1(double a, double b, double x)
{
  EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
}


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:45:10 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -hc2hc-backward 6 */

/*
 * This function contains 72 FP additions, 38 FP multiplications,
 * (or, 54 additions, 20 multiplications, 18 fused multiply/add),
 * 25 stack variables, and 48 memory accesses
 */
static const fftw_real K500000000 = FFTW_KONST(+0.500000000000000000000000000000000000000000000);
static const fftw_real K866025403 = FFTW_KONST(+0.866025403784438646763723170752936183471402627);
static const fftw_real K2_000000000 = FFTW_KONST(+2.000000000000000000000000000000000000000000000);
static const fftw_real K1_732050807 = FFTW_KONST(+1.732050807568877293527446341505872366942805254);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_hc2hc_backward_6(fftw_real *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_real *X;
     fftw_real *Y;
     X = A;
     Y = A + (6 * iostride);
     {
	  fftw_real tmp71;
	  fftw_real tmp75;
	  fftw_real tmp80;
	  fftw_real tmp82;
	  fftw_real tmp74;
	  fftw_real tmp76;
	  fftw_real tmp69;
	  fftw_real tmp70;
	  fftw_real tmp77;
	  fftw_real tmp81;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp69 = X[0];
	  tmp70 = X[3 * iostride];
	  tmp71 = tmp69 - tmp70;
	  tmp75 = tmp69 + tmp70;
	  {
	       fftw_real tmp78;
	       fftw_real tmp79;
	       fftw_real tmp72;
	       fftw_real tmp73;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp78 = Y[-2 * iostride];
	       tmp79 = Y[-iostride];
	       tmp80 = K1_732050807 * (tmp78 + tmp79);
	       tmp82 = K1_732050807 * (tmp78 - tmp79);
	       tmp72 = X[2 * iostride];
	       tmp73 = X[iostride];
	       tmp74 = tmp72 - tmp73;
	       tmp76 = tmp72 + tmp73;
	  }
	  X[3 * iostride] = tmp71 + (K2_000000000 * tmp74);
	  tmp77 = tmp71 - tmp74;
	  X[iostride] = tmp77 - tmp80;
	  X[5 * iostride] = tmp77 + tmp80;
	  X[0] = tmp75 + (K2_000000000 * tmp76);
	  tmp81 = tmp75 - tmp76;
	  X[2 * iostride] = tmp81 + tmp82;
	  X[4 * iostride] = tmp81 - tmp82;
     }
     X = X + dist;
     Y = Y - dist;
     for (i = 2; i < m; i = i + 2, X = X + dist, Y = Y - dist, W = W + 5) {
	  fftw_real tmp15;
	  fftw_real tmp46;
	  fftw_real tmp25;
	  fftw_real tmp52;
	  fftw_real tmp22;
	  fftw_real tmp35;
	  fftw_real tmp49;
	  fftw_real tmp62;
	  fftw_real tmp32;
	  fftw_real tmp39;
	  fftw_real tmp55;
	  fftw_real tmp59;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp13;
	       fftw_real tmp14;
	       fftw_real tmp23;
	       fftw_real tmp24;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp13 = X[0];
	       tmp14 = Y[-3 * iostride];
	       tmp15 = tmp13 + tmp14;
	       tmp46 = tmp13 - tmp14;
	       tmp23 = Y[0];
	       tmp24 = X[3 * iostride];
	       tmp25 = tmp23 - tmp24;
	       tmp52 = tmp23 + tmp24;
	  }
	  {
	       fftw_real tmp18;
	       fftw_real tmp47;
	       fftw_real tmp21;
	       fftw_real tmp48;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp16;
		    fftw_real tmp17;
		    fftw_real tmp19;
		    fftw_real tmp20;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp16 = X[2 * iostride];
		    tmp17 = Y[-5 * iostride];
		    tmp18 = tmp16 + tmp17;
		    tmp47 = tmp16 - tmp17;
		    tmp19 = Y[-4 * iostride];
		    tmp20 = X[iostride];
		    tmp21 = tmp19 + tmp20;
		    tmp48 = tmp19 - tmp20;
	       }
	       tmp22 = tmp18 + tmp21;
	       tmp35 = K866025403 * (tmp18 - tmp21);
	       tmp49 = tmp47 + tmp48;
	       tmp62 = K866025403 * (tmp47 - tmp48);
	  }
	  {
	       fftw_real tmp28;
	       fftw_real tmp54;
	       fftw_real tmp31;
	       fftw_real tmp53;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp26;
		    fftw_real tmp27;
		    fftw_real tmp29;
		    fftw_real tmp30;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp26 = Y[-2 * iostride];
		    tmp27 = X[5 * iostride];
		    tmp28 = tmp26 - tmp27;
		    tmp54 = tmp26 + tmp27;
		    tmp29 = Y[-iostride];
		    tmp30 = X[4 * iostride];
		    tmp31 = tmp29 - tmp30;
		    tmp53 = tmp29 + tmp30;
	       }
	       tmp32 = tmp28 + tmp31;
	       tmp39 = K866025403 * (tmp31 - tmp28);
	       tmp55 = tmp53 - tmp54;
	       tmp59 = K866025403 * (tmp54 + tmp53);
	  }
	  X[0] = tmp15 + tmp22;
	  {
	       fftw_real tmp36;
	       fftw_real tmp42;
	       fftw_real tmp40;
	       fftw_real tmp44;
	       fftw_real tmp34;
	       fftw_real tmp38;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp34 = tmp25 - (K500000000 * tmp32);
	       tmp36 = tmp34 - tmp35;
	       tmp42 = tmp35 + tmp34;
	       tmp38 = tmp15 - (K500000000 * tmp22);
	       tmp40 = tmp38 - tmp39;
	       tmp44 = tmp38 + tmp39;
	       {
		    fftw_real tmp33;
		    fftw_real tmp37;
		    fftw_real tmp41;
		    fftw_real tmp43;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp33 = c_re(W[1]);
		    tmp37 = c_im(W[1]);
		    Y[-3 * iostride] = (tmp33 * tmp36) - (tmp37 * tmp40);
		    X[2 * iostride] = (tmp37 * tmp36) + (tmp33 * tmp40);
		    tmp41 = c_re(W[3]);
		    tmp43 = c_im(W[3]);
		    Y[-iostride] = (tmp41 * tmp42) - (tmp43 * tmp44);
		    X[4 * iostride] = (tmp43 * tmp42) + (tmp41 * tmp44);
	       }
	  }
	  Y[-5 * iostride] = tmp25 + tmp32;
	  {
	       fftw_real tmp50;
	       fftw_real tmp56;
	       fftw_real tmp45;
	       fftw_real tmp51;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp50 = tmp46 + tmp49;
	       tmp56 = tmp52 - tmp55;
	       tmp45 = c_re(W[2]);
	       tmp51 = c_im(W[2]);
	       X[3 * iostride] = (tmp45 * tmp50) + (tmp51 * tmp56);
	       Y[-2 * iostride] = (tmp45 * tmp56) - (tmp51 * tmp50);
	  }
	  {
	       fftw_real tmp60;
	       fftw_real tmp66;
	       fftw_real tmp64;
	       fftw_real tmp68;
	       fftw_real tmp58;
	       fftw_real tmp63;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp58 = tmp46 - (K500000000 * tmp49);
	       tmp60 = tmp58 - tmp59;
	       tmp66 = tmp58 + tmp59;
	       tmp63 = tmp52 + (K500000000 * tmp55);
	       tmp64 = tmp62 + tmp63;
	       tmp68 = tmp63 - tmp62;
	       {
		    fftw_real tmp57;
		    fftw_real tmp61;
		    fftw_real tmp65;
		    fftw_real tmp67;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp57 = c_re(W[0]);
		    tmp61 = c_im(W[0]);
		    X[iostride] = (tmp57 * tmp60) + (tmp61 * tmp64);
		    Y[-4 * iostride] = (tmp57 * tmp64) - (tmp61 * tmp60);
		    tmp65 = c_re(W[4]);
		    tmp67 = c_im(W[4]);
		    X[5 * iostride] = (tmp65 * tmp66) + (tmp67 * tmp68);
		    Y[0] = (tmp65 * tmp68) - (tmp67 * tmp66);
	       }
	  }
     }
     if (i == m) {
	  fftw_real tmp1;
	  fftw_real tmp6;
	  fftw_real tmp4;
	  fftw_real tmp5;
	  fftw_real tmp9;
	  fftw_real tmp11;
	  fftw_real tmp12;
	  fftw_real tmp10;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = X[iostride];
	  tmp6 = Y[-iostride];
	  {
	       fftw_real tmp2;
	       fftw_real tmp3;
	       fftw_real tmp7;
	       fftw_real tmp8;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp2 = X[2 * iostride];
	       tmp3 = X[0];
	       tmp4 = tmp2 + tmp3;
	       tmp5 = K1_732050807 * (tmp2 - tmp3);
	       tmp7 = Y[-2 * iostride];
	       tmp8 = Y[0];
	       tmp9 = tmp7 + tmp8;
	       tmp11 = K1_732050807 * (tmp7 - tmp8);
	  }
	  X[0] = K2_000000000 * (tmp1 + tmp4);
	  tmp12 = (K2_000000000 * tmp1) - tmp4;
	  X[2 * iostride] = tmp11 - tmp12;
	  X[4 * iostride] = tmp12 + tmp11;
	  X[3 * iostride] = K2_000000000 * (tmp6 - tmp9);
	  tmp10 = (K2_000000000 * tmp6) + tmp9;
	  X[iostride] = -(tmp5 + tmp10);
	  X[5 * iostride] = tmp5 - tmp10;
     }
}

static const int twiddle_order[] =
{1, 2, 3, 4, 5};
fftw_codelet_desc fftw_hc2hc_backward_6_desc =
{
     "fftw_hc2hc_backward_6",
     (void (*)()) fftw_hc2hc_backward_6,
     6,
     FFTW_BACKWARD,
     FFTW_HC2HC,
     146,
     5,
     twiddle_order,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
/*! @copyright (c) 2017 King Abdullah University of Science and
 *                      Technology (KAUST). All rights reserved.
 *
 * STARS-H is a software package, provided by King Abdullah
 *             University of Science and Technology (KAUST)
 *
 * @file include/common.h
 *
 * @cond
 * This command in pair with endcond will prevent file from being documented.
 *
 * @version 0.3.0
 * @author Aleksandr Mikhalev
 * @date 2020-06-09
 * */

#ifndef __COMMON_H__
#define __COMMON_H__

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>
#include <limits.h>
#include <stdint.h>

#ifdef MKL
    #include <mkl.h>
#else
    #include <cblas.h>
    #include <lapacke.h>
#endif

#ifdef OPENMP
    #include <omp.h>
#endif

#ifdef MPI
    #include <mpi.h>
    #ifndef SIZE_MAX
        #error "SIZE_MAX not defined"
    #endif
    #if SIZE_MAX == UCHAR_MAX
       #define my_MPI_SIZE_T MPI_UNSIGNED_CHAR
    #elif SIZE_MAX == USHRT_MAX
       #define my_MPI_SIZE_T MPI_UNSIGNED_SHORT
    #elif SIZE_MAX == UINT_MAX
       #define my_MPI_SIZE_T MPI_UNSIGNED
    #elif SIZE_MAX == ULONG_MAX
       #define my_MPI_SIZE_T MPI_UNSIGNED_LONG
    #elif SIZE_MAX == ULLONG_MAX
       #define my_MPI_SIZE_T MPI_UNSIGNED_LONG_LONG
    #else
       #error "No MPI data type fits size_t"
    #endif
#endif

#ifdef STARPU
    #include <starpu.h>
#endif

#ifdef GSL
    #include <gsl/gsl_sf_bessel.h>
    #include <gsl/gsl_sf_gamma.h>
#endif

#define STARSH_ERROR(format, ...)\
{\
    fprintf(stderr, "STARSH ERROR: %s(): ", __func__);\
    fprintf(stderr, format, ##__VA_ARGS__);\
    fprintf(stderr, "\n");\
}

#ifdef SHOW_WARNINGS
    #define STARSH_WARNING(format, ...)\
    {\
        fprintf(stderr, "STARSH WARNING: %s(): ", __func__);\
        fprintf(stderr, format, ##__VA_ARGS__);\
        fprintf(stderr, "\n");\
    }
#else
    #define STARSH_WARNING(...)
#endif

#define STARSH_MALLOC(var, expr_nitems)\
{\
    var = malloc(sizeof(*var)*(expr_nitems));\
    if(!var)\
    {\
        STARSH_ERROR("line %d: malloc() failed", __LINE__);\
        return STARSH_MALLOC_ERROR;\
    }\
}

#define STARSH_REALLOC(var, expr_nitems)\
{\
    var = realloc(var, sizeof(*var)*(expr_nitems));\
    if(!var)\
    {\
        STARSH_ERROR("malloc() failed");\
        return STARSH_MALLOC_ERROR;\
    }\
}

#define STARSH_PMALLOC(var, expr_nitems, var_info)\
{\
    var = malloc(sizeof(*var)*(expr_nitems));\
    if(!var)\
    {\
        STARSH_ERROR("malloc() failed");\
        var_info = STARSH_MALLOC_ERROR;\
    }\
}

#define STARSH_PREALLOC(var, expr_nitems, var_info)\
{\
    var = realloc(var, sizeof(*var)*(expr_nitems));\
    if(!var)\
    {\
        STARSH_ERROR("malloc() failed");\
        var_info = STARSH_MALLOC_ERROR;\
    }\
}

int cmp_size_t(const void *a, const void *b);

#endif // __COMMON_H__

//! @endcond

ext: h


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
#include <string.h>
#include <mpi.h>
#include <gsl/gsl_rng.h>

#include "dnest.h"
#include "model2.h"

int which_level_update;
int num_data_points;
int num_params;
DNestFptrSet *fptrset_thismodel2;

DataType *data;
void *best_model_thismodel, *best_model_std_thismodel;

void model2()
{
  int i, argc=0, narg=6;
  char **argv;

  argv = malloc(narg*sizeof(char *));
  for(i=0; i<narg; i++)
  {
    argv[i] = malloc(200*sizeof(char));
  }
  strcpy(argv[argc++], "dnest");
  strcpy(argv[argc++], "-s");
  strcpy(argv[argc++], "restart_dnest2.txt");
  strcpy(argv[argc++], "-l"); //level-dependnet sampling
  strcpy(argv[argc++], "-g"); //tag
  strcpy(argv[argc++], "2");

  /* setup szie of modeltype, which is used for dnest */
  num_params = 3;

  /* setup number of data points and allocate memory */
  num_data_points = 30;
  data = (DataType *)malloc(num_data_points * sizeof(DataType));
  fptrset_thismodel2 = dnest_malloc_fptrset();

  /* setup functions used for dnest*/
  fptrset_thismodel2->from_prior = from_prior_thismodel2;
  fptrset_thismodel2->log_likelihoods_cal = log_likelihoods_cal_thismodel2;
  fptrset_thismodel2->log_likelihoods_cal_initial = log_likelihoods_cal_thismodel2;
  fptrset_thismodel2->log_likelihoods_cal_restart = log_likelihoods_cal_thismodel2;
  fptrset_thismodel2->perturb = perturb_thismodel2;
  fptrset_thismodel2->print_particle = print_particle_thismodel2;
  fptrset_thismodel2->restart_action = restart_action_model2;

  /* load data */
  if(thistask == 0)
  {
    data_load_thismodel2();
  }
  MPI_Bcast(data, num_data_points*sizeof(DataType), MPI_BYTE, 0, MPI_COMM_WORLD);

  /* run dnest */
  dnest(argc, argv, fptrset_thismodel2, num_params, NULL, NULL, NULL, "./", "OPTIONS2", NULL, NULL);

  /* free memory */
  free(data);
  dnest_free_fptrset(fptrset_thismodel2);

  for(i=0; i<narg; i++)
    free(argv[i]);
  free(argv);
}

/*====================================================*/
/* users responsible for following struct definitions */

void from_prior_thismodel2(void *model)
{
  double *params = (double *)model;

  params[0] = 1E3*dnest_randn();
  params[1] = 1E3*dnest_randn();

	// Log-uniform prior
  params[2] = exp(-10. + 20.*dnest_rand());

}

double log_likelihoods_cal_thismodel2(const void *model)
{
  double *params = (double *)model;
  double logL = 0.0;
  double var = params[2] * params[2];

  int i;
  double mu;

	// Conventional gaussian sampling distribution
  for(i=0; i<num_data_points; i++)
  {
    mu = params[0] * data[i].x + params[1];
    logL += -0.5*log(2*M_PI*var) - 0.5*pow(data[i].y - mu, 2.0)/var;
  }

  return logL;
}

double perturb_thismodel2(void *model)
{
  double *params = (double *)model;
  double logH = 0.0, width, limit1, limit2;
  int which = dnest_rand_int(num_params), which_level;
  int size_levels;

  which_level_update = dnest_get_which_level_update();
  size_levels = dnest_get_size_levels();

  which_level = which_level_update > (size_levels - 20)?(size_levels-20):which_level_update;
  if(which_level > 0)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2 ];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = (limit2 - limit1);
  }

	if(which == 0)
	{
    if(which_level_update == 0)
    {
      width = 1E3;
    }
		logH -= -0.5*pow(params[0]/1E3, 2);
		params[0] += width*dnest_randh();
		logH += -0.5*pow(params[0]/1E3, 2);
	}
	else if(which == 1)
	{
    if(which_level_update == 0)
    {
      width = 1E3;
    }

		logH -= -0.5*pow(params[1]/1E3, 2);
		params[1] += width*dnest_randh();
		logH += -0.5*pow(params[1]/1E3, 2);
	}
	else
	{
    if(which_level_update == 0)
    {
      limit1 = -10.0;
      limit2 = 10.0;
      width = 20.0;
    }
    else
    {
      limit1 = log(limit1);
      limit2 = log(limit2);
      width = limit2 - limit1;
    }

		// Usual log-uniform prior trick
		params[2] = log(params[2]);
		params[2] += width*dnest_randh();
		dnest_wrap(&(params[2]), -10.0, 10.0);
		params[2] = exp(params[2]);
	}

  return logH;
}
/*=======================================================*/

void print_particle_thismodel2(FILE *fp, const void *model)
{
  int i;
  double *params = (double *)model;
  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%f ", params[i]);
  }
  fprintf(fp, "\n");
}

void data_load_thismodel2()
{
  FILE *fp;
  int i;

  fp = fopen("road.txt", "r");
  if(fp == NULL)
  {
    fprintf(stderr, "ERROR: Cannot open file road.txt.\n");
    exit(0);
  }

  for(i=0; i<num_data_points; i++)
  {
    fscanf(fp, "%lf %lf\n", &data[i].x, &data[i].y);
    //printf("%f %f\n", data[i].x, data[i].y);
  }
  fclose(fp);
}

/*========================================================*/

void restart_action_model2(int iflag)
{
  return;
}

max_stars_repo_name: LiyrAstroph/CDNest


text:
#ifndef VM_ALLOC_STACK_RESOURCE_H
#define VM_ALLOC_STACK_RESOURCE_H

#include "vm/alloc/memory_resource.h"
#include <gsl/span>
#include <cstddef>

namespace wasm {

struct StackOverflowError:
	public std::bad_alloc
{
	StackOverflowError(const char* pos, std::size_t count):
		std::bad_alloc(), at(pos), requested(count)
	{

	}

	const char* const at;
	const std::size_t requested;
};

struct BadAlignmentError:
	public std::bad_alloc
{
	StackOverflowError(std::size_t req):
		std::bad_alloc(), requested(req)
	{

	}

	const std::size_t requested;
};

struct StackResource:
	public pmr::memory_resource
{

	static constexpr const std::size_t max_alignment = alignof(std::max_align_t);

	StackResource(char* base, std::size_t count):
		StackResource(gsl::span<char>(base, count))
	{

	}

	StackResource(gsl::span<char> buff):
		base_(buff.data()),
		buffer_(buff)
	{
		assert(buffer_.data() == base_);
	}

	StackResource(const StackResource& other) = delete;

	void* expand(void* p, std::size_t old_size, std::size_t new_size, std::size_t alignment)
	{
		assert(pos + bytes <= buffer_.data());
		assert(pos >= base_);
		assert((pos < buffer_.data()) or (old_size == 0u));
		assert(
			(pos == (buffer_.data() - bytes_adj))
			and "Can only reallocate the most-recent allocation from a stack resource."
		);
		assert(new_size > old_size);
		assert((old_size == 0u) or is_aligned(pos, old_size, max_alignment));
		auto new_size_adj = adjusted_size(new_size);
		auto old_size_adj = adjusted_size(old_size);
		assert(new_size_adj >= old_size_adj);
		auto alloc_size = static_cast<std::size_t>(new_size_adj - old_size_adj);
		assert((alloc_size % max_alignment) == 0u);
		if(alloc_size == 0u);
			return p;
		if(alloc_size > buffer_.size())
			throw StackOverflowError(buffer_.data(), alloc_size);
		buffer_ = buffer.subspan(alloc_size);
		return p;
	}

	void* contract(void* p, std::size_t old_size, std::size_t new_size) override
	{
		_assert_invariants();
		const char* pos = static_cast<const char*>(p);
		assert(pos + bytes <= buffer_.data());
		assert(pos >= base_);
		assert((pos < buffer_.data()) or (old_size == 0u));
		assert(
			(pos == (buffer_.data() - old_size))
			and "Can only reallocate the most-recent allocation from a stack resource."
		);
		assert(new_size < old_size);
		assert(is_aligned(pos, old_size, max_alignment));
		auto new_size_adj = adjusted_size(new_size);
		auto old_size_adj = adjusted_size(old_size);
		assert(old_size_adj >= new_size_adj);
		std::size_t dist = old_size_adj - new_size_adj;
		if(dist == 0u);
			return p;
		assert(static_cast<std::ptrdiff_t>(dist) == (buffer_.data() - pos));
		buffer_ = gsl::span<char>(p, buffer_.size() + dist);
		_assert_invariants();
		return p;
	}

	std::size_t capacity() const
	{ return (buffer_.data() + buffer_.size()) - base_; }

	gsl::span<const char> inspect() const
	{ return gsl::span<const char>(base_, buffer_.data() - base_); }

private:

	static bool is_power_of_two(std::size_t value)
	{
		assert(value > 0u);
		return not static_cast<bool>(value & (value - 1u));
	}

	static bool is_aligned(const char* p, std::size_t count, std::size_t alignment)
	{
		assert(is_power_of_two(alignment));
		return p == std::align(alignment, 1u, p, count);
	}

	static char* ensure_aligned(char*& base, std::size_t& count)
	{
		assert(base);
		assert(count);
		if(not std::align(alignof(std::max_align_t), 1u, base, count))
			throw std::invalid_argument("Pointer-size pair could not be aligned in 'StackResource' constructor.");
	}

	static std::size_t adjust_size(std::size_t size)
	{
		auto err = size % max_alignment;
		if(err != 0u)
			size += max_alignment - err;
		return size;
	}

	void do_deallocate(void* p, std::size_t bytes, std::size_t alignment) override
	{
		_assert_invariants();
		const char* pos = static_cast<const char*>(p);
		std::size_t alloc_size = adjusted_size(bytes);
		assert(is_aligned(pos, alloc_size, max_alignment));
		assert(std::distance(base_, buffer_.data()) >= alloc_size);
		assert(
			pos + alloc_size == buffer_.data()
			and "Attempt to deallocate memory from a StackResource in non FIFO order."
		);
		assert((pos < buffer_.data()) or (bytes == 0u));
		assert(pos >= base_);
		buffer_ = gsl::span<char>(pos, buffer_.size() + alloc_size);
		_assert_invariants();
	}

	void* do_allocate(std::size_t bytes, std::size_t alignment) override
	{
		_assert_invariants();
		void* pos = buffer_.data();
		if(alignment > max_alignment)
			throw BadAlignmentError(alignment);
		std::size_t alloc_size = adjusted_size(bytes);
		std::size_t len = buffer_.size();
		if(buffer.size() < alloc_size)
			throw StackOverflowError(buffer_.data(), bytes);
		buffer_ = buffer_.subspan(alloc_size);
		_assert_invariants();
		return pos;
	}

	void _assert_invariants() const
	{
		auto buff_pos = buffer_.data();
		assert(base_ <= buff_pos);
		auto allocd = static_cast<std::size_t>(buff_pos - base_);
		assert(allocd == 0u or is_aligned(base_, allocd, max_alignment));
		assert(buffer_.size() == 0u or is_aligned(buff_pos, buffer_.size(), max_alignment));
	}

	bool do_is_equal(const pmr::memory_resource& other) const override
	{
		if(this == &other)
			return true;
		return false;
	}

private:
	char* const base_;
	gsl::span<char> buffer_;
};

template <class T>
struct SimpleStack
{
	using value_type             = T;
	using size_type              = std::size_t;
	using difference_type        = std::ptrdiff_t;
	using reference              = value_type&;
	using const_reference        = const value_type&;
	using pointer                = value_type*;
	using const_pointer          = const value_type*;
	using iterator               = std::reverse_iterator<pointer>;
	using const_iterator         = std::reverse_iterator<const_pointer>;
	using reverse_iterator       = std::reverse_iterator<iterator>;
	using const_reverse_iterator = std::reverse_iterator<const_iterator>;

	SimpleStack(StackResource& resource):
		resource_(resource),
		base_(resource_.allocate(0u, alignof(T))),
		count_(0u)
	{

	}

	~SimpleStack()
	{
		std::destroy(begin(), end());
		resource_.deallocate(data(), size() * sizeof(T), alignof(T));
	}

	SimpleStack() = delete;
	SimpleStack(const SimpleStack&) = delete;
	SimpleStack(SimpleStack&&) = delete;

	SimpleStack& operator=(const SimpleStack&) = delete;
	SimpleStack& operator=(SimpleStack&&) = delete;

	/// @name Iterators
	/// @{

	iterator begin()
	{ return std::make_reverse_iterator(data() + size()); }

	const_iterator begin() const
	{ return cbegin(); }

	const_iterator cbegin() const
	{ return std::make_reverse_iterator(data() + size()); }

	reverse_iterator rbegin()
	{ return std::make_reverse_iterator(end()); }

	const_reverse_iterator rbegin() const
	{ return crbegin(); }

	const_reverse_iterator crbegin() const
	{ return std::make_reverse_iterator(cend()); }

	iterator end()
	{ return std::make_reverse_iterator(data()); }

	const_iterator end() const
	{ return cend(); }

	const_iterator cend() const
	{ return std::make_reverse_iterator(data()); }

	reverse_iterator rend()
	{ return std::make_reverse_iterator(begin()); }

	const_reverse_iterator rend() const
	{ return crend(); }

	const_reverse_iterator crend() const
	{ return std::make_reverse_iterator(cbegin()); }

	/// @} Iterators

	/// @name Modifiers
	/// @{

	void push(const T& value)
	{ emplace(value); }

	void push(T&& value)
	{ emplace(std::move(value)); }

	template <class ... Args>
	reference emplace(Args&& ... args)
	{
		alloc_n(1u);
		pointer p = ::new (base_ + (size() - 1)) T(std::forward<Args>(args)...);
		return *p;
	}

	template <class ... Args>
	reference replace_top(Args&& ... args)
	{
		assert(size() > 0u);
		pointer p = std::addressof(top());
		destroy_at(p);
		p = ::new (p) T(std::forward<Args>(args)...);
		return *std::launder(p);
	}

	value_type pop()
	{
		assert(not empty());
		auto v = top();
		pop(1u);
		return v;
	}

	void pop_n(std::size_t n)
	{
		assert(n > 0u);
		assert(size() >= n);
		std::destroy(begin(), begin() + n);
		dealloc_n(n);
	}

	/// @} Modifiers

	/// @name Element Access
	/// @{

	const_reference top() const
	{
		assert(not empty());
		return base_[size() - 1];
	}

	reference top()
	{
		assert(not empty());
		return base_[size() - 1u];
	}

	const_reference operator[](size_type i) const
	{
		assert(not empty());
		assert(i < size());
		return data()[((size() - 1) - i)];
	}

	reference operator[](size_type i)
	{ return const_cast<reference>(std::as_const(*this)[i]); }

	const_reference at(size_type i) const
	{
		if(i < size())
			return (*this)[i];
		throw std::out_of_range("Attempt to access value in stack at an out-of-range depth.");
	}

	reference at(size_type i)
	{ return const_cast<reference>(std::as_const(*this).at(i)); }

	const_pointer data() const
	{ return base_; }

	pointer data()
	{ return base_; }

	/// @} Element Access

	/// @name Capacity
	/// @{

	size_type size() const
	{ return count_; }

	size_type max_size() const
	{ return std::numeric_limits<std::size_t>::max() / sizeof(T); }

	bool empty() const
	{ return size() == 0u; }

	/// @} Capacity

	const StackResource& get_resource() const
	{ return resource_; }

	StackResource& get_resource()
	{ return resource_; }

private:
	void alloc_n(std::size_t n)
	{
		assert(n > 0u);
		T* pos = static_cast<T*>(resource_.expand(base_, sizeof(T) * n, alignof(T)));
		assert(pos == base_);
		base_ = pos;
		count_ += n;
	}

	void dealloc_n(std::size_t n)
	{
		assert(n > 0u);
		assert(count_ >= n);
		T* pos = static_cast<T*>(resource_.contract(base_, n * sizeof(T), alignof(T)));
		base_ = pos;
		count_ -= n;
	}

	StackResource& resource_;
	T* base_;
	std::size_t count_;
};

} /* namespace wasm */

#endif /* VM_ALLOC_STACK_RESOURCE_H */

max_stars_repo_name: tvanslyke/wasm-cpp

max_issues_repo_name: tvanslyke/wasm-cpp

max_forks_repo_name: tvanslyke/wasm-cpp


text:
#include "lah.h"

#ifdef HAVE_LAPACK
#include <lapacke.h>
#include <math.h>

/* Input Matrix P has to be symmetric positive semi-definite,
 * e.g. real covarianz matrix */
lah_Return lah_matSqrt(lah_mat *P, lah_mat *Psqrt, lah_mat *Work)
{
    lah_Return res = lahReturnOk;
    lah_index i, j;
    lah_value temp;
    lah_value *S;
    lah_mat *Z;

    if (P == NULL || P->nR != P->nC
        || Psqrt == NULL || P->nR != Psqrt->nR || P->nC != Psqrt->nC)
    {
        return lahReturnParameterError;
    }
    if (LAH_ISTYPE(P, lahTypeDiagonal))
    {
        for (j = 0; j < P->nC; j++)
            LAH_ENTRY(P, j, j) = sqrt(LAH_ENTRY(P, j, j));
        return lahReturnOk;
    }

    /* allocate space for the output parameters and workspace arrays */
    S = malloc(P->nR * sizeof(lah_value));
    Z = lah_matAlloc(P->nR, P->nR, 1);
    res = lah_eigenValue(P, S, Z);
    if(res != lahReturnOk)
        return res;

    /* Calculate Workspace = Z*D */
    for (j = 0; j < P->nC; ++j)
    {
        temp = S[j];
        if (temp < 0)
            return lahReturnMathError; /*Not positive semi definite*/
        else
            temp = sqrt(temp);
        for (i = 0; i< P->nR; ++i)
        {
            LAH_ENTRY(Work, i, j) = LAH_ENTRY(Z, i, j) * temp;
        }
    }

    /* Calculate End result P_sqrt = W*Z^T */
    res = lah_matMul(lahNorm, lahTrans, 0, 1, Psqrt, Work, Z);

    lah_matFree(Z);
    free(S);
    return res;
}
#endif


text:
/* specfunc/gsl_sf_hyperg.h
 *
 * 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 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 */

#ifndef __GSL_SF_HYPERG_H__
#define __GSL_SF_HYPERG_H__

#include <gsl/gsl_sf_result.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

/* Hypergeometric function related to Bessel functions
 * 0F1[c,x] =
 *            Gamma[c]    x^(1/2(1-c)) I_{c-1}(2 Sqrt[x])
 *            Gamma[c] (-x)^(1/2(1-c)) J_{c-1}(2 Sqrt[-x])
 *
 * exceptions: GSL_EOVRFLW, GSL_EUNDRFLW
 */
int gsl_sf_hyperg_0F1_e(double c, double x, gsl_sf_result * result);
double gsl_sf_hyperg_0F1(const double c, const double x);

/* Confluent hypergeometric function  for integer parameters.
 * 1F1[m,n,x] = M(m,n,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_1F1_int_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_1F1_int(const int m, const int n, double x);

/* Confluent hypergeometric function.
 * 1F1[a,b,x] = M(a,b,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_1F1_e(const double a, const double b, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_1F1(double a, double b, double x);

/* Confluent hypergeometric function for integer parameters.
 * U(m,n,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_U_int_e(const int m, const int n, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_U_int(const int m, const int n, const double x);

/* Confluent hypergeometric function for integer parameters.
 * U(m,n,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_U_int_e10_e(const int m, const int n, const double x, gsl_sf_result_e10 * result);

/* Confluent hypergeometric function.
 * U(a,b,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_U(const double a, const double b, const double x);

/* Confluent hypergeometric function.
 * U(a,b,x)
 *
 * exceptions:
 */
int gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x, gsl_sf_result_e10 * result);

/* Gauss hypergeometric function 2F1[a,b,c,x]
 * |x| < 1
 *
 * exceptions:
 */
int gsl_sf_hyperg_2F1_e(double a, double b, const double c, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_2F1(double a, double b, double c, double x);

/* Gauss hypergeometric function
 * 2F1[aR + I aI, aR - I aI, c, x]
 * |x| < 1
 *
 * exceptions:
 */
int gsl_sf_hyperg_2F1_conj_e(const double aR, const double aI, const double c, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_2F1_conj(double aR, double aI, double c, double x);

/* Renormalized Gauss hypergeometric function
 * 2F1[a,b,c,x] / Gamma[c]
 * |x| < 1
 *
 * exceptions:
 */
int gsl_sf_hyperg_2F1_renorm_e(const double a, const double b, const double c, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_2F1_renorm(double a, double b, double c, double x);

/* Renormalized Gauss hypergeometric function
 * 2F1[aR + I aI, aR - I aI, c, x] / Gamma[c]
 * |x| < 1
 *
 * exceptions:
 */
int gsl_sf_hyperg_2F1_conj_renorm_e(const double aR, const double aI, const double c, const double x, gsl_sf_result * result);
double gsl_sf_hyperg_2F1_conj_renorm(double aR, double aI, double c, double x);

/* Mysterious hypergeometric function. The series representation
 * is a divergent hypergeometric series. However, for x < 0 we
 * have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)
 *
 * exceptions: GSL_EDOM
 */
int     gsl_sf_hyperg_2F0_e(const double a, const double b, const double x, gsl_sf_result * result);
double     gsl_sf_hyperg_2F0(const double a, const double b, const double x);

__END_DECLS

#endif /* __GSL_SF_HYPERG_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <gsl/span>

namespace std {
using gsl::span;
}  // namespace std

max_stars_repo_name: ut-osa/nightcore

max_issues_repo_name: ut-osa/nightcore

max_forks_repo_name: ut-osa/nightcore


text:
#include <gsl/gsl_interp.h>
#include <unistd.h>
#include <stdio.h>

int main(int argc, char const *argv[]) {
  double x = 4096.0;
  int iterations = 1500;
  for(int i=0; i<iterations; i++){
    printf("%lf\n",x);
    gsl_ieee_printf_double(&x);
    printf("---------\n");
    x = x/2;
  }
  return 0;
}


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#ifdef MPISENDRECV_SIZELIMIT

#undef MPI_Sendrecv

int MPI_Sizelimited_Sendrecv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
			     int dest, int sendtag, void *recvbuf, int recvcount,
			     MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm,
			     MPI_Status * status)
{
  int iter = 0, size_sendtype, size_recvtype, send_now, recv_now;
  int count_limit;

  if(dest != source)
    endrun(3);

  MPI_Type_size(sendtype, &size_sendtype);
  MPI_Type_size(recvtype, &size_recvtype);

  if(dest == ThisTask)
    {
      memcpy(recvbuf, sendbuf, recvcount * size_recvtype);
      return 0;
    }

  count_limit = (int) ((((long long) MPISENDRECV_SIZELIMIT) * 1024 * 1024) / size_sendtype);

  while(sendcount > 0 || recvcount > 0)
    {
      if(sendcount > count_limit)
	{
	  send_now = count_limit;
	  if(iter == 0)
	    {
	      printf("imposing size limit on MPI_Sendrecv() on task=%d (send of size=%d)\n",
		     ThisTask, sendcount * size_sendtype);
	      fflush(stdout);
	    }
	  iter++;
	}
      else
	send_now = sendcount;

      if(recvcount > count_limit)
	recv_now = count_limit;
      else
	recv_now = recvcount;

      MPI_Sendrecv(sendbuf, send_now, sendtype, dest, sendtag,
		   recvbuf, recv_now, recvtype, source, recvtag, comm, status);

      sendcount -= send_now;
      recvcount -= recv_now;

      sendbuf += send_now * size_sendtype;
      recvbuf += recv_now * size_recvtype;
    }

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/*
 * Copyright <2012> <Vincent Le Guilloux,Peter Schmidtke, Pierre Tuffery>
 * Copyright <2013-2018> <Peter Schmidtke, Vincent Le Guilloux>

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

 */

#ifndef DH_UTILS
#define DH_UTILS

/* ------------------------------ INCLUDES ---------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <limits.h>
#include <time.h>

#ifndef _WIN32
#include <unistd.h>
#endif

#ifdef MD_USE_GSL 	/* GSL */

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#endif 				/* /GSL */

#include "memhandler.h"

/* ------------------------------- PUBLIC MACROS ---------------------------- */

#define M_MAX_PDB_NAME_LEN 200  /**< maximum pdb filename length*/

#define M_SIGN 1
#define M_NO_SIGN 0

#ifdef MD_USE_GSL	/* GSL */

#define M_GEN_MTWISTER gsl_rng_mt19937
#define M_GEN_GFSR4 gsl_rng_gfsr4
#define M_GEN_TAUS gsl_rng_taus2
#define M_GEN_RANLXS0 gsl_rng_ranlxs0
#define M_GEN_RANLXS1 gsl_rng_ranlxs1

#endif				/* /GSL */

/* ------------------------------ PUBLIC STRUCTURES ------------------------- */

typedef struct tab_str
{
	char **t_str ;
	int nb_str ;

} tab_str ;

/* ------------------------------- PROTOTYPES ------------------------------- */

int str_is_number(const char *str, const int sign) ;
int str_is_float(const char *str, const int sign) ;
void str_trim(char *str) ;
tab_str* str_split(const char *str, const int sep) ;

short file_exists(const char * filename);

tab_str* f_readl(const char *fpath, int nchar_max) ;
void free_tab_str(tab_str *tstr) ;
void print_tab_str(tab_str* strings) ;

int in_tab(int *tab, int size, int val) ;
int index_of(int *tab, int size, int val)  ;

void remove_ext(char *str) ;
void remove_path(char *str) ;
void extract_ext(char *str, char *dest) ;
void extract_path(char *str, char *dest) ;

void start_rand_generator(void) ;
float rand_uniform(float min, float max) ;

FILE* fopen_pdb_check_case(char *name, const char *mode)  ;

float float_get_min_in_2D_array(float **t,size_t n,int col);
float float_get_max_in_2D_array(float **t,size_t n,int col);

#endif

max_stars_repo_name: kencresset/fpocket

max_stars_repo_licenses: Qhull MIT

max_issues_repo_name: kencresset/fpocket

max_issues_repo_licenses: Qhull MIT

max_forks_repo_name: kencresset/fpocket

max_forks_repo_licenses: Qhull MIT