Repository: scipy/scipy
File: scipy/special/tests/test_orthogonal_eval.py (Python)

Code before:
        (2, [3.0, 2.16, 1.44, 0.84, 0.36, 0.0, -0.24, -0.36, -0.36, -0.24, 0.0]),
        (3, [-4.0, -1.98, -0.64, 0.14, 0.48, 0.5, 0.32, 0.06, -0.16, -0.22, 0.0]),
        (
            4,
            [5.0, 1.332, -0.288, -0.658, -0.408, 0.0, 0.272, 0.282, 0.072, -0.148, 0.0],
        ),
        (
            5,
            [
                -6.0,
                -0.43308,
                0.79104,
                0.36876,
                -0.21312,
                -0.375,
                -0.14208,
                0.15804,
                0.19776,
                -0.04812,
                0.0,
            ],
        ),
    ],
)
def test_jacobi_alpha_minus_one_beta_plus_one(n, expected):
    # gh-7001 - expected values were computed with mathematica
    x = np.linspace(-1.0, 1.0, 11)
    a, b = -1, 1  # alpha, beta
    assert_allclose(_ufuncs.eval_jacobi(n, a, b, x), expected)
    assert_allclose(_ufuncs.eval_jacobi(float(n), a, b, x), expected)


Review comment (suggestion):
Just wondering -- are we sure we should close the matching ticket here? I'm not saying that things are not getting better with your patch, but isn't the original point of contention that i.e., `scipy.special.eval_jacobi(2, -1, -1, 0)` produces `nan` (which it still does in this branch, and I think you also indicate as much at https://github.com/scipy/scipy/issues/7001#issuecomment-3647978502)?

That particular case/complaint in the original ticket still produces the allegedly-expected result with Wolfram Alpha,  with `JacobiP[2, -1, -1, 0]` returning `-0.25`. Do you agree that is the desired result?

Since you're labelling this an enhancement and not changing any old tests, it may very well be a step forward, I'm just thinking we might want to keep the cognate ticket open if one of the main concerns remains.

I did check locally that on the new test you've added here, many other `nan` return values/cases *do* appear to be getting fixed.

Code after:
        (2, [3.0, 2.16, 1.44, 0.84, 0.36, 0.0, -0.24, -0.36, -0.36, -0.24, 0.0]),
        (3, [-4.0, -1.98, -0.64, 0.14, 0.48, 0.5, 0.32, 0.06, -0.16, -0.22, 0.0]),
        (
            4,
            [5.0, 1.332, -0.288, -0.658, -0.408, 0.0, 0.272, 0.282, 0.072, -0.148, 0.0],
        ),
        (
            5,
            [
                -6.0,
                -0.43308,
                0.79104,
                0.36876,
                -0.21312,
                -0.375,
                -0.14208,
                0.15804,
                0.19776,
                -0.04812,
                0.0,
            ],
        ),
    ],
)
def test_jacobi_alpha_minus_one_beta_plus_one(n, expected):
    # gh-7001 - expected values were computed with mathematica
    x = np.linspace(-1.0, 1.0, 11)
    a, b = -1, 1  # alpha, beta
    assert_allclose(_ufuncs.eval_jacobi(n, a, b, x), expected, rtol=1e-10, atol=1e-14)
    assert_allclose(
        _ufuncs.eval_jacobi(float(n), a, b, x), expected, rtol=1e-10, atol=1e-14
    )


@pytest.mark.parametrize(
    "n, expected",
    [
        (0, [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]),
        (1, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
        (2, [0.0, -0.09, -0.16, -0.21, -0.24, -0.25, -0.24, -0.21, -0.16, -0.09, 0.0]),
        (
            3,
            [0.0, 0.144, 0.192, 0.168, 0.096, 0.0, -0.096, -0.168, -0.192, -0.144, 0.0],
        ),
        (
            4,
            [
                0.0,
                -0.1485,
                -0.096,
                0.0315,