Wrapper function to do polynomial fits with gsl
One issue is numerical stability. Although the GSL uses stable routines, the problem of fitting a polynomial in the monomial basis is so ill-conditioned (especially for large ranges in x and large (eg 50) degree), that you risk getting significant errors. Unless you really, really need the coefficients in this basis then I think you should use another basis, for example Legendre or Chebycheff polynomials. This will mean that any code that uses the coefficients will need to know how to do so, but the evaluation is in fact almost as efficient as evaluation in the polynomial basis.