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

   1 /* specfunc/ellint.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_precision.h>
  26 #include <gsl/gsl_sf_ellint.h>
  27
  28 #include "error.h"
  29
  30 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  31
  32 inline
  33 static double locMAX3(double x, double y, double z)
  34 {
  35   double xy = GSL_MAX(x, y);
  36   return GSL_MAX(xy, z);
  37 }
  38
  39 inline
  40 static double locMAX4(double x, double y, double z, double w)
  41 {
  42   double xy  = GSL_MAX(x,  y);
  43   double xyz = GSL_MAX(xy, z);
  44   return GSL_MAX(xyz, w);
  45 }
  46
  47
  48 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  49
  50
  51 /* based on Carlson's algorithms:
  52    [B. C. Carlson Numer. Math. 33, 1 (1979)]
  53
  54    see also:
  55    [B.C. Carlson, Special Functions of Applied Mathematics (1977)]
  56  */
  57
  58 /* According to Carlson's algorithm, the errtol parameter
  59    typically effects the relative error in the following way:
  60
  61    relative error < 16 errtol^6 / (1 - 2 errtol)
  62
  63      errtol     precision
  64      ------     ----------
  65      0.001       1.0e-17
  66      0.003       2.0e-14
  67      0.01        2.0e-11
  68      0.03        2.0e-8
  69      0.1         2.0e-5
  70 */
  71
  72
  73 int
  74 gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result)
  75 {
  76   const double lolim = 5.0 * GSL_DBL_MIN;
  77   const double uplim = 0.2 * GSL_DBL_MAX;
  78   const gsl_prec_t goal = GSL_MODE_PREC(mode);
  79   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
  80   const double prec   = gsl_prec_eps[goal];
  81
  82   if(x < 0.0 || y < 0.0 || x + y < lolim) {
  83     DOMAIN_ERROR(result);
  84   }
  85   else if(GSL_MAX(x, y) < uplim) {
  86     const double c1 = 1.0 / 7.0;
  87     const double c2 = 9.0 / 22.0;
  88     double xn = x;
  89     double yn = y;
  90     double mu, sn, lamda, s;
  91     while(1) {
  92       mu = (xn + yn + yn) / 3.0;
  93       sn = (yn + mu) / mu - 2.0;
  94       if (fabs(sn) < errtol) break;
  95       lamda = 2.0 * sqrt(xn) * sqrt(yn) + yn;
  96       xn = (xn + lamda) * 0.25;
  97       yn = (yn + lamda) * 0.25;
  98     }
  99     s = sn * sn * (0.3 + sn * (c1 + sn * (0.375 + sn * c2)));
 100     result->val = (1.0 + s) / sqrt(mu);
 101     result->err = prec * fabs(result->val);
 102     return GSL_SUCCESS;
 103   }
 104   else {
 105     DOMAIN_ERROR(result);
 106   }
 107 }
 108
 109
 110 int
 111 gsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
 112 {
 113   const gsl_prec_t goal = GSL_MODE_PREC(mode);
 114   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
 115   const double prec   = gsl_prec_eps[goal];
 116   const double lolim = 2.0/pow(GSL_DBL_MAX, 2.0/3.0);
 117   const double uplim = pow(0.1*errtol/GSL_DBL_MIN, 2.0/3.0);
 118
 119   if(GSL_MIN(x,y) < 0.0 || GSL_MIN(x+y,z) < lolim) {
 120     DOMAIN_ERROR(result);
 121   }
 122   else if(locMAX3(x,y,z) < uplim) {
 123     const double c1 = 3.0 / 14.0;
 124     const double c2 = 1.0 /  6.0;
 125     const double c3 = 9.0 / 22.0;
 126     const double c4 = 3.0 / 26.0;
 127     double xn = x;
 128     double yn = y;
 129     double zn = z;
 130     double sigma  = 0.0;
 131     double power4 = 1.0;
 132     double ea, eb, ec, ed, ef, s1, s2;
 133     double mu, xndev, yndev, zndev;
 134     while(1) {
 135       double xnroot, ynroot, znroot, lamda;
 136       double epslon;
 137       mu = (xn + yn + 3.0 * zn) * 0.2;
 138       xndev = (mu - xn) / mu;
 139       yndev = (mu - yn) / mu;
 140       zndev = (mu - zn) / mu;
 141       epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
 142       if (epslon < errtol) break;
 143       xnroot = sqrt(xn);
 144       ynroot = sqrt(yn);
 145       znroot = sqrt(zn);
 146       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
 147       sigma  += power4 / (znroot * (zn + lamda));
 148       power4 *= 0.25;
 149       xn = (xn + lamda) * 0.25;
 150       yn = (yn + lamda) * 0.25;
 151       zn = (zn + lamda) * 0.25;
 152     }
 153     ea = xndev * yndev;
 154     eb = zndev * zndev;
 155     ec = ea - eb;
 156     ed = ea - 6.0 * eb;
 157     ef = ed + ec + ec;
 158     s1 = ed * (- c1 + 0.25 * c3 * ed - 1.5 * c4 * zndev * ef);
 159     s2 = zndev * (c2 * ef + zndev * (- c3 * ec + zndev * c4 * ea));
 160     result->val = 3.0 * sigma + power4 * (1.0 + s1 + s2) / (mu * sqrt(mu));
 161     result->err = prec * fabs(result->val);
 162     return GSL_SUCCESS;
 163   }
 164   else {
 165     DOMAIN_ERROR(result);
 166   }
 167 }
 168
 169
 170 int
 171 gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
 172 {
 173   const double lolim = 5.0 * GSL_DBL_MIN;
 174   const double uplim = 0.2 * GSL_DBL_MAX;
 175   const gsl_prec_t goal = GSL_MODE_PREC(mode);
 176   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
 177   const double prec   = gsl_prec_eps[goal];
 178
 179   if(x < 0.0 || y < 0.0 || z < 0.0) {
 180     DOMAIN_ERROR(result);
 181   }
 182   else if(x+y < lolim || x+z < lolim || y+z < lolim) {
 183     DOMAIN_ERROR(result);
 184   }
 185   else if(locMAX3(x,y,z) < uplim) {
 186     const double c1 = 1.0 / 24.0;
 187     const double c2 = 3.0 / 44.0;
 188     const double c3 = 1.0 / 14.0;
 189     double xn = x;
 190     double yn = y;
 191     double zn = z;
 192     double mu, xndev, yndev, zndev, e2, e3, s;
 193     while(1) {
 194       double epslon, lamda;
 195       double xnroot, ynroot, znroot;
 196       mu = (xn + yn + zn) / 3.0;
 197       xndev = 2.0 - (mu + xn) / mu;
 198       yndev = 2.0 - (mu + yn) / mu;
 199       zndev = 2.0 - (mu + zn) / mu;
 200       epslon = locMAX3(fabs(xndev), fabs(yndev), fabs(zndev));
 201       if (epslon < errtol) break;
 202       xnroot = sqrt(xn);
 203       ynroot = sqrt(yn);
 204       znroot = sqrt(zn);
 205       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
 206       xn = (xn + lamda) * 0.25;
 207       yn = (yn + lamda) * 0.25;
 208       zn = (zn + lamda) * 0.25;
 209     }
 210     e2 = xndev * yndev - zndev * zndev;
 211     e3 = xndev * yndev * zndev;
 212     s = 1.0 + (c1 * e2 - 0.1 - c2 * e3) * e2 + c3 * e3;
 213     result->val = s / sqrt(mu);
 214     result->err = prec * fabs(result->val);
 215     return GSL_SUCCESS;
 216   }
 217   else {
 218     DOMAIN_ERROR(result);
 219   }
 220 }
 221
 222
 223 int
 224 gsl_sf_ellint_RJ_e(double x, double y, double z, double p, gsl_mode_t mode, gsl_sf_result * result)
 225 {
 226   const gsl_prec_t goal = GSL_MODE_PREC(mode);
 227   const double errtol = ( goal == GSL_PREC_DOUBLE ? 0.001 : 0.03 );
 228   const double prec   = gsl_prec_eps[goal];
 229   const double lolim =       pow(5.0 * GSL_DBL_MIN, 1.0/3.0);
 230   const double uplim = 0.3 * pow(0.2 * GSL_DBL_MAX, 1.0/3.0);
 231
 232   if(x < 0.0 || y < 0.0 || z < 0.0) {
 233     DOMAIN_ERROR(result);
 234   }
 235   else if(x + y < lolim || x + z < lolim || y + z < lolim || p < lolim) {
 236     DOMAIN_ERROR(result);
 237   }
 238   else if(locMAX4(x,y,z,p) < uplim) {
 239     const double c1 = 3.0 / 14.0;
 240     const double c2 = 1.0 /  3.0;
 241     const double c3 = 3.0 / 22.0;
 242     const double c4 = 3.0 / 26.0;
 243     double xn = x;
 244     double yn = y;
 245     double zn = z;
 246     double pn = p;
 247     double sigma = 0.0;
 248     double power4 = 1.0;
 249     double mu, xndev, yndev, zndev, pndev;
 250     double ea, eb, ec, e2, e3, s1, s2, s3;
 251     while(1) {
 252       double xnroot, ynroot, znroot;
 253       double lamda, alfa, beta;
 254       double epslon;
 255       gsl_sf_result rcresult;
 256       int rcstatus;
 257       mu = (xn + yn + zn + pn + pn) * 0.2;
 258       xndev = (mu - xn) / mu;
 259       yndev = (mu - yn) / mu;
 260       zndev = (mu - zn) / mu;
 261       pndev = (mu - pn) / mu;
 262       epslon = locMAX4(fabs(xndev), fabs(yndev), fabs(zndev), fabs(pndev));
 263       if(epslon < errtol) break;
 264       xnroot = sqrt(xn);
 265       ynroot = sqrt(yn);
 266       znroot = sqrt(zn);
 267       lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
 268       alfa = pn * (xnroot + ynroot + znroot) + xnroot * ynroot * znroot;
 269       alfa = alfa * alfa;
 270       beta = pn * (pn + lamda) * (pn + lamda);
 271       rcstatus = gsl_sf_ellint_RC_e(alfa, beta, mode, &rcresult);
 272       if(rcstatus != GSL_SUCCESS) {
 273         result->val = 0.0;
 274         result->err = 0.0;
 275         return rcstatus;
 276       }
 277       sigma  += power4 * rcresult.val;
 278       power4 *= 0.25;
 279       xn = (xn + lamda) * 0.25;
 280       yn = (yn + lamda) * 0.25;
 281       zn = (zn + lamda) * 0.25;
 282       pn = (pn + lamda) * 0.25;
 283     }
 284
 285     ea = xndev * (yndev + zndev) + yndev * zndev;
 286     eb = xndev * yndev * zndev;
 287     ec = pndev * pndev;
 288     e2 = ea - 3.0 * ec;
 289     e3 = eb + 2.0 * pndev * (ea - ec);
 290     s1 = 1.0 + e2 * (- c1 + 0.75 * c3 * e2 - 1.5 * c4 * e3);
 291     s2 = eb * (0.5 * c2 + pndev * (- c3 - c3 + pndev * c4));
 292     s3 = pndev * ea * (c2 - pndev * c3) - c2 * pndev * ec;
 293     result->val = 3.0 * sigma + power4 * (s1 + s2 + s3) / (mu * sqrt(mu));
 294     result->err = prec * fabs(result->val);
 295     return GSL_SUCCESS;
 296   }
 297   else {
 298     DOMAIN_ERROR(result);
 299   }
 300 }
 301
 302
 303 /* [Carlson, Numer. Math. 33 (1979) 1, (4.1)] */
 304 int
 305 gsl_sf_ellint_F_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
 306 {
 307   /* Angular reduction to -pi/2 < phi < pi/2 (we should really use an
 308      exact reduction but this will have to do for now) BJG */
 309
 310   double nc = floor(phi/M_PI + 0.5);
 311   double phi_red = phi - nc * M_PI;
 312   phi = phi_red;
 313
 314   {
 315     double sin_phi  = sin(phi);
 316     double sin2_phi = sin_phi*sin_phi;
 317     double x = 1.0 - sin2_phi;
 318     double y = 1.0 - k*k*sin2_phi;
 319     gsl_sf_result rf;
 320     int status = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
 321     result->val = sin_phi * rf.val;
 322     result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin_phi*rf.err);
 323     if (nc == 0) {
 324       return status;
 325     } else {
 326       gsl_sf_result rk;  /* add extra terms from periodicity */
 327       const int rkstatus = gsl_sf_ellint_Kcomp_e(k, mode, &rk);
 328       result->val += 2*nc*rk.val;
 329       result->err += 2*fabs(nc)*rk.err;
 330       return GSL_ERROR_SELECT_2(status, rkstatus);
 331     }
 332   }
 333 }
 334
 335
 336 /* [Carlson, Numer. Math. 33 (1979) 1, (4.2)] */
 337 int
 338 gsl_sf_ellint_E_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
 339 {
 340   /* Angular reduction to -pi/2 < phi < pi/2 (we should really use an
 341      exact reduction but this will have to do for now) BJG */
 342
 343   double nc = floor(phi/M_PI + 0.5);
 344   double phi_red = phi - nc * M_PI;
 345   phi = phi_red;
 346
 347   {
 348     const double sin_phi  = sin(phi);
 349     const double sin2_phi = sin_phi  * sin_phi;
 350     const double x = 1.0 - sin2_phi;
 351     const double y = 1.0 - k*k*sin2_phi;
 352
 353     if(x < GSL_DBL_EPSILON) {
 354       gsl_sf_result re;
 355       const int status = gsl_sf_ellint_Ecomp_e(k, mode, &re);
 356       /* could use A&S 17.4.14 to improve the value below */
 357       result->val = 2*nc*re.val + GSL_SIGN(sin_phi) * re.val;
 358       result->err = 2*fabs(nc)*re.err + re.err;
 359       return status;
 360     }
 361     else {
 362       gsl_sf_result rf, rd;
 363       const double sin3_phi = sin2_phi * sin_phi;
 364       const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
 365       const int rdstatus = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
 366       result->val  = sin_phi * rf.val - k*k/3.0 * sin3_phi * rd.val;
 367       result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
 368       result->err += fabs(sin_phi*rf.err);
 369       result->err += k*k/3.0 * GSL_DBL_EPSILON * fabs(sin3_phi * rd.val);
 370       result->err += k*k/3.0 * fabs(sin3_phi*rd.err);
 371       if (nc == 0) {
 372         return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
 373       } else {
 374         gsl_sf_result re;  /* add extra terms from periodicity */
 375         const int restatus = gsl_sf_ellint_Ecomp_e(k, mode, &re);
 376         result->val += 2*nc*re.val;
 377         result->err += 2*fabs(nc)*re.err;
 378         return GSL_ERROR_SELECT_3(rfstatus, rdstatus, restatus);
 379       }
 380     }
 381   }
 382 }
 383
 384
 385 /* [Carlson, Numer. Math. 33 (1979) 1, (4.3)] */
 386 int
 387 gsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
 388 {
 389   /* Angular reduction to -pi/2 < phi < pi/2 (we should really use an
 390      exact reduction but this will have to do for now) BJG */
 391
 392   double nc = floor(phi/M_PI + 0.5);
 393   double phi_red = phi - nc * M_PI;
 394   phi = phi_red;
 395
 396   /* FIXME: need to handle the case of small x, as for E,F */
 397
 398   {
 399     const double sin_phi  = sin(phi);
 400     const double sin2_phi = sin_phi  * sin_phi;
 401     const double sin3_phi = sin2_phi * sin_phi;
 402     const double x = 1.0 - sin2_phi;
 403     const double y = 1.0 - k*k*sin2_phi;
 404     gsl_sf_result rf;
 405     gsl_sf_result rj;
 406     const int rfstatus = gsl_sf_ellint_RF_e(x, y, 1.0, mode, &rf);
 407     const int rjstatus = gsl_sf_ellint_RJ_e(x, y, 1.0, 1.0 + n*sin2_phi, mode, &rj);
 408     result->val  = sin_phi * rf.val - n/3.0*sin3_phi * rj.val;
 409     result->err  = GSL_DBL_EPSILON * fabs(sin_phi * rf.val);
 410     result->err += fabs(sin_phi * rf.err);
 411     result->err += n/3.0 * GSL_DBL_EPSILON * fabs(sin3_phi*rj.val);
 412     result->err += n/3.0 * fabs(sin3_phi*rj.err);
 413     if (nc == 0) {
 414       return GSL_ERROR_SELECT_2(rfstatus, rjstatus);
 415     } else {
 416       gsl_sf_result rp;  /* add extra terms from periodicity */
 417       const int rpstatus = gsl_sf_ellint_Pcomp_e(k, n, mode, &rp);
 418       result->val += 2*nc*rp.val;
 419       result->err += 2*fabs(nc)*rp.err;
 420       return GSL_ERROR_SELECT_3(rfstatus, rjstatus, rpstatus);
 421     }
 422   }
 423 }
 424
 425
 426 /* [Carlson, Numer. Math. 33 (1979) 1, (4.4)] */
 427 int
 428 gsl_sf_ellint_D_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
 429 {
 430   /* Angular reduction to -pi/2 < phi < pi/2 (we should really use an
 431      exact reduction but this will have to do for now) BJG */
 432
 433   double nc = floor(phi/M_PI + 0.5);
 434   double phi_red = phi - nc * M_PI;
 435   phi = phi_red;
 436
 437   /* FIXME: need to handle the case of small x, as for E,F */
 438   {
 439     const double sin_phi  = sin(phi);
 440     const double sin2_phi = sin_phi  * sin_phi;
 441     const double sin3_phi = sin2_phi * sin_phi;
 442     const double x = 1.0 - sin2_phi;
 443     const double y = 1.0 - k*k*sin2_phi;
 444     gsl_sf_result rd;
 445     const int status = gsl_sf_ellint_RD_e(x, y, 1.0, mode, &rd);
 446     result->val = sin3_phi/3.0 * rd.val;
 447     result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(sin3_phi/3.0 * rd.err);
 448     if (nc == 0) {
 449       return status;
 450     } else {
 451       gsl_sf_result rd;  /* add extra terms from periodicity */
 452       const int rdstatus = gsl_sf_ellint_Dcomp_e(k, mode, &rd);
 453       result->val += 2*nc*rd.val;
 454       result->err += 2*fabs(nc)*rd.err;
 455       return GSL_ERROR_SELECT_2(status, rdstatus);
 456     }
 457   }
 458 }
 459
 460 int
 461 gsl_sf_ellint_Dcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
 462 {
 463   if(k*k >= 1.0) {
 464     DOMAIN_ERROR(result);
 465   } else {
 466     const double y = 1.0 - k*k;  /* FIXME: still need to handle k~=~1 */
 467     gsl_sf_result rd;
 468     const int status = gsl_sf_ellint_RD_e(0.0, y, 1.0, mode, &rd);
 469     result->val = (1.0/3.0) * rd.val;
 470     result->err = GSL_DBL_EPSILON * fabs(result->val) + fabs(1.0/3.0 * rd.err);
 471     return status;
 472   }
 473 }
 474
 475
 476 /* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
 477 int
 478 gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
 479 {
 480   if(k*k >= 1.0) {
 481     DOMAIN_ERROR(result);
 482   }
 483   else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
 484     /* [Abramowitz+Stegun, 17.3.33] */
 485     const double y = 1.0 - k*k;
 486     const double a[] = { 1.38629436112, 0.09666344259, 0.03590092383 };
 487     const double b[] = { 0.5, 0.12498593597, 0.06880248576 };
 488     const double ta = a[0] + y*(a[1] + y*a[2]);
 489     const double tb = -log(y) * (b[0] * y*(b[1] + y*b[2]));
 490     result->val = ta + tb;
 491     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
 492     return GSL_SUCCESS;
 493   }
 494   else {
 495     /* This was previously computed as,
 496
 497          return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);
 498
 499        but this underestimated the total error for small k, since the
 500        argument y=1-k^2 is not exact (there is an absolute error of
 501        GSL_DBL_EPSILON near y=0 due to cancellation in the subtraction).
 502        Taking the singular behavior of -log(y) above gives an error
 503        of 0.5*epsilon/y near y=0. (BJG) */
 504
 505     double y = 1.0 - k*k;
 506     int status = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, result);
 507     result->err += 0.5 * GSL_DBL_EPSILON / y;
 508     return status ;
 509   }
 510 }
 511
 512
 513 /* [Carlson, Numer. Math. 33 (1979) 1, (4.6)] */
 514 int
 515 gsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
 516 {
 517   if(k*k >= 1.0) {
 518     DOMAIN_ERROR(result);
 519   }
 520   else if(k*k >= 1.0 - GSL_SQRT_DBL_EPSILON) {
 521     /* [Abramowitz+Stegun, 17.3.36] */
 522     const double y = 1.0 - k*k;
 523     const double a[] = { 0.44325141463, 0.06260601220, 0.04757383546 };
 524     const double b[] = { 0.24998368310, 0.09200180037, 0.04069697526 };
 525     const double ta = 1.0 + y*(a[0] + y*(a[1] + a[2]*y));
 526     const double tb = -y*log(y) * (b[0] + y*(b[1] + b[2]*y));
 527     result->val = ta + tb;
 528     result->err = 2.0 * GSL_DBL_EPSILON * result->val;
 529     return GSL_SUCCESS;
 530   }
 531   else {
 532     gsl_sf_result rf;
 533     gsl_sf_result rd;
 534     const double y = 1.0 - k*k;
 535     const int rfstatus = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, &rf);
 536     const int rdstatus = gsl_sf_ellint_RD_e(0.0, y, 1.0, mode, &rd);
 537     result->val = rf.val - k*k/3.0 * rd.val;
 538     result->err = rf.err + k*k/3.0 * rd.err;
 539     return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
 540   }
 541 }
 542
 543 /* [Carlson, Numer. Math. 33 (1979) 1, (4.3) phi=pi/2] */
 544 int
 545 gsl_sf_ellint_Pcomp_e(double k, double n, gsl_mode_t mode, gsl_sf_result * result)
 546 {
 547   if(k*k >= 1.0) {
 548     DOMAIN_ERROR(result);
 549   }
 550   /* FIXME: need to handle k ~=~ 1  cancellations */
 551   else {
 552     gsl_sf_result rf;
 553     gsl_sf_result rj;
 554     const double y = 1.0 - k*k;
 555     const int rfstatus = gsl_sf_ellint_RF_e(0.0, y, 1.0, mode, &rf);
 556     const int rjstatus = gsl_sf_ellint_RJ_e(0.0, y, 1.0, 1.0 + n, mode, &rj);
 557     result->val = rf.val - (n/3.0) * rj.val;
 558     result->err = rf.err + fabs(n/3.0) * rj.err;
 559     return GSL_ERROR_SELECT_2(rfstatus, rjstatus);
 560   }
 561 }
 562
 563
 564
 565 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
 566
 567 #include "eval.h"
 568
 569 double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode)
 570 {
 571   EVAL_RESULT(gsl_sf_ellint_Kcomp_e(k, mode, &result));
 572 }
 573
 574 double gsl_sf_ellint_Ecomp(double k, gsl_mode_t mode)
 575 {
 576   EVAL_RESULT(gsl_sf_ellint_Ecomp_e(k, mode, &result));
 577 }
 578
 579 double gsl_sf_ellint_Pcomp(double k, double n, gsl_mode_t mode)
 580 {
 581   EVAL_RESULT(gsl_sf_ellint_Pcomp_e(k, n, mode, &result));
 582 }
 583
 584 double gsl_sf_ellint_Dcomp(double k, gsl_mode_t mode)
 585 {
 586   EVAL_RESULT(gsl_sf_ellint_Dcomp_e(k, mode, &result));
 587 }
 588
 589 double gsl_sf_ellint_F(double phi, double k, gsl_mode_t mode)
 590 {
 591   EVAL_RESULT(gsl_sf_ellint_F_e(phi, k, mode, &result));
 592 }
 593
 594 double gsl_sf_ellint_E(double phi, double k, gsl_mode_t mode)
 595 {
 596   EVAL_RESULT(gsl_sf_ellint_E_e(phi, k, mode, &result));
 597 }
 598
 599 double gsl_sf_ellint_P(double phi, double k, double n, gsl_mode_t mode)
 600 {
 601   EVAL_RESULT(gsl_sf_ellint_P_e(phi, k, n, mode, &result));
 602 }
 603
 604 double gsl_sf_ellint_D(double phi, double k, double n, gsl_mode_t mode)
 605 {
 606   EVAL_RESULT(gsl_sf_ellint_D_e(phi, k, n, mode, &result));
 607 }
 608
 609 double gsl_sf_ellint_RC(double x, double y, gsl_mode_t mode)
 610 {
 611   EVAL_RESULT(gsl_sf_ellint_RC_e(x, y, mode, &result));
 612 }
 613
 614 double gsl_sf_ellint_RD(double x, double y, double z, gsl_mode_t mode)
 615 {
 616   EVAL_RESULT(gsl_sf_ellint_RD_e(x, y, z, mode, &result));
 617 }
 618
 619 double gsl_sf_ellint_RF(double x, double y, double z, gsl_mode_t mode)
 620 {
 621   EVAL_RESULT(gsl_sf_ellint_RF_e(x, y, z, mode, &result));
 622 }
 623
 624 double gsl_sf_ellint_RJ(double x, double y, double z, double p, gsl_mode_t mode)
 625 {
 626   EVAL_RESULT(gsl_sf_ellint_RJ_e(x, y, z, p, mode, &result));
 627 }