Gnu Scientific Library poor preformance compared to libc

Question

I've implemented a simple program to measure sin calculation performance by gnu scientific library and libc. Here is a source code :

#include <iostream>
#include <vector>
#include <math.h>
#include <chrono>
#include <list>
#include "gsl/gsl_sf_trig.h"

int main()
{
    std::uint32_t numStepsToCalculate = 1000000;
    const double startX = 0.0;
    const double finishX = M_PI*2;
    const double stepX  = (finishX - startX)/numStepsToCalculate;

    double currentX = startX;
    std::list<double> res;

    auto startT = std::chrono::steady_clock::now();
    while ( currentX <= finishX ) {
        res.push_back( sin ( currentX ) );
        currentX += stepX;
    }
    auto endT = std::chrono::steady_clock::now();

    auto diffT = endT - startT;
    std::cout << "STD : " << std::chrono::duration <double, std::milli> (diffT).count() << " ms" << std::endl;
    std::cout << "STD res size " << res.size() << std::endl;

    std::list<double> resOpt;
    currentX = startX;
    auto startTopt = std::chrono::steady_clock::now();
    while ( currentX <= finishX ) {
        resOpt.push_back( gsl_sf_sin ( currentX ) );
        currentX += stepX;
    }

    auto endTopt = std::chrono::steady_clock::now();

    auto diffTopt = endTopt - startTopt;

    std::cout << "GSL : " << std::chrono::duration <double, std::milli> (diffTopt).count() << " ms" << std::endl;
    std::cout << "GSL res size " << resOpt.size() << std::endl;
    return 0;
}

Here is a result : STD : 57.8869 ms GSL : 106.787 ms

So is it OK for GSL to be slower than libc?

Answer 1

Unfortunately, I can't comment so I post this as a separate answer. I think that @geza's answer is correct, I would simply like to add a few details. According to GSL source codes, they have their own implementation of trigonometric functions:

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

/* I would have prefered just using the library sin() function.
 * But after some experimentation I decided that there was
 * no good way to understand the error; library sin() is just a black box.
 * So we have to roll our own.
 */
int
gsl_sf_sin_e(double x, gsl_sf_result * result)
{
/* CHECK_POINTER(result) */

{
const double P1 = 7.85398125648498535156e-1;
const double P2 = 3.77489470793079817668e-8;
const double P3 = 2.69515142907905952645e-15;

const double sgn_x = GSL_SIGN(x);
const double abs_x = fabs(x);

if(abs_x < GSL_ROOT4_DBL_EPSILON) {
  const double x2 = x*x;
  result->val = x * (1.0 - x2/6.0);
  result->err = fabs(x*x2*x2 / 100.0);
  return GSL_SUCCESS;
}
else {
  double sgn_result = sgn_x;
  double y = floor(abs_x/(0.25*M_PI));
  int octant = y - ldexp(floor(ldexp(y,-3)),3);
  int stat_cs;
  double z;

  if(GSL_IS_ODD(octant)) {
    octant += 1;
    octant &= 07;
    y += 1.0;
  }

  if(octant > 3) {
    octant -= 4;
    sgn_result = -sgn_result;
  }

  z = ((abs_x - y * P1) - y * P2) - y * P3;

  if(octant == 0) {
    gsl_sf_result sin_cs_result;
    const double t = 8.0*fabs(z)/M_PI - 1.0;
    stat_cs = cheb_eval_e(&sin_cs, t, &sin_cs_result);
    result->val = z * (1.0 + z*z * sin_cs_result.val);
  }
  else { /* octant == 2 */
    gsl_sf_result cos_cs_result;
    const double t = 8.0*fabs(z)/M_PI - 1.0;
    stat_cs = cheb_eval_e(&cos_cs, t, &cos_cs_result);
    result->val = 1.0 - 0.5*z*z * (1.0 - z*z * cos_cs_result.val);
  }

  result->val *= sgn_result;

  if(abs_x > 1.0/GSL_DBL_EPSILON) {
    result->err = fabs(result->val);
  }
  else if(abs_x > 100.0/GSL_SQRT_DBL_EPSILON) {
    result->err = 2.0 * abs_x * GSL_DBL_EPSILON * fabs(result->val);
  }
  else if(abs_x > 0.1/GSL_SQRT_DBL_EPSILON) {
    result->err = 2.0 * GSL_SQRT_DBL_EPSILON * fabs(result->val);
  }
  else {
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  }

  return stat_cs;
}
}
}

You should remember that GSL is a library created by scientists and for scientists. Therefore you need a way to control your results and their correctness (see that comment above and this answer for details), speed in this case is secondary.

At the same time you need to remember that system trigonometric functions were implemented by smart IBM guys who know how to make it fast (and absolutely unreadable).

Answer 2

GSL uses software sin (even when a CPU sin is available), and it can give you error estimates too: gsl_sf_sin is just a wrapper currently to gsl_sf_sin_e (which gives error estimates). This routine is not optimized too much for speed.

On the other hand, libc sin is usually optimized pretty well for speed.

Furthermore, your benchmark may be flawed, as sin might be optimized out by the compiler, if optimization is enabled.