-sqrt(-x^2 + 1)
sqrt(x^2 - 1)
5.30287619362453e-17 + 0.866025403784439*I
The result from eval_poly agrees with Mathematica, i.e.
LegendreP[1, 1, 0.5]
This was causing me difficulties because I was using gen_legendre_P(l,m,cos(theta)) to construct real spherical harmonics, but finding that for certain l, m combinations, the resulting spherical harmonics were not real.
Based on the above output, it seems to me that gen_legendre_P.eval_poly(1,1,cos(theta)) will always be real while gen_legendre_P(1,1,cos(theta)) will be complex (unless |cos(theta)| = 1), since cos(theta) is in the interval [-1,1].
Looking at the code for Func_associ_legendre_P._eval_special_values_, I suspect the culprit is the n == m case, which returns
factorial(2*m)/2**m/factorial(m) * (x**2-1)**(m/2)
Interestingly, it looks like this agrees with 14.7.15 at https://dlmf.nist.gov/14.7.E15, but differs from the definition on Wolfram Mathworld (http://mathworld.wolfram.com/AssociatedLegendrePolynomial.html, Eq. 10), where it is (1-x**2)**(m/2).
I am not an expert in orthogonal polynomials, so perhaps I am missing a subtlety here which means that we can expect this kind of discrepancy? I would be grateful if more experienced mathematicians or Sage Math developers could weigh in on this.
If it is indeed a bug, I am happy to file a report.