I have some special functions. I would like to deposit them into someone’s existing / new projects. An example:
function mylegenp3(n,m,z)# n,m integers , z > 1,associated Legendre
v=(z*z - big(1.))^(-m*big(1)/2)
v=v*factorial(big(n+m))/(((big(2))^n)* factorial(big(n)));
term = big(0.);
MXP=Int(trunc(((n+m)/2)))
for p=0:MXP
term= term +
((-big(1.))^p)*binomial(big(2*(n-p)),big(n-m))*
(z^(big(1)*(n + m -2*p)))*
binomial(big(n),big(p))
end
return term * v
end
Could you please give some references to your implementation? For instance: was any stability analysis done for this algorithm? term has contributions with alternating signs.
And for those like me who are not familiar with associated Legendre of type 3, what are they?
Note, however, that it looks like your code can be substantially improved. For on thing, the integer type used should depend on the type of the arguments—you shouldn’t unconditionally convert to BigInt. Also, when you’re computing polynomial series, you almost never want to compute each term in the series independently (with calls to ^, factorial, binomial etcetera). Instead, you want to compute each term as a recurrence from the previous term, analogous to Horner’s method. (This also helps to avoid overflow from the individual terms in ratios of factorials.)
And you should always look at the literature on specific special functions to see if there are better ways to compute them than just plugging in one of the definitions. For example, for associated Legendre polynomials there are recurrence relations that I’m guessing are a better method to evaluate them.
Another good thing to do is to compare to existing implementations, like the ones in GSL or SciPy. If your code is much slower then you are probably making a mistake.
References are Abramowitz and Stegun, Handbook of Mathematical Functions, also online Mathematica function reference, wikipedia, and python implementation.
I did unconditionally convert to BigInt because when dealing with factorials, things can get humungous. I used arbitrary precision for accuracy, not speed. Faster computation would use recurrence and your methods.
The trick is to arrange a recurrence carefully so that you can avoid overflow without resorting to bignum arithmetic. (Even though the individual factorial terms can be huge, the overall term in the series is usually easily within floating-point range.)
Performance is a big concern with special function implementations, because the whole point of using special functions is usually to gain performance over generic methods like quadrature etc. So it is really useful to compare your implementation to other implementations.
Matlab doesn’t calculate type 3 (z>1) function but does type 2 (-1<=x<.=1) function ; timing is approximately .0007sec compared to average .00014sec recNM3(n,m,z) function
function recNM3(n,m,z) #n,m positive integers, z number #reliable for values of n <= 17 for z=2.
# 0 <= m <= n
# z > 1.
M = (2m -1) # M must be odd
# (2m-1)!!= (2m)!/( m! 2^m) double factorial
dblfac=1
for j=1:M
if iseven(j)
continue
end
dblfac=jdblfac
end
if n == m
return( dblfac*(zz -1)^(m/2))
elseif n == m+1
return((2.m +1.)zdblfac(zz -1)^(m/2))
end
pj2=dblfac*(zz -1.)^(m/2)
pj1=z(2.*m+1.)pj2
for j = m+2 :n
pjj=(z(2.*j-1.)pj1 - pj2(j +m-1.)) /(j-m)
pj2=pj1
pj1=pjj
end
return pj1
end
I combined 8.6.6,.8.6.18, 8.2.5 (Abramaowitz & Stegun), used binomial expansion and differentiated to get a Finite sum.
Markdown :
A representation of P|σ|μ(z)(z2−1)|σ|/2 is used by combining 8.6.6,8.6.18,8.2.5 (Ref. 14),
P|σ|μ(z)=(μ+|σ|)!(z2−1)−|σ|/22μμ!(μ−|σ|)!dμ−|σ|dzμ−|σ|(z2−1)μ using binomial expansion and differentiating} ;
P|σ|μ(z)(z2−1)|σ|/2=(μ+|σ|)!2μμ!∑[μ+|σ|2]p=0(−)p(2μ−2pμ−|σ|)(μp)zμ+|σ|−2p
