From: johnmaddock_at_[hidden]
Date: 2007-06-05 05:25:54
Author: johnmaddock
Date: 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
New Revision: 4453
URL: http://svn.boost.org/trac/boost/changeset/4453
Log:
Added more tests: inverse beta and gamma tests, plus GSL and cephes comparisons.
Text files modified:
sandbox/math_toolkit/libs/math/performance/distributions.cpp | 143 ++++++++++++++++++++++++++++++++-------
sandbox/math_toolkit/libs/math/performance/main.cpp | 31 --------
sandbox/math_toolkit/libs/math/performance/test_erf.cpp | 58 +++++++++++++++
sandbox/math_toolkit/libs/math/performance/test_gamma.cpp | 130 ++++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/libs/math/performance/test_ibeta.cpp | 122 ++++++++++++++++++++++++++++++++++
sandbox/math_toolkit/libs/math/performance/test_igamma.cpp | 117 ++++++++++++++++++++++++++++++++
6 files changed, 543 insertions(+), 58 deletions(-)
Modified: sandbox/math_toolkit/libs/math/performance/distributions.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/distributions.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/distributions.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -192,12 +192,14 @@
BOOST_MATH_DISTRIBUTION2_TEST(lognormal, real_values, real_values, real_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(negative_binomial, int_values, probabilities, int_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(normal, real_values, real_values, real_values, probabilities)
-BOOST_MATH_DISTRIBUTION1_TEST(poisson, real_values, real_values, probabilities)
+BOOST_MATH_DISTRIBUTION1_TEST(poisson, real_values, int_values, probabilities)
BOOST_MATH_DISTRIBUTION1_TEST(students_t, int_values, real_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
#ifdef TEST_R
+#define MATHLIB_STANDALONE 1
+
extern "C" {
#include "Rmath.h"
}
@@ -216,7 +218,7 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_p##name (random_variable_table[k], param1_table[i], param2_table[j], 1, 0);\
+ result += p##name (random_variable_table[k], param1_table[i], param2_table[j], 1, 0);\
}\
}\
}\
@@ -237,7 +239,7 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_d##name (random_variable_table[k], param1_table[i], param2_table[j], 0);\
+ result += d##name (random_variable_table[k], param1_table[i], param2_table[j], 0);\
}\
}\
}\
@@ -258,7 +260,7 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_q##name (probability_table[k], param1_table[i], param2_table[j], 1, 0);\
+ result += q##name (probability_table[k], param1_table[i], param2_table[j], 1, 0);\
}\
}\
}\
@@ -278,7 +280,7 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_p##name (random_variable_table[k], param1_table[i], 1, 0);\
+ result += p##name (random_variable_table[k], param1_table[i], 1, 0);\
}\
}\
\
@@ -295,7 +297,7 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_d##name (random_variable_table[k], param1_table[i], 0);\
+ result += d##name (random_variable_table[k], param1_table[i], 0);\
}\
}\
\
@@ -312,32 +314,13 @@
{\
for(unsigned k = 0; k < c_size; ++k)\
{\
- result += Rf_q##name (probability_table[k], param1_table[i], 1, 0);\
+ result += q##name (probability_table[k], param1_table[i], 1, 0);\
}\
}\
\
consume_result(result);\
set_call_count(a_size * c_size);\
}
-//
-// The R lib actually gives these non-obvious names that we can't deduce
-// in our own macro code, so create unline forwarders and let R's own
-// macros have control over dispatch to the right algorithm:
-//
-inline double Rf_dnorm(double x, double mu, double sigma, int give_log)
-{
- return dnorm(x, mu, sigma, give_log);
-}
-
-inline double Rf_pnorm(double x, double mu, double sigma, int lower_tail, int give_log)
-{
- return pnorm(x, mu, sigma, lower_tail, give_log);
-}
-
-inline double Rf_qnorm(double p, double mu, double sigma, int lower_tail, int log_p)
-{
- return qnorm(p, mu, sigma, lower_tail, log_p);
-}
BOOST_MATH_R_DISTRIBUTION2_TEST(beta, probabilities, probabilities, probabilities, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(binom, int_values, probabilities, int_values, probabilities)
@@ -349,9 +332,115 @@
BOOST_MATH_R_DISTRIBUTION2_TEST(lnorm, real_values, real_values, real_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(nbinom, int_values, probabilities, int_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(norm, real_values, real_values, real_values, probabilities)
-BOOST_MATH_R_DISTRIBUTION1_TEST(pois, real_values, real_values, probabilities)
+BOOST_MATH_R_DISTRIBUTION1_TEST(pois, real_values, int_values, probabilities)
BOOST_MATH_R_DISTRIBUTION1_TEST(t, int_values, real_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
#endif
+#ifdef TEST_CEPHES
+
+extern "C"{
+
+double bdtr(int k, int n, double p);
+double bdtri(int k, int n, double p);
+
+double chdtr(double df, double x);
+double chdtri(double df, double p);
+
+double fdtr(int k, int n, double p);
+double fdtri(int k, int n, double p);
+
+double nbdtr(int k, int n, double p);
+double nbdtri(int k, int n, double p);
+
+}
+
+#define BOOST_MATH_CEPHES_DISTRIBUTION2_TEST(name, param1_table, param2_table, random_variable_table, probability_table) \
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" #name "-cephes-cdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ result += name##dtr (param1_table[i], param2_table[j], random_variable_table[k]);\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_quant, name), "dist-" #name "-cephes-quantile")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned b_size = sizeof(param2_table)/sizeof(param2_table[0]);\
+ unsigned c_size = sizeof(probability_table)/sizeof(probability_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned j = 0; j < b_size; ++j)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ result += name##dtri (param1_table[i], param2_table[j], probability_table[k]);\
+ }\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * b_size * c_size);\
+ }
+
+#define BOOST_MATH_CEPHES_DISTRIBUTION1_TEST(name, param1_table, random_variable_table, probability_table) \
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist, name), "dist-" #name "-cephes-cdf")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned c_size = sizeof(random_variable_table)/sizeof(random_variable_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ result += name##dtr (param1_table[i], random_variable_table[k]);\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * c_size);\
+ }\
+ BOOST_MATH_PERFORMANCE_TEST(BOOST_JOIN(dist_quant, name), "dist-" #name "-cephes-quantile")\
+ {\
+ double result = 0;\
+ unsigned a_size = sizeof(param1_table)/sizeof(param1_table[0]);\
+ unsigned c_size = sizeof(probability_table)/sizeof(probability_table[0]);\
+ \
+ for(unsigned i = 0; i < a_size; ++i)\
+ {\
+ for(unsigned k = 0; k < c_size; ++k)\
+ {\
+ result += name##dtri (param1_table[i], probability_table[k]);\
+ }\
+ }\
+ \
+ consume_result(result);\
+ set_call_count(a_size * c_size);\
+ }
+// Cephes inverse doesn't actually calculate the quantile!!!
+// BOOST_MATH_CEPHES_DISTRIBUTION2_TEST(b, int_values, int_values, probabilities, probabilities)
+BOOST_MATH_CEPHES_DISTRIBUTION1_TEST(ch, int_values, real_values, probabilities)
+BOOST_MATH_CEPHES_DISTRIBUTION2_TEST(f, int_values, int_values, real_values, probabilities)
+// Cephes inverse doesn't calculate the quantile!!!
+// BOOST_MATH_CEPHES_DISTRIBUTION2_TEST(nb, int_values, int_values, probabilities, probabilities)
+
+#endif
+
Modified: sandbox/math_toolkit/libs/math/performance/main.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/main.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/main.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -8,35 +8,6 @@
#include "performance_measure.hpp"
extern void reference_evaluate();
-/*
-extern void polynomial_evaluate();
-extern void polynomial_mixed_evaluate();
-extern void rational_evaluate();
-extern void rational_mixed_evaluate();
-extern void gamma_evaluate();
-extern void lgamma_evaluate();
-extern void erf_evaluate();
-extern void igamma_evaluate();
-extern void igamma_inv_evaluate();
-extern void ibeta_evaluate();
-extern void ibeta_inv_evaluate();
-
-
-
-test_info info[] = {
- { "polynomial", &polynomial_evaluate },
- { "polynomial_mixed", &polynomial_mixed_evaluate },
- { "rational", &rational_evaluate },
- { "rational_mixed", &rational_mixed_evaluate },
- { "gamma", &gamma_evaluate },
- { "lgamma", &lgamma_evaluate },
- { "erf", &erf_evaluate },
- { "igamma_inv", &igamma_inv_evaluate },
- { "igamma", &igamma_evaluate },
- { "ibeta_inv", &ibeta_inv_evaluate },
- { "ibeta", &ibeta_evaluate },
-};
-*/
std::map<std::string, double> times;
std::set<test_info> tests;
@@ -133,12 +104,12 @@
}
}
}
+ run_tests();
}
else
{
show_help();
}
- run_tests();
//
// This is just to confuse the optimisers:
//
Modified: sandbox/math_toolkit/libs/math/performance/test_erf.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/test_erf.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/test_erf.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -45,4 +45,60 @@
consume_result(result);
set_call_count((sizeof(erf_data) + sizeof(erf_large_data) + sizeof(erf_small_data)) / sizeof(erf_data[0]));
-}
\ No newline at end of file
+}
+
+#ifdef TEST_CEPHES
+
+extern "C" {
+
+double erf(double);
+
+}
+
+template <std::size_t N>
+double erf_evaluate_cephes(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += erf(data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(erf_test, "erf-cephes")
+{
+ double result = erf_evaluate_cephes(erf_data);
+ result += erf_evaluate_cephes(erf_large_data);
+ result += erf_evaluate_cephes(erf_small_data);
+
+ consume_result(result);
+ set_call_count((sizeof(erf_data) + sizeof(erf_large_data) + sizeof(erf_small_data)) / sizeof(erf_data[0]));
+}
+
+#endif
+
+#ifdef TEST_GSL
+
+#include <gsl/gsl_sf.h>
+
+template <std::size_t N>
+double erf_evaluate_gsl(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gsl_sf_erf (data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(erf_test, "erf-gsl")
+{
+ double result = erf_evaluate_gsl(erf_data);
+ result += erf_evaluate_gsl(erf_large_data);
+ result += erf_evaluate_gsl(erf_small_data);
+
+ consume_result(result);
+ set_call_count((sizeof(erf_data) + sizeof(erf_large_data) + sizeof(erf_small_data)) / sizeof(erf_data[0]));
+}
+
+#endif
+
+
Modified: sandbox/math_toolkit/libs/math/performance/test_gamma.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/test_gamma.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/test_gamma.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -62,3 +62,133 @@
+ sizeof(near_m10)
+ sizeof(near_m55)) / sizeof(factorials[0]));
}
+
+#ifdef TEST_CEPHES
+
+extern "C" {
+
+double gamma(double);
+double lgam(double);
+
+}
+
+template <std::size_t N>
+double gamma_evaluate_cephes(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gamma(data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(gamma_test, "gamma-cephes")
+{
+ double result = gamma_evaluate_cephes(factorials);
+ result += gamma_evaluate_cephes(near_1);
+ result += gamma_evaluate_cephes(near_2);
+ result += gamma_evaluate_cephes(near_0);
+ result += gamma_evaluate_cephes(near_m10);
+ result += gamma_evaluate_cephes(near_m55);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(factorials)
+ + sizeof(near_1)
+ + sizeof(near_2)
+ + sizeof(near_0)
+ + sizeof(near_m10)
+ + sizeof(near_m55)) / sizeof(factorials[0]));
+}
+
+template <std::size_t N>
+double lgamma_evaluate_cephes(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += lgam(data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(lgamma_test, "lgamma-cephes")
+{
+ double result = lgamma_evaluate_cephes(factorials);
+ result += lgamma_evaluate_cephes(near_1);
+ result += lgamma_evaluate_cephes(near_2);
+ result += lgamma_evaluate_cephes(near_0);
+ result += lgamma_evaluate_cephes(near_m10);
+ result += lgamma_evaluate_cephes(near_m55);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(factorials)
+ + sizeof(near_1)
+ + sizeof(near_2)
+ + sizeof(near_0)
+ + sizeof(near_m10)
+ + sizeof(near_m55)) / sizeof(factorials[0]));
+}
+
+#endif
+
+#ifdef TEST_GSL
+
+#include <gsl/gsl_sf_gamma.h>
+
+template <std::size_t N>
+double gamma_evaluate_gsl(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gsl_sf_gamma(data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(gamma_test, "gamma-gsl")
+{
+ double result = gamma_evaluate_gsl(factorials);
+ result += gamma_evaluate_gsl(near_1);
+ result += gamma_evaluate_gsl(near_2);
+ result += gamma_evaluate_gsl(near_0);
+ result += gamma_evaluate_gsl(near_m10);
+ result += gamma_evaluate_gsl(near_m55);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(factorials)
+ + sizeof(near_1)
+ + sizeof(near_2)
+ + sizeof(near_0)
+ + sizeof(near_m10)
+ + sizeof(near_m55)) / sizeof(factorials[0]));
+}
+
+template <std::size_t N>
+double lgamma_evaluate_gsl(const boost::array<boost::array<T, 3>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gsl_sf_lngamma(data[i][0]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(lgamma_test, "lgamma-gsl")
+{
+ double result = lgamma_evaluate_gsl(factorials);
+ result += lgamma_evaluate_gsl(near_1);
+ result += lgamma_evaluate_gsl(near_2);
+ result += lgamma_evaluate_gsl(near_0);
+ result += lgamma_evaluate_gsl(near_m10);
+ result += lgamma_evaluate_gsl(near_m55);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(factorials)
+ + sizeof(near_1)
+ + sizeof(near_2)
+ + sizeof(near_0)
+ + sizeof(near_m10)
+ + sizeof(near_m55)) / sizeof(factorials[0]));
+}
+
+#endif
+
Modified: sandbox/math_toolkit/libs/math/performance/test_ibeta.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/test_ibeta.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/test_ibeta.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -57,3 +57,125 @@
+ sizeof(ibeta_large_data)
+ sizeof(ibeta_small_data)) / sizeof(ibeta_data[0]));
}
+
+template <std::size_t N>
+double ibeta_invab_evaluate2(const boost::array<boost::array<T, 7>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ {
+ result += boost::math::ibeta_inva(data[i][1], data[i][2], data[i][5]);
+ result += boost::math::ibeta_invb(data[i][0], data[i][2], data[i][5]);
+ result += boost::math::ibetac_inva(data[i][1], data[i][2], data[i][6]);
+ result += boost::math::ibetac_invb(data[i][0], data[i][2], data[i][6]);
+ }
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(ibeta_test, "ibeta_invab")
+{
+ double result = ibeta_invab_evaluate2(ibeta_data);
+ result += ibeta_invab_evaluate2(ibeta_int_data);
+ result += ibeta_invab_evaluate2(ibeta_large_data);
+ result += ibeta_invab_evaluate2(ibeta_small_data);
+
+ consume_result(result);
+ set_call_count(
+ 4 * (sizeof(ibeta_data)
+ + sizeof(ibeta_int_data)
+ + sizeof(ibeta_large_data)
+ + sizeof(ibeta_small_data)) / sizeof(ibeta_data[0]));
+}
+
+#ifdef TEST_CEPHES
+
+extern "C" {
+
+double incbet(double a, double b, double x);
+double incbi(double a, double b, double y);
+
+}
+
+template <std::size_t N>
+double ibeta_evaluate_cephes(const boost::array<boost::array<T, 7>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += incbet(data[i][0], data[i][1], data[i][2]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(ibeta_test, "ibeta-cephes")
+{
+ double result = ibeta_evaluate_cephes(ibeta_data);
+ result += ibeta_evaluate_cephes(ibeta_int_data);
+ result += ibeta_evaluate_cephes(ibeta_large_data);
+ result += ibeta_evaluate_cephes(ibeta_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(ibeta_data)
+ + sizeof(ibeta_int_data)
+ + sizeof(ibeta_large_data)
+ + sizeof(ibeta_small_data)) / sizeof(ibeta_data[0]));
+}
+
+template <std::size_t N>
+double ibeta_inv_evaluate_cephes(const boost::array<boost::array<T, 7>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += incbi(data[i][0], data[i][1], data[i][5]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(ibeta_inv_test, "ibeta_inv-cephes")
+{
+ double result = ibeta_inv_evaluate_cephes(ibeta_data);
+ result += ibeta_inv_evaluate_cephes(ibeta_int_data);
+ result += ibeta_inv_evaluate_cephes(ibeta_large_data);
+ result += ibeta_inv_evaluate_cephes(ibeta_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(ibeta_data)
+ + sizeof(ibeta_int_data)
+ + sizeof(ibeta_large_data)
+ + sizeof(ibeta_small_data)) / sizeof(ibeta_data[0]));
+}
+
+#endif
+
+#ifdef TEST_GSL___
+//
+// This test segfaults inside GSL....
+//
+
+#include <gsl/gsl_sf.h>
+
+template <std::size_t N>
+double ibeta_evaluate_gsl(const boost::array<boost::array<T, 7>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gsl_sf_beta_inc(data[i][0], data[i][1], data[i][2]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(ibeta_test, "ibeta-gsl")
+{
+ double result = ibeta_evaluate_gsl(ibeta_data);
+ result += ibeta_evaluate_gsl(ibeta_int_data);
+ result += ibeta_evaluate_gsl(ibeta_large_data);
+ result += ibeta_evaluate_gsl(ibeta_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(ibeta_data)
+ + sizeof(ibeta_int_data)
+ + sizeof(ibeta_large_data)
+ + sizeof(ibeta_small_data)) / sizeof(ibeta_data[0]));
+}
+
+#endif
+
Modified: sandbox/math_toolkit/libs/math/performance/test_igamma.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/test_igamma.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/test_igamma.cpp 2007-06-05 05:25:54 EDT (Tue, 05 Jun 2007)
@@ -57,3 +57,120 @@
+ sizeof(igamma_med_data)
+ sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
}
+
+template <std::size_t N>
+double igamma_inva_evaluate2(const boost::array<boost::array<T, 6>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += boost::math::gamma_p_inva(data[i][1], data[i][5]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_inva_test, "igamma_inva")
+{
+ double result = igamma_inva_evaluate2(igamma_big_data);
+ result += igamma_inva_evaluate2(igamma_int_data);
+ result += igamma_inva_evaluate2(igamma_med_data);
+ result += igamma_inva_evaluate2(igamma_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(igamma_big_data)
+ + sizeof(igamma_int_data)
+ + sizeof(igamma_med_data)
+ + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+
+#ifdef TEST_CEPHES
+
+extern "C" {
+
+double igam(double, double);
+double igami(double, double);
+
+}
+
+template <std::size_t N>
+double igamma_evaluate_cephes(const boost::array<boost::array<T, 6>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += igam(data[i][0], data[i][1]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_test, "igamma-cephes")
+{
+ double result = igamma_evaluate_cephes(igamma_big_data);
+ result += igamma_evaluate_cephes(igamma_int_data);
+ result += igamma_evaluate_cephes(igamma_med_data);
+ result += igamma_evaluate_cephes(igamma_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(igamma_big_data)
+ + sizeof(igamma_int_data)
+ + sizeof(igamma_med_data)
+ + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+
+template <std::size_t N>
+double igamma_inv_evaluate_cephes(const boost::array<boost::array<T, 6>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += igami(data[i][0], data[i][3]); // note needs complement of probability!!
+ return result;
+}
+/*
+//
+// This test does not run to completion, gets stuck
+// in infinite loop inside cephes....
+//
+BOOST_MATH_PERFORMANCE_TEST(igamma_inv_test, "igamma_inv-cephes")
+{
+ double result = igamma_inv_evaluate_cephes(igamma_big_data);
+ result += igamma_inv_evaluate_cephes(igamma_int_data);
+ result += igamma_inv_evaluate_cephes(igamma_med_data);
+ result += igamma_inv_evaluate_cephes(igamma_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(igamma_big_data)
+ + sizeof(igamma_int_data)
+ + sizeof(igamma_med_data)
+ + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+*/
+#endif
+
+#ifdef TEST_GSL
+
+#include <gsl/gsl_sf_gamma.h>
+
+template <std::size_t N>
+double igamma_evaluate_gsl(const boost::array<boost::array<T, 6>, N>& data)
+{
+ double result = 0;
+ for(unsigned i = 0; i < N; ++i)
+ result += gsl_sf_gamma_inc_P(data[i][0], data[i][1]);
+ return result;
+}
+
+BOOST_MATH_PERFORMANCE_TEST(igamma_test, "igamma-gsl")
+{
+ double result = igamma_evaluate_gsl(igamma_big_data);
+ result += igamma_evaluate_gsl(igamma_int_data);
+ result += igamma_evaluate_gsl(igamma_med_data);
+ result += igamma_evaluate_gsl(igamma_small_data);
+
+ consume_result(result);
+ set_call_count(
+ (sizeof(igamma_big_data)
+ + sizeof(igamma_int_data)
+ + sizeof(igamma_med_data)
+ + sizeof(igamma_small_data)) / sizeof(igamma_big_data[0]));
+}
+
+#endif
\ No newline at end of file