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
default(fmt = :png)
#gr(); ENV["PLOTS_TEST"] = "true"
#clibrary(:colorcet)
#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

@show 3 + 1/(2*3)
@show 3 + 1/(2*3) - 1/(8*9*3)
@show 3 + 1/(2*3) - 1/(8*10*4)
@show √10;

x = symbols("x")
series((1+x)^(Sym(1)/2), n=5)

3/2/9, -3/8/9^2, 3/16/9^3

3 + 3/2/9 -3/8/9^2, 3/16/9^3

√(10)

factorial(10)

x = symbols("x")
sympy.init_printing(order="rev-lex")
for k in 0:5
    expand(diff(tan(x), x, k)) |> display
end
sympy.init_printing(order="lex") # default

x = symbols("x")
series(tan(x), n=10)

s = series(sin(x), n=7)
cinv = series(1/cos(x), n=6)
display(s)
display(cinv)
expand(s * cinv)

x = symbols("x")
sympy.init_printing(order="rev-lex")
for k in 0:5
    expand(diff(tanh(x), x, k)) |> display
end
sympy.init_printing(order="lex") # default

x = symbols("x")
series(tanh(x), n=10)

s = series(sinh(x), n=7)
cinv = series(1/cosh(x), n=6)
display(s)
display(cinv)
expand(s * cinv)

sum(1/factorial(k) for k in 0:9), exp(1)

f(x) = log((1+x)/(1-x))
g(n,x) = 2*sum(k->x^(2k+1)/(2k+1), 0:n-1)
R(n,x) = 2*x^(2n+1)/(2n+1)/(1-x^2)

sympy.init_printing(order="rev-lex")
x = symbols("x", positive=true)
g(2, x) |> display
sympy.init_printing(order="lex") # default

n = 2
y = 2
x = (y-1)/(y+1)
log(y), f(x), g(n,x), f(x) - g(n,x), R(n,x)

f(x) = log((1+x)/(1-x))
g(n,x) = 2*sum(k->x^(2k+1)/(2k+1), 0:n-1)
R(n,x) = 2*x^(2n+1)/(2n+1)/(1-x^2)

n = 3
y = 2
x = (y-1)/(y+1)
log(y), f(x), g(n,x), f(x) - g(n,x), R(n,x)

f(x) = log((1+x)/(1-x))
g(n,x) = 2*sum(k->x^(2k+1)/(2k+1), 0:n-1)
R(n,x) = 2*x^(2n+1)/(2n+1)/(1-x^2)

n = 2
y = 1.5
x = (y-1)/(y+1)
log(y), f(x), g(n,x), f(x) - g(n,x), R(n,x)

f(x) = log((1+x)/(1-x))
g(n,x) = 2*sum(k->x^(2k+1)/(2k+1), 0:n-1)
R(n,x) = 2*x^(2n+1)/(2n+1)/(1-x^2)

n = 3
y = 1.5
x = (y-1)/(y+1)
log(y), f(x), g(n,x), f(x) - g(n,x), R(n,x)

f(x) = x == 0 ? zero(x) : exp(-1/abs(x))
x = -2:0.001:2
plot(x, f.(x), label="y = f(x)", size=(400,250), legend=:top)

x = symbols("x")
series(e^x, n=10)

series(cos(x), n=10)

series(sin(x), n=10)

series(cosh(x), n=10)

series(sinh(x), n=10)

series(log(1+x), n=10)

series(-log(1-x), n=10)

series(asin(x), n=10)

k = symbols("k", integer=true)
doit(sympy.Sum(1/2^(2k)*factorial(2k)/(factorial(k)^2)*x^(2k+1)/(2k+1), (k,0,4)))

series(atan(x), n=10)

series(atanh(x), n=10)

sympy.Sum((-1)^k/(2k+1), (k,0,oo))

doit(sympy.Sum((-1)^k/(2k+1), (k,0,oo)))

f(n,k) = exp(lgamma(n+1)-lgamma(k+1)-lgamma(n-k+1)-n*log(2))
@time sum(k->f(2k,k)/(2k+1), 0:10^8-1), π/2

# Bernoulli数 B_{2k} の例

BernoulliNumber(n) = sympy.bernoulli(n)
B = [BernoulliNumber(2k) for k in 0:8]

# Euler数 E_{2k} の例

EulerNumber(n) = sympy.euler(n)
E = [EulerNumber(2k) for k in 0:8]

z = symbols("z")
series((z/2)*cot(z/2), z, n=10)

z = symbols("z")
series(cot(z), z, n=10)

z = symbols("z")
series(tan(z), z, n=10)

z = symbols("z")
series(csc(z), z, n=10)

z = symbols("z")
series(sec(z), z, n=10)

# 上のセルのべき級数の数値計算例

F(x; N=10^3) = sum(n->(-1)^n*x^(2.0^n), 0:N)
@show F(0.9)
@show F(0.99)
@show F(0.999)
@show F(0.9999)
@show F(0.99999)
@show F(0.999999);

# 上の f(x) が x→1 で収束しないことの数値的確認

F(x; N=10^3) = sum(n->(-1)^n*x^(2.0^n), 0:N)
x = 0:0.001:0.999
P1 = plot(x, F.(x), label="F(x)", ylims=(0,0.5), xlims=(0,1))
x = 0.999:0.000001:0.999999
P2 = plot(x, F.(x), label="F(x)", ylims=(0.497, 0.503), xlims=(0.999,1))
plot(P1, P2, size=(600, 240), legend=:topleft)

j, n = symbols("j n", integer=true)
doit(sympy.Sum((-1)^(j-1)*j^2, (j,1,n)))

aa(n) = n == floor(typeof(n), √n)^2 ? one(n) : zero(n)
ss(n) = sum(k->(-1)^k*aa(k), 0:n)
tt(n) = sum(k->ss(k), 0:n)
sigma(n) = tt(n)/(n+1)
N = 1000
n = 0:N
plot(size=(640, 180), xlim=(-0.005*N, 1.005*N))
@time plot!(n, sigma.(n), label="sigma_n")
hline!([1/2], label="1/2", ls=:dash)

f(x; N=2^5) = sum(j->(-1)^j*x^(typeof(x)(2)^j), 0:N)
@time [f(big"0.995"; N=2^m) for m in 0:5]

x = -1:0.005:1
y = -1:0.005:1
z = x' .+ im.*y
f(z, R) = (abs2(z)<1 && abs(1-z) ≤ R*(1-abs(z))) ? 0.0 : 1.0
t = 0:0.01:2π

PP = []
for RR in [1.0, 2.0, 5.0, 10.0]
    P = plot(aspectratio=1)
    plot!(xlim=(-1,1), ylim=(-1,1))
    plot!(title="R = $RR", titlefontsize=10)
    heatmap!(x, y, f.(z, RR), colorbar=false, color=:plasma)
    plot!(cos.(t), sin.(t), label="")
    push!(PP, P)
end

plot(PP..., size=(750, 170), layout=@layout([a b c d]))

f(z; N=39) = sum(k->(z^(2*3^k)-z^(3^k))/k, 1:N)
@show f(1)
@show f(0.9*e^(π*im/3^1))
@show f(0.9999*e^(π*im/3^2))
@show f(0.9999999*e^(π*im/3^3))
@show f(0.9999999999*e^(π*im/3^4))
@show f(0.9999999999999*e^(π*im/3^5));