module FastLegendre

# From
# O(1) Computation of Legendre polynomials and Gauss-Legendre nodes and weights for parallel computing
# Bogaert, Ignace and Michiels, Bart and Fostier, Jan
# SIAM JOURNAL ON SCIENTIFIC COMPUTING
# 34, 2012
# http://dx.doi.org/10.1137/110855442    
#
# as mentioned https://github.com/JuliaMath/SpecialFunctions.jl/issues/124
#
# Only Float64
#
# Gives a speed up over Clenshaw, which is expected, as this is O(1), not O(n):
#
# n= 1000 Clenshaw might be faster
# n= 10_000 10x faster
# n= 100_000 10x faster
# n= 1_000_000 500x faster
# n= 10_000_000 10,000x faster
# (measurement of basis(P,n)(x) compared to Basis(P,n)(x)

using FastGaussQuadrature
using SpecialFunctions


# 4 regions of evaluation
# specialized for Float64 values
function  fastlegendre(n, x)
    n == 0 && return one(x)
    n == 1 && return x
    n < 100 &&  return  direct_32(n,x)
    x < 0  &&  return  (iseven(n) ? 1 : -1) * fastlegendre(n,-x)
    θ = acos(x)
    (n+1)*sin(θ) > 25 && return asy_32(n,x)
    asy_33(n,x)
end

## ---- n < 100 case

# precompute first 99
const xᵢs = [FastGaussQuadrature.gausslegendre(n)[1][1:(n÷2)] for n in 1:99]
const xᵢs⁻¹ = [inv.(xᵢ) for xᵢ in xᵢs]
const xᵢs2⁻¹ = [inv.(1 .- xᵢ.^2) for xᵢ in xᵢs]

# precompute
const Cn = zeros(Float64, 99)
Cn[1] = 1.0
for i in 1:49
    n = 2i
    Cn[n]  = sqrt(pi) / gamma(i + 1)  / gamma((1-n)/2)
    Cn[n+1] = (n+1) * Cn[n]
end

# use cosine or sine
function  direct_32(n,x)
    if abs(x) > 0.707
        if iseven(n)
            direct_sine(Val(true), n,x)
        else
            direct_sine(Val(false), n,x)
        end
    else
        if iseven(n)
            direct_cosine(Val(true), n,x)
        else
            direct_cosine(Val(false), n,x)
        end
    end
end


#x = cos(theta) 1 - x^2 = sintheta0
function direct_sine(::Val{true}, n,x)
    tot = one(Float64)
    xs = xᵢs2⁻¹[n]
    for aᵢ in xs
        tot *= 1 - (1-x^2)*aᵢ
    end
    tot
end

direct_sine(::Val{false}, n,x) =  x * direct_sine(Val(true),n,x)


function direct_cosine(p::Val{true}, n,x)
    tot =  Cn[n]
    xs = xᵢs⁻¹[n]
    for aᵢ in  xs
        tot *= (1 - (x*aᵢ)^2)
    end
    tot
end
direct_cosine(p::Val{false}, n,x) = x * direct_cosine(Val(true), n, x)


## ----
const Cl0s = [gamma(n+1)/gamma(n + 3/2) for n in 1:9]
τ(x) = evalpoly(1/x, (1/1, 0, -1/64, 0, 21/8192, 0, -671/524288, 0, 180323/80323134, 0,
                    -20898423/8589934592, - 7426362705/1099511627776))

function Cl0(l)
    1 <= l < 10 && return Cl0s[l]
    1/sqrt(l+3/4)*τ(l+3/4)
end

function asy_32(n,x)
    M = 17 
    θ = acos(x) 
    val = sqrt((2/pi)*(1/sin(θ))) 
    tot = 0.0
    sθ = sin(θ)
    λ = Cl0(n)
    for m in 0:M-1
        α = (n+m+1/2)*θ - (m+1/2)*pi/2
        tot += λ * cos(α)
        λ *= (m+1/2)^2/2/(m+1)/(n+m+3/2) / sθ
    end
    val * tot
end

## ---
h(n,y) = y^n * besselj(n,y)

const cns = [[1/8, -1/12],
      [11/384, -7/160, 1/160],
      [173/15360, -101/3584, 671/80640, -61/120960],
      [22931/3440640,  -90497/3870720, 217/20480, -1261/967680, 1261/29030400],
      [1319183/247726080, -10918993/454164480, +1676287/113541120, -7034857/2554675200,
       +1501/8110080, -79/20275200],
      [233526463/43599790080, -1396004969/47233105920, +2323237523/101213798400,
       -72836747/12651724800, +3135577/5367398400, -1532789/61993451520, +66643/185980354560]
      ]

function fn(i,y)
    n = 2i
    cs = cns[i]
    tot = 0.0
    Δ = i-1
    for j in i:2i
        tot += cs[j-Δ]*h(j,y)
    end
    tot
end
    

function asy_33(n,x)
    v = n + 1/2
    θ = acos(x)
    y = v * θ

    # Eq (3.19)
    val = besselj(0,y)

    
    λ = 1/v^2
    for i in 1:6
        val += fn(i,y)*λ
        λ /= v^2
    end
    val

    
end

end
