Go to the documentation of this file.
10 static double didonato_SN(
double,
double,
unsigned,
double);
25 double a[4] = {0.213623493715853, 4.28342155967104,
26 11.6616720288968, 3.31125922108741};
27 double b[5] = {0.3611708101884203e-1, 1.27364489782223,
28 6.40691597760039, 6.61053765625462, 1};
57 double partial =
x / (
a + 1);
60 for(
i = 2;
i <=
N; ++
i) {
61 partial *=
x / (
a +
i);
63 if(partial < tolerance) {
97 if ((
b > 0.6) || ((
b >= 0.45) && (
a >= 0.3))) {
106 if((
b *
q > 1
e-8) && (
q > 1
e-5)) {
112 result = u / (1 - (u / (
a + 1)));
114 else if ((
a < 0.3) && (
b >= 0.35)) {
117 double u =
t *
exp(
t);
120 else if ((
b > 0.15) || (
a >= 0.3)) {
123 double u =
y - (1 -
a) *
log(
y);
129 double u =
y - (1 -
a) *
log(
y);
131 -
log((u * u + 2 * (3 -
a) * u + (2 -
a) * (3 -
a))
132 / (u * u + (5 -
a) * u + 2));
138 double c1_2 =
c1 *
c1;
139 double c1_3 = c1_2 *
c1;
140 double c1_4 = c1_2 * c1_2;
142 double a_3 = a_2 *
a;
144 double c2 = (
a - 1) * (1 +
c1);
145 double c3 = (
a - 1) * (-(c1_2 / 2)
148 double c4 = (
a - 1) * ((c1_3 / 3) - (3 *
a - 5) * c1_2 / 2
149 + (a_2 - 6 *
a + 7) *
c1
150 + (11 * a_2 - 46 *
a + 47) / 6);
151 double c5 = (
a - 1) * (-(c1_4 / 4)
152 + (11 *
a - 17) * c1_3 / 6
153 + (-3 * a_2 + 13 *
a -13) * c1_2
154 + (2 * a_3 - 25 * a_2 + 72 *
a - 61) *
c1 / 2
155 + (25 * a_3 - 195 * a_2 + 477 *
a - 379) / 12);
158 double y_3 = y_2 *
y;
159 double y_4 = y_2 * y_2;
160 result =
y +
c1 + (
c2 /
y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
168 double s_3 = s_2 *
s;
169 double s_4 = s_2 * s_2;
170 double s_5 = s_4 *
s;
173 double w =
a +
s * ra + (s_2 - 1) / 3;
174 w += (s_3 - 7 *
s) / (36 * ra);
175 w -= (3 * s_4 + 7 * s_2 - 16) / (810 *
a);
176 w += (9 * s_5 + 256 * s_3 - 433 *
s) / (38880 *
a * ra);
178 if ((
a >= 500) && (
fabs(1 -
w /
a) < 1
e-6)) {
186 double D =
fmax(2,
a * (
a - 1));
188 double lb =
log(
q) + lg;
193 double c1_2 =
c1 *
c1;
194 double c1_3 = c1_2 *
c1;
195 double c1_4 = c1_2 * c1_2;
197 double a_3 = a_2 *
a;
199 double c2 = (
a - 1) * (1 +
c1);
200 double c3 = (
a - 1) * (-(c1_2 / 2)
203 double c4 = (
a - 1) * ((c1_3 / 3)
204 - (3 *
a - 5) * c1_2 / 2
205 + (a_2 - 6 *
a + 7) *
c1
206 + (11 * a_2 - 46 *
a + 47) / 6);
207 double c5 = (
a - 1) * (-(c1_4 / 4)
208 + (11 *
a - 17) * c1_3 / 6
209 + (-3 * a_2 + 13 *
a -13) * c1_2
210 + (2 * a_3 - 25 * a_2 + 72 *
a - 61) *
c1 / 2
211 + (25 * a_3 - 195 * a_2 + 477 *
a - 379) / 12);
214 double y_3 = y_2 *
y;
215 double y_4 = y_2 * y_2;
216 result =
y +
c1 + (
c2 /
y) + (c3 / y_2) + (c4 / y_3) + (c5 / y_4);
220 double u = -lb + (
a - 1) *
log(
w) -
log(1 + (1 -
a) / (1 +
w));
229 if (
w < 0.15 * ap1) {
237 s =
log1p(
z / ap1 * (1 +
z / ap2 * (1 +
z / (
a + 3))));
241 if ((
z <= 0.01 * ap1) || (
z > 0.7 * ap1)) {
260 double x, fac, f_fp, fpp_fp;
265 else if ((
a < 0) || (
p < 0) || (
p > 1)) {
280 for (
i = 0;
i < 3;
i++) {
287 fpp_fp = -1.0 + (
a - 1) /
x;
293 x =
x - f_fp / (1.0 - 0.5 * f_fp * fpp_fp);
304 double x, fac, f_fp, fpp_fp;
309 else if ((
a < 0.0) || (
q < 0.0) || (
q > 1.0)) {
323 for (
i = 0;
i < 3;
i++) {
329 fpp_fp = -1.0 + (
a - 1) /
x;
334 x =
x - f_fp / (1.0 - 0.5 * f_fp * fpp_fp);
double igamc(double a, double x)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
static double didonato_SN(double, double, unsigned, double)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
double igam(double a, double x)
const EIGEN_DEVICE_FUNC LogReturnType log() const
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
static double find_inverse_s(double, double)
EIGEN_DEVICE_FUNC const Scalar & q
static double polevl(double x, const double coef[], int N)
double igam_fac(double a, double x)
void g(const string &key, int i)
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 fmax(const bfloat16 &a, const bfloat16 &b)
static double find_inverse_gamma(double, double, double)
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Array< int, Dynamic, 1 > v
Jet< T, N > sqrt(const Jet< T, N > &f)
double igamci(double a, double q)
double igami(double a, double p)
gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:44