Whi are chi-squared distributions in (C++) boost::random and C++ STL faster than in boost::math?

Question

I am trying to generate random chi-squared numbers in C++, according to some degree of freedom (which can be a float). Several libraries can be used for that purpose, among which the C++11 Standard Library with the header <random>, or two Boost libraries, boost::mathor boost::random.

These libraries don't exactly have the same user interface, but more surprisingly, they definitely don't produce numbers at the same speed. I have performed a test with the next toy code:

  #include <chrono>
  #include <random>
  #include <boost/math/distributions/chi_squared.hpp>
  #include <boost/random/chi_squared_distribution.hpp>

  unsigned int N = 1e6; 
  double df = 2.4;      //degree of freedom for chi-square

  //generator and distributions  
  std::mt19937 gen;  //random number engine, boost::mt19937 doesn't change anything
  std::chi_squared_distribution<> chi2STL(df); //STL implementation
  boost::random::chi_squared_distribution<> chi2BoostRandom(df); //boost::random 
  boost::math::chi_squared chi2BoostMath(df); //boost::math
  boost::uniform_real<> unif(0., 1.); //uniform distribution for boost::math
                                      //taking STL implementation does not change anything

  //time measure
  auto point1 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with STL
  for(unsigned int i = 0; i < N; ++i)
  {
    chi2STL(gen);
  }

  auto point2 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with boost::random
  for(unsigned int i = 0; i < N; ++i)
  {
    chi2BoostRandom(gen);
  }

  auto point3 = std::chrono::high_resolution_clock::now();

  //generate chi2 numbers with boost::math
  for(unsigned int i = 0; i < N; ++i)
  {
    boost::math::quantile(chi2BoostMath, unif(gen));
  }

  auto point4 = std::chrono::high_resolution_clock::now();

  //get time
  std::cerr << "Time STL loop: "   << std::chrono::duration_cast<std::chrono::microseconds>(point2 - point1).count()/1000000. << " s" << std::endl;
  std::cerr << "Time Boost Random loop: " << std::chrono::duration_cast<std::chrono::microseconds>(point3 - point2).count()/1000000. << " s" << std::endl;
  std::cerr << "Time Boost Math loop: "   << std::chrono::duration_cast<std::chrono::microseconds>(point4 - point3).count()/1000000. << " s" << std::endl;

With the next output:

Time STL loop: 0.474826 s
Time Boost Random loop: 0.374868 s
Time Boost Math loop: 5.02392 s

And it turns out that the boost::math method is way slower than the two other ones. I investigated a little more this surprising point, using a normal law instead of a chi-square law. In that case all three methods produce numbers at the same speed, which makes me think that the boost::quantile method is not the time-consuming one.

By the way, I checked the produced numbers, which all seem to be consistent with the chi-square law, at least for my needs.

It seems that the method implemented in boost::math is just less efficient than the two other ones. (Sadly, it is the only one for which I found implementation details online: it makes use of numerical calculations of the incomplete gamma function).

Is such a simplistic interpretation correct ? And in that case, shouldn't the boost::math implementation be deprecated ?