Hi there, I have created a test which times the performance of boost::math::quantile( ... ) when using a gamma distribution. I ran it against source code we use here at work, for ages. It can be found here: http://people.sc.fsu.edu/~burkardt/cpp_src/dcdflib/dcdflib.html The old source code is about 2x times faster than the boost version. MS Visual Studio 2005: boost: 35.4sec att_bell:19sec Intel 11.1: boost: 21.4sec att_bell: 11.2sec Question is if there is a way to incorporate such a function into boost::math? As far, as I can tell the results are almost identical. Regards, Christian Here the code: #include #include #include double min_mean = 2000; // 2,000 double max_mean = 500000000; //500,000,000 double min_std = 10000; // 10,000 double max_std = 100000000; // 100,000,000 double min_max = 600000000; // 600,000,000 double max_max = 1000000000; // 1,000,000,000 const std::size_t max_year = 5000000; // 5,000,000 const double right = 0.999; const double left = 0.001; inline double get_rand() { return static_cast< double >( std::rand() ) / static_cast< double >( RAND_MAX ); } inline void boost_( boost::math::gamma_distribution<>& d, double q ) { double value = boost::math::quantile( d, q ); } inline void att_bell( double alpha, double beta, double q ) { double q_Minus1 = 1 - q; double value = 0.0; double bound = 0.0; int which = 2; int status = 0; cdfgam( &which , &q , &q_Minus1 , &value , &alpha , &beta , &status , &bound ); } int main() { // boost { std::srand( 0 ); boost::timer timer; for( std::size_t y = 0; y < max_year; ++y ) { if(( y % 100000 ) == 0 ) std::cout << y << std::endl; double mean = get_rand() * ( max_mean - min_mean ) + min_mean; double std = get_rand() * ( max_std - min_std ) + min_std; double alpha = mean * mean / std / std; // shape parameter double beta = mean / alpha; // scale parameter boost::math::gamma_distribution<> d( alpha, beta ); boost_( d, right ); boost_( d, left ); } std::cout << "Boost - Time elapsed: " << timer.elapsed() << " sec" << std::endl; } // att bell { std::srand( 0 ); boost::timer timer; for( std::size_t y = 0; y < max_year; ++y ) { if(( y % 100000 ) == 0 ) std::cout << y << std::endl; double mean = get_rand() * ( max_mean - min_mean ) + min_mean; double std = get_rand() * ( max_std - min_std ) + min_std; double alpha = mean * mean / std / std; // shape parameter double beta = mean / alpha; // scale parameter att_bell( alpha, beta, right ); att_bell( alpha, beta, left ); } std::cout << "ATT Bell - Time elapsed: " << timer.elapsed() << " sec" << std::endl; } return 0; }