def jacobi(n, alpha, beta, monic=False):
if n < <NUM_LIT:0>:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>wfunc = lambda x: (<NUM_LIT:1> - x)**alpha * (<NUM_LIT:1> + x)**beta<EOL>if n == <NUM_LIT:0>:<EOL><INDENT>return orthopoly1d([], [], <NUM_LIT:1.0>, <NUM_LIT:1.0>, wfunc, (-<NUM_LIT:1>, <NUM_LIT:1>), monic,<EOL>eval_func=np.ones_like)<EOL><DEDENT>x, w, mu = j_roots(n, alpha, beta, mu=True)<EOL>ab1 = alpha + beta + <NUM_LIT:1.0><EOL>hn = <NUM_LIT:2>**ab1 / (<NUM_LIT:2> * n + ab1) * _gam(n + alpha + <NUM_LIT:1>)<EOL>hn *= _gam(n + beta + <NUM_LIT:1.0>) / _gam(n + <NUM_LIT:1>) / _gam(n + ab1)<EOL>kn = _gam(<NUM_LIT:2> * n + ab1) / <NUM_LIT>**n / _gam(n + <NUM_LIT:1>) / _gam(n + ab1)<EOL>p = orthopoly1d(x, w, hn, kn, wfunc, (-<NUM_LIT:1>, <NUM_LIT:1>), monic,<EOL>lambda x: eval_jacobi(n, alpha, beta, x))<EOL>return p<EOL>
Returns the nth order Jacobi polynomial, P^(alpha,beta)_n(x)
    orthogonal over [-1,1] with weighting function
    (1-x)**alpha (1+x)**beta with alpha,beta > -1.