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

   1 /* specfunc/legendre_H3d.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_exp.h>
  26 #include <gsl/gsl_sf_gamma.h>
  27 #include <gsl/gsl_sf_trig.h>
  28 #include <gsl/gsl_sf_legendre.h>
  29
  30 #include "error.h"
  31
  32 #include "legendre.h"
  33
  34 /* See [Abbott+Schaefer, Ap.J. 308, 546 (1986)] for
  35  * enough details to follow what is happening here.
  36  */
  37
  38
  39 /* Logarithm of normalization factor, Log[N(ell,lambda)].
  40  * N(ell,lambda) = Product[ lambda^2 + n^2, {n,0,ell} ]
  41  *               = |Gamma(ell + 1 + I lambda)|^2  lambda sinh(Pi lambda) / Pi
  42  * Assumes ell >= 0.
  43  */
  44 static
  45 int
  46 legendre_H3d_lnnorm(const int ell, const double lambda, double * result)
  47 {
  48   double abs_lam = fabs(lambda);
  49
  50   if(abs_lam == 0.0) {
  51     *result = 0.0;
  52     GSL_ERROR ("error", GSL_EDOM);
  53   }
  54   else if(lambda > (ell + 1.0)/GSL_ROOT3_DBL_EPSILON) {
  55     /* There is a cancellation between the sinh(Pi lambda)
  56      * term and the log(gamma(ell + 1 + i lambda) in the
  57      * result below, so we show some care and save some digits.
  58      * Note that the above guarantees that lambda is large,
  59      * since ell >= 0. We use Stirling and a simple expansion
  60      * of sinh.
  61      */
  62     double rat = (ell+1.0)/lambda;
  63     double ln_lam2ell2  = 2.0*log(lambda) + log(1.0 + rat*rat);
  64     double lg_corrected = -2.0*(ell+1.0) + M_LNPI + (ell+0.5)*ln_lam2ell2 + 1.0/(288.0*lambda*lambda);
  65     double angle_terms  = lambda * 2.0 * rat * (1.0 - rat*rat/3.0);
  66     *result = log(abs_lam) + lg_corrected + angle_terms - M_LNPI;
  67     return GSL_SUCCESS;
  68   }
  69   else {
  70     gsl_sf_result lg_r;
  71     gsl_sf_result lg_theta;
  72     gsl_sf_result ln_sinh;
  73     gsl_sf_lngamma_complex_e(ell+1.0, lambda, &lg_r, &lg_theta);
  74     gsl_sf_lnsinh_e(M_PI * abs_lam, &ln_sinh);
  75     *result = log(abs_lam) + ln_sinh.val + 2.0*lg_r.val - M_LNPI;
  76     return GSL_SUCCESS;
  77   }
  78 }
  79
  80
  81 /* Calculate series for small eta*lambda.
  82  * Assumes eta > 0, lambda != 0.
  83  *
  84  * This is just the defining hypergeometric for the Legendre function.
  85  *
  86  * P^{mu}_{-1/2 + I lam}(z) = 1/Gamma(l+3/2) ((z+1)/(z-1)^(mu/2)
  87  *                            2F1(1/2 - I lam, 1/2 + I lam; l+3/2; (1-z)/2)
  88  * We use
  89  *       z = cosh(eta)
  90  * (z-1)/2 = sinh^2(eta/2)
  91  *
  92  * And recall
  93  * H3d = sqrt(Pi Norm /(2 lam^2 sinh(eta))) P^{-l-1/2}_{-1/2 + I lam}(cosh(eta))
  94  */
  95 static
  96 int
  97 legendre_H3d_series(const int ell, const double lambda, const double eta,
  98                     gsl_sf_result * result)
  99 {
 100   const int nmax = 5000;
 101   const double shheta = sinh(0.5*eta);
 102   const double ln_zp1 = M_LN2 + log(1.0 + shheta*shheta);
 103   const double ln_zm1 = M_LN2 + 2.0*log(shheta);
 104   const double zeta = -shheta*shheta;
 105   gsl_sf_result lg_lp32;
 106   double term = 1.0;
 107   double sum  = 1.0;
 108   double sum_err = 0.0;
 109   gsl_sf_result lnsheta;
 110   double lnN;
 111   double lnpre_val, lnpre_err, lnprepow;
 112   int stat_e;
 113   int n;
 114
 115   gsl_sf_lngamma_e(ell + 3.0/2.0, &lg_lp32);
 116   gsl_sf_lnsinh_e(eta, &lnsheta);
 117   legendre_H3d_lnnorm(ell, lambda, &lnN);
 118   lnprepow = 0.5*(ell + 0.5) * (ln_zm1 - ln_zp1);
 119   lnpre_val  = lnprepow + 0.5*(lnN + M_LNPI - M_LN2 - lnsheta.val) - lg_lp32.val - log(fabs(lambda));
 120   lnpre_err  = lnsheta.err + lg_lp32.err + GSL_DBL_EPSILON * fabs(lnpre_val);
 121   lnpre_err += 2.0*GSL_DBL_EPSILON * (fabs(lnN) + M_LNPI + M_LN2);
 122   lnpre_err += 2.0*GSL_DBL_EPSILON * (0.5*(ell + 0.5) * (fabs(ln_zm1) + fabs(ln_zp1)));
 123   for(n=1; n<nmax; n++) {
 124     double aR = n - 0.5;
 125     term *= (aR*aR + lambda*lambda)*zeta/(ell + n + 0.5)/n;
 126     sum  += term;
 127     sum_err += 2.0*GSL_DBL_EPSILON*fabs(term);
 128     if(fabs(term/sum) < 2.0 * GSL_DBL_EPSILON) break;
 129   }
 130
 131   stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, sum, fabs(term)+sum_err, result);
 132   return GSL_ERROR_SELECT_2(stat_e, (n==nmax ? GSL_EMAXITER : GSL_SUCCESS));
 133 }
 134
 135
 136 /* Evaluate legendre_H3d(ell+1)/legendre_H3d(ell)
 137  * by continued fraction.
 138  */
 139 #if 0
 140 static
 141 int
 142 legendre_H3d_CF1(const int ell, const double lambda, const double coth_eta,
 143                  gsl_sf_result * result)
 144 {
 145   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
 146   const int maxiter = 5000;
 147   int n = 1;
 148   double Anm2 = 1.0;
 149   double Bnm2 = 0.0;
 150   double Anm1 = 0.0;
 151   double Bnm1 = 1.0;
 152   double a1 = hypot(lambda, ell+1.0);
 153   double b1 = (2.0*ell + 3.0) * coth_eta;
 154   double An = b1*Anm1 + a1*Anm2;
 155   double Bn = b1*Bnm1 + a1*Bnm2;
 156   double an, bn;
 157   double fn = An/Bn;
 158
 159   while(n < maxiter) {
 160     double old_fn;
 161     double del;
 162     n++;
 163     Anm2 = Anm1;
 164     Bnm2 = Bnm1;
 165     Anm1 = An;
 166     Bnm1 = Bn;
 167     an = -(lambda*lambda + ((double)ell + n)*((double)ell + n));
 168     bn = (2.0*ell + 2.0*n + 1.0) * coth_eta;
 169     An = bn*Anm1 + an*Anm2;
 170     Bn = bn*Bnm1 + an*Bnm2;
 171
 172     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
 173       An /= RECUR_BIG;
 174       Bn /= RECUR_BIG;
 175       Anm1 /= RECUR_BIG;
 176       Bnm1 /= RECUR_BIG;
 177       Anm2 /= RECUR_BIG;
 178       Bnm2 /= RECUR_BIG;
 179     }
 180
 181     old_fn = fn;
 182     fn = An/Bn;
 183     del = old_fn/fn;
 184
 185     if(fabs(del - 1.0) < 4.0*GSL_DBL_EPSILON) break;
 186   }
 187
 188   result->val = fn;
 189   result->err = 2.0 * GSL_DBL_EPSILON * (sqrt(n)+1.0) * fabs(fn);
 190
 191   if(n >= maxiter)
 192     GSL_ERROR ("error", GSL_EMAXITER);
 193   else
 194     return GSL_SUCCESS;
 195 }
 196 #endif /* 0 */
 197
 198
 199 /* Evaluate legendre_H3d(ell+1)/legendre_H3d(ell)
 200  * by continued fraction. Use the Gautschi (Euler)
 201  * equivalent series.
 202  */
 203  /* FIXME: Maybe we have to worry about this. The a_k are
 204   * not positive and there can be a blow-up. It happened
 205   * for J_nu once or twice. Then we should probably use
 206   * the method above.
 207   */
 208 static
 209 int
 210 legendre_H3d_CF1_ser(const int ell, const double lambda, const double coth_eta,
 211                      gsl_sf_result * result)
 212 {
 213   const double pre = hypot(lambda, ell+1.0)/((2.0*ell+3)*coth_eta);
 214   const int maxk = 20000;
 215   double tk   = 1.0;
 216   double sum  = 1.0;
 217   double rhok = 0.0;
 218   double sum_err = 0.0;
 219   int k;
 220
 221   for(k=1; k<maxk; k++) {
 222     double tlk = (2.0*ell + 1.0 + 2.0*k);
 223     double l1k = (ell + 1.0 + k);
 224     double ak = -(lambda*lambda + l1k*l1k)/(tlk*(tlk+2.0)*coth_eta*coth_eta);
 225     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
 226     tk  *= rhok;
 227     sum += tk;
 228     sum_err += 2.0 * GSL_DBL_EPSILON * k * fabs(tk);
 229     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
 230   }
 231
 232   result->val  = pre * sum;
 233   result->err  = fabs(pre * tk);
 234   result->err += fabs(pre * sum_err);
 235   result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
 236
 237   if(k >= maxk)
 238     GSL_ERROR ("error", GSL_EMAXITER);
 239   else
 240     return GSL_SUCCESS;
 241 }
 242
 243
 244
 245 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 246
 247 int
 248 gsl_sf_legendre_H3d_0_e(const double lambda, const double eta, gsl_sf_result * result)
 249 {
 250   /* CHECK_POINTER(result) */
 251
 252   if(eta < 0.0) {
 253     DOMAIN_ERROR(result);
 254   }
 255   else if(eta == 0.0 || lambda == 0.0) {
 256     result->val = 1.0;
 257     result->err = 0.0;
 258     return GSL_SUCCESS;
 259   }
 260   else {
 261     const double lam_eta = lambda * eta;
 262     gsl_sf_result s;
 263     gsl_sf_sin_err_e(lam_eta, 2.0*GSL_DBL_EPSILON * fabs(lam_eta), &s);
 264     if(eta > -0.5*GSL_LOG_DBL_EPSILON) {
 265       double f = 2.0 / lambda * exp(-eta);
 266       result->val  = f * s.val;
 267       result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
 268       result->err += fabs(f) * s.err;
 269       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 270     }
 271     else {
 272       double f = 1.0/(lambda*sinh(eta));
 273       result->val  = f * s.val;
 274       result->err  = fabs(f * s.val) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
 275       result->err += fabs(f) * s.err;
 276       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 277     }
 278     return GSL_SUCCESS;
 279   }
 280 }
 281
 282
 283 int
 284 gsl_sf_legendre_H3d_1_e(const double lambda, const double eta, gsl_sf_result * result)
 285 {
 286   const double xi    = fabs(eta*lambda);
 287   const double lsq   = lambda*lambda;
 288   const double lsqp1 = lsq + 1.0;
 289
 290   /* CHECK_POINTER(result) */
 291
 292   if(eta < 0.0) {
 293     DOMAIN_ERROR(result);
 294   }
 295   else if(eta == 0.0 || lambda == 0.0) {
 296     result->val = 0.0;
 297     result->err = 0.0;
 298     return GSL_SUCCESS;
 299   }
 300   else if(xi < GSL_ROOT5_DBL_EPSILON && eta < GSL_ROOT5_DBL_EPSILON) {
 301     double etasq = eta*eta;
 302     double xisq  = xi*xi;
 303     double term1 = (etasq + xisq)/3.0;
 304     double term2 = -(2.0*etasq*etasq + 5.0*etasq*xisq + 3.0*xisq*xisq)/90.0;
 305     double sinh_term = 1.0 - eta*eta/6.0 * (1.0 - 7.0/60.0*eta*eta);
 306     double pre = sinh_term/sqrt(lsqp1) / eta;
 307     result->val  = pre * (term1 + term2);
 308     result->err  = pre * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2));
 309     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 310     return GSL_SUCCESS;
 311   }
 312   else {
 313     double sin_term;     /*  Sin(xi)/xi     */
 314     double cos_term;     /*  Cos(xi)        */
 315     double coth_term;    /*  eta/Tanh(eta)  */
 316     double sinh_term;    /*  eta/Sinh(eta)  */
 317     double sin_term_err;
 318     double cos_term_err;
 319     double t1;
 320     double pre_val;
 321     double pre_err;
 322     double term1;
 323     double term2;
 324     if(xi < GSL_ROOT5_DBL_EPSILON) {
 325       sin_term = 1.0 - xi*xi/6.0 * (1.0 - xi*xi/20.0);
 326       cos_term = 1.0 - 0.5*xi*xi * (1.0 - xi*xi/12.0);
 327       sin_term_err = GSL_DBL_EPSILON;
 328       cos_term_err = GSL_DBL_EPSILON;
 329     }
 330     else {
 331       gsl_sf_result sin_xi_result;
 332       gsl_sf_result cos_xi_result;
 333       gsl_sf_sin_e(xi, &sin_xi_result);
 334       gsl_sf_cos_e(xi, &cos_xi_result);
 335       sin_term = sin_xi_result.val/xi;
 336       cos_term = cos_xi_result.val;
 337       sin_term_err = sin_xi_result.err/fabs(xi);
 338       cos_term_err = cos_xi_result.err;
 339     }
 340     if(eta < GSL_ROOT5_DBL_EPSILON) {
 341       coth_term = 1.0 + eta*eta/3.0 * (1.0 - eta*eta/15.0);
 342       sinh_term = 1.0 - eta*eta/6.0 * (1.0 - 7.0/60.0*eta*eta);
 343     }
 344     else {
 345       coth_term = eta/tanh(eta);
 346       sinh_term = eta/sinh(eta);
 347     }
 348     t1 = sqrt(lsqp1) * eta;
 349     pre_val = sinh_term/t1;
 350     pre_err = 2.0 * GSL_DBL_EPSILON * fabs(pre_val);
 351     term1 = sin_term*coth_term;
 352     term2 = cos_term;
 353     result->val  = pre_val * (term1 - term2);
 354     result->err  = pre_err * fabs(term1 - term2);
 355     result->err += pre_val * (sin_term_err * coth_term + cos_term_err);
 356     result->err += pre_val * fabs(term1-term2) * (fabs(eta) + 1.0) * GSL_DBL_EPSILON;
 357     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 358     return GSL_SUCCESS;
 359   }
 360 }
 361
 362
 363 int
 364 gsl_sf_legendre_H3d_e(const int ell, const double lambda, const double eta,
 365                          gsl_sf_result * result)
 366 {
 367   const double abs_lam = fabs(lambda);
 368   const double lsq     = abs_lam*abs_lam;
 369   const double xi      = abs_lam * eta;
 370   const double cosh_eta = cosh(eta);
 371
 372   /* CHECK_POINTER(result) */
 373
 374   if(eta < 0.0) {
 375     DOMAIN_ERROR(result);
 376   }
 377   else if(eta > GSL_LOG_DBL_MAX) {
 378     /* cosh(eta) is too big. */
 379     OVERFLOW_ERROR(result);
 380   }
 381   else if(ell == 0) {
 382     return gsl_sf_legendre_H3d_0_e(lambda, eta, result);
 383   }
 384   else if(ell == 1) {
 385     return gsl_sf_legendre_H3d_1_e(lambda, eta, result);
 386   }
 387   else if(eta == 0.0) {
 388     result->val = 0.0;
 389     result->err = 0.0;
 390     return GSL_SUCCESS;
 391   }
 392   else if(xi < 1.0) {
 393     return legendre_H3d_series(ell, lambda, eta, result);
 394   }
 395   else if((ell*ell+lsq)/sqrt(1.0+lsq)/(cosh_eta*cosh_eta) < 5.0*GSL_ROOT3_DBL_EPSILON) {
 396     /* Large argument.
 397      */
 398     gsl_sf_result P;
 399     double lm;
 400     int stat_P = gsl_sf_conicalP_large_x_e(-ell-0.5, lambda, cosh_eta, &P, &lm);
 401     if(P.val == 0.0) {
 402       result->val = 0.0;
 403       result->err = 0.0;
 404       return stat_P;
 405     }
 406     else {
 407       double lnN;
 408       gsl_sf_result lnsh;
 409       double ln_abslam;
 410       double lnpre_val, lnpre_err;
 411       int stat_e;
 412       gsl_sf_lnsinh_e(eta, &lnsh);
 413       legendre_H3d_lnnorm(ell, lambda, &lnN);
 414       ln_abslam = log(abs_lam);
 415       lnpre_val  = 0.5*(M_LNPI + lnN - M_LN2 - lnsh.val) - ln_abslam;
 416       lnpre_err  = lnsh.err;
 417       lnpre_err += 2.0 * GSL_DBL_EPSILON * (0.5*(M_LNPI + M_LN2 + fabs(lnN)) + fabs(ln_abslam));
 418       lnpre_err += 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
 419       stat_e = gsl_sf_exp_mult_err_e(lnpre_val + lm, lnpre_err, P.val, P.err, result);
 420       return GSL_ERROR_SELECT_2(stat_e, stat_P);
 421     }
 422   }
 423   else if(abs_lam > 1000.0*ell*ell) {
 424     /* Large degree.
 425      */
 426     gsl_sf_result P;
 427     double lm;
 428     int stat_P = gsl_sf_conicalP_xgt1_neg_mu_largetau_e(ell+0.5,
 429                                                            lambda,
 430                                                            cosh_eta, eta,
 431                                                            &P, &lm);
 432     if(P.val == 0.0) {
 433       result->val = 0.0;
 434       result->err = 0.0;
 435       return stat_P;
 436     }
 437     else {
 438       double lnN;
 439       gsl_sf_result lnsh;
 440       double ln_abslam;
 441       double lnpre_val, lnpre_err;
 442       int stat_e;
 443       gsl_sf_lnsinh_e(eta, &lnsh);
 444       legendre_H3d_lnnorm(ell, lambda, &lnN);
 445       ln_abslam = log(abs_lam);
 446       lnpre_val  = 0.5*(M_LNPI + lnN - M_LN2 - lnsh.val) - ln_abslam;
 447       lnpre_err  = lnsh.err;
 448       lnpre_err += GSL_DBL_EPSILON * (0.5*(M_LNPI + M_LN2 + fabs(lnN)) + fabs(ln_abslam));
 449       lnpre_err += 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val);
 450       stat_e = gsl_sf_exp_mult_err_e(lnpre_val + lm, lnpre_err, P.val, P.err, result);
 451       return GSL_ERROR_SELECT_2(stat_e, stat_P);
 452     }
 453   }
 454   else {
 455     /* Backward recurrence.
 456      */
 457     const double coth_eta = 1.0/tanh(eta);
 458     const double coth_err_mult = fabs(eta) + 1.0;
 459     gsl_sf_result rH;
 460     int stat_CF1 = legendre_H3d_CF1_ser(ell, lambda, coth_eta, &rH);
 461     double Hlm1;
 462     double Hl    = GSL_SQRT_DBL_MIN;
 463     double Hlp1  = rH.val * Hl;
 464     int lp;
 465     for(lp=ell; lp>0; lp--) {
 466       double root_term_0 = hypot(lambda,lp);
 467       double root_term_1 = hypot(lambda,lp+1.0);
 468       Hlm1 = ((2.0*lp + 1.0)*coth_eta*Hl - root_term_1 * Hlp1)/root_term_0;
 469       Hlp1 = Hl;
 470       Hl   = Hlm1;
 471     }
 472
 473     if(fabs(Hl) > fabs(Hlp1)) {
 474       gsl_sf_result H0;
 475       int stat_H0 = gsl_sf_legendre_H3d_0_e(lambda, eta, &H0);
 476       result->val  = GSL_SQRT_DBL_MIN/Hl * H0.val;
 477       result->err  = GSL_SQRT_DBL_MIN/fabs(Hl) * H0.err;
 478       result->err += fabs(rH.err/rH.val) * (ell+1.0) * coth_err_mult * fabs(result->val);
 479       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 480       return GSL_ERROR_SELECT_2(stat_H0, stat_CF1);
 481     }
 482     else {
 483       gsl_sf_result H1;
 484       int stat_H1 = gsl_sf_legendre_H3d_1_e(lambda, eta, &H1);
 485       result->val  = GSL_SQRT_DBL_MIN/Hlp1 * H1.val;
 486       result->err  = GSL_SQRT_DBL_MIN/fabs(Hlp1) * H1.err;
 487       result->err += fabs(rH.err/rH.val) * (ell+1.0) * coth_err_mult * fabs(result->val);
 488       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 489       return GSL_ERROR_SELECT_2(stat_H1, stat_CF1);
 490     }
 491   }
 492 }
 493
 494
 495 int
 496 gsl_sf_legendre_H3d_array(const int lmax, const double lambda, const double eta, double * result_array)
 497 {
 498   /* CHECK_POINTER(result_array) */
 499
 500  if(eta < 0.0 || lmax < 0) {
 501     int ell;
 502     for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;
 503     GSL_ERROR ("domain error", GSL_EDOM);
 504   }
 505   else if(eta > GSL_LOG_DBL_MAX) {
 506     /* cosh(eta) is too big. */
 507     int ell;
 508     for(ell=0; ell<=lmax; ell++) result_array[ell] = 0.0;
 509     GSL_ERROR ("overflow", GSL_EOVRFLW);
 510   }
 511   else if(lmax == 0) {
 512     gsl_sf_result H0;
 513     int stat = gsl_sf_legendre_H3d_e(0, lambda, eta, &H0);
 514     result_array[0] = H0.val;
 515     return stat;
 516   }
 517   else {
 518     /* Not the most efficient method. But what the hell... it's simple.
 519      */
 520     gsl_sf_result r_Hlp1;
 521     gsl_sf_result r_Hl;
 522     int stat_lmax   = gsl_sf_legendre_H3d_e(lmax,   lambda, eta, &r_Hlp1);
 523     int stat_lmaxm1 = gsl_sf_legendre_H3d_e(lmax-1, lambda, eta, &r_Hl);
 524     int stat_max = GSL_ERROR_SELECT_2(stat_lmax, stat_lmaxm1);
 525
 526     const double coth_eta = 1.0/tanh(eta);
 527     int stat_recursion = GSL_SUCCESS;
 528     double Hlp1 = r_Hlp1.val;
 529     double Hl   = r_Hl.val;
 530     double Hlm1;
 531     int ell;
 532
 533     result_array[lmax]   = Hlp1;
 534     result_array[lmax-1] = Hl;
 535
 536     for(ell=lmax-1; ell>0; ell--) {
 537       double root_term_0 = hypot(lambda,ell);
 538       double root_term_1 = hypot(lambda,ell+1.0);
 539       Hlm1 = ((2.0*ell + 1.0)*coth_eta*Hl - root_term_1 * Hlp1)/root_term_0;
 540       result_array[ell-1] = Hlm1;
 541       if(!(Hlm1 < GSL_DBL_MAX)) stat_recursion = GSL_EOVRFLW;
 542       Hlp1 = Hl;
 543       Hl   = Hlm1;
 544     }
 545
 546     return GSL_ERROR_SELECT_2(stat_recursion, stat_max);
 547   }
 548 }
 549
 550
 551 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 552
 553 #include "eval.h"
 554
 555 double gsl_sf_legendre_H3d_0(const double lambda, const double eta)
 556 {
 557   EVAL_RESULT(gsl_sf_legendre_H3d_0_e(lambda, eta, &result));
 558 }
 559
 560 double gsl_sf_legendre_H3d_1(const double lambda, const double eta)
 561 {
 562   EVAL_RESULT(gsl_sf_legendre_H3d_1_e(lambda, eta, &result));
 563 }
 564
 565 double gsl_sf_legendre_H3d(const int l, const double lambda, const double eta)
 566 {
 567   EVAL_RESULT(gsl_sf_legendre_H3d_e(l, lambda, eta, &result));
 568 }