Trac #33962: wrong sign for value of Legendre polynomial at 0

As discussed in [https://groups.google.com/g/sage-
devel/c/ES00LNWv6DQ/m/YLPkmAxRAgAJ this sage-devel thread],
`legendre_P(n, 0)` should be negative when `n` is congruent to 2 modulo
4, but sagemath returns a positive value:
{{{
sage: [legendre_P(n, 0) for n in range(0, 10, 2)]
[1, 1/2, 3/8, 5/16, 35/128]
}}}
The correct values are `[1, -1/2, 3/8, -5/16, 35/128]`. (The signs
should alternate when restricted to even values of `n`.)

This is a pynac bug. It only arises in the code branch where `n` is an
integer variable, so we get the correct values when `n` is real:
{{{
sage: [QQ(legendre_P(RR(n), 0)) for n in range(0, 10, 2)]
[1, -1/2, 3/8, -5/16, 35/128]
}}}

URL: https://trac.sagemath.org/33962
Reported by: gh-DaveWitteMorris
Ticket author(s): Dave Morris
Reviewer(s): Travis Scrimshaw

src/sage/functions/orthogonal_polys.py

@@ -1290,6 +1290,14 @@ class Func_legendre_P(GinacFunction):
         Traceback (most recent call last):
         ...
         RuntimeError: derivative w.r.t. to the index is not supported yet
+
+    TESTS:
+
+    Verify that :trac:`33962` is fixed::
+
+        sage: [legendre_P(n, 0) for n in range(-10, 10)]
+        [0, 35/128, 0, -5/16, 0, 3/8, 0, -1/2, 0, 1,
+         1, 0, -1/2, 0, 3/8, 0, -5/16, 0, 35/128, 0]
     """
     def __init__(self):
         r"""

src/sage/symbolic/ginac/inifcns_orthopoly.cpp

@@ -185,10 +185,17 @@ static ex legp_eval(const ex& n_, const ex& x)
                                 if (n.info(info_flags::even)) {
                                         if (is_exactly_a<numeric>(n)) {
                                                 const numeric& numn = ex_to<numeric>(n);
-                                                return (numn+*_num_1_p).factorial() / numn.mul(*_num1_2_p).factorial().pow_intexp(2) * numn / _num2_p->pow_intexp(numn.to_int());
+                                                return (numn + *_num_1_p).factorial()
+                                                        / numn.mul(*_num1_2_p).factorial().pow_intexp(2)
+                                                        * numn
+                                                        / _num2_p->pow_intexp(numn.to_int())
+                                                        * _num_1_p->pow_intexp(numn.mul(*_num1_2_p).to_int());
                                         }
-                                        return gamma(n) / pow(gamma(n/_ex2),
-                                                        _ex2) / pow(_ex2, n-_ex2) / n;
+                                        return gamma(n)
+                                                / pow(gamma(n / _ex2), _ex2)
+                                                / pow(_ex2, n - _ex2)
+                                                / n
+                                                * pow(_ex_1, n / _ex2);
                                 }
                         }
                 if (is_exactly_a<numeric>(n)