14#ifndef GEMMI_BESSEL_HPP_
15#define GEMMI_BESSEL_HPP_
23 static_assert(
N > 1,
"");
25 for (
int i =
N-2;
i >= 0; --
i)
33 static const double P1[8];
34 static const double Q1[9];
35 static const double P2[5];
36 static const double Q2[5];
37 static const double Q3[3];
39template<
class Dummy>
const double BesselTables_<Dummy>::P1[8] = {
49template<
class Dummy>
const double BesselTables_<Dummy>::Q1[9] = {
50 1.00000003928615375e+00,
51 2.49999576572179639e-01,
52 2.77785268558399407e-02,
53 1.73560257755821695e-03,
54 6.96166518788906424e-05,
55 1.89645733877137904e-06,
56 4.29455004657565361e-08,
57 3.90565476357034480e-10,
58 1.48095934745267240e-11
60template<
class Dummy>
const double BesselTables_<Dummy>::P2[5] = {
61 3.98942115977513013e-01,
62 -1.49581264836620262e-01,
63 -4.76475741878486795e-02,
64 -2.65157315524784407e-02,
65 -1.47148600683672014e-01
67template<
class Dummy>
const double BesselTables_<Dummy>::Q2[5] = {
68 3.98942651588301770e-01,
69 4.98327234176892844e-02,
70 2.91866904423115499e-02,
71 1.35614940793742178e-02,
72 1.31409251787866793e-01
74template<
class Dummy>
const double BesselTables_<Dummy>::Q3[3] = {
75 3.98942391532752700e-01,
76 4.98455950638200020e-02,
77 2.94835666900682535e-02
105 double a = x * x / 4;
110 double ex = std::exp(x / 2);
119 double a = x * x / 4;
124 return x + std::log(
q / std::sqrt(x));
double evaluate_polynomial(const double(&poly)[N], double x)
double bessel_i0(double x)
double bessel_i1_over_i0(double x)
double log_bessel_i0(double x)