https://gcc.gnu.org/bugzilla/show_bug.cgi?id=66689
--- Comment #1 from John Maddock <john at johnmaddock dot co.uk> --- I think I've found the problem - make the nu parameter negative and you get the correct answers. These functions appear to use the same definition for ellint_3 as GSL uses (same code?), but TR1 uses a different one which matches A&S: note the change of sign in the nu parameter in the integral: GSL has Pi(\phi,k,n) = \int_0^\phi dt 1/((1 + n \sin^2(t))\sqrt(1 - k^2 \sin^2(t))) TR1 has Pi(\phi,k,n) = \int_0^\phi dt 1/((1 - n \sin^2(t))\sqrt(1 - k^2 \sin^2(t)))