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

   1 /* specfunc/legendre_con.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_poly.h>
  26 #include <gsl/gsl_sf_exp.h>
  27 #include <gsl/gsl_sf_trig.h>
  28 #include <gsl/gsl_sf_gamma.h>
  29 #include <gsl/gsl_sf_ellint.h>
  30 #include <gsl/gsl_sf_pow_int.h>
  31 #include <gsl/gsl_sf_bessel.h>
  32 #include <gsl/gsl_sf_hyperg.h>
  33 #include <gsl/gsl_sf_legendre.h>
  34
  35 #include "error.h"
  36 #include "legendre.h"
  37
  38 #define Root_2OverPi_  0.797884560802865355879892
  39 #define locEPS         (1000.0*GSL_DBL_EPSILON)
  40
  41
  42 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  43
  44
  45 #define RECURSE_LARGE  (1.0e-5*GSL_DBL_MAX)
  46 #define RECURSE_SMALL  (1.0e+5*GSL_DBL_MIN)
  47
  48
  49 /* Continued fraction for f_{ell+1}/f_ell
  50  * f_ell := P^{-mu-ell}_{-1/2 + I tau}(x),  x < 1.0
  51  *
  52  * Uses standard CF method from Temme's book.
  53  */
  54 static
  55 int
  56 conicalP_negmu_xlt1_CF1(const double mu, const int ell, const double tau,
  57                         const double x, gsl_sf_result * result)
  58 {
  59   const double RECUR_BIG = GSL_SQRT_DBL_MAX;
  60   const int maxiter = 5000;
  61   int n = 1;
  62   double xi = x/(sqrt(1.0-x)*sqrt(1.0+x));
  63   double Anm2 = 1.0;
  64   double Bnm2 = 0.0;
  65   double Anm1 = 0.0;
  66   double Bnm1 = 1.0;
  67   double a1 = 1.0;
  68   double b1 = 2.0*(mu + ell + 1.0) * xi;
  69   double An = b1*Anm1 + a1*Anm2;
  70   double Bn = b1*Bnm1 + a1*Bnm2;
  71   double an, bn;
  72   double fn = An/Bn;
  73
  74   while(n < maxiter) {
  75     double old_fn;
  76     double del;
  77     n++;
  78     Anm2 = Anm1;
  79     Bnm2 = Bnm1;
  80     Anm1 = An;
  81     Bnm1 = Bn;
  82     an = tau*tau + (mu - 0.5 + ell + n)*(mu - 0.5 + ell + n);
  83     bn = 2.0*(ell + mu + n) * xi;
  84     An = bn*Anm1 + an*Anm2;
  85     Bn = bn*Bnm1 + an*Bnm2;
  86
  87     if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
  88       An /= RECUR_BIG;
  89       Bn /= RECUR_BIG;
  90       Anm1 /= RECUR_BIG;
  91       Bnm1 /= RECUR_BIG;
  92       Anm2 /= RECUR_BIG;
  93       Bnm2 /= RECUR_BIG;
  94     }
  95
  96     old_fn = fn;
  97     fn = An/Bn;
  98     del = old_fn/fn;
  99
 100     if(fabs(del - 1.0) < 2.0*GSL_DBL_EPSILON) break;
 101   }
 102
 103   result->val = fn;
 104   result->err = 4.0 * GSL_DBL_EPSILON * (sqrt(n) + 1.0) * fabs(fn);
 105
 106   if(n >= maxiter)
 107     GSL_ERROR ("error", GSL_EMAXITER);
 108   else
 109     return GSL_SUCCESS;
 110 }
 111
 112
 113 /* Continued fraction for f_{ell+1}/f_ell
 114  * f_ell := P^{-mu-ell}_{-1/2 + I tau}(x),  x >= 1.0
 115  *
 116  * Uses Gautschi (Euler) equivalent series.
 117  */
 118 static
 119 int
 120 conicalP_negmu_xgt1_CF1(const double mu, const int ell, const double tau,
 121                         const double x, gsl_sf_result * result)
 122 {
 123   const int maxk = 20000;
 124   const double gamma = 1.0-1.0/(x*x);
 125   const double pre = sqrt(x-1.0)*sqrt(x+1.0) / (x*(2.0*(ell+mu+1.0)));
 126   double tk   = 1.0;
 127   double sum  = 1.0;
 128   double rhok = 0.0;
 129   int k;
 130
 131   for(k=1; k<maxk; k++) {
 132     double tlk = 2.0*(ell + mu + k);
 133     double l1k = (ell + mu - 0.5 + 1.0 + k);
 134     double ak = -(tau*tau + l1k*l1k)/(tlk*(tlk+2.0)) * gamma;
 135     rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0 + rhok));
 136     tk  *= rhok;
 137     sum += tk;
 138     if(fabs(tk/sum) < GSL_DBL_EPSILON) break;
 139   }
 140
 141   result->val  = pre * sum;
 142   result->err  = fabs(pre * tk);
 143   result->err += 2.0 * GSL_DBL_EPSILON * (sqrt(k) + 1.0) * fabs(pre*sum);
 144
 145   if(k >= maxk)
 146     GSL_ERROR ("error", GSL_EMAXITER);
 147   else
 148     return GSL_SUCCESS;
 149 }
 150
 151
 152 /* Implementation of large negative mu asymptotic
 153  * [Dunster, Proc. Roy. Soc. Edinburgh 119A, 311 (1991), p. 326]
 154  */
 155
 156 inline
 157 static double olver_U1(double beta2, double p)
 158 {
 159   return (p-1.0)/(24.0*(1.0+beta2)) * (3.0 + beta2*(2.0 + 5.0*p*(1.0+p)));
 160 }
 161
 162 inline
 163 static double olver_U2(double beta2, double p)
 164 {
 165   double beta4 = beta2*beta2;
 166   double p2    = p*p;
 167   double poly1 =  4.0*beta4 + 84.0*beta2 - 63.0;
 168   double poly2 = 16.0*beta4 + 90.0*beta2 - 81.0;
 169   double poly3 = beta2*p2*(97.0*beta2 - 432.0 + 77.0*p*(beta2-6.0) - 385.0*beta2*p2*(1.0 + p));
 170   return (1.0-p)/(1152.0*(1.0+beta2)) * (poly1 + poly2 + poly3);
 171 }
 172
 173 static const double U3c1[] = {   -1307.0,   -1647.0,    3375.0,    3675.0 };
 174 static const double U3c2[] = {   29366.0,   35835.0, -252360.0, -272630.0,
 175                                 276810.0,  290499.0 };
 176 static const double U3c3[] = {  -29748.0,   -8840.0, 1725295.0, 1767025.0,
 177                               -7313470.0, -754778.0, 6309875.0, 6480045.0 };
 178 static const double U3c4[] = {    2696.0,    -16740.0,   -524250.0,  -183975.0,
 179                               14670540.0,  14172939.0, -48206730.0, -48461985.0,
 180                               36756720.0,  37182145.0 };
 181 static const double U3c5[] = {       9136.0,      22480.0,     12760.0,
 182                                   -252480.0,   -4662165.0,   -1705341.0,
 183                                  92370135.0,   86244015.0, -263678415.0,
 184                                -260275015.0, 185910725.0,  185910725.0 };
 185
 186 #if 0
 187 static double olver_U3(double beta2, double p)
 188 {
 189   double beta4 = beta2*beta2;
 190   double beta6 = beta4*beta2;
 191   double opb2s = (1.0+beta2)*(1.0+beta2);
 192   double den   = 39813120.0 * opb2s*opb2s;
 193   double poly1 = gsl_poly_eval(U3c1, 4, p);
 194   double poly2 = gsl_poly_eval(U3c2, 6, p);
 195   double poly3 = gsl_poly_eval(U3c3, 8, p);
 196   double poly4 = gsl_poly_eval(U3c4, 10, p);
 197   double poly5 = gsl_poly_eval(U3c5, 12, p);
 198
 199   return (p-1.0)*(     1215.0*poly1 + 324.0*beta2*poly2
 200                  + 54.0*beta4*poly3 +  12.0*beta6*poly4
 201                  + beta4*beta4*poly5
 202                  ) / den;
 203 }
 204 #endif /* 0 */
 205
 206
 207 /* Large negative mu asymptotic
 208  * P^{-mu}_{-1/2 + I tau}, mu -> Inf
 209  * |x| < 1
 210  *
 211  * [Dunster, Proc. Roy. Soc. Edinburgh 119A, 311 (1991), p. 326]
 212  */
 213 int
 214 gsl_sf_conicalP_xlt1_large_neg_mu_e(double mu, double tau, double x,
 215                                        gsl_sf_result * result, double * ln_multiplier)
 216 {
 217   double beta  = tau/mu;
 218   double beta2 = beta*beta;
 219   double S     = beta * acos((1.0-beta2)/(1.0+beta2));
 220   double p     = x/sqrt(beta2*(1.0-x*x) + 1.0);
 221   gsl_sf_result lg_mup1;
 222   int lg_stat = gsl_sf_lngamma_e(mu+1.0, &lg_mup1);
 223   double ln_pre_1 =  0.5*mu*(S - log(1.0+beta2) + log((1.0-p)/(1.0+p))) - lg_mup1.val;
 224   double ln_pre_2 = -0.25 * log(1.0 + beta2*(1.0-x));
 225   double ln_pre_3 = -tau * atan(p*beta);
 226   double ln_pre = ln_pre_1 + ln_pre_2 + ln_pre_3;
 227   double sum   = 1.0 - olver_U1(beta2, p)/mu + olver_U2(beta2, p)/(mu*mu);
 228
 229   if(sum == 0.0) {
 230     result->val = 0.0;
 231     result->err = 0.0;
 232     *ln_multiplier = 0.0;
 233     return GSL_SUCCESS;
 234   }
 235   else {
 236     int stat_e = gsl_sf_exp_mult_e(ln_pre, sum, result);
 237     if(stat_e != GSL_SUCCESS) {
 238       result->val = sum;
 239       result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
 240       *ln_multiplier = ln_pre;
 241     }
 242     else {
 243       *ln_multiplier = 0.0;
 244     }
 245     return lg_stat;
 246   }
 247 }
 248
 249
 250 /* Implementation of large tau asymptotic
 251  *
 252  * A_n^{-mu}, B_n^{-mu}  [Olver, p.465, 469]
 253  */
 254
 255 inline
 256 static double olver_B0_xi(double mu, double xi)
 257 {
 258   return (1.0 - 4.0*mu*mu)/(8.0*xi) * (1.0/tanh(xi) - 1.0/xi);
 259 }
 260
 261 static double olver_A1_xi(double mu, double xi, double x)
 262 {
 263   double B = olver_B0_xi(mu, xi);
 264   double psi;
 265   if(fabs(x - 1.0) < GSL_ROOT4_DBL_EPSILON) {
 266     double y = x - 1.0;
 267     double s = -1.0/3.0 + y*(2.0/15.0 - y *(61.0/945.0 - 452.0/14175.0*y));
 268     psi = (4.0*mu*mu - 1.0)/16.0 * s;
 269   }
 270   else {
 271     psi = (4.0*mu*mu - 1.0)/16.0 * (1.0/(x*x-1.0) - 1.0/(xi*xi));
 272   }
 273   return 0.5*xi*xi*B*B + (mu+0.5)*B - psi + mu/6.0*(0.25 - mu*mu);
 274 }
 275
 276 inline
 277 static double olver_B0_th(double mu, double theta)
 278 {
 279   return -(1.0 - 4.0*mu*mu)/(8.0*theta) * (1.0/tan(theta) - 1.0/theta);
 280 }
 281
 282 static double olver_A1_th(double mu, double theta, double x)
 283 {
 284   double B = olver_B0_th(mu, theta);
 285   double psi;
 286   if(fabs(x - 1.0) < GSL_ROOT4_DBL_EPSILON) {
 287     double y = 1.0 - x;
 288     double s = -1.0/3.0 + y*(2.0/15.0 - y *(61.0/945.0 - 452.0/14175.0*y));
 289     psi = (4.0*mu*mu - 1.0)/16.0 * s;
 290   }
 291   else {
 292     psi = (4.0*mu*mu - 1.0)/16.0 * (1.0/(x*x-1.0) + 1.0/(theta*theta));
 293   }
 294   return -0.5*theta*theta*B*B + (mu+0.5)*B - psi + mu/6.0*(0.25 - mu*mu);
 295 }
 296
 297
 298 /* Large tau uniform asymptotics
 299  * P^{-mu}_{-1/2 + I tau}
 300  * 1 < x
 301  * tau -> Inf
 302  * [Olver, p. 469]
 303  */
 304 int
 305 gsl_sf_conicalP_xgt1_neg_mu_largetau_e(const double mu, const double tau,
 306                                           const double x, double acosh_x,
 307                                           gsl_sf_result * result, double * ln_multiplier)
 308 {
 309   double xi = acosh_x;
 310   double ln_xi_pre;
 311   double ln_pre;
 312   double sumA, sumB, sum;
 313   double arg;
 314   gsl_sf_result J_mup1;
 315   gsl_sf_result J_mu;
 316   double J_mum1;
 317
 318   if(xi < GSL_ROOT4_DBL_EPSILON) {
 319     ln_xi_pre = -xi*xi/6.0;           /* log(1.0 - xi*xi/6.0) */
 320   }
 321   else {
 322     gsl_sf_result lnshxi;
 323     gsl_sf_lnsinh_e(xi, &lnshxi);
 324     ln_xi_pre = log(xi) - lnshxi.val;     /* log(xi/sinh(xi) */
 325   }
 326
 327   ln_pre = 0.5*ln_xi_pre - mu*log(tau);
 328
 329   arg = tau*xi;
 330
 331   gsl_sf_bessel_Jnu_e(mu + 1.0,   arg, &J_mup1);
 332   gsl_sf_bessel_Jnu_e(mu,         arg, &J_mu);
 333   J_mum1 = -J_mup1.val + 2.0*mu/arg*J_mu.val;      /* careful of mu < 1 */
 334
 335   sumA = 1.0 - olver_A1_xi(-mu, xi, x)/(tau*tau);
 336   sumB = olver_B0_xi(-mu, xi);
 337   sum  = J_mu.val * sumA - xi/tau * J_mum1 * sumB;
 338
 339   if(sum == 0.0) {
 340     result->val = 0.0;
 341     result->err = 0.0;
 342     *ln_multiplier = 0.0;
 343     return GSL_SUCCESS;
 344   }
 345   else {
 346     int stat_e = gsl_sf_exp_mult_e(ln_pre, sum, result);
 347     if(stat_e != GSL_SUCCESS) {
 348       result->val = sum;
 349       result->err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
 350       *ln_multiplier = ln_pre;
 351     }
 352     else {
 353       *ln_multiplier = 0.0;
 354     }
 355     return GSL_SUCCESS;
 356   }
 357 }
 358
 359
 360 /* Large tau uniform asymptotics
 361  * P^{-mu}_{-1/2 + I tau}
 362  * -1 < x < 1
 363  * tau -> Inf
 364  * [Olver, p. 473]
 365  */
 366 int
 367 gsl_sf_conicalP_xlt1_neg_mu_largetau_e(const double mu, const double tau,
 368                                           const double x, const double acos_x,
 369                                           gsl_sf_result * result, double * ln_multiplier)
 370 {
 371   double theta = acos_x;
 372   double ln_th_pre;
 373   double ln_pre;
 374   double sumA, sumB, sum, sumerr;
 375   double arg;
 376   gsl_sf_result I_mup1, I_mu;
 377   double I_mum1;
 378
 379   if(theta < GSL_ROOT4_DBL_EPSILON) {
 380     ln_th_pre = theta*theta/6.0;   /* log(1.0 + theta*theta/6.0) */
 381   }
 382   else {
 383     ln_th_pre = log(theta/sin(theta));
 384   }
 385
 386   ln_pre = 0.5 * ln_th_pre - mu * log(tau);
 387
 388   arg = tau*theta;
 389   gsl_sf_bessel_Inu_e(mu + 1.0,   arg, &I_mup1);
 390   gsl_sf_bessel_Inu_e(mu,         arg, &I_mu);
 391   I_mum1 = I_mup1.val + 2.0*mu/arg * I_mu.val; /* careful of mu < 1 */
 392
 393   sumA = 1.0 - olver_A1_th(-mu, theta, x)/(tau*tau);
 394   sumB = olver_B0_th(-mu, theta);
 395   sum  = I_mu.val * sumA - theta/tau * I_mum1 * sumB;
 396   sumerr  = fabs(I_mu.err * sumA);
 397   sumerr += fabs(I_mup1.err * theta/tau * sumB);
 398   sumerr += fabs(I_mu.err   * theta/tau * sumB * 2.0 * mu/arg);
 399
 400   if(sum == 0.0) {
 401     result->val = 0.0;
 402     result->err = 0.0;
 403     *ln_multiplier = 0.0;
 404     return GSL_SUCCESS;
 405   }
 406   else {
 407     int stat_e = gsl_sf_exp_mult_e(ln_pre, sum, result);
 408     if(stat_e != GSL_SUCCESS) {
 409       result->val  = sum;
 410       result->err  = sumerr;
 411       result->err += GSL_DBL_EPSILON * fabs(sum);
 412       *ln_multiplier = ln_pre;
 413     }
 414     else {
 415       *ln_multiplier = 0.0;
 416     }
 417     return GSL_SUCCESS;
 418   }
 419 }
 420
 421
 422 /* Hypergeometric function which appears in the
 423  * large x expansion below:
 424  *
 425  *   2F1(1/4 - mu/2 - I tau/2, 3/4 - mu/2 - I tau/2, 1 - I tau, y)
 426  *
 427  * Note that for the usage below y = 1/x^2;
 428  */
 429 static
 430 int
 431 conicalP_hyperg_large_x(const double mu, const double tau, const double y,
 432                         double * reF, double * imF)
 433 {
 434   const int kmax = 1000;
 435   const double re_a = 0.25 - 0.5*mu;
 436   const double re_b = 0.75 - 0.5*mu;
 437   const double re_c = 1.0;
 438   const double im_a = -0.5*tau;
 439   const double im_b = -0.5*tau;
 440   const double im_c = -tau;
 441
 442   double re_sum = 1.0;
 443   double im_sum = 0.0;
 444   double re_term = 1.0;
 445   double im_term = 0.0;
 446   int k;
 447
 448   for(k=1; k<=kmax; k++) {
 449     double re_ak = re_a + k - 1.0;
 450     double re_bk = re_b + k - 1.0;
 451     double re_ck = re_c + k - 1.0;
 452     double im_ak = im_a;
 453     double im_bk = im_b;
 454     double im_ck = im_c;
 455     double den = re_ck*re_ck + im_ck*im_ck;
 456     double re_multiplier = ((re_ak*re_bk - im_ak*im_bk)*re_ck + im_ck*(im_ak*re_bk + re_ak*im_bk)) / den;
 457     double im_multiplier = ((im_ak*re_bk + re_ak*im_bk)*re_ck - im_ck*(re_ak*re_bk - im_ak*im_bk)) / den;
 458     double re_tmp = re_multiplier*re_term - im_multiplier*im_term;
 459     double im_tmp = im_multiplier*re_term + re_multiplier*im_term;
 460     double asum = fabs(re_sum) + fabs(im_sum);
 461     re_term = y/k * re_tmp;
 462     im_term = y/k * im_tmp;
 463     if(fabs(re_term/asum) < GSL_DBL_EPSILON && fabs(im_term/asum) < GSL_DBL_EPSILON) break;
 464     re_sum += re_term;
 465     im_sum += im_term;
 466   }
 467
 468   *reF = re_sum;
 469   *imF = im_sum;
 470
 471   if(k == kmax)
 472     GSL_ERROR ("error", GSL_EMAXITER);
 473   else
 474     return GSL_SUCCESS;
 475 }
 476
 477
 478 /* P^{mu}_{-1/2 + I tau}
 479  * x->Inf
 480  */
 481 int
 482 gsl_sf_conicalP_large_x_e(const double mu, const double tau, const double x,
 483                              gsl_sf_result * result, double * ln_multiplier)
 484 {
 485   /* 2F1 term
 486    */
 487   double y = ( x < 0.5*GSL_SQRT_DBL_MAX ? 1.0/(x*x) : 0.0 );
 488   double reF, imF;
 489   int stat_F = conicalP_hyperg_large_x(mu, tau, y, &reF, &imF);
 490
 491   /* f = Gamma(+i tau)/Gamma(1/2 - mu + i tau)
 492    * FIXME: shift so it's better for tau-> 0
 493    */
 494   gsl_sf_result lgr_num, lgth_num;
 495   gsl_sf_result lgr_den, lgth_den;
 496   int stat_gn = gsl_sf_lngamma_complex_e(0.0,tau,&lgr_num,&lgth_num);
 497   int stat_gd = gsl_sf_lngamma_complex_e(0.5-mu,tau,&lgr_den,&lgth_den);
 498
 499   double angle = lgth_num.val - lgth_den.val + atan2(imF,reF);
 500
 501   double lnx   = log(x);
 502   double lnxp1 = log(x+1.0);
 503   double lnxm1 = log(x-1.0);
 504   double lnpre_const = 0.5*M_LN2 - 0.5*M_LNPI;
 505   double lnpre_comm = (mu-0.5)*lnx - 0.5*mu*(lnxp1 + lnxm1);
 506   double lnpre_err  =   GSL_DBL_EPSILON * (0.5*M_LN2 + 0.5*M_LNPI)
 507                       + GSL_DBL_EPSILON * fabs((mu-0.5)*lnx)
 508                       + GSL_DBL_EPSILON * fabs(0.5*mu)*(fabs(lnxp1)+fabs(lnxm1));
 509
 510   /*  result = pre*|F|*|f| * cos(angle - tau * (log(x)+M_LN2))
 511    */
 512   gsl_sf_result cos_result;
 513   int stat_cos = gsl_sf_cos_e(angle + tau*(log(x) + M_LN2), &cos_result);
 514   int status = GSL_ERROR_SELECT_4(stat_cos, stat_gd, stat_gn, stat_F);
 515   if(cos_result.val == 0.0) {
 516     result->val = 0.0;
 517     result->err = 0.0;
 518     return status;
 519   }
 520   else {
 521     double lnFf_val = 0.5*log(reF*reF+imF*imF) + lgr_num.val - lgr_den.val;
 522     double lnFf_err = lgr_num.err + lgr_den.err + GSL_DBL_EPSILON * fabs(lnFf_val);
 523     double lnnoc_val = lnpre_const + lnpre_comm + lnFf_val;
 524     double lnnoc_err = lnpre_err + lnFf_err + GSL_DBL_EPSILON * fabs(lnnoc_val);
 525     int stat_e = gsl_sf_exp_mult_err_e(lnnoc_val, lnnoc_err,
 526                                           cos_result.val, cos_result.err,
 527                                           result);
 528     if(stat_e == GSL_SUCCESS) {
 529       *ln_multiplier = 0.0;
 530     }
 531     else {
 532       result->val  = cos_result.val;
 533       result->err  = cos_result.err;
 534       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 535       *ln_multiplier = lnnoc_val;
 536     }
 537     return status;
 538   }
 539 }
 540
 541
 542 /* P^{mu}_{-1/2 + I tau}  first hypergeometric representation
 543  * -1 < x < 1
 544  * This is more effective for |x| small, however it will work w/o
 545  * reservation for any x < 0 because everything is positive
 546  * definite in that case.
 547  *
 548  * [Kolbig,   (3)] (note typo in args of gamma functions)
 549  * [Bateman, (22)] (correct form)
 550  */
 551 static
 552 int
 553 conicalP_xlt1_hyperg_A(double mu, double tau, double x, gsl_sf_result * result)
 554 {
 555   double x2 = x*x;
 556   double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
 557   double pre_val = M_SQRTPI / pow(0.5*sqrt(1-x2), mu);
 558   double pre_err = err_amp * GSL_DBL_EPSILON * (fabs(mu)+1.0) * fabs(pre_val) ;
 559   gsl_sf_result ln_g1, ln_g2, arg_g1, arg_g2;
 560   gsl_sf_result F1, F2;
 561   gsl_sf_result pre1, pre2;
 562   double t1_val, t1_err;
 563   double t2_val, t2_err;
 564
 565   int stat_F1 = gsl_sf_hyperg_2F1_conj_e(0.25 - 0.5*mu, 0.5*tau, 0.5, x2, &F1);
 566   int stat_F2 = gsl_sf_hyperg_2F1_conj_e(0.75 - 0.5*mu, 0.5*tau, 1.5, x2, &F2);
 567   int status = GSL_ERROR_SELECT_2(stat_F1, stat_F2);
 568
 569   gsl_sf_lngamma_complex_e(0.75 - 0.5*mu, -0.5*tau, &ln_g1, &arg_g1);
 570   gsl_sf_lngamma_complex_e(0.25 - 0.5*mu, -0.5*tau, &ln_g2, &arg_g2);
 571
 572   gsl_sf_exp_err_e(-2.0*ln_g1.val, 2.0*ln_g1.err, &pre1);
 573   gsl_sf_exp_err_e(-2.0*ln_g2.val, 2.0*ln_g2.err, &pre2);
 574   pre2.val *= -2.0*x;
 575   pre2.err *=  2.0*fabs(x);
 576   pre2.err +=  GSL_DBL_EPSILON * fabs(pre2.val);
 577
 578   t1_val = pre1.val * F1.val;
 579   t1_err = fabs(pre1.val) * F1.err + pre1.err * fabs(F1.val);
 580   t2_val = pre2.val * F2.val;
 581   t2_err = fabs(pre2.val) * F2.err + pre2.err * fabs(F2.val);
 582
 583   result->val  = pre_val * (t1_val + t2_val);
 584   result->err  = pre_val * (t1_err + t2_err);
 585   result->err += pre_err * fabs(t1_val + t2_val);
 586   result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 587
 588   return status;
 589 }
 590
 591
 592 /* P^{mu}_{-1/2 + I tau}
 593  * defining hypergeometric representation
 594  * [Abramowitz+Stegun, 8.1.2]
 595  * 1 < x < 3
 596  * effective for x near 1
 597  *
 598  */
 599 #if 0
 600 static
 601 int
 602 conicalP_def_hyperg(double mu, double tau, double x, double * result)
 603 {
 604   double F;
 605   int stat_F = gsl_sf_hyperg_2F1_conj_renorm_e(0.5, tau, 1.0-mu, 0.5*(1.0-x), &F);
 606   *result = pow((x+1.0)/(x-1.0), 0.5*mu) * F;
 607   return stat_F;
 608 }
 609 #endif /* 0 */
 610
 611
 612 /* P^{mu}_{-1/2 + I tau}  second hypergeometric representation
 613  * [Zhurina+Karmazina, (3.1)]
 614  * -1 < x < 3
 615  * effective for x near 1
 616  *
 617  */
 618 #if 0
 619 static
 620 int
 621 conicalP_xnear1_hyperg_C(double mu, double tau, double x, double * result)
 622 {
 623   double ln_pre, arg_pre;
 624   double ln_g1, arg_g1;
 625   double ln_g2, arg_g2;
 626   double F;
 627
 628   int stat_F = gsl_sf_hyperg_2F1_conj_renorm_e(0.5+mu, tau, 1.0+mu, 0.5*(1.0-x), &F);
 629
 630   gsl_sf_lngamma_complex_e(0.5+mu, tau, &ln_g1, &arg_g1);
 631   gsl_sf_lngamma_complex_e(0.5-mu, tau, &ln_g2, &arg_g2);
 632
 633   ln_pre  = mu*M_LN2 - 0.5*mu*log(fabs(x*x-1.0)) + ln_g1 - ln_g2;
 634   arg_pre = arg_g1 - arg_g2;
 635
 636   *result = exp(ln_pre) * F;
 637   return stat_F;
 638 }
 639 #endif /* 0 */
 640
 641
 642 /* V0, V1 from Kolbig, m = 0
 643  */
 644 static
 645 int
 646 conicalP_0_V(const double t, const double f, const double tau, const double sgn,
 647              double * V0, double * V1)
 648 {
 649   double C[8];
 650   double T[8];
 651   double H[8];
 652   double V[12];
 653   int i;
 654   T[0] = 1.0;
 655   H[0] = 1.0;
 656   V[0] = 1.0;
 657   for(i=1; i<=7; i++) {
 658     T[i] = T[i-1] * t;
 659     H[i] = H[i-1] * (t*f);
 660   }
 661   for(i=1; i<=11; i++) {
 662     V[i] = V[i-1] * tau;
 663   }
 664
 665   C[0] = 1.0;
 666   C[1] = (H[1]-1.0)/(8.0*T[1]);
 667   C[2] = (9.0*H[2] + 6.0*H[1] - 15.0 - sgn*8.0*T[2])/(128.0*T[2]);
 668   C[3] = 5.0*(15.0*H[3] + 27.0*H[2] + 21.0*H[1] - 63.0 - sgn*T[2]*(16.0*H[1]+24.0))/(1024.0*T[3]);
 669   C[4] = 7.0*(525.0*H[4] + 1500.0*H[3] + 2430.0*H[2] + 1980.0*H[1] - 6435.0
 670               + 192.0*T[4] - sgn*T[2]*(720.0*H[2]+1600.0*H[1]+2160.0)
 671               ) / (32768.0*T[4]);
 672   C[5] = 21.0*(2835.0*H[5] + 11025.0*H[4] + 24750.0*H[3] + 38610.0*H[2]
 673                + 32175.0*H[1] - 109395.0 + T[4]*(1984.0*H[1]+4032.0)
 674                - sgn*T[2]*(4800.0*H[3]+15120.0*H[2]+26400.0*H[1]+34320.0)
 675                ) / (262144.0*T[5]);
 676   C[6] = 11.0*(218295.0*H[6] + 1071630.0*H[5] + 3009825.0*H[4] + 6142500.0*H[3]
 677                + 9398025.0*H[2] + 7936110.0*H[1] - 27776385.0
 678                + T[4]*(254016.0*H[2]+749952.0*H[1]+1100736.0)
 679                - sgn*T[2]*(441000.0*H[4] + 1814400.0*H[3] + 4127760.0*H[2]
 680                          + 6552000.0*H[1] + 8353800.0 + 31232.0*T[4]
 681                          )
 682                ) / (4194304.0*T[6]);
 683
 684   *V0 = C[0] + (-4.0*C[3]/T[1]+C[4])/V[4]
 685              + (-192.0*C[5]/T[3]+144.0*C[6]/T[2])/V[8]
 686              + sgn * (-C[2]/V[2]
 687                       + (-24.0*C[4]/T[2]+12.0*C[5]/T[1]-C[6])/V[6]
 688                       + (-1920.0*C[6]/T[4])/V[10]
 689                       );
 690   *V1 = C[1]/V[1] + (8.0*(C[3]/T[2]-C[4]/T[1])+C[5])/V[5]
 691                   + (384.0*C[5]/T[4] - 768.0*C[6]/T[3])/V[9]
 692                   + sgn * ((2.0*C[2]/T[1]-C[3])/V[3]
 693                            + (48.0*C[4]/T[3]-72.0*C[5]/T[2] + 18.0*C[6]/T[1])/V[7]
 694                            + (3840.0*C[6]/T[5])/V[11]
 695                            );
 696
 697   return GSL_SUCCESS;
 698 }
 699
 700
 701 /* V0, V1 from Kolbig, m = 1
 702  */
 703 static
 704 int
 705 conicalP_1_V(const double t, const double f, const double tau, const double sgn,
 706              double * V0, double * V1)
 707 {
 708   double Cm1;
 709   double C[8];
 710   double T[8];
 711   double H[8];
 712   double V[12];
 713   int i;
 714   T[0] = 1.0;
 715   H[0] = 1.0;
 716   V[0] = 1.0;
 717   for(i=1; i<=7; i++) {
 718     T[i] = T[i-1] * t;
 719     H[i] = H[i-1] * (t*f);
 720   }
 721   for(i=1; i<=11; i++) {
 722     V[i] = V[i-1] * tau;
 723   }
 724
 725   Cm1  = -1.0;
 726   C[0] = 3.0*(1.0-H[1])/(8.0*T[1]);
 727   C[1] = (-15.0*H[2]+6.0*H[1]+9.0+sgn*8.0*T[2])/(128.0*T[2]);
 728   C[2] = 3.0*(-35.0*H[3] - 15.0*H[2] + 15.0*H[1] + 35.0 + sgn*T[2]*(32.0*H[1]+8.0))/(1024.0*T[3]);
 729   C[3] = (-4725.0*H[4] - 6300.0*H[3] - 3150.0*H[2] + 3780.0*H[1] + 10395.0
 730           -1216.0*T[4] + sgn*T[2]*(6000.0*H[2]+5760.0*H[1]+1680.0)) / (32768.0*T[4]);
 731   C[4] = 7.0*(-10395.0*H[5] - 23625.0*H[4] - 28350.0*H[3] - 14850.0*H[2]
 732               +19305.0*H[1] + 57915.0 - T[4]*(6336.0*H[1]+6080.0)
 733               + sgn*T[2]*(16800.0*H[3] + 30000.0*H[2] + 25920.0*H[1] + 7920.0)
 734               ) / (262144.0*T[5]);
 735   C[5] = (-2837835.0*H[6] - 9168390.0*H[5] - 16372125.0*H[4] - 18918900*H[3]
 736           -10135125.0*H[2] + 13783770.0*H[1] + 43648605.0
 737           -T[4]*(3044160.0*H[2] + 5588352.0*H[1] + 4213440.0)
 738           +sgn*T[2]*(5556600.0*H[4] + 14817600.0*H[3] + 20790000.0*H[2]
 739                      + 17297280.0*H[1] + 5405400.0 + 323072.0*T[4]
 740                      )
 741           ) / (4194304.0*T[6]);
 742   C[6] = 0.0;
 743
 744   *V0 = C[0] + (-4.0*C[3]/T[1]+C[4])/V[4]
 745              + (-192.0*C[5]/T[3]+144.0*C[6]/T[2])/V[8]
 746              + sgn * (-C[2]/V[2]
 747                       + (-24.0*C[4]/T[2]+12.0*C[5]/T[1]-C[6])/V[6]
 748                       );
 749   *V1 = C[1]/V[1] + (8.0*(C[3]/T[2]-C[4]/T[1])+C[5])/V[5]
 750                   + (384.0*C[5]/T[4] - 768.0*C[6]/T[3])/V[9]
 751                   + sgn * (Cm1*V[1] + (2.0*C[2]/T[1]-C[3])/V[3]
 752                            + (48.0*C[4]/T[3]-72.0*C[5]/T[2] + 18.0*C[6]/T[1])/V[7]
 753                            );
 754
 755   return GSL_SUCCESS;
 756 }
 757
 758
 759
 760 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
 761
 762 /* P^0_{-1/2 + I lambda}
 763  */
 764 int
 765 gsl_sf_conicalP_0_e(const double lambda, const double x, gsl_sf_result * result)
 766 {
 767   /* CHECK_POINTER(result) */
 768
 769   if(x <= -1.0) {
 770     DOMAIN_ERROR(result);
 771   }
 772   else if(x == 1.0) {
 773     result->val = 1.0;
 774     result->err = 0.0;
 775     return GSL_SUCCESS;
 776   }
 777   else if(lambda == 0.0) {
 778     gsl_sf_result K;
 779     int stat_K;
 780     if(x < 1.0) {
 781       const double th = acos(x);
 782       const double s  = sin(0.5*th);
 783       stat_K = gsl_sf_ellint_Kcomp_e(s, GSL_MODE_DEFAULT, &K);
 784       result->val  = 2.0/M_PI * K.val;
 785       result->err  = 2.0/M_PI * K.err;
 786       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 787       return stat_K;
 788     }
 789     else {
 790       const double xi = acosh(x);
 791       const double c  = cosh(0.5*xi);
 792       const double t  = tanh(0.5*xi);
 793       stat_K = gsl_sf_ellint_Kcomp_e(t, GSL_MODE_DEFAULT, &K);
 794       result->val  = 2.0/M_PI / c * K.val;
 795       result->err  = 2.0/M_PI / c * K.err;
 796       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 797       return stat_K;
 798     }
 799   }
 800   else if(   (x <= 0.0 && lambda < 1000.0)
 801           || (x <  0.1 && lambda < 17.0)
 802           || (x <  0.2 && lambda < 5.0 )
 803     ) {
 804     return conicalP_xlt1_hyperg_A(0.0, lambda, x, result);
 805   }
 806   else if(   (x <= 0.2 && lambda < 17.0)
 807           || (x <= 1.5 && lambda < 20.0)
 808     ) {
 809     return gsl_sf_hyperg_2F1_conj_e(0.5, lambda, 1.0, (1.0-x)/2, result);
 810   }
 811   else if(1.5 < x && lambda < GSL_MAX(x,20.0)) {
 812     gsl_sf_result P;
 813     double lm;
 814     int stat_P = gsl_sf_conicalP_large_x_e(0.0, lambda, x,
 815                                               &P, &lm
 816                                               );
 817     int stat_e = gsl_sf_exp_mult_err_e(lm, 2.0*GSL_DBL_EPSILON * fabs(lm),
 818                                           P.val, P.err,
 819                                           result);
 820     return GSL_ERROR_SELECT_2(stat_e, stat_P);
 821   }
 822   else {
 823     double V0, V1;
 824     if(x < 1.0) {
 825       double th  = acos(x);
 826       double sth = sqrt(1.0-x*x);  /* sin(th) */
 827       gsl_sf_result I0, I1;
 828       int stat_I0 = gsl_sf_bessel_I0_scaled_e(th * lambda, &I0);
 829       int stat_I1 = gsl_sf_bessel_I1_scaled_e(th * lambda, &I1);
 830       int stat_I  = GSL_ERROR_SELECT_2(stat_I0, stat_I1);
 831       int stat_V  = conicalP_0_V(th, x/sth, lambda, -1.0, &V0, &V1);
 832       double bessterm = V0 * I0.val + V1 * I1.val;
 833       double besserr  = fabs(V0) * I0.err + fabs(V1) * I1.err;
 834       double arg1 = th*lambda;
 835       double sqts = sqrt(th/sth);
 836       int stat_e = gsl_sf_exp_mult_err_e(arg1, 4.0 * GSL_DBL_EPSILON * fabs(arg1),
 837                                             sqts * bessterm, sqts * besserr,
 838                                             result);
 839       return GSL_ERROR_SELECT_3(stat_e, stat_V, stat_I);
 840     }
 841     else {
 842       double sh = sqrt(x-1.0)*sqrt(x+1.0);  /* sinh(xi)      */
 843       double xi = log(x + sh);              /* xi = acosh(x) */
 844       gsl_sf_result J0, J1;
 845       int stat_J0 = gsl_sf_bessel_J0_e(xi * lambda, &J0);
 846       int stat_J1 = gsl_sf_bessel_J1_e(xi * lambda, &J1);
 847       int stat_J  = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
 848       int stat_V  = conicalP_0_V(xi, x/sh, lambda, 1.0, &V0, &V1);
 849       double bessterm = V0 * J0.val + V1 * J1.val;
 850       double besserr  = fabs(V0) * J0.err + fabs(V1) * J1.err;
 851       double pre_val = sqrt(xi/sh);
 852       double pre_err = 2.0 * fabs(pre_val);
 853       result->val  = pre_val * bessterm;
 854       result->err  = pre_val * besserr;
 855       result->err += pre_err * fabs(bessterm);
 856       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 857       return GSL_ERROR_SELECT_2(stat_V, stat_J);
 858     }
 859   }
 860 }
 861
 862
 863 /* P^1_{-1/2 + I lambda}
 864  */
 865 int
 866 gsl_sf_conicalP_1_e(const double lambda, const double x, gsl_sf_result * result)
 867 {
 868   /* CHECK_POINTER(result) */
 869
 870   if(x <= -1.0) {
 871     DOMAIN_ERROR(result);
 872   }
 873   else if(lambda == 0.0) {
 874     gsl_sf_result K, E;
 875     int stat_K, stat_E;
 876     if(x == 1.0) {
 877       result->val = 0.0;
 878       result->err = 0.0;
 879       return GSL_SUCCESS;
 880     }
 881     else if(x < 1.0) {
 882       if(1.0-x < GSL_SQRT_DBL_EPSILON) {
 883         double err_amp = GSL_MAX_DBL(1.0, 1.0/(GSL_DBL_EPSILON + fabs(1.0-x)));
 884         result->val = 0.25/M_SQRT2 * sqrt(1.0-x) * (1.0 + 5.0/16.0 * (1.0-x));
 885         result->err = err_amp * 3.0 * GSL_DBL_EPSILON * fabs(result->val);
 886         return GSL_SUCCESS;
 887       }
 888       else {
 889         const double th = acos(x);
 890         const double s  = sin(0.5*th);
 891         const double c2 = 1.0 - s*s;
 892         const double sth = sin(th);
 893         const double pre = 2.0/(M_PI*sth);
 894         stat_K = gsl_sf_ellint_Kcomp_e(s, GSL_MODE_DEFAULT, &K);
 895         stat_E = gsl_sf_ellint_Ecomp_e(s, GSL_MODE_DEFAULT, &E);
 896         result->val  = pre * (E.val - c2 * K.val);
 897         result->err  = pre * (E.err + fabs(c2) * K.err);
 898         result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 899         return stat_K;
 900       }
 901     }
 902     else {
 903       if(x-1.0 < GSL_SQRT_DBL_EPSILON) {
 904         double err_amp = GSL_MAX_DBL(1.0, 1.0/(GSL_DBL_EPSILON + fabs(1.0-x)));
 905         result->val = -0.25/M_SQRT2 * sqrt(x-1.0) * (1.0 - 5.0/16.0 * (x-1.0));
 906         result->err = err_amp * 3.0 * GSL_DBL_EPSILON * fabs(result->val);
 907         return GSL_SUCCESS;
 908       }
 909       else {
 910         const double xi = acosh(x);
 911         const double c  = cosh(0.5*xi);
 912         const double t  = tanh(0.5*xi);
 913         const double sxi = sinh(xi);
 914         const double pre = 2.0/(M_PI*sxi) * c;
 915         stat_K = gsl_sf_ellint_Kcomp_e(t, GSL_MODE_DEFAULT, &K);
 916         stat_E = gsl_sf_ellint_Ecomp_e(t, GSL_MODE_DEFAULT, &E);
 917         result->val  = pre * (E.val - K.val);
 918         result->err  = pre * (E.err + K.err);
 919         result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 920         return stat_K;
 921       }
 922     }
 923   }
 924   else if(   (x <= 0.0 && lambda < 1000.0)
 925           || (x <  0.1 && lambda < 17.0)
 926           || (x <  0.2 && lambda < 5.0 )
 927     ) {
 928     return conicalP_xlt1_hyperg_A(1.0, lambda, x, result);
 929   }
 930   else if(   (x <= 0.2 && lambda < 17.0)
 931           || (x <  1.5 && lambda < 20.0)
 932     ) {
 933     const double arg = fabs(x*x - 1.0);
 934     const double sgn = GSL_SIGN(1.0 - x);
 935     const double pre = 0.5*(lambda*lambda + 0.25) * sgn * sqrt(arg);
 936     gsl_sf_result F;
 937     int stat_F = gsl_sf_hyperg_2F1_conj_e(1.5, lambda, 2.0, (1.0-x)/2, &F);
 938     result->val  = pre * F.val;
 939     result->err  = fabs(pre) * F.err;
 940     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
 941     return stat_F;
 942   }
 943   else if(1.5 <= x && lambda < GSL_MAX(x,20.0)) {
 944     gsl_sf_result P;
 945     double lm;
 946     int stat_P = gsl_sf_conicalP_large_x_e(1.0, lambda, x,
 947                                               &P, &lm
 948                                               );
 949     int stat_e = gsl_sf_exp_mult_err_e(lm, 2.0 * GSL_DBL_EPSILON * fabs(lm),
 950                                           P.val, P.err,
 951                                           result);
 952     return GSL_ERROR_SELECT_2(stat_e, stat_P);
 953   }
 954   else {
 955     double V0, V1;
 956     if(x < 1.0) {
 957       const double sqrt_1mx = sqrt(1.0 - x);
 958       const double sqrt_1px = sqrt(1.0 + x);
 959       const double th  = acos(x);
 960       const double sth = sqrt_1mx * sqrt_1px;  /* sin(th) */
 961       gsl_sf_result I0, I1;
 962       int stat_I0 = gsl_sf_bessel_I0_scaled_e(th * lambda, &I0);
 963       int stat_I1 = gsl_sf_bessel_I1_scaled_e(th * lambda, &I1);
 964       int stat_I  = GSL_ERROR_SELECT_2(stat_I0, stat_I1);
 965       int stat_V  = conicalP_1_V(th, x/sth, lambda, -1.0, &V0, &V1);
 966       double bessterm = V0 * I0.val + V1 * I1.val;
 967       double besserr  =  fabs(V0) * I0.err + fabs(V1) * I1.err
 968                        + 2.0 * GSL_DBL_EPSILON * fabs(V0 * I0.val)
 969                        + 2.0 * GSL_DBL_EPSILON * fabs(V1 * I1.val);
 970       double arg1 = th * lambda;
 971       double sqts = sqrt(th/sth);
 972       int stat_e = gsl_sf_exp_mult_err_e(arg1, 2.0 * GSL_DBL_EPSILON * fabs(arg1),
 973                                             sqts * bessterm, sqts * besserr,
 974                                             result);
 975       result->err *= 1.0/sqrt_1mx;
 976       return GSL_ERROR_SELECT_3(stat_e, stat_V, stat_I);
 977     }
 978     else {
 979       const double sqrt_xm1 = sqrt(x - 1.0);
 980       const double sqrt_xp1 = sqrt(x + 1.0);
 981       const double sh = sqrt_xm1 * sqrt_xp1;  /* sinh(xi)      */
 982       const double xi = log(x + sh);          /* xi = acosh(x) */
 983       const double xi_lam = xi * lambda;
 984       gsl_sf_result J0, J1;
 985       const int stat_J0 = gsl_sf_bessel_J0_e(xi_lam, &J0);
 986       const int stat_J1 = gsl_sf_bessel_J1_e(xi_lam, &J1);
 987       const int stat_J  = GSL_ERROR_SELECT_2(stat_J0, stat_J1);
 988       const int stat_V  = conicalP_1_V(xi, x/sh, lambda, 1.0, &V0, &V1);
 989       const double bessterm = V0 * J0.val + V1 * J1.val;
 990       const double besserr  = fabs(V0) * J0.err + fabs(V1) * J1.err
 991                        + 512.0 * 2.0 * GSL_DBL_EPSILON * fabs(V0 * J0.val)
 992                        + 512.0 * 2.0 * GSL_DBL_EPSILON * fabs(V1 * J1.val)
 993                        + GSL_DBL_EPSILON * fabs(xi_lam * V0 * J1.val)
 994                        + GSL_DBL_EPSILON * fabs(xi_lam * V1 * J0.val);
 995       const double pre = sqrt(xi/sh);
 996       result->val  = pre * bessterm;
 997       result->err  = pre * besserr * sqrt_xp1 / sqrt_xm1;
 998       result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
 999       return GSL_ERROR_SELECT_2(stat_V, stat_J);
1000     }
1001   }
1002 }
1003
1004
1005 /* P^{1/2}_{-1/2 + I lambda} (x)
1006  * [Abramowitz+Stegun 8.6.8, 8.6.12]
1007  * checked OK [GJ] Fri May  8 12:24:36 MDT 1998
1008  */
1009 int gsl_sf_conicalP_half_e(const double lambda, const double x,
1010                               gsl_sf_result * result
1011                               )
1012 {
1013   /* CHECK_POINTER(result) */
1014
1015   if(x <= -1.0) {
1016     DOMAIN_ERROR(result);
1017   }
1018   else if(x < 1.0) {
1019     double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
1020     double ac  = acos(x);
1021     double den = sqrt(sqrt(1.0-x)*sqrt(1.0+x));
1022     result->val  = Root_2OverPi_ / den * cosh(ac * lambda);
1023     result->err  = err_amp * 3.0 * GSL_DBL_EPSILON * fabs(result->val);
1024     result->err *= fabs(ac * lambda) + 1.0;
1025     return GSL_SUCCESS;
1026   }
1027   else if(x == 1.0) {
1028     result->val = 0.0;
1029     result->err = 0.0;
1030     return GSL_SUCCESS;
1031   }
1032   else {
1033     /* x > 1 */
1034     double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
1035     double sq_term = sqrt(x-1.0)*sqrt(x+1.0);
1036     double ln_term = log(x + sq_term);
1037     double den = sqrt(sq_term);
1038     double carg_val = lambda * ln_term;
1039     double carg_err = 2.0 * GSL_DBL_EPSILON * fabs(carg_val);
1040     gsl_sf_result cos_result;
1041     int stat_cos = gsl_sf_cos_err_e(carg_val, carg_err, &cos_result);
1042     result->val  = Root_2OverPi_ / den * cos_result.val;
1043     result->err  = err_amp * Root_2OverPi_ / den * cos_result.err;
1044     result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val);
1045     return stat_cos;
1046   }
1047 }
1048
1049
1050 /* P^{-1/2}_{-1/2 + I lambda} (x)
1051  * [Abramowitz+Stegun 8.6.9, 8.6.14]
1052  * checked OK [GJ] Fri May  8 12:24:43 MDT 1998
1053  */
1054 int gsl_sf_conicalP_mhalf_e(const double lambda, const double x, gsl_sf_result * result)
1055 {
1056   /* CHECK_POINTER(result) */
1057
1058   if(x <= -1.0) {
1059     DOMAIN_ERROR(result);
1060   }
1061   else if(x < 1.0) {
1062     double ac  = acos(x);
1063     double den = sqrt(sqrt(1.0-x)*sqrt(1.0+x));
1064     double arg = ac * lambda;
1065     double err_amp = 1.0 + 1.0/(GSL_DBL_EPSILON + fabs(1.0-fabs(x)));
1066     if(fabs(arg) < GSL_SQRT_DBL_EPSILON) {
1067       result->val  = Root_2OverPi_ / den * ac;
1068       result->err  = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1069       result->err *= err_amp;
1070     }
1071     else {
1072       result->val  = Root_2OverPi_ / (den*lambda) * sinh(arg);
1073       result->err  = GSL_DBL_EPSILON * (fabs(arg)+1.0) * fabs(result->val);
1074       result->err *= err_amp;
1075       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1076     }
1077     return GSL_SUCCESS;
1078   }
1079   else if(x == 1.0) {
1080     result->val = 0.0;
1081     result->err = 0.0;
1082     return GSL_SUCCESS;
1083   }
1084   else {
1085     /* x > 1 */
1086     double sq_term = sqrt(x-1.0)*sqrt(x+1.0);
1087     double ln_term = log(x + sq_term);
1088     double den = sqrt(sq_term);
1089     double arg_val = lambda * ln_term;
1090     double arg_err = 2.0 * GSL_DBL_EPSILON * fabs(arg_val);
1091     if(arg_val < GSL_SQRT_DBL_EPSILON) {
1092       result->val = Root_2OverPi_ / den * ln_term;
1093       result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1094       return GSL_SUCCESS;
1095     }
1096     else {
1097       gsl_sf_result sin_result;
1098       int stat_sin = gsl_sf_sin_err_e(arg_val, arg_err, &sin_result);
1099       result->val  = Root_2OverPi_ / (den*lambda) * sin_result.val;
1100       result->err  = Root_2OverPi_ / fabs(den*lambda) * sin_result.err;
1101       result->err += 3.0 * GSL_DBL_EPSILON * fabs(result->val);
1102       return stat_sin;
1103     }
1104   }
1105 }
1106
1107
1108 int gsl_sf_conicalP_sph_reg_e(const int l, const double lambda,
1109                                  const double x,
1110                                  gsl_sf_result * result
1111                                  )
1112 {
1113   /* CHECK_POINTER(result) */
1114
1115   if(x <= -1.0 || l < -1) {
1116     DOMAIN_ERROR(result);
1117   }
1118   else if(l == -1) {
1119     return gsl_sf_conicalP_half_e(lambda, x, result);
1120   }
1121   else if(l == 0) {
1122     return gsl_sf_conicalP_mhalf_e(lambda, x, result);
1123   }
1124   else if(x == 1.0) {
1125     result->val = 0.0;
1126     result->err = 0.0;
1127     return GSL_SUCCESS;
1128   }
1129   else if(x < 0.0) {
1130     double c = 1.0/sqrt(1.0-x*x);
1131     gsl_sf_result r_Pellm1;
1132     gsl_sf_result r_Pell;
1133     int stat_0 = gsl_sf_conicalP_half_e(lambda, x, &r_Pellm1);  /* P^( 1/2) */
1134     int stat_1 = gsl_sf_conicalP_mhalf_e(lambda, x, &r_Pell);   /* P^(-1/2) */
1135     int stat_P = GSL_ERROR_SELECT_2(stat_0, stat_1);
1136     double Pellm1 = r_Pellm1.val;
1137     double Pell   = r_Pell.val;
1138     double Pellp1;
1139     int ell;
1140
1141     for(ell=0; ell<l; ell++) {
1142       double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
1143       Pellp1 = (Pellm1 - (2.0*ell+1.0)*c*x * Pell) / d;
1144       Pellm1 = Pell;
1145       Pell   = Pellp1;
1146     }
1147
1148     result->val  = Pell;
1149     result->err  = (0.5*l + 1.0) * GSL_DBL_EPSILON * fabs(Pell);
1150     result->err += GSL_DBL_EPSILON * l * fabs(result->val);
1151     return stat_P;
1152   }
1153   else if(x < 1.0) {
1154     const double xi = x/(sqrt(1.0-x)*sqrt(1.0+x));
1155     gsl_sf_result rat;
1156     gsl_sf_result Phf;
1157     int stat_CF1 = conicalP_negmu_xlt1_CF1(0.5, l, lambda, x, &rat);
1158     int stat_Phf = gsl_sf_conicalP_half_e(lambda, x, &Phf);
1159     double Pellp1 = rat.val * GSL_SQRT_DBL_MIN;
1160     double Pell   = GSL_SQRT_DBL_MIN;
1161     double Pellm1;
1162     int ell;
1163
1164     for(ell=l; ell>=0; ell--) {
1165       double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
1166       Pellm1 = (2.0*ell+1.0)*xi * Pell + d * Pellp1;
1167       Pellp1 = Pell;
1168       Pell   = Pellm1;
1169     }
1170
1171     result->val  = GSL_SQRT_DBL_MIN * Phf.val / Pell;
1172     result->err  = GSL_SQRT_DBL_MIN * Phf.err / fabs(Pell);
1173     result->err += fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
1174     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1175
1176     return GSL_ERROR_SELECT_2(stat_Phf, stat_CF1);
1177   }
1178   else if(x == 1.0) {
1179     result->val = 0.0;
1180     result->err = 0.0;
1181     return GSL_SUCCESS;
1182   }
1183   else {
1184     /* x > 1.0 */
1185
1186     const double xi = x/sqrt((x-1.0)*(x+1.0));
1187     gsl_sf_result rat;
1188     int stat_CF1 = conicalP_negmu_xgt1_CF1(0.5, l, lambda, x, &rat);
1189     int stat_P;
1190     double Pellp1 = rat.val * GSL_SQRT_DBL_MIN;
1191     double Pell   = GSL_SQRT_DBL_MIN;
1192     double Pellm1;
1193     int ell;
1194
1195     for(ell=l; ell>=0; ell--) {
1196       double d = (ell+1.0)*(ell+1.0) + lambda*lambda;
1197       Pellm1 = (2.0*ell+1.0)*xi * Pell - d * Pellp1;
1198       Pellp1 = Pell;
1199       Pell   = Pellm1;
1200     }
1201
1202     if(fabs(Pell) > fabs(Pellp1)){
1203       gsl_sf_result Phf;
1204       stat_P = gsl_sf_conicalP_half_e(lambda, x, &Phf);
1205       result->val  =       GSL_SQRT_DBL_MIN * Phf.val / Pell;
1206       result->err  = 2.0 * GSL_SQRT_DBL_MIN * Phf.err / fabs(Pell);
1207       result->err += 2.0 * fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
1208       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1209     }
1210     else {
1211       gsl_sf_result Pmhf;
1212       stat_P = gsl_sf_conicalP_mhalf_e(lambda, x, &Pmhf);
1213       result->val  =       GSL_SQRT_DBL_MIN * Pmhf.val / Pellp1;
1214       result->err  = 2.0 * GSL_SQRT_DBL_MIN * Pmhf.err / fabs(Pellp1);
1215       result->err += 2.0 * fabs(rat.err/rat.val) * (l + 1.0) * fabs(result->val);
1216       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1217     }
1218
1219     return GSL_ERROR_SELECT_2(stat_P, stat_CF1);
1220   }
1221 }
1222
1223
1224 int gsl_sf_conicalP_cyl_reg_e(const int m, const double lambda,
1225                                  const double x,
1226                                  gsl_sf_result * result
1227                                  )
1228 {
1229   /* CHECK_POINTER(result) */
1230
1231   if(x <= -1.0 || m < -1) {
1232     DOMAIN_ERROR(result);
1233   }
1234   else if(m == -1) {
1235     return gsl_sf_conicalP_1_e(lambda, x, result);
1236   }
1237   else if(m == 0) {
1238     return gsl_sf_conicalP_0_e(lambda, x, result);
1239   }
1240   else if(x == 1.0) {
1241     result->val = 0.0;
1242     result->err = 0.0;
1243     return GSL_SUCCESS;
1244   }
1245   else if(x < 0.0) {
1246     double c = 1.0/sqrt(1.0-x*x);
1247     gsl_sf_result r_Pkm1;
1248     gsl_sf_result r_Pk;
1249     int stat_0 = gsl_sf_conicalP_1_e(lambda, x, &r_Pkm1);  /* P^1 */
1250     int stat_1 = gsl_sf_conicalP_0_e(lambda, x, &r_Pk);    /* P^0 */
1251     int stat_P = GSL_ERROR_SELECT_2(stat_0, stat_1);
1252     double Pkm1 = r_Pkm1.val;
1253     double Pk   = r_Pk.val;
1254     double Pkp1;
1255     int k;
1256
1257     for(k=0; k<m; k++) {
1258       double d = (k+0.5)*(k+0.5) + lambda*lambda;
1259       Pkp1 = (Pkm1 - 2.0*k*c*x * Pk) / d;
1260       Pkm1 = Pk;
1261       Pk   = Pkp1;
1262     }
1263
1264     result->val  = Pk;
1265     result->err  = (m + 2.0) * GSL_DBL_EPSILON * fabs(Pk);
1266     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1267
1268     return stat_P;
1269   }
1270   else if(x < 1.0) {
1271     const double xi = x/(sqrt(1.0-x)*sqrt(1.0+x));
1272     gsl_sf_result rat;
1273     gsl_sf_result P0;
1274     int stat_CF1 = conicalP_negmu_xlt1_CF1(0.0, m, lambda, x, &rat);
1275     int stat_P0  = gsl_sf_conicalP_0_e(lambda, x, &P0);
1276     double Pkp1 = rat.val * GSL_SQRT_DBL_MIN;
1277     double Pk   = GSL_SQRT_DBL_MIN;
1278     double Pkm1;
1279     int k;
1280
1281     for(k=m; k>0; k--) {
1282       double d = (k+0.5)*(k+0.5) + lambda*lambda;
1283       Pkm1 = 2.0*k*xi * Pk + d * Pkp1;
1284       Pkp1 = Pk;
1285       Pk   = Pkm1;
1286     }
1287
1288     result->val  = GSL_SQRT_DBL_MIN * P0.val / Pk;
1289     result->err  = 2.0 * GSL_SQRT_DBL_MIN * P0.err / fabs(Pk);
1290     result->err += 2.0 * fabs(rat.err/rat.val) * (m + 1.0) * fabs(result->val);
1291     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1292
1293     return GSL_ERROR_SELECT_2(stat_P0, stat_CF1);
1294   }
1295   else if(x == 1.0) {
1296     result->val = 0.0;
1297     result->err = 0.0;
1298     return GSL_SUCCESS;
1299   }
1300   else {
1301     /* x > 1.0 */
1302
1303     const double xi = x/sqrt((x-1.0)*(x+1.0));
1304     gsl_sf_result rat;
1305     int stat_CF1 = conicalP_negmu_xgt1_CF1(0.0, m, lambda, x, &rat);
1306     int stat_P;
1307     double Pkp1 = rat.val * GSL_SQRT_DBL_MIN;
1308     double Pk   = GSL_SQRT_DBL_MIN;
1309     double Pkm1;
1310     int k;
1311
1312     for(k=m; k>-1; k--) {
1313       double d = (k+0.5)*(k+0.5) + lambda*lambda;
1314       Pkm1 = 2.0*k*xi * Pk - d * Pkp1;
1315       Pkp1 = Pk;
1316       Pk   = Pkm1;
1317     }
1318
1319     if(fabs(Pk) > fabs(Pkp1)){
1320       gsl_sf_result P1;
1321       stat_P = gsl_sf_conicalP_1_e(lambda, x, &P1);
1322       result->val  = GSL_SQRT_DBL_MIN * P1.val / Pk;
1323       result->err  = 2.0 * GSL_SQRT_DBL_MIN * P1.err / fabs(Pk);
1324       result->err += 2.0 * fabs(rat.err/rat.val) * (m+2.0) * fabs(result->val);
1325       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1326     }
1327     else {
1328       gsl_sf_result P0;
1329       stat_P = gsl_sf_conicalP_0_e(lambda, x, &P0);
1330       result->val  = GSL_SQRT_DBL_MIN * P0.val / Pkp1;
1331       result->err  = 2.0 * GSL_SQRT_DBL_MIN * P0.err / fabs(Pkp1);
1332       result->err += 2.0 * fabs(rat.err/rat.val) * (m+2.0) * fabs(result->val);
1333       result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
1334     }
1335
1336     return GSL_ERROR_SELECT_2(stat_P, stat_CF1);
1337   }
1338 }
1339
1340
1341 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
1342
1343 #include "eval.h"
1344
1345 double gsl_sf_conicalP_0(const double lambda, const double x)
1346 {
1347   EVAL_RESULT(gsl_sf_conicalP_0_e(lambda, x, &result));
1348 }
1349
1350 double gsl_sf_conicalP_1(const double lambda, const double x)
1351 {
1352   EVAL_RESULT(gsl_sf_conicalP_1_e(lambda, x, &result));
1353 }
1354
1355 double gsl_sf_conicalP_half(const double lambda, const double x)
1356 {
1357   EVAL_RESULT(gsl_sf_conicalP_half_e(lambda, x, &result));
1358 }
1359
1360 double gsl_sf_conicalP_mhalf(const double lambda, const double x)
1361 {
1362   EVAL_RESULT(gsl_sf_conicalP_mhalf_e(lambda, x, &result));
1363 }
1364
1365 double gsl_sf_conicalP_sph_reg(const int l, const double lambda, const double x)
1366 {
1367   EVAL_RESULT(gsl_sf_conicalP_sph_reg_e(l, lambda, x, &result));
1368 }
1369
1370 double gsl_sf_conicalP_cyl_reg(const int m, const double lambda, const double x)
1371 {
1372   EVAL_RESULT(gsl_sf_conicalP_cyl_reg_e(m, lambda, x, &result));
1373 }