Thermal-FIST 1.5
Package for hadron resonance gas model applications
Loading...
Searching...
No Matches
xMath.cpp
Go to the documentation of this file.
1/*
2 * Thermal-FIST package
3 *
4 * Copyright (c) 2014-2018 Volodymyr Vovchenko
5 *
6 * GNU General Public License (GPLv3 or later)
7 */
8#include "HRGBase/xMath.h"
9
10#include <cmath>
11#include <sstream>
12#include <stdexcept>
13#include <cfloat>
14
15using namespace std;
16
17namespace thermalfist {
18
19 // Implementation of the special functions below adapted from
20 // CERN-ROOT package: https://root.cern.ch/
21
22 //______________________________________________________________________________
23 double xMath::BesselI0(double x)
24 {
25 // Compute the modified Bessel function I_0(x) for any real x.
26 //
27 //--- NvE 12-mar-2000 UU-SAP Utrecht
28
29 // Parameters of the polynomial approximation
30 const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424,
31 p4 = 1.2067492, p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
32
33 const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
34 q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
35 q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
36
37 const double k1 = 3.75;
38 double ax = fabs(x);
39
40 double y = 0, result = 0;
41
42 if (ax < k1) {
43 double xx = x / k1;
44 y = xx * xx;
45 result = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))));
46 }
47 else {
48 y = k1 / ax;
49 result = (exp(ax) / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
50 }
51 return result;
52 }
53
54 //______________________________________________________________________________
55 double xMath::BesselK0(double x)
56 {
57 // Compute the modified Bessel function K_0(x) for positive real x.
58 //
59 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
60 // Applied Mathematics Series vol. 55 (1964), Washington.
61 //
62 //--- NvE 12-mar-2000 UU-SAP Utrecht
63
64 // Parameters of the polynomial approximation
65 const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
66 p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
67
68 const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
69 q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3, q7 = 5.3208e-4;
70
71 if (x <= 0) {
72 std::stringstream os;
73 os << "Invalid input argument " << x << ". Argument must be positive.";
74 throw std::invalid_argument(os.str());
75 }
76
77 double y = 0, result = 0;
78
79 if (x <= 2) {
80 y = x * x / 4;
81 result = (-log(x / 2.)*xMath::BesselI0(x)) + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
82 }
83 else {
84 y = 2 / x;
85 result = (exp(-x) / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
86 }
87 return result;
88 }
89
90 //______________________________________________________________________________
91 double xMath::BesselI1(double x)
92 {
93 // Compute the modified Bessel function I_1(x) for any real x.
94 //
95 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
96 // Applied Mathematics Series vol. 55 (1964), Washington.
97 //
98 //--- NvE 12-mar-2000 UU-SAP Utrecht
99
100 // Parameters of the polynomial approximation
101 const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869,
102 p4 = 0.15084934, p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
103
104 const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
105 q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
106 q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
107
108 const double k1 = 3.75;
109 double ax = fabs(x);
110
111 double y = 0, result = 0;
112
113 if (ax < k1) {
114 double xx = x / k1;
115 y = xx * xx;
116 result = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
117 }
118 else {
119 y = k1 / ax;
120 result = (exp(ax) / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
121 if (x < 0) result = -result;
122 }
123 return result;
124 }
125
126 //______________________________________________________________________________
127 double xMath::BesselK1(double x)
128 {
129 // Compute the modified Bessel function K_1(x) for positive real x.
130 //
131 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
132 // Applied Mathematics Series vol. 55 (1964), Washington.
133 //
134 //--- NvE 12-mar-2000 UU-SAP Utrecht
135
136 // Parameters of the polynomial approximation
137 const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579,
138 p4 = -0.18156897, p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
139
140 const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
141 q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3, q7 = -6.8245e-4;
142 if (x <= 0) {
143 std::stringstream os;
144 os << "Invalid argument x = " << x << ". Argument must be positive.";
145 throw std::invalid_argument(os.str());
146 }
147
148 double y = 0, result = 0;
149
150 if (x <= 2) {
151 y = x * x / 4;
152 result = (log(x / 2.)*xMath::BesselI1(x)) + (1. / x)*(p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))));
153 }
154 else {
155 y = 2 / x;
156 result = (exp(-x) / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
157 }
158 return result;
159 }
160
161 //______________________________________________________________________________
162 double xMath::BesselK(int n, double x)
163 {
164 if (n < 0) n = -n;
165 // Compute the Integer Order Modified Bessel function K_n(x)
166 // for n=0,1,2,... and positive real x.
167 //
168 //--- NvE 12-mar-2000 UU-SAP Utrecht
169 if (x <= 0 || n < 0) {
170 std::stringstream os;
171 os << "Invalid input argument(s) (n,x) = (" << n << ", " << x << "). Argument x must be positive.";
172 throw std::invalid_argument(os.str());
173 }
174
175 if (n == 0) return xMath::BesselK0(x);
176 if (n == 1) return xMath::BesselK1(x);
177
178 // Perform upward recurrence for all x
179 double tox = 2 / x;
180 double bkm = xMath::BesselK0(x);
181 double bk = xMath::BesselK1(x);
182 double bkp = 0;
183 for (int j = 1; j < n; j++) {
184 bkp = bkm + double(j)*tox*bk;
185 bkm = bk;
186 bk = bkp;
187 }
188 return bk;
189 }
190
191 //______________________________________________________________________________
192 double xMath::BesselI(int n, double x)
193 {
194 // Compute the Integer Order Modified Bessel function I_n(x)
195 // for n=0,1,2,... and any real x.
196 //
197 //--- NvE 12-mar-2000 UU-SAP Utrecht
198 if (n < 0) n = -n;
199
200 int iacc = 40; // Increase to enhance accuracy
201 const double kBigPositive = 1.e10;
202 const double kBigNegative = 1.e-10;
203
204 if (n == 0) return xMath::BesselI0(x);
205 if (n == 1) return xMath::BesselI1(x);
206
207 if (x == 0) return 0;
208 if (fabs(x) > kBigPositive) return 0;
209
210 double tox = 2 / fabs(x);
211 double bip = 0, bim = 0;
212 double bi = 1;
213 double result = 0;
214 int m = 2 * ((n + int(sqrt(float(iacc*n)))));
215 for (int j = m; j >= 1; j--) {
216 bim = bip + double(j)*tox*bi;
217 bip = bi;
218 bi = bim;
219 // Renormalise to prevent overflows
220 if (fabs(bi) > kBigPositive) {
221 result *= kBigNegative;
222 bi *= kBigNegative;
223 bip *= kBigNegative;
224 }
225 if (j == n) result = bip;
226 }
227
228 result *= xMath::BesselI0(x) / bi; // Normalise with BesselI0(x)
229 if ((x < 0) && (n % 2 == 1)) result = -result;
230
231 return result;
232 }
233
234 //______________________________________________________________________________
235 double xMath::BesselJ0(double x)
236 {
237 // Returns the Bessel function J0(x) for any real x.
238
239 double ax, z;
240 double xx, y, result, result1, result2;
241 const double p1 = 57568490574.0, p2 = -13362590354.0, p3 = 651619640.7;
242 const double p4 = -11214424.18, p5 = 77392.33017, p6 = -184.9052456;
243 const double p7 = 57568490411.0, p8 = 1029532985.0, p9 = 9494680.718;
244 const double p10 = 59272.64853, p11 = 267.8532712;
245
246 const double q1 = 0.785398164;
247 const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
248 const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
249 const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
250 const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
251 const double q10 = 0.934935152e-7, q11 = 0.636619772;
252
253 if ((ax = fabs(x)) < 8) {
254 y = x * x;
255 result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
256 result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
257 result = result1 / result2;
258 }
259 else {
260 z = 8 / ax;
261 y = z * z;
262 xx = ax - q1;
263 result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
264 result2 = q6 + y * (q7 + y * (q8 + y * (q9 - y * q10)));
265 result = sqrt(q11 / ax)*(cos(xx)*result1 - z * sin(xx)*result2);
266 }
267 return result;
268 }
269
270 //______________________________________________________________________________
271 double xMath::BesselJ1(double x)
272 {
273 // Returns the Bessel function J1(x) for any real x.
274
275 double ax, z;
276 double xx, y, result, result1, result2;
277 const double p1 = 72362614232.0, p2 = -7895059235.0, p3 = 242396853.1;
278 const double p4 = -2972611.439, p5 = 15704.48260, p6 = -30.16036606;
279 const double p7 = 144725228442.0, p8 = 2300535178.0, p9 = 18583304.74;
280 const double p10 = 99447.43394, p11 = 376.9991397;
281
282 const double q1 = 2.356194491;
283 const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
284 const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
285 const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
286 const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
287 const double q10 = 0.105787412e-6, q11 = 0.636619772;
288
289 if ((ax = fabs(x)) < 8) {
290 y = x * x;
291 result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
292 result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
293 result = result1 / result2;
294 }
295 else {
296 z = 8 / ax;
297 y = z * z;
298 xx = ax - q1;
299 result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
300 result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
301 result = sqrt(q11 / ax)*(cos(xx)*result1 - z * sin(xx)*result2);
302 if (x < 0) result = -result;
303 }
304 return result;
305 }
306
307 //______________________________________________________________________________
308 double xMath::BesselY0(double x)
309 {
310 // Returns the Bessel function Y0(x) for positive x.
311
312 double z, xx, y, result, result1, result2;
313 const double p1 = -2957821389., p2 = 7062834065.0, p3 = -512359803.6;
314 const double p4 = 10879881.29, p5 = -86327.92757, p6 = 228.4622733;
315 const double p7 = 40076544269., p8 = 745249964.8, p9 = 7189466.438;
316 const double p10 = 47447.26470, p11 = 226.1030244, p12 = 0.636619772;
317
318 const double q1 = 0.785398164;
319 const double q2 = -0.1098628627e-2, q3 = 0.2734510407e-4;
320 const double q4 = -0.2073370639e-5, q5 = 0.2093887211e-6;
321 const double q6 = -0.1562499995e-1, q7 = 0.1430488765e-3;
322 const double q8 = -0.6911147651e-5, q9 = 0.7621095161e-6;
323 const double q10 = -0.934945152e-7, q11 = 0.636619772;
324
325 if (x < 8) {
326 y = x * x;
327 result1 = p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6))));
328 result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y))));
329 result = (result1 / result2) + p12 * xMath::BesselJ0(x)*log(x);
330 }
331 else {
332 z = 8 / x;
333 y = z * z;
334 xx = x - q1;
335 result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
336 result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
337 result = sqrt(q11 / x)*(sin(xx)*result1 + z * cos(xx)*result2);
338 }
339 return result;
340 }
341
342 //______________________________________________________________________________
343 double xMath::BesselY1(double x)
344 {
345 // Returns the Bessel function Y1(x) for positive x.
346
347 double z, xx, y, result, result1, result2;
348 const double p1 = -0.4900604943e13, p2 = 0.1275274390e13;
349 const double p3 = -0.5153438139e11, p4 = 0.7349264551e9;
350 const double p5 = -0.4237922726e7, p6 = 0.8511937935e4;
351 const double p7 = 0.2499580570e14, p8 = 0.4244419664e12;
352 const double p9 = 0.3733650367e10, p10 = 0.2245904002e8;
353 const double p11 = 0.1020426050e6, p12 = 0.3549632885e3;
354 const double p13 = 0.636619772;
355 const double q1 = 2.356194491;
356 const double q2 = 0.183105e-2, q3 = -0.3516396496e-4;
357 const double q4 = 0.2457520174e-5, q5 = -0.240337019e-6;
358 const double q6 = 0.04687499995, q7 = -0.2002690873e-3;
359 const double q8 = 0.8449199096e-5, q9 = -0.88228987e-6;
360 const double q10 = 0.105787412e-6, q11 = 0.636619772;
361
362 if (x < 8) {
363 y = x * x;
364 result1 = x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * p6)))));
365 result2 = p7 + y * (p8 + y * (p9 + y * (p10 + y * (p11 + y * (p12 + y)))));
366 result = (result1 / result2) + p13 * (xMath::BesselJ1(x)*log(x) - 1 / x);
367 }
368 else {
369 z = 8 / x;
370 y = z * z;
371 xx = x - q1;
372 result1 = 1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
373 result2 = q6 + y * (q7 + y * (q8 + y * (q9 + y * q10)));
374 result = sqrt(q11 / x)*(sin(xx)*result1 + z * cos(xx)*result2);
375 }
376 return result;
377 }
378
379 //______________________________________________________________________________
380 double xMath::StruveH0(double x)
381 {
382 // Struve Functions of Order 0
383 //
384 // Converted from CERNLIB M342 by Rene Brun.
385
386 const int n1 = 15;
387 const int n2 = 25;
388 const double c1[16] = { 1.00215845609911981, -1.63969292681309147,
389 1.50236939618292819, -.72485115302121872,
390 .18955327371093136, -.03067052022988,
391 .00337561447375194, -2.6965014312602e-4,
392 1.637461692612e-5, -7.8244408508e-7,
393 3.021593188e-8, -9.6326645e-10,
394 2.579337e-11, -5.8854e-13,
395 1.158e-14, -2e-16 };
396 const double c2[26] = { .99283727576423943, -.00696891281138625,
397 1.8205103787037e-4, -1.063258252844e-5,
398 9.8198294287e-7, -1.2250645445e-7,
399 1.894083312e-8, -3.44358226e-9,
400 7.1119102e-10, -1.6288744e-10,
401 4.065681e-11, -1.091505e-11,
402 3.12005e-12, -9.4202e-13,
403 2.9848e-13, -9.872e-14,
404 3.394e-14, -1.208e-14,
405 4.44e-15, -1.68e-15,
406 6.5e-16, -2.6e-16,
407 1.1e-16, -4e-17,
408 2e-17, -1e-17 };
409
410 const double c0 = 2 / xMath::Pi();
411
412 int i;
413 double alfa, h, r, y, b0, b1, b2;
414 double v = fabs(x);
415
416 v = fabs(x);
417 if (v < 8) {
418 y = v / 8;
419 h = 2 * y*y - 1;
420 alfa = h + h;
421 b0 = 0;
422 b1 = 0;
423 b2 = 0;
424 for (i = n1; i >= 0; --i) {
425 b0 = c1[i] + alfa * b1 - b2;
426 b2 = b1;
427 b1 = b0;
428 }
429 h = y * (b0 - h * b2);
430 }
431 else {
432 r = 1 / v;
433 h = 128 * r*r - 1;
434 alfa = h + h;
435 b0 = 0;
436 b1 = 0;
437 b2 = 0;
438 for (i = n2; i >= 0; --i) {
439 b0 = c2[i] + alfa * b1 - b2;
440 b2 = b1;
441 b1 = b0;
442 }
443 h = xMath::BesselY0(v) + r * c0*(b0 - h * b2);
444 }
445 if (x < 0) h = -h;
446 return h;
447 }
448
449 //______________________________________________________________________________
450 double xMath::StruveH1(double x)
451 {
452 // Struve Functions of Order 1
453 //
454 // Converted from CERNLIB M342 by Rene Brun.
455
456 const int n3 = 16;
457 const int n4 = 22;
458 const double c3[17] = { .5578891446481605, -.11188325726569816,
459 -.16337958125200939, .32256932072405902,
460 -.14581632367244242, .03292677399374035,
461 -.00460372142093573, 4.434706163314e-4,
462 -3.142099529341e-5, 1.7123719938e-6,
463 -7.416987005e-8, 2.61837671e-9,
464 -7.685839e-11, 1.9067e-12,
465 -4.052e-14, 7.5e-16,
466 -1e-17 };
467 const double c4[23] = { 1.00757647293865641, .00750316051248257,
468 -7.043933264519e-5, 2.66205393382e-6,
469 -1.8841157753e-7, 1.949014958e-8,
470 -2.6126199e-9, 4.236269e-10,
471 -7.955156e-11, 1.679973e-11,
472 -3.9072e-12, 9.8543e-13,
473 -2.6636e-13, 7.645e-14,
474 -2.313e-14, 7.33e-15,
475 -2.42e-15, 8.3e-16,
476 -3e-16, 1.1e-16,
477 -4e-17, 2e-17,-1e-17 };
478
479 const double c0 = 2 / xMath::Pi();
480 const double cc = 2 / (3 * xMath::Pi());
481
482 int i, i1;
483 double alfa, h, r, y, b0, b1, b2;
484 double v = fabs(x);
485
486 if (v == 0) {
487 h = 0;
488 }
489 else if (v <= 0.3) {
490 y = v * v;
491 r = 1;
492 h = 1;
493 i1 = (int)(-8. / log10(v));
494 for (i = 1; i <= i1; ++i) {
495 h = -h * y / ((2 * i + 1)*(2 * i + 3));
496 r += h;
497 }
498 h = cc * y*r;
499 }
500 else if (v < 8) {
501 h = v * v / 32 - 1;
502 alfa = h + h;
503 b0 = 0;
504 b1 = 0;
505 b2 = 0;
506 for (i = n3; i >= 0; --i) {
507 b0 = c3[i] + alfa * b1 - b2;
508 b2 = b1;
509 b1 = b0;
510 }
511 h = b0 - h * b2;
512 }
513 else {
514 h = 128 / (v*v) - 1;
515 alfa = h + h;
516 b0 = 0;
517 b1 = 0;
518 b2 = 0;
519 for (i = n4; i >= 0; --i) {
520 b0 = c4[i] + alfa * b1 - b2;
521 b2 = b1;
522 b1 = b0;
523 }
524 h = xMath::BesselY1(v) + c0 * (b0 - h * b2);
525 }
526 return h;
527 }
528
529
530 //______________________________________________________________________________
531 double xMath::StruveL0(double x)
532 {
533 // Modified Struve Function of Order 0.
534 // By Kirill Filimonov.
535
536 const double pi = xMath::Pi();
537
538 double s = 1.0;
539 double r = 1.0;
540
541 double a0, sl0, a1, bi0;
542
543 int km, i;
544
545 if (x <= 20.) {
546 a0 = 2.0*x / pi;
547 for (i = 1; i <= 60; i++) {
548 r *= (x / (2 * i + 1))*(x / (2 * i + 1));
549 s += r;
550 if (fabs(r / s) < 1.e-12) break;
551 }
552 sl0 = a0 * s;
553 }
554 else {
555 km = int(5 * (x + 1.0));
556 if (x >= 50.0)km = 25;
557 for (i = 1; i <= km; i++) {
558 r *= (2 * i - 1)*(2 * i - 1) / x / x;
559 s += r;
560 if (fabs(r / s) < 1.0e-12) break;
561 }
562 a1 = exp(x) / sqrt(2 * pi*x);
563 r = 1.0;
564 bi0 = 1.0;
565 for (i = 1; i <= 16; i++) {
566 r = 0.125*r*(2.0*i - 1.0)*(2.0*i - 1.0) / (i*x);
567 bi0 += r;
568 if (fabs(r / bi0) < 1.0e-12) break;
569 }
570
571 bi0 = a1 * bi0;
572 sl0 = -2.0 / (pi*x)*s + bi0;
573 }
574 return sl0;
575 }
576
577 //______________________________________________________________________________
578 double xMath::StruveL1(double x)
579 {
580 // Modified Struve Function of Order 1.
581 // By Kirill Filimonov.
582
583 const double pi = xMath::Pi();
584 double a1, sl1, bi1, s;
585 double r = 1.0;
586 int km, i;
587
588 if (x <= 20.) {
589 s = 0.0;
590 for (i = 1; i <= 60; i++) {
591 r *= x * x / (4.0*i*i - 1.0);
592 s += r;
593 if (fabs(r) < fabs(s)*1.e-12) break;
594 }
595 sl1 = 2.0 / pi * s;
596 }
597 else {
598 s = 1.0;
599 km = int(0.5*x);
600 if (x > 50.0)km = 25;
601 for (i = 1; i <= km; i++) {
602 r *= (2 * i + 3)*(2 * i + 1) / x / x;
603 s += r;
604 if (fabs(r / s) < 1.0e-12) break;
605 }
606 sl1 = 2.0 / pi * (-1.0 + 1.0 / (x*x) + 3.0*s / (x*x*x*x));
607 a1 = exp(x) / sqrt(2 * pi*x);
608 r = 1.0;
609 bi1 = 1.0;
610 for (i = 1; i <= 16; i++) {
611 r = -0.125*r*(4.0 - (2.0*i - 1.0)*(2.0*i - 1.0)) / (i*x);
612 bi1 += r;
613 if (fabs(r / bi1) < 1.0e-12) break;
614 }
615 sl1 += a1 * bi1;
616 }
617 return sl1;
618 }
619
620
621
622 //______________________________________________________________________________
623 double xMath::BesselK0exp(double x)
624 {
625 // Compute the modified Bessel function K_0(x) for positive real x.
626 //
627 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
628 // Applied Mathematics Series vol. 55 (1964), Washington.
629 //
630 //--- NvE 12-mar-2000 UU-SAP Utrecht
631
632 // Parameters of the polynomial approximation
633 const double p1 = -0.57721566, p2 = 0.42278420, p3 = 0.23069756,
634 p4 = 3.488590e-2, p5 = 2.62698e-3, p6 = 1.0750e-4, p7 = 7.4e-6;
635
636 const double q1 = 1.25331414, q2 = -7.832358e-2, q3 = 2.189568e-2,
637 q4 = -1.062446e-2, q5 = 5.87872e-3, q6 = -2.51540e-3, q7 = 5.3208e-4;
638
639 if (x <= 0) {
640 std::stringstream os;
641 os << "Invalid input argument " << x << ". Argument must be positive.";
642 throw std::invalid_argument(os.str());
643 }
644 double y = 0, result = 0;
645 if (x <= 2) {
646 y = x * x / 4;
647 result = exp(x)*((-log(x / 2.)*xMath::BesselI0(x)) + (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
648 }
649 else {
650 y = 2 / x;
651 result = (1. / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
652 }
653 return result;
654 }
655
656 //______________________________________________________________________________
657 double xMath::BesselK1exp(double x)
658 {
659 // Compute the modified Bessel function K_1(x) for positive real x.
660 //
661 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
662 // Applied Mathematics Series vol. 55 (1964), Washington.
663 //
664 //--- NvE 12-mar-2000 UU-SAP Utrecht
665
666 // Parameters of the polynomial approximation
667 const double p1 = 1., p2 = 0.15443144, p3 = -0.67278579,
668 p4 = -0.18156897, p5 = -1.919402e-2, p6 = -1.10404e-3, p7 = -4.686e-5;
669
670 const double q1 = 1.25331414, q2 = 0.23498619, q3 = -3.655620e-2,
671 q4 = 1.504268e-2, q5 = -7.80353e-3, q6 = 3.25614e-3, q7 = -6.8245e-4;
672
673 if (x <= 0) {
674 std::stringstream os;
675 os << "Invalid input argument " << x << ". Argument must be positive.";
676 throw std::invalid_argument(os.str());
677 return 0;
678 }
679
680 double y = 0, result = 0;
681 if (x <= 2) {
682 y = x * x / 4;
683 result = exp(x) * ((log(x / 2.)*xMath::BesselI1(x)) + (1. / x)*(p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
684 }
685 else {
686 y = 2 / x;
687 result = (1. / sqrt(x))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * q7))))));
688 }
689 return result;
690 }
691
692 //______________________________________________________________________________
693 double xMath::BesselKexp(int n, double x)
694 {
695 if (n < 0) n = -n;
696 // Compute the Integer Order Modified Bessel function K_n(x)
697 // for n=0,1,2,... and positive real x.
698 //
699 //--- NvE 12-mar-2000 UU-SAP Utrecht
700 if (x <= 0 || n < 0) {
701 std::stringstream os;
702 os << "Invalid input argument(s) (n,x) = (" << n << ", " << x << "). Argument x must be positive.";
703 throw std::invalid_argument(os.str());
704 }
705
706 if (n == 0) return xMath::BesselK0exp(x);
707 if (n == 1) return xMath::BesselK1exp(x);
708
709 // Perform upward recurrence for all x
710 double tox = 2 / x;
711 double bkm = xMath::BesselK0exp(x);
712 double bk = xMath::BesselK1exp(x);
713 double bkp = 0;
714 for (int j = 1; j < n; j++) {
715 bkp = bkm + double(j)*tox*bk;
716 bkm = bk;
717 bk = bkp;
718 }
719 return bk;
720 }
721
722 //______________________________________________________________________________
723 double xMath::BesselI0exp(double x)
724 {
725 // Compute the modified Bessel function I_0(x) for any real x.
726 //
727 //--- NvE 12-mar-2000 UU-SAP Utrecht
728
729 // Parameters of the polynomial approximation
730 const double p1 = 1.0, p2 = 3.5156229, p3 = 3.0899424,
731 p4 = 1.2067492, p5 = 0.2659732, p6 = 3.60768e-2, p7 = 4.5813e-3;
732
733 const double q1 = 0.39894228, q2 = 1.328592e-2, q3 = 2.25319e-3,
734 q4 = -1.57565e-3, q5 = 9.16281e-3, q6 = -2.057706e-2,
735 q7 = 2.635537e-2, q8 = -1.647633e-2, q9 = 3.92377e-3;
736
737 const double k1 = 3.75;
738 double ax = fabs(x);
739
740 double y = 0, result = 0;
741
742 if (ax < k1) {
743 double xx = x / k1;
744 y = xx * xx;
745 result = exp(-ax) * ( p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7))))) );
746 }
747 else {
748 y = k1 / ax;
749 result = (1. / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
750 }
751 return result;
752 }
753
754 //______________________________________________________________________________
755 double xMath::BesselI1exp(double x)
756 {
757 // Compute the modified Bessel function I_1(x) for any real x.
758 //
759 // M.Abramowitz and I.A.Stegun, Handbook of Mathematical Functions,
760 // Applied Mathematics Series vol. 55 (1964), Washington.
761 //
762 //--- NvE 12-mar-2000 UU-SAP Utrecht
763
764 // Parameters of the polynomial approximation
765 const double p1 = 0.5, p2 = 0.87890594, p3 = 0.51498869,
766 p4 = 0.15084934, p5 = 2.658733e-2, p6 = 3.01532e-3, p7 = 3.2411e-4;
767
768 const double q1 = 0.39894228, q2 = -3.988024e-2, q3 = -3.62018e-3,
769 q4 = 1.63801e-3, q5 = -1.031555e-2, q6 = 2.282967e-2,
770 q7 = -2.895312e-2, q8 = 1.787654e-2, q9 = -4.20059e-3;
771
772 const double k1 = 3.75;
773 double ax = fabs(x);
774
775 double y = 0, result = 0;
776
777 if (ax < k1) {
778 double xx = x / k1;
779 y = xx * xx;
780 result = exp(-ax) * (x * (p1 + y * (p2 + y * (p3 + y * (p4 + y * (p5 + y * (p6 + y * p7)))))));
781 }
782 else {
783 y = k1 / ax;
784 result = (1. / sqrt(ax))*(q1 + y * (q2 + y * (q3 + y * (q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))))))));
785 if (x < 0) result = -result;
786 }
787 return result;
788 }
789
790 //______________________________________________________________________________
791 double xMath::BesselIexp(int n, double x)
792 {
793 // Compute the Integer Order Modified Bessel function I_n(x)
794 // for n=0,1,2,... and any real x.
795 //
796 //--- NvE 12-mar-2000 UU-SAP Utrecht
797 if (n < 0) n = -n;
798
799 int iacc = 40; // Increase to enhance accuracy
800 const double kBigPositive = 1.e10;
801 const double kBigNegative = 1.e-10;
802
803 if (n == 0) return xMath::BesselI0exp(x);
804 if (n == 1) return xMath::BesselI1exp(x);
805
806 if (x == 0) return 0;
807 if (fabs(x) > kBigPositive) return 0;
808
809 double tox = 2 / fabs(x);
810 double bip = 0, bim = 0;
811 double bi = 1;
812 double result = 0;
813 int m = 2 * ((n + int(sqrt(float(iacc*n)))));
814 for (int j = m; j >= 1; j--) {
815 bim = bip + double(j)*tox*bi;
816 bip = bi;
817 bi = bim;
818 // Renormalise to prevent overflows
819 if (fabs(bi) > kBigPositive) {
820 result *= kBigNegative;
821 bi *= kBigNegative;
822 bip *= kBigNegative;
823 }
824 if (j == n) result = bip;
825 }
826
827 result *= xMath::BesselI0exp(x) / bi; // Normalise with BesselI0exp(x)
828 if ((x < 0) && (n % 2 == 1)) result = -result;
829
830 return result;
831 }
832
833
834 double xMath::Gamma
835 (
836 double x // We require x > 0
837 )
838 {
839 if (x <= 0.0)
840 {
841 std::stringstream os;
842 os << "Invalid input argument " << x << ". Argument must be positive.";
843 throw std::invalid_argument(os.str());
844 }
845
846 // Split the function domain into three intervals:
847 // (0, 0.001), [0.001, 12), and (12, infinity)
848
850 // First interval: (0, 0.001)
851 //
852 // For small x, 1/Gamma(x) has power series x + gamma x^2 - ...
853 // So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
854 // The relative error over this interval is less than 6e-7.
855
856 const double gamma = 0.577215664901532860606512090; // Euler's gamma constant
857
858 if (x < 0.001)
859 return 1.0 / (x*(1.0 + gamma * x));
860
862 // Second interval: [0.001, 12)
863
864 if (x < 12.0)
865 {
866 // The algorithm directly approximates gamma over (1,2) and uses
867 // reduction identities to reduce other arguments to this interval.
868
869 double y = x;
870 int n = 0;
871 bool arg_was_less_than_one = (y < 1.0);
872
873 // Add or subtract integers as necessary to bring y into (1,2)
874 // Will correct for this below
875 if (arg_was_less_than_one)
876 {
877 y += 1.0;
878 }
879 else
880 {
881 n = static_cast<int> (floor(y)) - 1; // will use n later
882 y -= n;
883 }
884
885 // numerator coefficients for approximation over the interval (1,2)
886 static const double p[] =
887 {
888 -1.71618513886549492533811E+0,
889 2.47656508055759199108314E+1,
890 -3.79804256470945635097577E+2,
891 6.29331155312818442661052E+2,
892 8.66966202790413211295064E+2,
893 -3.14512729688483675254357E+4,
894 -3.61444134186911729807069E+4,
895 6.64561438202405440627855E+4
896 };
897
898 // denominator coefficients for approximation over the interval (1,2)
899 static const double q[] =
900 {
901 -3.08402300119738975254353E+1,
902 3.15350626979604161529144E+2,
903 -1.01515636749021914166146E+3,
904 -3.10777167157231109440444E+3,
905 2.25381184209801510330112E+4,
906 4.75584627752788110767815E+3,
907 -1.34659959864969306392456E+5,
908 -1.15132259675553483497211E+5
909 };
910
911 double num = 0.0;
912 double den = 1.0;
913 int i;
914
915 double z = y - 1;
916 for (i = 0; i < 8; i++)
917 {
918 num = (num + p[i])*z;
919 den = den * z + q[i];
920 }
921 double result = num / den + 1.0;
922
923 // Apply correction if argument was not initially in (1,2)
924 if (arg_was_less_than_one)
925 {
926 // Use identity gamma(z) = gamma(z+1)/z
927 // The variable "result" now holds gamma of the original y + 1
928 // Thus we use y-1 to get back the orginal y.
929 result /= (y - 1.0);
930 }
931 else
932 {
933 // Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
934 for (i = 0; i < n; i++)
935 result *= y++;
936 }
937
938 return result;
939 }
940
942 // Third interval: [12, infinity)
943
944 if (x > 171.624)
945 {
946 // Correct answer too large to display. Force +infinity.
947 double temp = DBL_MAX;
948 return temp * 2.0;
949 }
950
951 return exp(LogGamma(x));
952 }
953
954 double xMath::LogGamma
955 (
956 double x // x must be positive
957 )
958 {
959 if (x <= 0.0)
960 {
961 std::stringstream os;
962 os << "Invalid input argument " << x << ". Argument must be positive.";
963 throw std::invalid_argument(os.str());
964 }
965
966 if (x < 12.0)
967 {
968 return log(fabs(Gamma(x)));
969 }
970
971 // Abramowitz and Stegun 6.1.41
972 // Asymptotic series should be good to at least 11 or 12 figures
973 // For error analysis, see Whittiker and Watson
974 // A Course in Modern Analysis (1927), page 252
975
976 static const double c[8] =
977 {
978 1.0 / 12.0,
979 -1.0 / 360.0,
980 1.0 / 1260.0,
981 -1.0 / 1680.0,
982 1.0 / 1188.0,
983 -691.0 / 360360.0,
984 1.0 / 156.0,
985 -3617.0 / 122400.0
986 };
987 double z = 1.0 / (x*x);
988 double sum = c[7];
989 for (int i = 6; i >= 0; i--)
990 {
991 sum *= z;
992 sum += c[i];
993 }
994 double series = sum / x;
995
996 static const double halfLogTwoPi = 0.91893853320467274178032973640562;
997 double logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;
998 return logGamma;
999 }
1000
1001} // namespace thermalfist
thermalfist::xMath::BesselK0exp
double BesselK0exp(double x)
Modified Bessel function K_0(x), divided by exponential factor.
Definition xMath.cpp:623
thermalfist::xMath::BesselK
double BesselK(int n, double x)
Integer order modified Bessel function K_n(x)
Definition xMath.cpp:162
thermalfist::xMath::BesselI1exp
double BesselI1exp(double x)
Modified Bessel function I_1(x), divided by exponential factor.
Definition xMath.cpp:755
thermalfist::xMath::BesselI0
double BesselI0(double x)
Modified Bessel function I_0(x)
Definition xMath.cpp:23
thermalfist::xMath::BesselY0
double BesselY0(double x)
Bessel function Y0(x) for positive x.
Definition xMath.cpp:308
thermalfist::xMath::Pi
constexpr double Pi()
Pi constant.
Definition xMath.h:23
thermalfist::xMath::BesselI0exp
double BesselI0exp(double x)
Modified Bessel function I_0(x), divided by exponential factor.
Definition xMath.cpp:723
thermalfist::xMath::BesselY1
double BesselY1(double x)
Bessel function Y1(x) for positive x.
Definition xMath.cpp:343
thermalfist::xMath::StruveL0
double StruveL0(double x)
Modified Struve function of order 0.
Definition xMath.cpp:531
thermalfist::xMath::BesselJ0
double BesselJ0(double x)
Bessel function J0(x) for any real x.
Definition xMath.cpp:235
thermalfist::xMath::BesselKexp
double BesselKexp(int n, double x)
Modified Bessel function K_n(x), divided by exponential factor.
Definition xMath.cpp:693
thermalfist::xMath::StruveL1
double StruveL1(double x)
Modified Struve function of order 1.
Definition xMath.cpp:578
thermalfist::xMath::Gamma
double Gamma(double x)
Computes the Gamma function.
Definition xMath.cpp:835
thermalfist::xMath::BesselJ1
double BesselJ1(double x)
Bessel function J1(x) for any real x.
Definition xMath.cpp:271
thermalfist::xMath::BesselK1exp
double BesselK1exp(double x)
Modified Bessel function K_1(x), divided by exponential factor.
Definition xMath.cpp:657
thermalfist::xMath::BesselIexp
double BesselIexp(int n, double x)
Modified Bessel function I_n(x), divided by exponential factor.
Definition xMath.cpp:791
thermalfist::xMath::BesselK1
double BesselK1(double x)
Modified Bessel function K_1(x)
Definition xMath.cpp:127
thermalfist::xMath::LogGamma
double LogGamma(double x)
Computes the logarithm of the Gamma function.
Definition xMath.cpp:955
thermalfist::xMath::StruveH0
double StruveH0(double x)
Struve function of order 0.
Definition xMath.cpp:380
thermalfist::xMath::BesselI
double BesselI(int n, double x)
Integer order modified Bessel function I_n(x)
Definition xMath.cpp:192
thermalfist::xMath::BesselI1
double BesselI1(double x)
Modified Bessel function I_1(x)
Definition xMath.cpp:91
thermalfist::xMath::BesselK0
double BesselK0(double x)
Modified Bessel function K_0(x)
Definition xMath.cpp:55
thermalfist::xMath::StruveH1
double StruveH1(double x)
Struve function of order 1.
Definition xMath.cpp:450
thermalfist
The main namespace where all classes and functions of the Thermal-FIST library reside.
Definition CosmicEoS.h:9
xMath.h
Contains some extra mathematical functions used in the code.
Generated by doxygen 1.13.2