[adrian_b]

If any of those libraries uses the Taylor expansion for approximation, that is a big mistake, because the approximation error becomes large at the upper end of the argument interval, even if it is small close to zero.

What is much more likely is that if you will carefully compare the polynomial coefficients with those of the Taylor series, you will see that the last decimals are different and the difference from the Taylor series increases towards the coefficients corresponding to higher degrees.

Towards zero, any approximation polynomial begins to resemble the Taylor series in the low-degree terms, because the high-degree terms become negligible and the limit of the Taylor series and of the approximation polynomial is the same.

So when looking at the polynomial coefficients, they should resemble those of the Taylor series in the low-degree coefficients, but an accurate coefficient computation should demonstrate that the polynomials are different.

[xyzzyz]

I did some testing, and you are correct that they are slightly different:

  >>> print("{:.30f}".format(sin(pi/4, 0, 0, *standard_coeff)))
  0.707106781186567889818661569734
  >>> print("{:.30f}".format(sin(pi/4, 0, 0, *musl_coeff)))
  0.707106781186547461715008466854

The difference at the very edge of the interval occurs at the 14th digit of decimal expansion, and it's at the edge of accuracy of double, at 16th digit: after ...6547, the exact value starts with ...6547_5244, instead of 4617. I wouldn't exactly call it a big mistake, as the difference would not be relevant in almost all practical uses, but that would be a mistake nevertheless, and I'm sure someone would be bitten by this. Thanks, I learned something new today!