Hi all, As I alluded to on my thread regarding polynomial evaluation mpmath appears to have some issues evaluating Jacobi polynomials:
In [3]: for i in range(5):
...: for j in range(5):
...: for k in range(5):
...: try:
...: mp.jacobi(i, j, k, 0)
...: except:
...: bad += 1
...:
In [4]: bad
Out[4]: 12
The implementation in mpmath currently treats Jacobi polynomials as a
special case of 2F1. However, in the case of an integer n there exists
a rather nice recurrence relation. Using this I have hacked up the
following (for integer n):
def jacobi(n, a, b, z):
half = mp.mpf(0.5)
if n == 0:
return 1
elif n == 1:
return half*((a + b + 2)*z + a - b)
else:
jqm1, jqm2 = half*((a + b + 2)*z + a - b), 1
apb, bbmaa = a + b, b*b - a*a
for q in xrange(2, n + 1):
qapbpq, apbp2q = q*(apb + q), apb + 2*q
apbp2qm1, apbp2qm2 = apbp2q - 1, apbp2q - 2
aq = half*apbp2q*apbp2qm1/qapbpq
bq = half*apbp2qm1*bbmaa/(qapbpq*apbp2qm2)
cq = apbp2q*(a + q - 1)*(b + q - 1)/(qapbpq*apbp2qm2)
# Update
jqm1, jqm2 = (aq*z - bq)*jqm1 - cq*jqm2, jqm1
return jqm1
which in my testing between [-1,1] for -1 < (n,a,b) < 5 agrees with
mp.jacobi (except for when mp.jacobi fails to converge). Longer term I
would like to get this as part of upstream mpmath. However, the issue
at the moment is how to best handle the parameters a, b and z. These
can be: integers, floats, complex floats, mp.mpf, mp.mpc. How should I
cast/treat them such that the above always does 'the right thing'. For
example, when a and b are integers the above will fail on Python 2 as cq
is integer/integer => integer (truncation). What is the usual guidance
here, should I just cast everything to an mp.mpc (albeit at a
performance cost and the need to possibly chop off an imaginary part)?
Regards, Freddie.
signature.asc
Description: OpenPGP digital signature
