using Base.MathConstants
using Base64
using Printf
using Statistics
const e = ℯ
endof(a) = lastindex(a)
linspace(start, stop, length) = range(start, stop, length=length)

using Plots
#gr(); ENV["PLOTS_TEST"] = "true"
#clibrary(:colorcet)
#clibrary(:misc)
default(fmt=:png)

function pngplot(P...; kwargs...)
    sleep(0.1)
    pngfile = tempname() * ".png"
    savefig(plot(P...; kwargs...), pngfile)
    showimg("image/png", pngfile)
end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))

showimg(mime, fn) = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<img src="data:$mime;base64,$base64">""")
end

using SymPy
#sympy.init_printing(order="lex") # default
#sympy.init_printing(order="rev-lex")

using SpecialFunctions
using QuadGK

# Override the Base.show definition of SymPy.jl:
# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/types.jl#L87-L105

@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
    print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_frac=false)))
end
@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
    function toeqnarray(x::Vector{Sym})
        a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")
        """\\left[ \\begin{array}{r}$a\\end{array} \\right]"""
    end
    function toeqnarray(x::AbstractArray{Sym,2})
        sz = size(x)
        a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz[1]], "\\\\")
        "\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
    end
    print(io, as_markdown(toeqnarray(x)))
end

[(n, factorial(n)*(1/10)^(n+1)) for n in 7:13]

F₀(x) = exp(1/x)*quadgk(u->exp(-1/(x*u))/u, 0, 1)[1]
F₀_ae(x, n=10) = sum(k->(-1)^k*factorial(k)*x^(k+1), 0:n-1)

@show Y = F₀(1/10)
@show Y_ae9 = F₀_ae(1/10, 9)
@show Y_ae10 = F₀_ae(1/10, 10)

n = 1:20
plot(size=(500,350))
plot!(ylims=(-0.0001,0.0001))
plot!(legend=:top)
plot!(n, F₀_ae.(1/10, n) .- Y, label="error")

function G_cf(z; n::Int=2)
    cf = 1 + (n+1)/z
    for i = n:-1:1
        cf = z + (1+i)/cf
        cf = 1 + i/cf
    end
    return 1 / (z + 1/cf)
end

G(z) = exp(z)*quadgk(u->exp(-z/u)/u, 0, 1)[1]

sympy.init_printing(order="lex")
z = symbols("z")
for n in 1:3
    cf = G_cf(z, n=n)
    display(cf)
end

sympy.init_printing(order="rev-lex")
display(series(sum((-1)^n*factorial(n)*z^(n+1) for n in 0:9), n=11))
for n in 1:3
    display(series(G_cf(1/z, n=n), n=11))
end

G_cf(4), G(4)

x = 0.05:0.05:4
@time y = G.(x)
@time y_cf1 = G_cf.(x, n=1)
@time y_cf2 = G_cf.(x, n=2)
@time y_cf3 = G_cf.(x, n=3)
plot(size=(500,350))
plot!(title="y = G(x) = e^x E_1(x)", titlefontsize=11)
plot!(xlims=(0,maximum(x)), ylims=(0,2.3))
plot!(x, y, label="numetical integration")
plot!(x, y_cf1, label="continued fractional approximation n=1", ls=:dash)
plot!(x, y_cf2, label="continued fractional approximation n=2", ls=:dash)
plot!(x, y_cf3, label="continued fractional approximation n=3", ls=:dash)