Using version 2.6 of GSL from ftp://ftp.gnu.org/gnu/gsl/gsl-2.6.tar.gz
OS is Ubuntu 18.04.2 LTS Hardware is a Dell workstation with a 3.60 GHz Intel i7-4790 and 32 GB of RAM Compiler is gcc 7.4.0 with the -w option When the second argument is 0, mpmath and scipy's implementations of the gegenbauer polynomial agree that the output should be uniformly 0, but gsl_sf_gegenpoly_n consistently returns nonzero values. This happens regardless of the value of the first argument, though in particular when the first argument is 1 and the second is 0, the function returns twice the third argument... #include <gsl/gsl_sf.h> #include <stdio.h> #include <gsl/gsl_math.h> #include <float.h> int main (void){ double out; out = gsl_sf_gegenpoly_n(1, 0, 5); printf("%.*f\n", DBL_DIG-1, out); printf("%.21e\n", out); return 0; }