Let sphericalharmonicy call FastTransforms.sphevaluatepi?

[MikaelSlevinsky]

I wonder how accurate and efficient the normalization is

HarmonicOrthogonalPolynomials.jl/src/HarmonicOrthogonalPolynomials.jl

Lines 175 to 178 in 4049e8f

	function sphericalharmonicy(ℓ, m, θ, φ)
	m̃ = abs(m)
	exp((lgamma(ℓ+m̃+1)+lgamma(ℓ-m̃+1)-2lgamma(ℓ+1))/2)*sqrt((2ℓ+1)/(4π)) * exp(im*m*φ) * sin(θ/2)^m̃ * cos(θ/2)^m̃ * jacobip(ℓ-m̃,m̃,m̃,cos(θ))
	end

compared to recurrence relations with normalization built-in

julia> θ, φ, ℓ, m = 0.1, 0.2, 2, 1
(0.1, 0.2, 2, 1)

julia> @time FastTransforms.sphevaluatepi(θ/π, ℓ, m)*exp(im*m*φ)/sqrt(2π)
  0.000011 seconds (7 allocations: 176 bytes)
0.07521112971423365 + 0.015246050775019672im

julia> θ, φ, ℓ, m = π/2, 0.2, 10_000, 5000
(1.5707963267948966, 0.2, 10000, 5000)

julia> @time FastTransforms.sphevaluatepi(θ/π, ℓ, m)*exp(im*m*φ)/sqrt(2π)
  0.000060 seconds (7 allocations: 176 bytes)
0.19235804128133974 + 0.2828286746396144im

[TSGut]

Here are some basic tests.

julia> using BenchmarkTools, FastTransforms

julia> θ, φ, ℓ, m = 0.1, 0.2, 2, 1
(0.1, 0.2, 2, 1)

julia> @btime FastTransforms.sphevaluatepi(θ/π, ℓ, m)*exp(im*m*φ)/sqrt(2π)
  150.510 ns (8 allocations: 208 bytes)
0.07521112971423365 + 0.015246050775019672im

julia> @btime sphericalharmonicy(ℓ, m, θ, φ)
  136.160 ns (2 allocations: 128 bytes)
0.07521112971423363 + 0.015246050775019674im

julia> abs(FastTransforms.sphevaluatepi(θ/π, ℓ, m)*exp(im*m*φ)/sqrt(2π)-sphericalharmonicy(ℓ, m, θ, φ))
1.3985787785349405e-17

There's an issue for high values of m, though:

julia> θ, φ, ℓ, m = π/2, 0.2, 10000, 5000
(1.5707963267948966, 0.2, 10000, 5000)

julia> @btime FastTransforms.sphevaluatepi(θ/π, ℓ, m)*exp(im*m*φ)/sqrt(2π)
  51.562 μs (8 allocations: 208 bytes)
0.19235804128133974 + 0.2828286746396144im

julia> # the following clearly doesn't work

julia> sphericalharmonicy(ℓ, m, θ, φ)
NaN + NaN*im

julia> # with BigFloat and BigInt it works but is slow

julia> sphericalharmonicy(BigInt(ℓ), m, BigFloat(θ), φ)
0.1923580412813398318764142897415155812528907393225858348966156677614988766433293 + 0.2828286746396146063051417153359899381451149677419309920950959124851106430082544im

julia> # this originates in the call to jacobip 
julia> θ, φ, ℓ, m = π/2, 0.2, 10000, 1000
(1.5707963267948966, 0.2, 10000, 1000)

julia> jacobip(ℓ-abs(m),abs(m),abs(m),cos(θ))
1.5242151618568348e277

julia> θ, φ, ℓ, m = π/2, 0.2, 10000, 5000
(1.5707963267948966, 0.2, 10000, 5000)

julia> jacobip(ℓ-abs(m),abs(m),abs(m),cos(θ))
NaN