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)
default(fmt=:png)
#clibrary(:misc)

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

# 1/1+1/2+...+1/n と log n の比較
# nを大きくするとそれらの差はEuler定数 γ = 0.5772… に近付く
# ゆえに, 1/1+1/2+...+1/n は log n + γ でよく近似される.
# 以下のプロットでもそれらのグラフはほとんど一致している.

f(x) = sum(k->1/k, 1:floor(Int, x))
n = 1:100
x = 1:0.1:100
plot(size=(500,350), legend=:bottomright)
plot!(n, f.(n), label="1/1+1/2+...+1/n", lw=2.5)
plot!(x, log.(x), label="log(x)", ls=:dashdot, color=:orange)
plot!(x, log.(x.+1), label="log(x+1)", ls=:dashdot, color=:darkgreen)
plot!(x, log.(x).+eulergamma, label="log(x) + Euler's gamma", lw=2, ls=:dash, color=:red)

?eulergamma

a = 1.01
l = 100
n = 1:10^3:2*10^5
f(n) = l*log(n)
g(n) = n*log(a)
plot(size=(500, 350), legend=:topleft, xlabel="n")
plot!(n, f.(n), label="$l*log(n)")
plot!(n, g.(n), label="n*log($a)")

f(x) = x*log(x)
x = 0.001:0.001:1
plot(x, f.(x), xlim=(-0.1,1), ylim=(-0.4,0), label="y = x log x", legend=:top)

lfact(n) = lgamma(n+1) # = log n!
f1(n) = n*log(n)
f2(n) = n*log(n) - n
f3(n) = n*log(n) - n + 0.5*log(n)
f4(n) = n*log(n) - n + 0.5*log(n) + 0.5*log(2π)

function plot_logfactorial(n)
    plot(legend=:topleft)
    plot!(n, lfact.(n), label="log n!", lw=1.5)
    plot!(n, f1.(n), label="n log n", ls=:dash)
    plot!(n, f2.(n), label="n log n - n", ls=:dash)
    plot!(n, f3.(n), label="n log n - n + 0.5 log n", ls=:dash)
    plot!(n, f4.(n), label="n log n - n + (1/2)log n + (1/2)log 2pi", color=:red, ls=:dash)
end

plot_logfactorial(1:20)

plot_logfactorial(10:10:1000)

plot_logfactorial(10^16:10^16:10^18)

nmax = 100
n = 1:nmax
plot(n, @.(10/n + (-1)^n), label="a_n", ylims=(-1.5,10), marker=:o, markersize=3, markerstrokealpha=0)
plot!(n, fill(1, size(n)), label="1")

nmax = 200
n = 1:nmax
plot(n, @.(10/n + sin(4*n)), label="b_n", ylims=(-1.5, 10), marker=:o, markersize=2, markerstrokealpha=0)
plot!(n, fill(1, size(n)), label="1")