11 #include "../Eigen/SpecialFunctions" 13 template<
typename X,
typename Y>
16 for(
Index i=0; i<x.size(); ++i)
19 VERIFY_IS_APPROX( x(i), y(i) );
23 VERIFY_IS_EQUAL( x(i), y(i) );
31 typedef typename ArrayType::Scalar Scalar;
34 Scalar plusinf = std::numeric_limits<Scalar>::infinity();
35 Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
37 Index rows = internal::random<Index>(1,30);
42 ArrayType m1 = ArrayType::Random(rows,cols);
43 #if EIGEN_HAS_C99_MATH 44 VERIFY_IS_APPROX(m1.lgamma(),
lgamma(m1));
45 VERIFY_IS_APPROX(m1.digamma(),
digamma(m1));
46 VERIFY_IS_APPROX(m1.erf(),
erf(m1));
47 VERIFY_IS_APPROX(m1.erfc(),
erfc(m1));
48 #endif // EIGEN_HAS_C99_MATH 52 #if EIGEN_HAS_C99_MATH 58 ArrayType m1 = ArrayType::Random(rows,cols);
59 ArrayType m2 = ArrayType::Random(rows,cols);
67 ArrayType a = m1.abs() + 2;
68 ArrayType x = m2.abs() + 2;
69 ArrayType zero = ArrayType::Zero(rows, cols);
70 ArrayType one = ArrayType::Constant(rows, cols, Scalar(1.0));
71 ArrayType a_m1 = a - one;
73 ArrayType Gamma_a_m1_x =
Eigen::igammac(a_m1, x) * a_m1.lgamma().exp();
75 ArrayType gamma_a_m1_x =
Eigen::igamma(a_m1, x) * a_m1.lgamma().exp();
81 VERIFY_IS_APPROX(Gamma_a_x + gamma_a_x, a.lgamma().exp());
84 VERIFY_IS_APPROX(Gamma_a_x, (a - 1) * Gamma_a_m1_x + x.pow(a-1) * (-x).
exp());
87 VERIFY_IS_APPROX(gamma_a_x, (a - 1) * gamma_a_m1_x - x.pow(a-1) * (-x).
exp());
92 Scalar a_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
93 Scalar x_s[] = {Scalar(0), Scalar(1), Scalar(1.5), Scalar(4), Scalar(0.0001), Scalar(1000.5)};
96 Scalar igamma_s[][6] = {{0.0, nan, nan, nan, nan, nan},
97 {0.0, 0.6321205588285578, 0.7768698398515702,
98 0.9816843611112658, 9.999500016666262e-05, 1.0},
99 {0.0, 0.4275932955291202, 0.608374823728911,
100 0.9539882943107686, 7.522076445089201e-07, 1.0},
101 {0.0, 0.01898815687615381, 0.06564245437845008,
102 0.5665298796332909, 4.166333347221828e-18, 1.0},
103 {0.0, 0.9999780593618628, 0.9999899967080838,
104 0.9999996219837988, 0.9991370418689945, 1.0},
105 {0.0, 0.0, 0.0, 0.0, 0.0, 0.5042041932513908}};
106 Scalar igammac_s[][6] = {{nan, nan, nan, nan, nan, nan},
107 {1.0, 0.36787944117144233, 0.22313016014842982,
108 0.018315638888734182, 0.9999000049998333, 0.0},
109 {1.0, 0.5724067044708798, 0.3916251762710878,
110 0.04601170568923136, 0.9999992477923555, 0.0},
111 {1.0, 0.9810118431238462, 0.9343575456215499,
112 0.4334701203667089, 1.0, 0.0},
113 {1.0, 2.1940638138146658e-05, 1.0003291916285e-05,
114 3.7801620118431334e-07, 0.0008629581310054535,
116 {1.0, 1.0, 1.0, 1.0, 1.0, 0.49579580674813944}};
117 for (
int i = 0; i < 6; ++i) {
118 for (
int j = 0; j < 6; ++j) {
134 #endif // EIGEN_HAS_C99_MATH 138 ArrayType x(7),
q(7), res(7), ref(7);
139 x << 1.5, 4, 10.5, 10000.5, 3, 1, 0.9;
140 q << 2, 1.5, 3, 1.0001, -2.5, 1.2345, 1.2345;
141 ref << 1.61237534869, 0.234848505667, 1.03086757337e-5, 0.367879440865, 0.054102025820864097, plusinf, nan;
149 ArrayType x(7), res(7), ref(7);
150 x << 1, 1.5, 4, -10.5, 10000.5, 0, -1;
151 ref << -0.5772156649015329, 0.03648997397857645, 1.2561176684318, 2.398239129535781, 9.210340372392849, plusinf, plusinf;
159 #if EIGEN_HAS_C99_MATH 161 ArrayType n(11), x(11), res(11), ref(11);
162 n << 1, 1, 1, 1.5, 17, 31, 28, 8, 42, 147, 170;
163 x << 2, 3, 25.5, 1.5, 4.7, 11.8, 17.7, 30.2, 15.8, 54.1, 64;
164 ref << 0.644934066848, 0.394934066848, 0.0399946696496, nan, 293.334565435, 0.445487887616, -2.47810300902e-07, -8.29668781082e-09, -0.434562276666, 0.567742190178, -0.0108615497927;
167 if(
sizeof(RealScalar)>=8) {
179 #if EIGEN_HAS_C99_MATH 196 a << 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
197 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
198 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
199 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
200 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
201 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
202 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
203 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
204 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
205 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
206 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999,
207 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999,
208 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999, 0.999,
209 31.62177660168379, 31.62177660168379, 31.62177660168379,
210 31.62177660168379, 31.62177660168379, 31.62177660168379,
211 31.62177660168379, 31.62177660168379, 31.62177660168379,
212 31.62177660168379, 31.62177660168379, 31.62177660168379,
213 31.62177660168379, 31.62177660168379, 31.62177660168379,
214 31.62177660168379, 31.62177660168379, 31.62177660168379,
215 31.62177660168379, 31.62177660168379, 31.62177660168379,
216 31.62177660168379, 31.62177660168379, 31.62177660168379,
217 31.62177660168379, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999,
218 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999,
219 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999, 999.999,
220 999.999, 999.999, 999.999;
222 b << 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379, 0.03062277660168379,
223 0.03062277660168379, 0.03062277660168379, 0.03062277660168379, 0.999,
224 0.999, 0.999, 0.999, 0.999, 31.62177660168379, 31.62177660168379,
225 31.62177660168379, 31.62177660168379, 31.62177660168379, 999.999,
226 999.999, 999.999, 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0,
227 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
228 0.03062277660168379, 0.03062277660168379, 0.999, 0.999, 0.999, 0.999,
229 0.999, 31.62177660168379, 31.62177660168379, 31.62177660168379,
230 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999,
231 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379,
232 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
233 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999,
234 31.62177660168379, 31.62177660168379, 31.62177660168379,
235 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999,
236 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379,
237 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
238 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999,
239 31.62177660168379, 31.62177660168379, 31.62177660168379,
240 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999,
241 999.999, 999.999, 0.0, 0.0, 0.0, 0.0, 0.0, 0.03062277660168379,
242 0.03062277660168379, 0.03062277660168379, 0.03062277660168379,
243 0.03062277660168379, 0.999, 0.999, 0.999, 0.999, 0.999,
244 31.62177660168379, 31.62177660168379, 31.62177660168379,
245 31.62177660168379, 31.62177660168379, 999.999, 999.999, 999.999,
248 x << -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5,
249 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2,
250 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1,
251 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1,
252 -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8,
253 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5,
254 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2,
255 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1,
256 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5, 0.8, 1.1, -0.1, 0.2, 0.5,
259 v << nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
260 nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
261 nan, nan, nan, 0.47972119876364683, 0.5, 0.5202788012363533, nan, nan,
262 0.9518683957740043, 0.9789663010413743, 0.9931729188073435, nan, nan,
263 0.999995949033062, 0.9999999999993698, 0.9999999999999999, nan, nan,
264 0.9999999999999999, 0.9999999999999999, 0.9999999999999999, nan, nan,
265 nan, nan, nan, nan, nan, 0.006827081192655869, 0.0210336989586256,
266 0.04813160422599567, nan, nan, 0.20014344256217678, 0.5000000000000001,
267 0.7998565574378232, nan, nan, 0.9991401428435834, 0.999999999698403,
268 0.9999999999999999, nan, nan, 0.9999999999999999, 0.9999999999999999,
269 0.9999999999999999, nan, nan, nan, nan, nan, nan, nan,
270 1.0646600232370887e-25, 6.301722877826246e-13, 4.050966937974938e-06,
271 nan, nan, 7.864342668429763e-23, 3.015969667594166e-10,
272 0.0008598571564165444, nan, nan, 6.031987710123844e-08,
273 0.5000000000000007, 0.9999999396801229, nan, nan, 0.9999999999999999,
274 0.9999999999999999, 0.9999999999999999, nan, nan, nan, nan, nan, nan,
275 nan, 0.0, 7.029920380986636e-306, 2.2450728208591345e-101, nan, nan,
276 0.0, 9.275871147869727e-302, 1.2232913026152827e-97, nan, nan, 0.0,
277 3.0891393081932924e-252, 2.9303043666183996e-60, nan, nan,
278 2.248913486879199e-196, 0.5000000000004947, 0.9999999999999999, nan;
280 CALL_SUBTEST(res =
betainc(a, b, x);
286 ArrayType m1 = ArrayType::Random(32);
287 ArrayType m2 = ArrayType::Random(32);
288 ArrayType m3 = ArrayType::Random(32);
289 ArrayType one = ArrayType::Constant(32, Scalar(1.0));
290 const Scalar eps = std::numeric_limits<Scalar>::epsilon();
291 ArrayType a = (m1 * 4.0).exp();
292 ArrayType
b = (m2 * 4.0).exp();
293 ArrayType x = m3.abs();
298 ArrayType expected = x.pow(a);
304 ArrayType expected = one - (one - x).
pow(b);
310 ArrayType expected = one;
315 ArrayType num = x.pow(a) * (one - x).
pow(b);
316 ArrayType denom = a * (a.lgamma() + b.lgamma() - (a +
b).
lgamma()).
exp();
319 ArrayType expected =
betainc(a, b, x) - num / denom + eps;
321 if (
sizeof(Scalar) >= 8) {
332 ArrayType num = x.pow(a) * (one - x).
pow(b);
333 ArrayType denom = b * (a.lgamma() + b.lgamma() - (a +
b).
lgamma()).
exp();
334 ArrayType expected =
betainc(a, b, x) + num / denom + eps;
335 ArrayType test =
betainc(a, b + one, x) + eps;
343 CALL_SUBTEST_1(array_special_functions<ArrayXf>());
344 CALL_SUBTEST_2(array_special_functions<ArrayXd>());
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool() isfinite(const half &a)
const CwiseBinaryOp< internal::scalar_zeta_op< Scalar >, const Derived, const DerivedQ > zeta(const EIGEN_CURRENT_STORAGE_BASE_CLASS< DerivedQ > &q) const
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
EIGEN_DEVICE_FUNC const ErfReturnType erf() const
EIGEN_DEVICE_FUNC const SqrtReturnType sqrt() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const AbsReturnType abs() const
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igammac_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igammac(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_igamma_op< typename Derived::Scalar >, const Derived, const ExponentDerived > igamma(const Eigen::ArrayBase< Derived > &a, const Eigen::ArrayBase< ExponentDerived > &x)
void verify_component_wise(const X &x, const Y &y)
void array_special_functions()
internal::enable_if< !(internal::is_same< typename Derived::Scalar, ScalarExponent >::value)&&EIGEN_SCALAR_BINARY_SUPPORTED(pow, typename Derived::Scalar, ScalarExponent), const EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(Derived, ScalarExponent, pow) >::type pow(const Eigen::ArrayBase< Derived > &x, const ScalarExponent &exponent)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
EIGEN_DEVICE_FUNC const Scalar & q
void test_special_functions()
EIGEN_DEVICE_FUNC const ErfcReturnType erfc() const
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_polygamma_op< typename DerivedX::Scalar >, const DerivedN, const DerivedX > polygamma(const Eigen::ArrayBase< DerivedN > &n, const Eigen::ArrayBase< DerivedX > &x)
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bool() isnan(const half &a)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const TensorCwiseTernaryOp< internal::scalar_betainc_op< typename XDerived::Scalar >, const ADerived, const BDerived, const XDerived > betainc(const ADerived &a, const BDerived &b, const XDerived &x)
EIGEN_DEVICE_FUNC const Scalar & b
EIGEN_DEVICE_FUNC const DigammaReturnType digamma() const
EIGEN_DEVICE_FUNC const LgammaReturnType lgamma() const