htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/hyperg_1F1.c

   1 /* specfunc/hyperg_1F1.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_elementary.h>
  26 #include <gsl/gsl_sf_exp.h>
  27 #include <gsl/gsl_sf_bessel.h>
  28 #include <gsl/gsl_sf_gamma.h>
  29 #include <gsl/gsl_sf_laguerre.h>
  30 #include <gsl/gsl_sf_hyperg.h>
  31
  32 #include "error.h"
  33 #include "hyperg.h"
  34
  35 #define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON)
  36
  37
  38 /* Asymptotic result for 1F1(a, b, x)  x -> -Infinity.
  39  * Assumes b-a != neg integer and b != neg integer.
  40  */
  41 static
  42 int
  43 hyperg_1F1_asymp_negx(const double a, const double b, const double x,
  44                      gsl_sf_result * result)
  45 {
  46   gsl_sf_result lg_b;
  47   gsl_sf_result lg_bma;
  48   double sgn_b;
  49   double sgn_bma;
  50
  51   int stat_b   = gsl_sf_lngamma_sgn_e(b,   &lg_b,   &sgn_b);
  52   int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma);
  53
  54   if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) {
  55     gsl_sf_result F;
  56     int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F);
  57     if(F.val != 0) {
  58       double ln_term_val = a*log(-x);
  59       double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val));
  60       double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val;
  61       double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err;
  62       int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
  63                                             sgn_bma*sgn_b*F.val, F.err,
  64                                             result);
  65       return GSL_ERROR_SELECT_2(stat_e, stat_F);
  66     }
  67     else {
  68       result->val = 0.0;
  69       result->err = 0.0;
  70       return stat_F;
  71     }
  72   }
  73   else {
  74     DOMAIN_ERROR(result);
  75   }
  76 }
  77
  78
  79 /* Asymptotic result for 1F1(a, b, x)  x -> +Infinity
  80  * Assumes b != neg integer and a != neg integer
  81  */
  82 static
  83 int
  84 hyperg_1F1_asymp_posx(const double a, const double b, const double x,
  85                       gsl_sf_result * result)
  86 {
  87   gsl_sf_result lg_b;
  88   gsl_sf_result lg_a;
  89   double sgn_b;
  90   double sgn_a;
  91
  92   int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b);
  93   int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a);
  94
  95   if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) {
  96     gsl_sf_result F;
  97     int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F);
  98     if(stat_F == GSL_SUCCESS && F.val != 0) {
  99       double lnx = log(x);
 100       double ln_term_val = (a-b)*lnx;
 101       double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx)
 102                          + 2.0 * GSL_DBL_EPSILON * fabs(a-b);
 103       double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x;
 104       double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x);
 105       int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
 106                                             sgn_a*sgn_b*F.val, F.err,
 107                                             result);
 108       return GSL_ERROR_SELECT_2(stat_e, stat_F);
 109     }
 110     else {
 111       result->val = 0.0;
 112       result->err = 0.0;
 113       return stat_F;
 114     }
 115   }
 116   else {
 117     DOMAIN_ERROR(result);
 118   }
 119 }
 120
 121 /* Asymptotic result from Slater 4.3.7
 122  *
 123  * To get the general series, write M(a,b,x) as
 124  *
 125  *  M(a,b,x)=sum ((a)_n/(b)_n) (x^n / n!)
 126  *
 127  * and expand (b)_n in inverse powers of b as follows
 128  *
 129  * -log(1/(b)_n) = sum_(k=0)^(n-1) log(b+k)
 130  *             = n log(b) + sum_(k=0)^(n-1) log(1+k/b)
 131  *
 132  * Do a taylor expansion of the log in 1/b and sum the resulting terms
 133  * using the standard algebraic formulas for finite sums of powers of
 134  * k.  This should then give
 135  *
 136  * M(a,b,x) = sum_(n=0)^(inf) (a_n/n!) (x/b)^n * (1 - n(n-1)/(2b)
 137  *                          + (n-1)n(n+1)(3n-2)/(24b^2) + ...
 138  *
 139  * which can be summed explicitly. The trick for summing it is to take
 140  * derivatives of sum_(i=0)^(inf) a_n*y^n/n! = (1-y)^(-a);
 141  *
 142  * [BJG 16/01/2007]
 143  */
 144
 145 static
 146 int
 147 hyperg_1F1_largebx(const double a, const double b, const double x, gsl_sf_result * result)
 148 {
 149   double y = x/b;
 150   double f = exp(-a*log1p(-y));
 151   double t1 = -((a*(a+1.0))/(2*b))*pow((y/(1.0-y)),2.0);
 152   double t2 = (1/(24*b*b))*((a*(a+1)*y*y)/pow(1-y,4))*(12+8*(2*a+1)*y+(3*a*a-a-2)*y*y);
 153   double t3 = (-1/(48*b*b*b*pow(1-y,6)))*a*((a + 1)*((y*((a + 1)*(a*(y*(y*((y*(a - 2) + 16)*(a - 1)) + 72)) + 96)) + 24)*pow(y, 2)));
 154   result->val = f * (1 + t1 + t2 + t3);
 155   result->err = 2*fabs(f*t3) + 2*GSL_DBL_EPSILON*fabs(result->val);
 156   return GSL_SUCCESS;
 157 }
 158
 159 /* Asymptotic result for x < 2b-4a, 2b-4a large.
 160  * [Abramowitz+Stegun, 13.5.21]
 161  *
 162  * assumes 0 <= x/(2b-4a) <= 1
 163  */
 164 static
 165 int
 166 hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result)
 167 {
 168   double eta    = 2.0*b - 4.0*a;
 169   double cos2th = x/eta;
 170   double sin2th = 1.0 - cos2th;
 171   double th = acos(sqrt(cos2th));
 172   double pre_h  = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th;
 173   gsl_sf_result lg_b;
 174   int stat_lg = gsl_sf_lngamma_e(b, &lg_b);
 175   double t1 = 0.5*(1.0-b)*log(0.25*x*eta);
 176   double t2 = 0.25*log(pre_h);
 177   double lnpre_val = lg_b.val + 0.5*x + t1 - t2;
 178   double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2));
 179 #if SMALL_ANGLE
 180   const double eps = asin(sqrt(cos2th));  /* theta = pi/2 - eps */
 181   double s1 = (fmod(a, 1.0) == 0.0) ? 0.0 : sin(a*M_PI);
 182   double eta_reduc = (fmod(eta + 1, 4.0) == 0.0) ? 0.0 : fmod(eta + 1, 8.0);
 183   double phi1 = 0.25*eta_reduc*M_PI;
 184   double phi2 = 0.25*eta*(2*eps + sin(2.0*eps));
 185   double s2 = sin(phi1 - phi2);
 186 #else
 187   double s1 = sin(a*M_PI);
 188   double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI);
 189 #endif
 190   double ser_val = s1 + s2;
 191   double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2));
 192   int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
 193                                         ser_val, ser_err,
 194                                         result);
 195   return GSL_ERROR_SELECT_2(stat_e, stat_lg);
 196 }
 197
 198
 199 /* Luke's rational approximation.
 200  * See [Luke, Algorithms for the Computation of Mathematical Functions, p.182]
 201  *
 202  * Like the case of the 2F1 rational approximations, these are
 203  * probably guaranteed to converge for x < 0, barring gross
 204  * numerical instability in the pre-asymptotic regime.
 205  */
 206 static
 207 int
 208 hyperg_1F1_luke(const double a, const double c, const double xin,
 209                 gsl_sf_result * result)
 210 {
 211   const double RECUR_BIG = 1.0e+50;
 212   const int nmax = 5000;
 213   int n = 3;
 214   const double x  = -xin;
 215   const double x3 = x*x*x;
 216   const double t0 = a/c;
 217   const double t1 = (a+1.0)/(2.0*c);
 218   const double t2 = (a+2.0)/(2.0*(c+1.0));
 219   double F = 1.0;
 220   double prec;
 221
 222   double Bnm3 = 1.0;                                  /* B0 */
 223   double Bnm2 = 1.0 + t1 * x;                         /* B1 */
 224   double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x);    /* B2 */
 225
 226   double Anm3 = 1.0;                                                      /* A0 */
 227   double Anm2 = Bnm2 - t0 * x;                                            /* A1 */
 228   double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x;   /* A2 */
 229
 230   while(1) {
 231     double npam1 = n + a - 1;
 232     double npcm1 = n + c - 1;
 233     double npam2 = n + a - 2;
 234     double npcm2 = n + c - 2;
 235     double tnm1  = 2*n - 1;
 236     double tnm3  = 2*n - 3;
 237     double tnm5  = 2*n - 5;
 238     double F1 =  (n-a-2) / (2*tnm3*npcm1);
 239     double F2 =  (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1);
 240     double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
 241     double E  = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1);
 242
 243     double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
 244     double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
 245     double r = An/Bn;
 246
 247     prec = fabs((F - r)/F);
 248     F = r;
 249
 250     if(prec < GSL_DBL_EPSILON || n > nmax) break;
 251
 252     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 253       An   /= RECUR_BIG;
 254       Bn   /= RECUR_BIG;
 255       Anm1 /= RECUR_BIG;
 256       Bnm1 /= RECUR_BIG;
 257       Anm2 /= RECUR_BIG;
 258       Bnm2 /= RECUR_BIG;
 259       Anm3 /= RECUR_BIG;
 260       Bnm3 /= RECUR_BIG;
 261     }
 262     else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
 263       An   *= RECUR_BIG;
 264       Bn   *= RECUR_BIG;
 265       Anm1 *= RECUR_BIG;
 266       Bnm1 *= RECUR_BIG;
 267       Anm2 *= RECUR_BIG;
 268       Bnm2 *= RECUR_BIG;
 269       Anm3 *= RECUR_BIG;
 270       Bnm3 *= RECUR_BIG;
 271     }
 272
 273     n++;
 274     Bnm3 = Bnm2;
 275     Bnm2 = Bnm1;
 276     Bnm1 = Bn;
 277     Anm3 = Anm2;
 278     Anm2 = Anm1;
 279     Anm1 = An;
 280   }
 281
 282   result->val  = F;
 283   result->err  = 2.0 * fabs(F * prec);
 284   result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F);
 285
 286   return GSL_SUCCESS;
 287 }
 288
 289 /* Series for 1F1(1,b,x)
 290  * b > 0
 291  */
 292 static
 293 int
 294 hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result)
 295 {
 296   double sum_val = 1.0;
 297   double sum_err = 0.0;
 298   double term = 1.0;
 299   double n    = 1.0;
 300   while(fabs(term/sum_val) > 0.25*GSL_DBL_EPSILON) {
 301     term *= x/(b+n-1);
 302     sum_val += term;
 303     sum_err += 8.0*GSL_DBL_EPSILON*fabs(term) + GSL_DBL_EPSILON*fabs(sum_val);
 304     n += 1.0;
 305   }
 306   result->val  = sum_val;
 307   result->err  = sum_err;
 308   result->err += 2.0 *  fabs(term);
 309   return GSL_SUCCESS;
 310 }
 311
 312
 313 /* 1F1(1,b,x)
 314  * b >= 1, b integer
 315  */
 316 static
 317 int
 318 hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result)
 319 {
 320   if(b < 1) {
 321     DOMAIN_ERROR(result);
 322   }
 323   else if(b == 1) {
 324     return gsl_sf_exp_e(x, result);
 325   }
 326   else if(b == 2) {
 327     return gsl_sf_exprel_e(x, result);
 328   }
 329   else if(b == 3) {
 330     return gsl_sf_exprel_2_e(x, result);
 331   }
 332   else {
 333     return gsl_sf_exprel_n_e(b-1, x, result);
 334   }
 335 }
 336
 337
 338 /* 1F1(1,b,x)
 339  * b >=1, b real
 340  *
 341  * checked OK: [GJ] Thu Oct  1 16:46:35 MDT 1998
 342  */
 343 static
 344 int
 345 hyperg_1F1_1(const double b, const double x, gsl_sf_result * result)
 346 {
 347   double ax = fabs(x);
 348   double ib = floor(b + 0.1);
 349
 350   if(b < 1.0) {
 351     DOMAIN_ERROR(result);
 352   }
 353   else if(b == 1.0) {
 354     return gsl_sf_exp_e(x, result);
 355   }
 356   else if(b >= 1.4*ax) {
 357     return hyperg_1F1_1_series(b, x, result);
 358   }
 359   else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) {
 360     return hyperg_1F1_1_int((int)ib, x, result);
 361   }
 362   else if(x > 0.0) {
 363     if(x > 100.0 && b < 0.75*x) {
 364       return hyperg_1F1_asymp_posx(1.0, b, x, result);
 365     }
 366     else if(b < 1.0e+05) {
 367       /* Recurse backward on b, from a
 368        * chosen offset point. For x > 0,
 369        * which holds here, this should
 370        * be a stable direction.
 371        */
 372       const double off = ceil(1.4*x-b) + 1.0;
 373       double bp = b + off;
 374       gsl_sf_result M;
 375       int stat_s = hyperg_1F1_1_series(bp, x, &M);
 376       const double err_rat = M.err / fabs(M.val);
 377       while(bp > b+0.1) {
 378         /* M(1,b-1) = x/(b-1) M(1,b) + 1 */
 379         bp -= 1.0;
 380         M.val  = 1.0 + x/bp * M.val;
 381       }
 382       result->val  = M.val;
 383       result->err  = err_rat * fabs(M.val);
 384       result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val);
 385       return stat_s;
 386     } else if (fabs(x) < fabs(b) && fabs(x) < sqrt(fabs(b)) * fabs(b-x)) {
 387       return hyperg_1F1_largebx(1.0, b, x, result);
 388     } else if (fabs(x) > fabs(b)) {
 389       return hyperg_1F1_1_series(b, x, result);
 390     } else {
 391       return hyperg_1F1_large2bm4a(1.0, b, x, result);
 392     }
 393   }
 394   else {
 395     /* x <= 0 and b not large compared to |x|
 396      */
 397     if(ax < 10.0 && b < 10.0) {
 398       return hyperg_1F1_1_series(b, x, result);
 399     }
 400     else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) {
 401       return hyperg_1F1_asymp_negx(1.0, b, x, result);
 402     }
 403     else {
 404       return hyperg_1F1_luke(1.0, b, x, result);
 405     }
 406   }
 407 }
 408
 409
 410 /* 1F1(a,b,x)/Gamma(b) for b->0
 411  * [limit of Abramowitz+Stegun 13.3.7]
 412  */
 413 static
 414 int
 415 hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result)
 416 {
 417   double eta = a*x;
 418   if(eta > 0.0) {
 419     double root_eta = sqrt(eta);
 420     gsl_sf_result I1_scaled;
 421     int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled);
 422     if(I1_scaled.val <= 0.0) {
 423       result->val = 0.0;
 424       result->err = 0.0;
 425       return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM);
 426     }
 427     else {
 428       /* Note that 13.3.7 contains higher terms which are zeroth order
 429          in b.  These make a non-negligible contribution to the sum.
 430          With the first correction term, the I1 above is replaced by
 431          I1 + (2/3)*a*(x/(4a))**(3/2)*I2(2*root_eta).  We will add
 432          this as part of the result and error estimate. */
 433
 434       const double corr1 =(2.0/3.0)*a*pow(x/(4.0*a),1.5)*gsl_sf_bessel_In_scaled(2, 2.0*root_eta)
 435  ;
 436       const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(2.0*root_eta) + log(I1_scaled.val+corr1);
 437       const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs((I1_scaled.err+corr1)/I1_scaled.val);
 438       return gsl_sf_exp_err_e(lnr_val, lnr_err, result);
 439     }
 440   }
 441   else if(eta == 0.0) {
 442     result->val = 0.0;
 443     result->err = 0.0;
 444     return GSL_SUCCESS;
 445   }
 446   else {
 447     /* eta < 0 */
 448     double root_eta = sqrt(-eta);
 449     gsl_sf_result J1;
 450     int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1);
 451     if(J1.val <= 0.0) {
 452       result->val = 0.0;
 453       result->err = 0.0;
 454       return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM);
 455     }
 456     else {
 457       const double t1 = 0.5*x;
 458       const double t2 = 0.5*log(-eta);
 459       const double t3 = fabs(x);
 460       const double t4 = log(J1.val);
 461       const double lnr_val = t1 + t2 + t3 + t4;
 462       const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val);
 463       gsl_sf_result ex;
 464       int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex);
 465       result->val = -ex.val;
 466       result->err =  ex.err;
 467       return stat_e;
 468     }
 469   }
 470
 471 }
 472
 473
 474 /* 1F1'(a,b,x)/1F1(a,b,x)
 475  * Uses Gautschi's version of the CF.
 476  * [Gautschi, Math. Comp. 31, 994 (1977)]
 477  *
 478  * Supposedly this suffers from the "anomalous convergence"
 479  * problem when b < x. I have seen anomalous convergence
 480  * in several of the continued fractions associated with
 481  * 1F1(a,b,x). This particular CF formulation seems stable
 482  * for b > x. However, it does display a painful artifact
 483  * of the anomalous convergence; the convergence plateaus
 484  * unless b >>> x. For example, even for b=1000, x=1, this
 485  * method locks onto a ratio which is only good to about
 486  * 4 digits. Apparently the rest of the digits are hiding
 487  * way out on the plateau, but finite-precision lossage
 488  * means you will never get them.
 489  */
 490 #if 0
 491 static
 492 int
 493 hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result)
 494 {
 495   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
 496   const int maxiter = 5000;
 497   int n = 1;
 498   double Anm2 = 1.0;
 499   double Bnm2 = 0.0;
 500   double Anm1 = 0.0;
 501   double Bnm1 = 1.0;
 502   double a1 = 1.0;
 503   double b1 = 1.0;
 504   double An = b1*Anm1 + a1*Anm2;
 505   double Bn = b1*Bnm1 + a1*Bnm2;
 506   double an, bn;
 507   double fn = An/Bn;
 508
 509   while(n < maxiter) {
 510     double old_fn;
 511     double del;
 512     n++;
 513     Anm2 = Anm1;
 514     Bnm2 = Bnm1;
 515     Anm1 = An;
 516     Bnm1 = Bn;
 517     an = (a+n)*x/((b-x+n-1)*(b-x+n));
 518     bn = 1.0;
 519     An = bn*Anm1 + an*Anm2;
 520     Bn = bn*Bnm1 + an*Bnm2;
 521
 522     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 523       An /= RECUR_BIG;
 524       Bn /= RECUR_BIG;
 525       Anm1 /= RECUR_BIG;
 526       Bnm1 /= RECUR_BIG;
 527       Anm2 /= RECUR_BIG;
 528       Bnm2 /= RECUR_BIG;
 529     }
 530
 531     old_fn = fn;
 532     fn = An/Bn;
 533     del = old_fn/fn;
 534
 535     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
 536   }
 537
 538   *result = a/(b-x) * fn;
 539
 540   if(n == maxiter)
 541     GSL_ERROR ("error", GSL_EMAXITER);
 542   else
 543     return GSL_SUCCESS;
 544 }
 545 #endif /* 0 */
 546
 547
 548 /* 1F1'(a,b,x)/1F1(a,b,x)
 549  * Uses Gautschi's series transformation of the
 550  * continued fraction. This is apparently the best
 551  * method for getting this ratio in the stable region.
 552  * The convergence is monotone and supergeometric
 553  * when b > x.
 554  * Assumes a >= -1.
 555  */
 556 static
 557 int
 558 hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result)
 559 {
 560   if(a == 0.0) {
 561     *result = 0.0;
 562     return GSL_SUCCESS;
 563   }
 564   else {
 565     const int maxiter = 5000;
 566     double sum  = 1.0;
 567     double pk   = 1.0;
 568     double rhok = 0.0;
 569     int k;
 570     for(k=1; k<maxiter; k++) {
 571       double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k));
 572       rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok));
 573       pk  *= rhok;
 574       sum += pk;
 575       if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break;
 576     }
 577     *result = a/(b-x) * sum;
 578     if(k == maxiter)
 579       GSL_ERROR ("error", GSL_EMAXITER);
 580     else
 581       return GSL_SUCCESS;
 582   }
 583 }
 584
 585
 586 /* 1F1(a+1,b,x)/1F1(a,b,x)
 587  *
 588  * I think this suffers from typical "anomalous convergence".
 589  * I could not find a region where it was truly useful.
 590  */
 591 #if 0
 592 static
 593 int
 594 hyperg_1F1_CF1(const double a, const double b, const double x, double * result)
 595 {
 596   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
 597   const int maxiter = 5000;
 598   int n = 1;
 599   double Anm2 = 1.0;
 600   double Bnm2 = 0.0;
 601   double Anm1 = 0.0;
 602   double Bnm1 = 1.0;
 603   double a1 = b - a - 1.0;
 604   double b1 = b - x - 2.0*(a+1.0);
 605   double An = b1*Anm1 + a1*Anm2;
 606   double Bn = b1*Bnm1 + a1*Bnm2;
 607   double an, bn;
 608   double fn = An/Bn;
 609
 610   while(n < maxiter) {
 611     double old_fn;
 612     double del;
 613     n++;
 614     Anm2 = Anm1;
 615     Bnm2 = Bnm1;
 616     Anm1 = An;
 617     Bnm1 = Bn;
 618     an = (a + n - 1.0) * (b - a - n);
 619     bn = b - x - 2.0*(a+n);
 620     An = bn*Anm1 + an*Anm2;
 621     Bn = bn*Bnm1 + an*Bnm2;
 622
 623     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 624       An /= RECUR_BIG;
 625       Bn /= RECUR_BIG;
 626       Anm1 /= RECUR_BIG;
 627       Bnm1 /= RECUR_BIG;
 628       Anm2 /= RECUR_BIG;
 629       Bnm2 /= RECUR_BIG;
 630     }
 631
 632     old_fn = fn;
 633     fn = An/Bn;
 634     del = old_fn/fn;
 635
 636     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
 637   }
 638
 639   *result = fn;
 640   if(n == maxiter)
 641     GSL_ERROR ("error", GSL_EMAXITER);
 642   else
 643     return GSL_SUCCESS;
 644 }
 645 #endif /* 0 */
 646
 647
 648 /* 1F1(a,b+1,x)/1F1(a,b,x)
 649  *
 650  * This seemed to suffer from "anomalous convergence".
 651  * However, I have no theory for this recurrence.
 652  */
 653 #if 0
 654 static
 655 int
 656 hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result)
 657 {
 658   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
 659   const int maxiter = 5000;
 660   int n = 1;
 661   double Anm2 = 1.0;
 662   double Bnm2 = 0.0;
 663   double Anm1 = 0.0;
 664   double Bnm1 = 1.0;
 665   double a1 = b + 1.0;
 666   double b1 = (b + 1.0) * (b - x);
 667   double An = b1*Anm1 + a1*Anm2;
 668   double Bn = b1*Bnm1 + a1*Bnm2;
 669   double an, bn;
 670   double fn = An/Bn;
 671
 672   while(n < maxiter) {
 673     double old_fn;
 674     double del;
 675     n++;
 676     Anm2 = Anm1;
 677     Bnm2 = Bnm1;
 678     Anm1 = An;
 679     Bnm1 = Bn;
 680     an = (b + n) * (b + n - 1.0 - a) * x;
 681     bn = (b + n) * (b + n - 1.0 - x);
 682     An = bn*Anm1 + an*Anm2;
 683     Bn = bn*Bnm1 + an*Bnm2;
 684
 685     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 686       An /= RECUR_BIG;
 687       Bn /= RECUR_BIG;
 688       Anm1 /= RECUR_BIG;
 689       Bnm1 /= RECUR_BIG;
 690       Anm2 /= RECUR_BIG;
 691       Bnm2 /= RECUR_BIG;
 692     }
 693
 694     old_fn = fn;
 695     fn = An/Bn;
 696     del = old_fn/fn;
 697
 698     if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
 699   }
 700
 701   *result = fn;
 702   if(n == maxiter)
 703     GSL_ERROR ("error", GSL_EMAXITER);
 704   else
 705     return GSL_SUCCESS;
 706 }
 707 #endif /* 0 */
 708
 709
 710 /* 1F1(a,b,x)
 711  * |a| <= 1, b > 0
 712  */
 713 static
 714 int
 715 hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result)
 716 {
 717   const double bma = b-a;
 718   const double oma = 1.0-a;
 719   const double ap1mb = 1.0+a-b;
 720   const double abs_bma = fabs(bma);
 721   const double abs_oma = fabs(oma);
 722   const double abs_ap1mb = fabs(ap1mb);
 723
 724   const double ax = fabs(x);
 725
 726   if(a == 0.0) {
 727     result->val = 1.0;
 728     result->err = 0.0;
 729     return GSL_SUCCESS;
 730   }
 731   else if(a == 1.0 && b >= 1.0) {
 732     return hyperg_1F1_1(b, x, result);
 733   }
 734   else if(a == -1.0) {
 735     result->val  = 1.0 + a/b * x;
 736     result->err  = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x));
 737     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 738     return GSL_SUCCESS;
 739   }
 740   else if(b >= 1.4*ax) {
 741     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
 742   }
 743   else if(x > 0.0) {
 744     if(x > 100.0 && abs_bma*abs_oma < 0.5*x) {
 745       return hyperg_1F1_asymp_posx(a, b, x, result);
 746     }
 747     else if(b < 5.0e+06) {
 748       /* Recurse backward on b from
 749        * a suitably high point.
 750        */
 751       const double b_del = ceil(1.4*x-b) + 1.0;
 752       double bp = b + b_del;
 753       gsl_sf_result r_Mbp1;
 754       gsl_sf_result r_Mb;
 755       double Mbp1;
 756       double Mb;
 757       double Mbm1;
 758       int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1);
 759       int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp,     x, &r_Mb);
 760       const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val);
 761       Mbp1 = r_Mbp1.val;
 762       Mb   = r_Mb.val;
 763       while(bp > b+0.1) {
 764         /* Do backward recursion. */
 765         Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0);
 766         bp -= 1.0;
 767         Mbp1 = Mb;
 768         Mb   = Mbm1;
 769       }
 770       result->val  = Mb;
 771       result->err  = err_rat * (fabs(b_del)+1.0) * fabs(Mb);
 772       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb);
 773       return GSL_ERROR_SELECT_2(stat_0, stat_1);
 774     }
 775     else if (fabs(x) < fabs(b) && fabs(a*x) < sqrt(fabs(b)) * fabs(b-x)) {
 776       return hyperg_1F1_largebx(a, b, x, result);
 777     } else {
 778       return hyperg_1F1_large2bm4a(a, b, x, result);
 779     }
 780   }
 781   else {
 782     /* x < 0 and b not large compared to |x|
 783      */
 784     if(ax < 10.0 && b < 10.0) {
 785       return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
 786     }
 787     else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) {
 788       return hyperg_1F1_asymp_negx(a, b, x, result);
 789     }
 790     else {
 791       return hyperg_1F1_luke(a, b, x, result);
 792     }
 793   }
 794 }
 795
 796
 797 /* 1F1(b+eps,b,x)
 798  * |eps|<=1, b > 0
 799  */
 800 static
 801 int
 802 hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result)
 803 {
 804   if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) {
 805     /* If b-a is very small and x/b is not too large we can
 806      * use this explicit approximation.
 807      *
 808      * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...)
 809      *
 810      *   v2 = a/(2b^2(b+1))
 811      *   v3 = a(b-2a)/(3b^3(b+1)(b+2))
 812      *   ...
 813      *
 814      * See [Luke, Mathematical Functions and Their Approximations, p.292]
 815      *
 816      * This cannot be used for b near a negative integer or zero.
 817      * Also, if x/b is large the deviation from exp(x) behaviour grows.
 818      */
 819     double a = b + eps;
 820     gsl_sf_result exab;
 821     int stat_e = gsl_sf_exp_e(a*x/b, &exab);
 822     double v2 = a/(2.0*b*b*(b+1.0));
 823     double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0));
 824     double v  = v2 + v3 * x;
 825     double f  = (1.0 - eps*x*x*v);
 826     result->val  = exab.val * f;
 827     result->err  = exab.err * fabs(f);
 828     result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v));
 829     result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
 830     return stat_e;
 831   }
 832   else {
 833     /* Otherwise use a Kummer transformation to reduce
 834      * it to the small a case.
 835      */
 836     gsl_sf_result Kummer_1F1;
 837     int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1);
 838     if(Kummer_1F1.val != 0.0) {
 839       int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x),
 840                                             Kummer_1F1.val, Kummer_1F1.err,
 841                                             result);
 842       return GSL_ERROR_SELECT_2(stat_e, stat_K);
 843     }
 844     else {
 845       result->val = 0.0;
 846       result->err = 0.0;
 847       return stat_K;
 848     }
 849   }
 850 }
 851
 852
 853 /* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2)
 854  *
 855  * E(x) = exp(x) x > 0
 856  *      = 1      x < 0
 857  *
 858  * a >= 1/2
 859  */
 860 static
 861 int
 862 hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result)
 863 {
 864   if(x == 0.0) {
 865     result->val = 1.0;
 866     result->err = 0.0;
 867     return GSL_SUCCESS;
 868   }
 869   else {
 870     gsl_sf_result I;
 871     int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I);
 872     gsl_sf_result lg;
 873     int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg);
 874     double ln_term   = (0.5-a)*log(0.25*fabs(x));
 875     double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term;
 876     double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x));
 877     int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
 878                                           I.val, I.err,
 879                                           result);
 880     return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I);
 881   }
 882 }
 883
 884
 885 /* Determine middle parts of diagonal recursion along b=2a
 886  * from two endpoints, i.e.
 887  *
 888  * given:  M(a,b)      and  M(a+1,b+2)
 889  * get:    M(a+1,b+1)  and  M(a,b+1)
 890  */
 891 #if 0
 892 inline
 893 static
 894 int
 895 hyperg_1F1_diag_step(const double a, const double b, const double x,
 896                      const double Mab, const double Map1bp2,
 897                      double * Map1bp1, double * Mabp1)
 898 {
 899   if(a == b) {
 900     *Map1bp1 = Mab;
 901     *Mabp1   = Mab - x/(b+1.0) * Map1bp2;
 902   }
 903   else {
 904     *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
 905     *Mabp1   = (a * *Map1bp1 - b * Mab)/(a-b);
 906   }
 907   return GSL_SUCCESS;
 908 }
 909 #endif /* 0 */
 910
 911
 912 /* Determine endpoint of diagonal recursion.
 913  *
 914  * given:  M(a,b)    and  M(a+1,b+2)
 915  * get:    M(a+1,b)  and  M(a+1,b+1)
 916  */
 917 #if 0
 918 inline
 919 static
 920 int
 921 hyperg_1F1_diag_end_step(const double a, const double b, const double x,
 922                          const double Mab, const double Map1bp2,
 923                          double * Map1b, double * Map1bp1)
 924 {
 925   *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2;
 926   *Map1b   = Mab + x/b * *Map1bp1;
 927   return GSL_SUCCESS;
 928 }
 929 #endif /* 0 */
 930
 931
 932 /* Handle the case of a and b both positive integers.
 933  * Assumes a > 0 and b > 0.
 934  */
 935 static
 936 int
 937 hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result)
 938 {
 939   double ax = fabs(x);
 940
 941   if(a == b) {
 942     return gsl_sf_exp_e(x, result);             /* 1F1(a,a,x) */
 943   }
 944   else if(a == 1) {
 945     return gsl_sf_exprel_n_e(b-1, x, result);   /* 1F1(1,b,x) */
 946   }
 947   else if(b == a + 1) {
 948     gsl_sf_result K;
 949     int stat_K = gsl_sf_exprel_n_e(a, -x, &K);  /* 1F1(1,1+a,-x) */
 950     int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
 951                                           K.val, K.err,
 952                                           result);
 953     return GSL_ERROR_SELECT_2(stat_e, stat_K);
 954   }
 955   else if(a == b + 1) {
 956     gsl_sf_result ex;
 957     int stat_e = gsl_sf_exp_e(x, &ex);
 958     result->val  = ex.val * (1.0 + x/b);
 959     result->err  = ex.err * (1.0 + x/b);
 960     result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b));
 961     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 962     return stat_e;
 963   }
 964   else if(a == b + 2) {
 965     gsl_sf_result ex;
 966     int stat_e = gsl_sf_exp_e(x, &ex);
 967     double poly  = (1.0 + x/b*(2.0 + x/(b+1.0)));
 968     result->val  = ex.val * poly;
 969     result->err  = ex.err * fabs(poly);
 970     result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0))));
 971     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 972     return stat_e;
 973   }
 974   else if(b == 2*a) {
 975     return hyperg_1F1_beq2a_pos(a, x, result);  /* 1F1(a,2a,x) */
 976   }
 977   else if(   ( b < 10 && a < 10 && ax < 5.0 )
 978           || ( b > a*ax )
 979           || ( b > a && ax < 5.0 )
 980     ) {
 981     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
 982   }
 983   else if(b > a && b >= 2*a + x) {
 984     /* Use the Gautschi CF series, then
 985      * recurse backward to a=0 for normalization.
 986      * This will work for either sign of x.
 987      */
 988     double rap;
 989     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
 990     double ra = 1.0 + x/a * rap;
 991     double Ma   = GSL_SQRT_DBL_MIN;
 992     double Map1 = ra * Ma;
 993     double Mnp1 = Map1;
 994     double Mn   = Ma;
 995     double Mnm1;
 996     int n;
 997     for(n=a; n>0; n--) {
 998       Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n);
 999       Mnp1 = Mn;
1000       Mn   = Mnm1;
1001     }
1002     result->val = Ma/Mn;
1003     result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn);
1004     return stat_CF1;
1005   }
1006   else if(b > a && b < 2*a + x && b > x) {
1007     /* Use the Gautschi series representation of
1008      * the continued fraction. Then recurse forward
1009      * to the a=b line for normalization. This will
1010      * work for either sign of x, although we do need
1011      * to check for b > x, for when x is positive.
1012      */
1013     double rap;
1014     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
1015     double ra = 1.0 + x/a * rap;
1016     gsl_sf_result ex;
1017     int stat_ex;
1018
1019     double Ma   = GSL_SQRT_DBL_MIN;
1020     double Map1 = ra * Ma;
1021     double Mnm1 = Ma;
1022     double Mn   = Map1;
1023     double Mnp1;
1024     int n;
1025     for(n=a+1; n<b; n++) {
1026       Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1027       Mnm1 = Mn;
1028       Mn   = Mnp1;
1029     }
1030
1031     stat_ex = gsl_sf_exp_e(x, &ex);  /* 1F1(b,b,x) */
1032     result->val  = ex.val * Ma/Mn;
1033     result->err  = ex.err * fabs(Ma/Mn);
1034     result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
1035     return GSL_ERROR_SELECT_2(stat_ex, stat_CF1);
1036   }
1037   else if(x >= 0.0) {
1038
1039     if(b < a) {
1040       /* The point b,b is below the b=2a+x line.
1041        * Forward recursion on a from b,b+1 is possible.
1042        * Note that a > b + 1 as well, since we already tried a = b + 1.
1043        */
1044       if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) {
1045         double ex = exp(x);
1046         int n;
1047         double Mnm1 = ex;                 /* 1F1(b,b,x)   */
1048         double Mn   = ex * (1.0 + x/b);   /* 1F1(b+1,b,x) */
1049         double Mnp1;
1050         for(n=b+1; n<a; n++) {
1051           Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1052           Mnm1 = Mn;
1053           Mn   = Mnp1;
1054         }
1055         result->val  = Mn;
1056         result->err  = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn);
1057         result->err *= fabs(a-b)+1.0;
1058         return GSL_SUCCESS;
1059       }
1060       else {
1061         OVERFLOW_ERROR(result);
1062       }
1063     }
1064     else {
1065       /* b > a
1066        * b < 2a + x
1067        * b <= x (otherwise we would have finished above)
1068        *
1069        * Gautschi anomalous convergence region. However, we can
1070        * recurse forward all the way from a=0,1 because we are
1071        * always underneath the b=2a+x line.
1072        */
1073       gsl_sf_result r_Mn;
1074       double Mnm1 = 1.0;    /* 1F1(0,b,x) */
1075       double Mn;            /* 1F1(1,b,x)  */
1076       double Mnp1;
1077       int n;
1078       gsl_sf_exprel_n_e(b-1, x, &r_Mn);
1079       Mn = r_Mn.val;
1080       for(n=1; n<a; n++) {
1081         Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1082         Mnm1 = Mn;
1083         Mn   = Mnp1;
1084       }
1085       result->val  = Mn;
1086       result->err  = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val);
1087       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
1088       return GSL_SUCCESS;
1089     }
1090   }
1091   else {
1092     /* x < 0
1093      * b < a (otherwise we would have tripped one of the above)
1094      */
1095
1096     if(a <= 0.5*(b-x) || a >= -x) {
1097       /* Gautschi continued fraction is in the anomalous region,
1098        * so we must find another way. We recurse down in b,
1099        * from the a=b line.
1100        */
1101       double ex = exp(x);
1102       double Manp1 = ex;
1103       double Man   = ex * (1.0 + x/(a-1.0));
1104       double Manm1;
1105       int n;
1106       for(n=a-1; n>b; n--) {
1107         Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
1108         Manp1 = Man;
1109         Man = Manm1;
1110       }
1111       result->val  = Man;
1112       result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man);
1113       result->err *= fabs(b-a)+1.0;
1114       return GSL_SUCCESS;
1115     }
1116     else {
1117       /* Pick a0 such that b ~= 2a0 + x, then
1118        * recurse down in b from a0,a0 to determine
1119        * the values near the line b=2a+x. Then recurse
1120        * forward on a from a0.
1121        */
1122       int a0 = ceil(0.5*(b-x));
1123       double Ma0b;    /* M(a0,b)   */
1124       double Ma0bp1;  /* M(a0,b+1) */
1125       double Ma0p1b;  /* M(a0+1,b) */
1126       double Mnm1;
1127       double Mn;
1128       double Mnp1;
1129       int n;
1130       {
1131         double ex = exp(x);
1132         double Ma0np1 = ex;
1133         double Ma0n   = ex * (1.0 + x/(a0-1.0));
1134         double Ma0nm1;
1135         for(n=a0-1; n>b; n--) {
1136           Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
1137           Ma0np1 = Ma0n;
1138           Ma0n = Ma0nm1;
1139         }
1140         Ma0bp1 = Ma0np1;
1141         Ma0b   = Ma0n;
1142         Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b);
1143       }
1144
1145       /* Initialise the recurrence correctly BJG */
1146
1147       if (a0 >= a)
1148         {
1149           Mn = Ma0b;
1150         }
1151       else if (a0 + 1>= a)
1152         {
1153           Mn = Ma0p1b;
1154         }
1155       else
1156         {
1157           Mnm1 = Ma0b;
1158           Mn   = Ma0p1b;
1159
1160           for(n=a0+1; n<a; n++) {
1161             Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1162             Mnm1 = Mn;
1163             Mn   = Mnp1;
1164           }
1165         }
1166
1167       result->val  = Mn;
1168       result->err  = (fabs(x) + 1.0) * GSL_DBL_EPSILON *  fabs(Mn);
1169       result->err *= fabs(b-a)+1.0;
1170       return GSL_SUCCESS;
1171     }
1172   }
1173 }
1174
1175
1176 /* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner)
1177  * When the terms are all positive, this
1178  * must work. We will assume this here.
1179  */
1180 static
1181 int
1182 hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result)
1183 {
1184   if(a == 0) {
1185     result->val = 1.0;
1186     result->err = 1.0;
1187     return GSL_SUCCESS;
1188   }
1189   else {
1190     int N = -a;
1191     double poly = 1.0;
1192     int k;
1193     for(k=N-1; k>=0; k--) {
1194       double t = (a+k)/(b+k) * (x/(k+1));
1195       double r = t + 1.0/poly;
1196       if(r > 0.9*GSL_DBL_MAX/poly) {
1197         OVERFLOW_ERROR(result);
1198       }
1199       else {
1200         poly *= r;  /* P_n = 1 + t_n P_{n-1} */
1201       }
1202     }
1203     result->val = poly;
1204     result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly);
1205     return GSL_SUCCESS;
1206   }
1207 }
1208
1209
1210 /* Evaluate negative integer a case by relation
1211  * to Laguerre polynomials. This is more general than
1212  * the direct polynomial evaluation, but is safe
1213  * for all values of x.
1214  *
1215  * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x]
1216  *             = n B(b,n) Laguerre[n,b-1,x]
1217  *
1218  * assumes b is not a negative integer
1219  */
1220 static
1221 int
1222 hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result)
1223 {
1224   const int n = -a;
1225
1226   gsl_sf_result lag;
1227   const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag);
1228   if(b < 0.0) {
1229     gsl_sf_result lnfact;
1230     gsl_sf_result lng1;
1231     gsl_sf_result lng2;
1232     double s1, s2;
1233     const int stat_f  = gsl_sf_lnfact_e(n, &lnfact);
1234     const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1);
1235     const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2);
1236     const double lnpre_val = lnfact.val - (lng1.val - lng2.val);
1237     const double lnpre_err = lnfact.err + lng1.err + lng2.err
1238       + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
1239     const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err,
1240                                                 s1*s2*lag.val, lag.err,
1241                                                 result);
1242     return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f);
1243   }
1244   else {
1245     gsl_sf_result lnbeta;
1246     gsl_sf_lnbeta_e(b, n, &lnbeta);
1247     if(fabs(lnbeta.val) < 0.1) {
1248       /* As we have noted, when B(x,y) is near 1,
1249        * evaluating log(B(x,y)) is not accurate.
1250        * Instead we evaluate B(x,y) directly.
1251        */
1252       const double ln_term_val = log(1.25*n);
1253       const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val;
1254       gsl_sf_result beta;
1255       int stat_b = gsl_sf_beta_e(b, n, &beta);
1256       int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
1257                                             lag.val, lag.err,
1258                                             result);
1259       result->val *= beta.val/1.25;
1260       result->err *= beta.val/1.25;
1261       return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b);
1262     }
1263     else {
1264       /* B(x,y) was not near 1, so it is safe to use
1265        * the logarithmic values.
1266        */
1267       const double ln_n = log(n);
1268       const double ln_term_val = lnbeta.val + ln_n;
1269       const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n);
1270       int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err,
1271                                             lag.val, lag.err,
1272                                             result);
1273       return GSL_ERROR_SELECT_2(stat_e, stat_l);
1274     }
1275   }
1276 }
1277
1278
1279 /* Handle negative integer a case for x > 0 and
1280  * generic b.
1281  *
1282  * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27]
1283  * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x)
1284  */
1285 #if 0
1286 static
1287 int
1288 hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result)
1289 {
1290   const int n = -a;
1291   const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 );
1292   double sgpoch;
1293   gsl_sf_result lnpoch;
1294   gsl_sf_result U;
1295   const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch);
1296   const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U);
1297   const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err,
1298                                               sgn * sgpoch * U.val, U.err,
1299                                               result);
1300   return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p);
1301 }
1302 #endif
1303
1304
1305 /* Assumes a <= -1,  b <= -1, and b <= a.
1306  */
1307 static
1308 int
1309 hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result)
1310 {
1311   if(x == 0.0) {
1312     result->val = 1.0;
1313     result->err = 0.0;
1314     return GSL_SUCCESS;
1315   }
1316   else if(x > 0.0) {
1317     return hyperg_1F1_a_negint_poly(a, b, x, result);
1318   }
1319   else {
1320     /* Apply a Kummer transformation to make x > 0 so
1321      * we can evaluate the polynomial safely. Of course,
1322      * this assumes b <= a, which must be true for
1323      * a<0 and b<0, since otherwise the thing is undefined.
1324      */
1325     gsl_sf_result K;
1326     int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K);
1327     int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x),
1328                                           K.val, K.err,
1329                                           result);
1330     return GSL_ERROR_SELECT_2(stat_e, stat_K);
1331   }
1332 }
1333
1334
1335 /* [Abramowitz+Stegun, 13.1.3]
1336  *
1337  * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) *
1338  *            { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) }
1339  *
1340  * b not an integer >= 2
1341  * a-b not a negative integer
1342  */
1343 static
1344 int
1345 hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result)
1346 {
1347   const double bp = 2.0 - b;
1348   const double ap = a - b + 1.0;
1349
1350   gsl_sf_result lg_ap, lg_bp;
1351   double sg_ap;
1352   int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap);
1353   int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp);
1354   int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1);
1355   double t1 = (bp-1.0) * log(x);
1356   double lnpre_val = lg_ap.val - lg_bp.val + t1;
1357   double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1);
1358
1359   gsl_sf_result lg_2mbp, lg_1papmbp;
1360   double sg_2mbp, sg_1papmbp;
1361   int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp,    &lg_2mbp,    &sg_2mbp);
1362   int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp);
1363   int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4);
1364   double lnc1_val = lg_2mbp.val - lg_1papmbp.val;
1365   double lnc1_err = lg_2mbp.err + lg_1papmbp.err
1366                     + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val));
1367
1368   gsl_sf_result M;
1369   gsl_sf_result_e10 U;
1370   int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M);
1371   int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U);
1372   int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U);
1373
1374   gsl_sf_result_e10 term_M;
1375   int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err,
1376                                              sg_2mbp*sg_1papmbp*M.val, M.err,
1377                                              &term_M);
1378
1379   const double ombp = 1.0 - bp;
1380   const double Uee_val = U.e10*M_LN10;
1381   const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val);
1382   const double Mee_val = term_M.e10*M_LN10;
1383   const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val);
1384   int stat_e1;
1385
1386   /* Do a little dance with the exponential prefactors
1387    * to avoid overflows in intermediate results.
1388    */
1389   if(Uee_val > Mee_val) {
1390     const double factorM_val = exp(Mee_val-Uee_val);
1391     const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val;
1392     const double inner_val = term_M.val*factorM_val - ombp*U.val;
1393     const double inner_err =
1394         term_M.err*factorM_val + fabs(ombp) * U.err
1395       + fabs(term_M.val) * factorM_err
1396       + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val));
1397     stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err,
1398                                        sg_ap*inner_val, inner_err,
1399                                        result);
1400   }
1401   else {
1402     const double factorU_val = exp(Uee_val - Mee_val);
1403     const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val;
1404     const double inner_val = term_M.val - ombp*factorU_val*U.val;
1405     const double inner_err =
1406         term_M.err + fabs(ombp*factorU_val*U.err)
1407       + fabs(ombp*factorU_err*U.val)
1408       + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val));
1409     stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err,
1410                                        sg_ap*inner_val, inner_err,
1411                                        result);
1412   }
1413
1414   return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2);
1415 }
1416
1417
1418 /* Handle case of generic positive a, b.
1419  * Assumes b-a is not a negative integer.
1420  */
1421 static
1422 int
1423 hyperg_1F1_ab_pos(const double a, const double b,
1424                   const double x,
1425                   gsl_sf_result * result)
1426 {
1427   const double ax = fabs(x);
1428
1429   if(   ( b < 10.0 && a < 10.0 && ax < 5.0 )
1430      || ( b > a*ax )
1431      || ( b > a && ax < 5.0 )
1432     ) {
1433     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
1434   }
1435   else if(   x < -100.0
1436           && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x)
1437     ) {
1438     /* Large negative x asymptotic.
1439      */
1440     return hyperg_1F1_asymp_negx(a, b, x, result);
1441   }
1442   else if(   x > 100.0
1443           && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x)
1444     ) {
1445     /* Large positive x asymptotic.
1446      */
1447     return hyperg_1F1_asymp_posx(a, b, x, result);
1448   }
1449   else if(fabs(b-a) <= 1.0) {
1450     /* Directly handle b near a.
1451      */
1452     return hyperg_1F1_beps_bgt0(a-b, b, x, result);  /* a = b + eps */
1453   }
1454
1455   else if(b > a && b >= 2*a + x) {
1456     /* Use the Gautschi CF series, then
1457      * recurse backward to a near 0 for normalization.
1458      * This will work for either sign of x.
1459      */
1460     double rap;
1461     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
1462     double ra = 1.0 + x/a * rap;
1463
1464     double Ma   = GSL_SQRT_DBL_MIN;
1465     double Map1 = ra * Ma;
1466     double Mnp1 = Map1;
1467     double Mn   = Ma;
1468     double Mnm1;
1469     gsl_sf_result Mn_true;
1470     int stat_Mt;
1471     double n;
1472     for(n=a; n>0.5; n -= 1.0) {
1473       Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n);
1474       Mnp1 = Mn;
1475       Mn   = Mnm1;
1476     }
1477
1478     stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true);
1479
1480     result->val  = (Ma/Mn) * Mn_true.val;
1481     result->err  = fabs(Ma/Mn) * Mn_true.err;
1482     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val);
1483     return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
1484   }
1485   else if(b > a && b < 2*a + x && b > x) {
1486     /* Use the Gautschi series representation of
1487      * the continued fraction. Then recurse forward
1488      * to near the a=b line for normalization. This will
1489      * work for either sign of x, although we do need
1490      * to check for b > x, which is relevant when x is positive.
1491      */
1492     gsl_sf_result Mn_true;
1493     int stat_Mt;
1494     double rap;
1495     int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap);
1496     double ra = 1.0 + x/a * rap;
1497     double Ma   = GSL_SQRT_DBL_MIN;
1498     double Mnm1 = Ma;
1499     double Mn   = ra * Mnm1;
1500     double Mnp1;
1501     double n;
1502     for(n=a+1.0; n<b-0.5; n += 1.0) {
1503       Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1504       Mnm1 = Mn;
1505       Mn   = Mnp1;
1506     }
1507     stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true);
1508     result->val  = Ma/Mn * Mn_true.val;
1509     result->err  = fabs(Ma/Mn) * Mn_true.err;
1510     result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val);
1511     return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1);
1512   }
1513   else if(x >= 0.0) {
1514
1515     if(b < a) {
1516       /* Forward recursion on a from a=b+eps-1,b+eps.
1517        */
1518       double N   = floor(a-b);
1519       double eps = a - b - N;
1520       gsl_sf_result r_M0;
1521       gsl_sf_result r_M1;
1522       int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0);
1523       int stat_1 = hyperg_1F1_beps_bgt0(eps,     b, x, &r_M1);
1524       double M0 = r_M0.val;
1525       double M1 = r_M1.val;
1526
1527       double Mam1 = M0;
1528       double Ma   = M1;
1529       double Map1;
1530       double ap;
1531       double start_pair = fabs(M0) + fabs(M1);
1532       double minim_pair = GSL_DBL_MAX;
1533       double pair_ratio;
1534       double rat_0 = fabs(r_M0.err/r_M0.val);
1535       double rat_1 = fabs(r_M1.err/r_M1.val);
1536       for(ap=b+eps; ap<a-0.1; ap += 1.0) {
1537         Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap;
1538         Mam1 = Ma;
1539         Ma   = Map1;
1540         minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair);
1541       }
1542       pair_ratio = start_pair/minim_pair;
1543       result->val  = Ma;
1544       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma);
1545       result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma);
1546       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma);
1547       return GSL_ERROR_SELECT_2(stat_0, stat_1);
1548     }
1549     else {
1550       /* b > a
1551        * b < 2a + x
1552        * b <= x
1553        *
1554        * Recurse forward on a from a=eps,eps+1.
1555        */
1556       double eps = a - floor(a);
1557       gsl_sf_result r_Mnm1;
1558       gsl_sf_result r_Mn;
1559       int stat_0 = hyperg_1F1_small_a_bgt0(eps,     b, x, &r_Mnm1);
1560       int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn);
1561       double Mnm1 = r_Mnm1.val;
1562       double Mn   = r_Mn.val;
1563       double Mnp1;
1564
1565       double n;
1566       double start_pair = fabs(Mn) + fabs(Mnm1);
1567       double minim_pair = GSL_DBL_MAX;
1568       double pair_ratio;
1569       double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val);
1570       double rat_1 = fabs(r_Mn.err/r_Mn.val);
1571       for(n=eps+1.0; n<a-0.1; n++) {
1572         Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1573         Mnm1 = Mn;
1574         Mn   = Mnp1;
1575         minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair);
1576       }
1577       pair_ratio = start_pair/minim_pair;
1578       result->val  = Mn;
1579       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn);
1580       result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn);
1581       result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn);
1582       return GSL_ERROR_SELECT_2(stat_0, stat_1);
1583     }
1584   }
1585   else {
1586     /* x < 0
1587      * b < a
1588      */
1589
1590     if(a <= 0.5*(b-x) || a >= -x) {
1591       /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1.
1592        */
1593       double N   = floor(a - b);
1594       double eps = 1.0 + N - a + b;
1595       gsl_sf_result r_Manp1;
1596       gsl_sf_result r_Man;
1597       int stat_0 = hyperg_1F1_beps_bgt0(-eps,    a+eps,     x, &r_Manp1);
1598       int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man);
1599       double Manp1 = r_Manp1.val;
1600       double Man   = r_Man.val;
1601       double Manm1;
1602
1603       double n;
1604       double start_pair = fabs(Manp1) + fabs(Man);
1605       double minim_pair = GSL_DBL_MAX;
1606       double pair_ratio;
1607       double rat_0 = fabs(r_Manp1.err/r_Manp1.val);
1608       double rat_1 = fabs(r_Man.err/r_Man.val);
1609       for(n=a+eps-1.0; n>b+0.1; n -= 1.0) {
1610         Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0));
1611         Manp1 = Man;
1612         Man = Manm1;
1613         minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair);
1614       }
1615
1616       /* FIXME: this is a nasty little hack; there is some
1617          (transient?) instability in this recurrence for some
1618          values. I can tell when it happens, which is when
1619          this pair_ratio is large. But I do not know how to
1620          measure the error in terms of it. I guessed quadratic
1621          below, but it is probably worse than that.
1622          */
1623       pair_ratio = start_pair/minim_pair;
1624       result->val  = Man;
1625       result->err  = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man);
1626       result->err *= pair_ratio*pair_ratio + 1.0;
1627       return GSL_ERROR_SELECT_2(stat_0, stat_1);
1628     }
1629     else {
1630       /* Pick a0 such that b ~= 2a0 + x, then
1631        * recurse down in b from a0,a0 to determine
1632        * the values near the line b=2a+x. Then recurse
1633        * forward on a from a0.
1634        */
1635       double epsa = a - floor(a);
1636       double a0   = floor(0.5*(b-x)) + epsa;
1637       double N    = floor(a0 - b);
1638       double epsb = 1.0 + N - a0 + b;
1639       double Ma0b;
1640       double Ma0bp1;
1641       double Ma0p1b;
1642       int stat_a0;
1643       double Mnm1;
1644       double Mn;
1645       double Mnp1;
1646       double n;
1647       double err_rat;
1648       {
1649         gsl_sf_result r_Ma0np1;
1650         gsl_sf_result r_Ma0n;
1651         int stat_0 = hyperg_1F1_beps_bgt0(-epsb,    a0+epsb,     x, &r_Ma0np1);
1652         int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n);
1653         double Ma0np1 = r_Ma0np1.val;
1654         double Ma0n   = r_Ma0n.val;
1655         double Ma0nm1;
1656
1657         err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val);
1658
1659         for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) {
1660           Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0));
1661           Ma0np1 = Ma0n;
1662           Ma0n = Ma0nm1;
1663         }
1664         Ma0bp1 = Ma0np1;
1665         Ma0b   = Ma0n;
1666         Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */
1667         stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1);
1668       }
1669
1670
1671       /* Initialise the recurrence correctly BJG */
1672
1673       if (a0 >= a - 0.1)
1674         {
1675           Mn = Ma0b;
1676         }
1677       else if (a0 + 1>= a - 0.1)
1678         {
1679           Mn = Ma0p1b;
1680         }
1681       else
1682         {
1683           Mnm1 = Ma0b;
1684           Mn   = Ma0p1b;
1685
1686           for(n=a0+1.0; n<a-0.1; n += 1.0) {
1687             Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n;
1688             Mnm1 = Mn;
1689             Mn   = Mnp1;
1690           }
1691         }
1692
1693       result->val = Mn;
1694       result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn);
1695       return stat_a0;
1696     }
1697   }
1698 }
1699
1700
1701 /* Assumes b != integer
1702  * Assumes a != integer when x > 0
1703  * Assumes b-a != neg integer when x < 0
1704  */
1705 static
1706 int
1707 hyperg_1F1_ab_neg(const double a, const double b, const double x,
1708                   gsl_sf_result * result)
1709 {
1710   const double bma = b - a;
1711   const double abs_x = fabs(x);
1712   const double abs_a = fabs(a);
1713   const double abs_b = fabs(b);
1714   const double size_a = GSL_MAX(abs_a, 1.0);
1715   const double size_b = GSL_MAX(abs_b, 1.0);
1716   const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD );
1717
1718   if(   (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0)
1719      || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x))
1720     ) {
1721     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
1722   }
1723   else if(   x > 0.0
1724           && size_b > size_a
1725           && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0
1726     ) {
1727     /* Series terms are positive definite up until
1728      * there is a sign change. But by then the
1729      * terms are small due to the last condition.
1730      */
1731     return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
1732   }
1733   else if(   (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0)
1734           || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x)
1735     ) {
1736     /* Use Kummer transformation to render series safe.
1737      */
1738     gsl_sf_result Kummer_1F1;
1739     int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1);
1740     int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1741                                       Kummer_1F1.val, Kummer_1F1.err,
1742                                       result);
1743     return GSL_ERROR_SELECT_2(stat_e, stat_K);
1744   }
1745   else if(   x < -30.0
1746           && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)
1747     ) {
1748     /* Large negative x asymptotic.
1749      * Note that we do not check if b-a is a negative integer.
1750      */
1751     return hyperg_1F1_asymp_negx(a, b, x, result);
1752   }
1753   else if(   x > 100.0
1754           && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x)
1755     ) {
1756     /* Large positive x asymptotic.
1757      * Note that we do not check if a is a negative integer.
1758      */
1759     return hyperg_1F1_asymp_posx(a, b, x, result);
1760   }
1761   else if(x > 0.0 && !(bma_integer && bma > 0.0)) {
1762     return hyperg_1F1_U(a, b, x, result);
1763   }
1764   else {
1765     /* FIXME:  if all else fails, try the series... BJG */
1766     if (x < 0.0) {
1767       /* Apply Kummer Transformation */
1768       int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result);
1769       double K_factor = exp(x);
1770       result->val *= K_factor;
1771       result->err *= K_factor;
1772       return status;
1773     } else {
1774       int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result);
1775       return status;
1776     }
1777
1778     /* Sadness... */
1779     /* result->val = 0.0; */
1780     /* result->err = 0.0; */
1781     /* GSL_ERROR ("error", GSL_EUNIMPL); */
1782   }
1783 }
1784
1785
1786 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
1787
1788 int
1789 gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result)
1790 {
1791   /* CHECK_POINTER(result) */
1792
1793   if(x == 0.0) {
1794     result->val = 1.0;
1795     result->err = 0.0;
1796     return GSL_SUCCESS;
1797   }
1798   else if(a == b) {
1799     return gsl_sf_exp_e(x, result);
1800   }
1801   else if(b == 0) {
1802     DOMAIN_ERROR(result);
1803   }
1804   else if(a == 0) {
1805     result->val = 1.0;
1806     result->err = 0.0;
1807     return GSL_SUCCESS;
1808   }
1809   else if(b < 0 && (a < b || a > 0)) {
1810     /* Standard domain error due to singularity. */
1811     DOMAIN_ERROR(result);
1812   }
1813   else if(x > 100.0  && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) {
1814     /* x -> +Inf asymptotic */
1815     return hyperg_1F1_asymp_posx(a, b, x, result);
1816   }
1817   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)) {
1818     /* x -> -Inf asymptotic */
1819     return hyperg_1F1_asymp_negx(a, b, x, result);
1820   }
1821   else if(a < 0 && b < 0) {
1822     return hyperg_1F1_ab_negint(a, b, x, result);
1823   }
1824   else if(a < 0 && b > 0) {
1825     /* Use Kummer to reduce it to the positive integer case.
1826      * Note that b > a, strictly, since we already trapped b = a.
1827      */
1828     gsl_sf_result Kummer_1F1;
1829     int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1);
1830     int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1831                                       Kummer_1F1.val, Kummer_1F1.err,
1832                                       result);
1833     return GSL_ERROR_SELECT_2(stat_e, stat_K);
1834   }
1835   else {
1836     /* a > 0 and b > 0 */
1837     return hyperg_1F1_ab_posint(a, b, x, result);
1838   }
1839 }
1840
1841
1842 int
1843 gsl_sf_hyperg_1F1_e(const double a, const double b, const double x,
1844                        gsl_sf_result * result
1845                        )
1846 {
1847   const double bma = b - a;
1848   const double rinta = floor(a + 0.5);
1849   const double rintb = floor(b + 0.5);
1850   const double rintbma = floor(bma + 0.5);
1851   const int a_integer   = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX );
1852   const int b_integer   = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX );
1853   const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX );
1854   const int b_neg_integer   = ( b < -0.1 && b_integer );
1855   const int a_neg_integer   = ( a < -0.1 && a_integer );
1856   const int bma_neg_integer = ( bma < -0.1 &&  bma_integer );
1857
1858   /* CHECK_POINTER(result) */
1859
1860   if(x == 0.0) {
1861     /* Testing for this before testing a and b
1862      * is somewhat arbitrary. The result is that
1863      * we have 1F1(a,0,0) = 1.
1864      */
1865     result->val = 1.0;
1866     result->err = 0.0;
1867     return GSL_SUCCESS;
1868   }
1869   else if(b == 0.0) {
1870     DOMAIN_ERROR(result);
1871   }
1872   else if(a == 0.0) {
1873     result->val = 1.0;
1874     result->err = 0.0;
1875     return GSL_SUCCESS;
1876   }
1877   else if(a == b) {
1878     /* case: a==b; exp(x)
1879      * It's good to test exact equality now.
1880      * We also test approximate equality later.
1881      */
1882     return gsl_sf_exp_e(x, result);
1883   } else if(fabs(b) < _1F1_INT_THRESHOLD && fabs(a) < _1F1_INT_THRESHOLD) {
1884     /* a and b near zero: 1 + a/b (exp(x)-1)
1885      */
1886
1887     /* Note that neither a nor b is zero, since
1888      * we eliminated that with the above tests.
1889      */
1890
1891     gsl_sf_result exm1;
1892     int stat_e = gsl_sf_expm1_e(x, &exm1);
1893     double sa = ( a > 0.0 ? 1.0 : -1.0 );
1894     double sb = ( b > 0.0 ? 1.0 : -1.0 );
1895     double lnab = log(fabs(a/b)); /* safe */
1896     gsl_sf_result hx;
1897     int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab),
1898                                         sa * sb * exm1.val, exm1.err,
1899                                         &hx);
1900     result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val);  /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/
1901     result->err = hx.err;
1902     return GSL_ERROR_SELECT_2(stat_hx, stat_e);
1903   } else if (fabs(b) < _1F1_INT_THRESHOLD && fabs(x*a) < 1) {
1904       /* b near zero and a not near zero
1905        */
1906       const double m_arg = 1.0/(0.5*b);
1907       gsl_sf_result F_renorm;
1908       int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm);
1909       int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg,
1910                                             0.5*F_renorm.val, 0.5*F_renorm.err,
1911                                             result);
1912       return GSL_ERROR_SELECT_2(stat_m, stat_F);
1913   }
1914   else if(a_integer && b_integer) {
1915     /* Check for reduction to the integer case.
1916      * Relies on the arbitrary "near an integer" test.
1917      */
1918     return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result);
1919   }
1920   else if(b_neg_integer && !(a_neg_integer && a > b)) {
1921     /* Standard domain error due to
1922      * uncancelled singularity.
1923      */
1924     DOMAIN_ERROR(result);
1925   }
1926   else if(a_neg_integer) {
1927     return hyperg_1F1_a_negint_lag((int)rinta, b, x, result);
1928   }
1929   else if(b > 0.0) {
1930     if(-1.0 <= a && a <= 1.0) {
1931       /* Handle small a explicitly.
1932        */
1933       return hyperg_1F1_small_a_bgt0(a, b, x, result);
1934     }
1935     else if(bma_neg_integer) {
1936       /* Catch this now, to avoid problems in the
1937        * generic evaluation code.
1938        */
1939       gsl_sf_result Kummer_1F1;
1940       int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1);
1941       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1942                                             Kummer_1F1.val, Kummer_1F1.err,
1943                                             result);
1944       return GSL_ERROR_SELECT_2(stat_e, stat_K);
1945     }
1946     else if(a < 0.0 && fabs(x) < 100.0) {
1947       /* Use Kummer to reduce it to the generic positive case.
1948        * Note that b > a, strictly, since we already trapped b = a.
1949        * Also b-(b-a)=a, and a is not a negative integer here,
1950        * so the generic evaluation is safe.
1951        */
1952       gsl_sf_result Kummer_1F1;
1953       int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1);
1954       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1955                                             Kummer_1F1.val, Kummer_1F1.err,
1956                                             result);
1957       return GSL_ERROR_SELECT_2(stat_e, stat_K);
1958     }
1959     else if (a > 0) {
1960       /* a > 0.0 */
1961       return hyperg_1F1_ab_pos(a, b, x, result);
1962     } else {
1963       return gsl_sf_hyperg_1F1_series_e(a, b, x, result);
1964     }
1965   }
1966   else {
1967     /* b < 0.0 */
1968
1969     if(bma_neg_integer && x < 0.0) {
1970       /* Handle this now to prevent problems
1971        * in the generic evaluation.
1972        */
1973       gsl_sf_result K;
1974       int stat_K;
1975       int stat_e;
1976       if(a < 0.0) {
1977         /* Kummer transformed version of safe polynomial.
1978          * The condition a < 0 is equivalent to b < b-a,
1979          * which is the condition required for the series
1980          * to be positive definite here.
1981          */
1982         stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K);
1983       }
1984       else {
1985         /* Generic eval for negative integer a. */
1986         stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K);
1987       }
1988       stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1989                                         K.val, K.err,
1990                                         result);
1991       return GSL_ERROR_SELECT_2(stat_e, stat_K);
1992     }
1993     else if(a > 0.0) {
1994       /* Use Kummer to reduce it to the generic negative case.
1995        */
1996       gsl_sf_result K;
1997       int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K);
1998       int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x),
1999                                             K.val, K.err,
2000                                             result);
2001       return GSL_ERROR_SELECT_2(stat_e, stat_K);
2002     }
2003     else {
2004       return hyperg_1F1_ab_neg(a, b, x, result);
2005     }
2006   }
2007 }
2008
2009
2010
2011 #if 0
2012     /* Luke in the canonical case.
2013    */
2014   if(x < 0.0 && !a_neg_integer && !bma_neg_integer) {
2015     double prec;
2016     return hyperg_1F1_luke(a, b, x, result, &prec);
2017   }
2018
2019
2020   /* Luke with Kummer transformation.
2021    */
2022   if(x > 0.0 && !a_neg_integer && !bma_neg_integer) {
2023     double prec;
2024     double Kummer_1F1;
2025     double ex;
2026     int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec);
2027     int stat_e = gsl_sf_exp_e(x, &ex);
2028     if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) {
2029       double lnr = log(fabs(Kummer_1F1)) + x;
2030       if(lnr < GSL_LOG_DBL_MAX) {
2031         *result = ex * Kummer_1F1;
2032         return GSL_SUCCESS;
2033       }
2034       else {
2035         *result = GSL_POSINF;
2036         GSL_ERROR ("overflow", GSL_EOVRFLW);
2037       }
2038     }
2039     else if(stat_F != GSL_SUCCESS) {
2040       *result = 0.0;
2041       return stat_F;
2042     }
2043     else {
2044       *result = 0.0;
2045       return stat_e;
2046     }
2047   }
2048 #endif
2049
2050
2051
2052 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
2053
2054 #include "eval.h"
2055
2056 double gsl_sf_hyperg_1F1_int(const int m, const int n, double x)
2057 {
2058   EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result));
2059 }
2060
2061 double gsl_sf_hyperg_1F1(double a, double b, double x)
2062 {
2063   EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result));
2064 }