13 const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424, p4 = 1.2067492,
14 p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
16 const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
17 q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
18 q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
20 const double k1 = 3.75;
23 double y = 0, result = 0;
28 result = p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7)))));
46 +
y * (q8 +
y * q9))))))));
61 const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
62 p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
64 const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
65 q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3,
73 double y = 0, result = 0;
79 + (p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7))))));
84 * (q1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * (q5 +
y * (q6 +
y * q7))))));
99 const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869, p4 = 0.15084934,
100 p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
102 const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
103 q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
104 q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
106 const double k1 = 3.75;
109 double y = 0, result = 0;
115 x * (p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7))))));
133 +
y * (q8 +
y * q9))))))));
134 if (
x < 0) result = -result;
149 const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579, p4 = -0.18156897,
150 p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
152 const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
153 q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3,
161 double y = 0, result = 0;
169 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7))))));
174 * (q1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * (q5 +
y * (q6 +
y * q7))))));
186 if (
x <= 0 || n < 0) {
199 for (
int j = 1; j < n; j++) {
200 bkp = bkm + double(j) * tox * bk;
215 const double kBigPositive = 1.e10;
216 const double kBigNegative = 1.e-10;
226 if (
x == 0)
return 0;
227 if (
fabs(
x) > kBigPositive)
return 0;
229 double tox = 2 /
fabs(
x);
230 double bip = 0, bim = 0;
233 int m = 2 * ((n + int(
sqrt(
float(iacc * n)))));
234 for (
int j =
m; j >= 1; j--) {
235 bim = bip + double(j) * tox * bi;
239 if (
fabs(bi) > kBigPositive) {
240 result *= kBigNegative;
244 if (j == n) result = bip;
248 if ((
x < 0) && (n % 2 == 1)) result = -result;
258 double xx,
y, result, result1, result2;
259 const double p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
260 const double p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
261 const double p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
262 const double p10 = 59272.64853, p11 = 267.8532712;
264 const double q1 = 0.785398164;
265 const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
266 const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
267 const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
268 const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
269 const double q10 = 0.934935152e-7, q11 = 0.636619772;
271 if ((ax =
fabs(
x)) < 8) {
273 result1 = p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * p6))));
274 result2 = p7 +
y * (p8 +
y * (p9 +
y * (p10 +
y * (p11 +
y))));
275 result = result1 / result2;
280 result1 = 1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * q5)));
281 result2 = q6 +
y * (q7 +
y * (q8 +
y * (q9 -
y * q10)));
282 result =
sqrt(q11 / ax) * (
cos(xx) * result1 - z *
sin(xx) * result2);
292 double xx,
y, result, result1, result2;
293 const double p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
294 const double p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
295 const double p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
296 const double p10 = 99447.43394, p11 = 376.9991397;
298 const double q1 = 2.356194491;
299 const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
300 const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
301 const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
302 const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
303 const double q10 = 0.105787412e-6, q11 = 0.636619772;
305 if ((ax =
fabs(
x)) < 8) {
307 result1 =
x * (p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * p6)))));
308 result2 = p7 +
y * (p8 +
y * (p9 +
y * (p10 +
y * (p11 +
y))));
309 result = result1 / result2;
314 result1 = 1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * q5)));
315 result2 = q6 +
y * (q7 +
y * (q8 +
y * (q9 +
y * q10)));
316 result =
sqrt(q11 / ax) * (
cos(xx) * result1 - z *
sin(xx) * result2);
317 if (
x < 0) result = -result;
326 double z, xx,
y, result, result1, result2;
327 const double p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
328 const double p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
329 const double p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
330 const double p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
332 const double q1 = 0.785398164;
333 const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
334 const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
335 const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
336 const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
337 const double q10 = -0.934945152e-7, q11 = 0.636619772;
341 result1 = p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * p6))));
342 result2 = p7 +
y * (p8 +
y * (p9 +
y * (p10 +
y * (p11 +
y))));
348 result1 = 1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * q5)));
349 result2 = q6 +
y * (q7 +
y * (q8 +
y * (q9 +
y * q10)));
350 result =
sqrt(q11 /
x) * (
sin(xx) * result1 + z *
cos(xx) * result2);
359 double z, xx,
y, result, result1, result2;
360 const double p1 = -0.4900604943e13, p2 = 0.1275274390e13;
361 const double p3 = -0.5153438139e11, p4 = 0.7349264551e9;
362 const double p5 = -0.4237922726e7, p6 = 0.8511937935e4;
363 const double p7 = 0.2499580570e14, p8 = 0.4244419664e12;
364 const double p9 = 0.3733650367e10, p10 = 0.2245904002e8;
365 const double p11 = 0.1020426050e6, p12 = 0.3549632885e3;
366 const double p13 = 0.636619772;
367 const double q1 = 2.356194491;
368 const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
369 const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
370 const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
371 const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
372 const double q10 = 0.105787412e-6, q11 = 0.636619772;
376 result1 =
x * (p1 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * p6)))));
377 result2 = p7 +
y * (p8 +
y * (p9 +
y * (p10 +
y * (p11 +
y * (p12 +
y)))));
383 result1 = 1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * q5)));
384 result2 = q6 +
y * (q7 +
y * (q8 +
y * (q9 +
y * q10)));
385 result =
sqrt(q11 /
x) * (
sin(xx) * result1 + z *
cos(xx) * result2);
398 const double c1[16] = {1.00215845609911981,
399 -1.63969292681309147,
414 const double c2[26] = {.99283727576423943,
444 double alfa,
h, r,
y, b0, b1, b2;
455 for (
i = n1;
i >= 0; --
i) {
456 b0 = c1[
i] + alfa * b1 - b2;
460 h =
y * (b0 -
h * b2);
468 for (
i = n2;
i >= 0; --
i) {
469 b0 = c2[
i] + alfa * b1 - b2;
487 const double c3[17] = {.5578891446481605,
504 const double c4[23] = {1.00757647293865641,
532 double alfa,
h, r,
y, b0, b1, b2;
537 }
else if (
v <= 0.3) {
541 i1 =
static_cast<int>((-8. / log10(
v)));
542 for (
i = 1;
i <= i1; ++
i) {
543 h = -
h *
y / ((2 *
i + 1) * (2 *
i + 3));
553 for (
i = n3;
i >= 0; --
i) {
554 b0 = c3[
i] + alfa * b1 - b2;
560 h = 128 / (
v *
v) - 1;
565 for (
i = n4;
i >= 0; --
i) {
566 b0 = c4[
i] + alfa * b1 - b2;
586 double a0, sl0, a1, bi0;
592 for (
i = 1;
i <= 60;
i++) {
593 r *= (
x / (2 *
i + 1)) * (
x / (2 *
i + 1));
595 if (
fabs(r / s) < 1.e-12)
break;
599 km = int(5 * (
x + 1.0));
600 if (
x >= 50.0) km = 25;
601 for (
i = 1;
i <= km;
i++) {
602 r *= (2 *
i - 1) * (2 *
i - 1) /
x /
x;
604 if (
fabs(r / s) < 1.0e-12)
break;
609 for (
i = 1;
i <= 16;
i++) {
610 r = 0.125 * r * (2.0 *
i - 1.0) * (2.0 *
i - 1.0) / (
i *
x);
612 if (
fabs(r / bi0) < 1.0e-12)
break;
616 sl0 = -2.0 / (pi *
x) * s + bi0;
627 double a1, sl1, bi1, s;
633 for (
i = 1;
i <= 60;
i++) {
634 r *=
x *
x / (4.0 *
i *
i - 1.0);
636 if (
fabs(r) <
fabs(s) * 1.e-12)
break;
642 if (
x > 50.0) km = 25;
643 for (
i = 1;
i <= km;
i++) {
644 r *= (2 *
i + 3) * (2 *
i + 1) /
x /
x;
646 if (
fabs(r / s) < 1.0e-12)
break;
648 sl1 = 2.0 / pi * (-1.0 + 1.0 / (
x *
x) + 3.0 * s / (
x *
x *
x *
x));
652 for (
i = 1;
i <= 16;
i++) {
653 r = -0.125 * r * (4.0 - (2.0 *
i - 1.0) * (2.0 *
i - 1.0)) / (
i *
x);
655 if (
fabs(r / bi1) < 1.0e-12)
break;
673 const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
674 p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
676 const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
677 q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3,
685 double y = 0, result = 0;
693 +
y * (p2 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7)))))));
698 * (q1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * (q5 +
y * (q6 +
y * q7))))));
713 const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579, p4 = -0.18156897,
714 p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
716 const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
717 q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3,
725 double y = 0, result = 0;
736 +
y * (p3 +
y * (p4 +
y * (p5 +
y * (p6 +
y * p7)))))));
741 * (q1 +
y * (q2 +
y * (q3 +
y * (q4 +
y * (q5 +
y * (q6 +
y * q7))))));
753 if (
x <= 0 || n < 0) {
766 for (
int j = 1; j < n; j++) {
767 bkp = bkm + double(j) * tox * bk;