3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
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.
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.
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.
20 /* Author: G. Jungman */
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>
30 /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
33 static double locMAX3(double x, double y, double z)
35 double xy = GSL_MAX(x, y);
36 return GSL_MAX(xy, z);
40 static double locMAX4(double x, double y, double z, double w)
42 double xy = GSL_MAX(x, y);
43 double xyz = GSL_MAX(xy, z);
44 return GSL_MAX(xyz, w);
48 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
51 /* based on Carlson's algorithms:
52 [B. C. Carlson Numer. Math. 33, 1 (1979)]
55 [B.C. Carlson, Special Functions of Applied Mathematics (1977)]
58 /* According to Carlson's algorithm, the errtol parameter
59 typically effects the relative error in the following way:
61 relative error < 16 errtol^6 / (1 - 2 errtol)
74 gsl_sf_ellint_RC_e(double x, double y, gsl_mode_t mode, gsl_sf_result * result)
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];
82 if(x < 0.0 || y < 0.0 || x + y < lolim) {
85 else if(GSL_MAX(x, y) < uplim) {
86 const double c1 = 1.0 / 7.0;
87 const double c2 = 9.0 / 22.0;
90 double mu, sn, lamda, s;
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;
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);
105 DOMAIN_ERROR(result);
111 gsl_sf_ellint_RD_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
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);
119 if(GSL_MIN(x,y) < 0.0 || GSL_MIN(x+y,z) < lolim) {
120 DOMAIN_ERROR(result);
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;
132 double ea, eb, ec, ed, ef, s1, s2;
133 double mu, xndev, yndev, zndev;
135 double xnroot, ynroot, znroot, lamda;
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;
146 lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
147 sigma += power4 / (znroot * (zn + lamda));
149 xn = (xn + lamda) * 0.25;
150 yn = (yn + lamda) * 0.25;
151 zn = (zn + lamda) * 0.25;
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);
165 DOMAIN_ERROR(result);
171 gsl_sf_ellint_RF_e(double x, double y, double z, gsl_mode_t mode, gsl_sf_result * result)
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];
179 if(x < 0.0 || y < 0.0 || z < 0.0) {
180 DOMAIN_ERROR(result);
182 else if(x+y < lolim || x+z < lolim || y+z < lolim) {
183 DOMAIN_ERROR(result);
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;
192 double mu, xndev, yndev, zndev, e2, e3, s;
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;
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;
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);
218 DOMAIN_ERROR(result);
224 gsl_sf_ellint_RJ_e(double x, double y, double z, double p, gsl_mode_t mode, gsl_sf_result * result)
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);
232 if(x < 0.0 || y < 0.0 || z < 0.0) {
233 DOMAIN_ERROR(result);
235 else if(x + y < lolim || x + z < lolim || y + z < lolim || p < lolim) {
236 DOMAIN_ERROR(result);
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;
249 double mu, xndev, yndev, zndev, pndev;
250 double ea, eb, ec, e2, e3, s1, s2, s3;
252 double xnroot, ynroot, znroot;
253 double lamda, alfa, beta;
255 gsl_sf_result rcresult;
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;
267 lamda = xnroot * (ynroot + znroot) + ynroot * znroot;
268 alfa = pn * (xnroot + ynroot + znroot) + xnroot * ynroot * znroot;
270 beta = pn * (pn + lamda) * (pn + lamda);
271 rcstatus = gsl_sf_ellint_RC_e(alfa, beta, mode, &rcresult);
272 if(rcstatus != GSL_SUCCESS) {
277 sigma += power4 * rcresult.val;
279 xn = (xn + lamda) * 0.25;
280 yn = (yn + lamda) * 0.25;
281 zn = (zn + lamda) * 0.25;
282 pn = (pn + lamda) * 0.25;
285 ea = xndev * (yndev + zndev) + yndev * zndev;
286 eb = xndev * yndev * zndev;
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);
298 DOMAIN_ERROR(result);
303 /* [Carlson, Numer. Math. 33 (1979) 1, (4.1)] */
305 gsl_sf_ellint_F_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
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 */
310 double nc = floor(phi/M_PI + 0.5);
311 double phi_red = phi - nc * M_PI;
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;
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);
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);
336 /* [Carlson, Numer. Math. 33 (1979) 1, (4.2)] */
338 gsl_sf_ellint_E_e(double phi, double k, gsl_mode_t mode, gsl_sf_result * result)
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 */
343 double nc = floor(phi/M_PI + 0.5);
344 double phi_red = phi - nc * M_PI;
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;
353 if(x < GSL_DBL_EPSILON) {
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;
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);
372 return GSL_ERROR_SELECT_2(rfstatus, rdstatus);
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);
385 /* [Carlson, Numer. Math. 33 (1979) 1, (4.3)] */
387 gsl_sf_ellint_P_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
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 */
392 double nc = floor(phi/M_PI + 0.5);
393 double phi_red = phi - nc * M_PI;
396 /* FIXME: need to handle the case of small x, as for E,F */
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;
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);
414 return GSL_ERROR_SELECT_2(rfstatus, rjstatus);
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);
426 /* [Carlson, Numer. Math. 33 (1979) 1, (4.4)] */
428 gsl_sf_ellint_D_e(double phi, double k, double n, gsl_mode_t mode, gsl_sf_result * result)
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 */
433 double nc = floor(phi/M_PI + 0.5);
434 double phi_red = phi - nc * M_PI;
437 /* FIXME: need to handle the case of small x, as for E,F */
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;
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);
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);
461 gsl_sf_ellint_Dcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
464 DOMAIN_ERROR(result);
466 const double y = 1.0 - k*k; /* FIXME: still need to handle k~=~1 */
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);
476 /* [Carlson, Numer. Math. 33 (1979) 1, (4.5)] */
478 gsl_sf_ellint_Kcomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
481 DOMAIN_ERROR(result);
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;
495 /* This was previously computed as,
497 return gsl_sf_ellint_RF_e(0.0, 1.0 - k*k, 1.0, mode, result);
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) */
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;
513 /* [Carlson, Numer. Math. 33 (1979) 1, (4.6)] */
515 gsl_sf_ellint_Ecomp_e(double k, gsl_mode_t mode, gsl_sf_result * result)
518 DOMAIN_ERROR(result);
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;
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);
543 /* [Carlson, Numer. Math. 33 (1979) 1, (4.3) phi=pi/2] */
545 gsl_sf_ellint_Pcomp_e(double k, double n, gsl_mode_t mode, gsl_sf_result * result)
548 DOMAIN_ERROR(result);
550 /* FIXME: need to handle k ~=~ 1 cancellations */
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);
565 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
569 double gsl_sf_ellint_Kcomp(double k, gsl_mode_t mode)
571 EVAL_RESULT(gsl_sf_ellint_Kcomp_e(k, mode, &result));
574 double gsl_sf_ellint_Ecomp(double k, gsl_mode_t mode)
576 EVAL_RESULT(gsl_sf_ellint_Ecomp_e(k, mode, &result));
579 double gsl_sf_ellint_Pcomp(double k, double n, gsl_mode_t mode)
581 EVAL_RESULT(gsl_sf_ellint_Pcomp_e(k, n, mode, &result));
584 double gsl_sf_ellint_Dcomp(double k, gsl_mode_t mode)
586 EVAL_RESULT(gsl_sf_ellint_Dcomp_e(k, mode, &result));
589 double gsl_sf_ellint_F(double phi, double k, gsl_mode_t mode)
591 EVAL_RESULT(gsl_sf_ellint_F_e(phi, k, mode, &result));
594 double gsl_sf_ellint_E(double phi, double k, gsl_mode_t mode)
596 EVAL_RESULT(gsl_sf_ellint_E_e(phi, k, mode, &result));
599 double gsl_sf_ellint_P(double phi, double k, double n, gsl_mode_t mode)
601 EVAL_RESULT(gsl_sf_ellint_P_e(phi, k, n, mode, &result));
604 double gsl_sf_ellint_D(double phi, double k, double n, gsl_mode_t mode)
606 EVAL_RESULT(gsl_sf_ellint_D_e(phi, k, n, mode, &result));
609 double gsl_sf_ellint_RC(double x, double y, gsl_mode_t mode)
611 EVAL_RESULT(gsl_sf_ellint_RC_e(x, y, mode, &result));
614 double gsl_sf_ellint_RD(double x, double y, double z, gsl_mode_t mode)
616 EVAL_RESULT(gsl_sf_ellint_RD_e(x, y, z, mode, &result));
619 double gsl_sf_ellint_RF(double x, double y, double z, gsl_mode_t mode)
621 EVAL_RESULT(gsl_sf_ellint_RF_e(x, y, z, mode, &result));
624 double gsl_sf_ellint_RJ(double x, double y, double z, double p, gsl_mode_t mode)
626 EVAL_RESULT(gsl_sf_ellint_RJ_e(x, y, z, p, mode, &result));