id,summary,reporter,owner,description,type,status,milestone,component,version,severity,resolution,keywords,cc
3408,gamma distribution's quantile performance,chhenning@…,John Maddock,"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.

Here the code:

#include <dcdflib.h>

#include <boost/math/distributions/gamma.hpp>

#include <boost/timer.hpp>

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;
}
",Tasks,assigned,Boost 1.41.0,math,Boost 1.40.0,Optimization,,,