module lib_legendre
!函数默认逐元
Implicit None
private
public::legendre,associated_legendre,mod_legendre!求模函数输出为双精度
interface legendre
module procedure legendre_sp!勒让德多项式，输入输出均为单精度
module procedure legendre_dp!双精度
end interface
interface associated_legendre!未除模进行归一化，值会比较大
module procedure associated_legendre_sp!单精度连带勒让德
module procedure associated_legendre_dp!双精度
end interface
contains
elemental function combination(m,n)!组合数
!C(m,n)=m!/n!(m-n)!=exp(SIGMA(i=m-n+1,m)(ln(m_i))-SIGMA(ln(n)))
integer,intent(in)::m,n
integer::combination
integer::counter
real(kind=8),allocatable::top(:),bottom(:)!在逐元函数中，虚参不能用于声明变量，除非是查询函数
allocate(top(n),bottom(n))!只能这样手动配置
top=[(counter,counter=m-n+1,m,1)]
bottom=[(counter,counter=1,n,1)]
combination=nint(exp(sum(dlog(top))-sum(dlog(bottom))))
end function combination
pure function mod_legendre(l,m) result(ans)!勒让德多项式的模
integer,intent(in)::l
integer,intent(in),optional::m
integer::counter
real(kind=8)::ans
if (present(m)) then
ans=sqrt(product([(counter,counter=l-m+1,l+m,1)])*2/real(2*l+1))
else
ans=sqrt(2/real(2*l+1))
end if
end function mod_legendre
elemental function legendre_sp(x,n) result(ans)!P_{n}(x)!n阶勒让德函数，输入-1<x<1
real(kind=4),intent(in)::x
integer,intent(in)::n
integer::limit,counter
real(kind=4),allocatable::x_vector(:),l_vector(:)
integer,allocatable::m(:)
real(kind=4)::ans
limit=floor(real(n)/2)
allocate(x_vector(0:limit),l_vector(0:limit),m(0:limit))
m=[(counter,counter=0,limit,1)] !求和指标
l_vector=real(((-1)**m)*combination(n,m)*combination(2*(n-m),n))/real(2**n)!系数，注意要用real函数把整型变实型
x_vector=x**(n-2*m)!自变量
ans=dot_product(l_vector,x_vector)
end function legendre_sp
elemental function associated_legendre_sp(x,n,m) result(ans) !P_{n}^{m}(x)=-1^{m}(1-x^{2})^{m/2}\frac{d^{m}}{dx^{m}}P_{n}(x)
real::ans !n阶连带勒让德
real , intent(in)::x
integer,intent(in)::n,m
integer::limit,counter
real(kind=4),allocatable::x_vector(:),l_vector(:)
integer,allocatable::k(:)
if (m == 0) then
ans=legendre_sp(x,n)
else
limit=floor(real(n)/2)
allocate(x_vector(0:limit),l_vector(0:limit),k(0:limit))
k=[(counter,counter=0,limit,1)] !求和指标
l_vector=real(((-1)**k)*combination(n,k)*combination(2*(n-k),n))/real(2**n)!系数
x_vector=product(reshape([(n-2*k-counter,counter=0,m-1,1)],[limit+1,m]),dim=2)*(x**(n-2*k-m))!自变量
where (n-2*k-m < 0) x_vector=0
ans=dot_product(l_vector,x_vector)*((1-x**2)**(real(m)/2))*((-1)**m)
end if
end function associated_legendre_sp
elemental function legendre_dp(x,n) result(ans)!P_{n}(x)!n阶勒让德函数，输入-1<x<1
real(kind=8),intent(in)::x
integer,intent(in)::n
integer::limit,counter
real(kind=8),allocatable::x_vector(:),l_vector(:)
integer,allocatable::m(:)
real(kind=8)::ans
limit=floor(real(n)/2)
allocate(x_vector(0:limit),l_vector(0:limit),m(0:limit))
m=[(counter,counter=0,limit,1)] !求和指标
l_vector=real(((-1)**m)*combination(n,m)*combination(2*(n-m),n))/real(2**n)!系数，注意要用real函数把整型变实型
x_vector=x**(n-2*m)!自变量
ans=dot_product(l_vector,x_vector)
end function legendre_dp
elemental function associated_legendre_dp(x,n,m) result(ans) !P_{n}^{m}(x)=-1^{m}(1-x^{2})^{m/2}\frac{d^{m}}{dx^{m}}P_{n}(x)
real(kind=8)::ans !n阶连带勒让德
real(kind=8),intent(in)::x
integer,intent(in)::n,m
integer::limit,counter
real(kind=8),allocatable::x_vector(:),l_vector(:)
integer,allocatable::k(:)
if ( m == 0 ) then
ans=legendre_dp(x,n)
else
limit=floor(real(n)/2)
allocate(x_vector(0:limit),l_vector(0:limit),k(0:limit))
k=[(counter,counter=0,limit,1)] !求和指标
l_vector=real(((-1)**k)*combination(n,k)*combination(2*(n-k),n))/real(2**n)!系数
x_vector=product(reshape([(n-2*k-counter,counter=0,m-1,1)],[limit+1,m]),dim=2)*(x**(n-2*k-m))!自变量
where (n-2*k-m < 0) x_vector=0
ans=dot_product(l_vector,x_vector)*((1-x**2)**(real(m)/2))*((-1)**m)
end if
end function associated_legendre_dp
end module lib_legendre
