53 for (
unsigned int i = 1;
i < 25;
i += 2)
55 Sum += Sign * Power / Fact;
57 Fact *= (
i + 1) * (
i + 2);
73 for (
unsigned int i = 2;
i < 25;
i += 2)
78 Sum += Sign * Power / Fact;
107 E = 2.71828182845905;
122 L = (P / (
exp(N - 1.0)));
141 if (gamma_in > 1.0 + ep)
143 return ((gamma_in - 1) *
gamma(gamma_in - 1));
148 if ((gamma_in - ep < 1.0) && (gamma_in + ep > 1.0))
150 else if ((gamma_in - ep < 0.5) && (gamma_in + ep > 0.5))
153 return std::numeric_limits<Mdouble>::quiet_NaN();
166 return std::exp(std::lgamma(z) + std::lgamma(w) - std::lgamma(z + w));
176 Mdouble mainfactor = pow(x, k / 2.0 - 1) *
exp(x / -2.0);
178 return mainfactor / prefactor;
192 const int num_steps_per_unit = 100;
195 long int num_steps =
static_cast<int>(num_steps_per_unit * x_max);
197 for (
int i = 0;
i < num_steps;
i++)
199 x = x_max / num_steps * (
i + 0.5);
202 return 1.0 - sum * x_max / num_steps;
209 if (std::abs(max - min) < endCondition)
211 return 0.5 * (min + max);
213 logger(
INFO,
"Min=% Max=% diff= %", min, max, max - min);
214 Mdouble resphi = 2 - 0.5 * (1 + std::sqrt(5));
216 if (max - cur > cur - min)
218 x = cur + resphi * (max - cur);
222 x = cur - resphi * (cur - min);
224 if (std::isnan(curVal))
225 curVal =
function(cur);
229 if (max - cur > cur - min)
240 if (max - cur > cur - min)
253 return std::abs(v1 - v2) <= absError;
300 logger.assert_debug(i > 0,
"i is greater than 0");
305 b0 = x * b1 - b2 + *p++;
308 return (0.5 * (b0 - b2));
316 -4.415341646479339379501E-18,
317 3.330794518822238097831E-17,
318 -2.431279846547954693591E-16,
319 1.715391285555133030611E-15,
320 -1.168533287799345168081E-14,
321 7.676185498604935616881E-14,
322 -4.856446783111929460901E-13,
323 2.955052663129639834611E-12,
324 -1.726826291441555707231E-11,
325 9.675809035373236912241E-11,
326 -5.189795601635262906661E-10,
327 2.659823724682386650351E-9,
328 -1.300025009986248042121E-8,
329 6.046995022541918949321E-8,
330 -2.670793853940611733911E-7,
331 1.117387539120103718151E-6,
332 -4.416738358458750563591E-6,
333 1.644844807072889708931E-5,
334 -5.754195010082103703981E-5,
335 1.885028850958416557291E-4,
336 -5.763755745385823658851E-4,
337 1.639475616941335798421E-3,
338 -4.324309995050575944301E-3,
339 1.054646039459499831831E-2,
340 -2.373741480589946881561E-2,
341 4.930528423967070848781E-2,
342 -9.490109704804764442101E-2,
343 1.716209015222087753491E-1,
344 -3.046826723431983986831E-1,
345 6.767952744094760849951E-1
351 -7.233180487874753954561E-18,
352 -4.830504485944182071261E-18,
353 4.465621420296759999011E-17,
354 3.461222867697461093101E-17,
355 -2.827623980516583484941E-16,
356 -3.425485619677219134621E-16,
357 1.772560133056526383601E-15,
358 3.811680669352622420751E-15,
359 -9.554846698828307648701E-15,
360 -4.150569347287222086631E-14,
361 1.540086217521409826911E-14,
362 3.852778382742142701141E-13,
363 7.180124451383666233671E-13,
364 -1.794178531506806117781E-12,
365 -1.321581184044771311881E-11,
366 -3.149916527963241364541E-11,
367 1.188914710784643834241E-11,
368 4.940602388224969589101E-10,
369 3.396232025708386345151E-9,
370 2.266668990498178064591E-8,
371 2.048918589469063741831E-7,
372 2.891370520834756482971E-6,
373 6.889758346916823984261E-5,
374 3.369116478255694089901E-3,
375 8.044904110141088316081E-1
387 return (
chebyshev(32.0 / x - 2.0, B, 25) / sqrt(x));
405 std::size_t nTerms = 0.5 * (n + 1) * (n + 2);
408 size_t location_current;
409 size_t location_previous;
413 for (
int l = 1; l <= n; l++)
416 location_current = 0.5 * (l + 1) * (l + 2) - 1;
417 location_previous = location_current - (l + 1);
418 polynomials(location_current) = -(2.0 * (l - 1.0) + 1.0) * std::sqrt(1.0 - x * x) *
419 polynomials(location_previous);
422 polynomials(location_current - 1) =
423 x * (2.0 * (l - 1) + 1) * polynomials(location_previous);
427 for (
int l = 2; l <= n; l++)
429 for (
int m = (l - 2); m >= 0; m--)
431 location_current = (0.5 * (l + 1) * (l + 2) - 1) - l + m;
432 temp = polynomials(location_current + 2) +
433 2.0 * (m + 1) * x / sqrt(1 - x * x) * polynomials(location_current + 1);
434 polynomials(location_current) = temp / ((m - l) * (l + m + 1));
447 std::size_t nTerms = 0.5 * (p + 1) * (2 * p + 2);
452 for (
int n = 0; n <= p; n++)
454 for (
int mt = -n; mt <= n; mt++)
458 std::size_t location_current = n * n + (m + n);
459 std::size_t location_polynomial = 0.5 * n * (n + 1) + m_abs;
462 Mdouble fact = 1.0 * fact1 / fact2;
463 std::complex<Mdouble> value =
465 Y(location_current) = value;
475 std::size_t nTerms = 0.5 * (p + 1) * (2 * p + 2);
478 for (
int n = 0; n <= p; n++)
480 for (
int m = -n; m <= n; m++)
482 std::size_t location = n * n + (m + n);
486 squaredFactorials(location) = std::pow(-1.0, n) / std::sqrt(fact);
490 return squaredFactorials;
Implementation of a 3D quaternion (by Vitaliy).
Mdouble X
the vector components
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
Definition of different loggers with certain modules. A user can define its own custom logger here...
Mdouble goldenSectionSearch(Mdouble(*function)(const Mdouble), Mdouble min, Mdouble cur, Mdouble max, Mdouble endCondition, Mdouble curVal=std::numeric_limits< Mdouble >::quiet_NaN())
This function performs a golden section search to find the location of the minimum of a function...
Mdouble getComponent(int index) const
Returns the requested component of this Quaternion.
LL< Log::INFO > INFO
Info log level.
Mdouble exp(Mdouble Exponent)
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
const std::complex< Mdouble > i
Mdouble log(Mdouble Power)
bool isEqual(Mdouble v1, Mdouble v2, Mdouble absError)
Compares the difference of two Mdouble with an absolute error, useful in UnitTests.
Mdouble chi_squared_prob(Mdouble x, unsigned int k)
This is the function which actually gives the probability back using a chi squared test...
constexpr T factorial(const T t)
factorial function
Mdouble chi_squared(Mdouble x, unsigned int k)
This is a chi_squared function return the value x and degrees of freedom k.
Mdouble I0_exp(Mdouble x)
Mdouble chebyshev(Mdouble x, const Mdouble coef[], int N)
Namespace for evaluating the zeroth modified Bessel function of the first kind, I0(x), required in StatisticsPoint.hcc.
NumericalVector< std::complex< Mdouble > > sphericalHarmonics(int p, Mdouble theta, Mdouble phi)
Mdouble XX
all nine matrix elements
Mdouble gamma(Mdouble gamma_in)
This is the gamma function returns the true value for the half integer value.
NumericalVector associatedLegendrePolynomials(int n, Mdouble x)
Implementation of a 3D matrix.
Implementation of a 3D symmetric matrix.
NumericalVector computeSquaredFactorialValues(int p)
Mdouble XX
The six distinctive matrix elements.