(连带)勒让德多项式计算代码

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