using Logging; disable_logging(Logging.Warn)
using Printf
using Base64

using Plots
pythonplot(fmt=:png)

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

using SymPy
using LaTeXStrings
using SpecialFunctions
using QuadGK
using Elliptic.Jacobi: cd, sn
using LinearAlgebra: det

# 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

showimg("image/jpeg", "images/trigonometric1.jpg", scale="40%")

showimg("image/jpeg", "images/trigonometric2.jpg", scale="40%")

showimg("image/jpeg", "images/trigonometric3.jpg", scale="80%")

showimg("image/jpeg", "images/trigonometric4.jpg", scale="80%")

# 1の原始3乗根

x, y, z = symbols("x y z")
ω₃ = factor((-1+√Sym(-3))/2)
latexstring(raw"\ds ω =", sympy.latex(ω₃))

# 因数分解の公式の確認

ω = symbols("ω")
f3_factored = (x+y+z)*(x+ω*y+ω^2*z)*(x+ω^2*y+ω*z)
f3 = simplify(f3_factored(ω=>ω₃))
latexstring(sympy.latex(f3_factored), "=", sympy.latex(f3))

pp = y*z
qq = y^3+z^3
latexstring("p=", sympy.latex(pp)) |> display
latexstring("q=", sympy.latex(qq)) |> display

p, q = symbols("p q")
latexstring(sympy.latex(x^3 - 3p*x + q), "=", sympy.latex((x^3 - 3pp*x + qq)))

# (m-y^3)(m-z^3) の係数

m = symbols("m")
f2_org = (m-y^3)*(m-z^3)
f2 = expand((m-y^3)*(m-z^3))
latexstring(sympy.latex(f2_org), "=", sympy.latex(f2)) |> display
c2 = sympy.Poly(f2, m).all_coeffs()
latexstring(raw"\text{coefficients:}\ ", sympy.latex(c2))

# 2次方程式の解

p, q = symbols("p q")
equ = m^2-q*m+p^3
sol = factor.(solve(equ, m))
latexstring(sympy.latex(equ), "=0") |> display
latexstring(raw"\ds m = ", sympy.latex(sol[1]), ",", sympy.latex(sol[2]))

# 3次方程式の解が得られていることの確認

y = sol[2]^(Sym(1)/3)
z = p/y
X = -[
    y+z
    ω*y+ω^2*z
    ω^2*y+ω*z
]
equ = @. X^3 - 3p*X + q
res = @.(simplify((f->f(ω=>ω₃))(equ)))
for i in 1:3
    latexstring("x_{$i}^3-3px_{$i}+q=", sympy.latex(res[i])) |> display
end
latexstring(raw"\ds x_1 = ", sympy.latex(X[1])) |> display
latexstring(raw"\ds x_2 = ", sympy.latex(X[2])) |> display
latexstring(raw"\ds x_3 = ", sympy.latex(X[3])) |> display

# 3次方程式の解 (SymPyのsolve函数による直接計算)

x, p, q = symbols("x p q")
equ = x^3-3p*x+q
sol = solve(equ, x)
latexstring("x^3-3px+q=0") |> display
display(L"x = x_1,x_2,x_3")
latexstring(raw"\ds x_1 = ", sympy.latex(sol[1])) |> display
latexstring(raw"\ds x_2 = ", sympy.latex(sol[2])) |> display
latexstring(raw"\ds x_3 = ", sympy.latex(sol[3])) |> display

# 展開結果の確認

w, x, y, z = symbols("w x y z")
f4_org = (w+x+y+z)*(w+x-y-z)*(w-x+y-z)*(w-x-y+z)
f4 = expand(f4_org)
latexstring(sympy.latex(f4_org)) |> display
latexstring("=", sympy.latex(f4)) |> display

# 展開結果の 0,1,2,3,4 次の係数の確認.
# 3次の係数は0になっている.

display(L"\text{coefficients:}\ ")
c4 = sympy.Poly(f4, w).all_coeffs()

p = c4[3]/(-2)
latexstring("p=", sympy.latex(p))

q = c4[4]/8
latexstring("q=", sympy.latex(q))

r = x^2*y^2 + x^2*z^2 + y^2*z^2
@show simplify(c4[5] - (p^2-4r))
latexstring("r=", sympy.latex(r))

# p,q,rを使った公式の確認

ff4 = w^4 - 2p*w^2 + 8q*w + (p^2-4r)
simplify(f4 - ff4) |> display
latexstring("w^4-2pw^2+8qw+(p^2-4r)") |> display
latexstring("=", sympy.latex(f4)) |> display

m = symbols("m")
equ = m^3 - p*m^2 + r*m - q^2
sol = simplify(factor(equ))
latexstring(sympy.latex(equ), "=", sympy.latex(sol))

# 補助的な3次方程式の形の確認

p, q, r, m, M = symbols("p q r m M")
(f3 = m^3 - p*m^2 + r*m - q^2) |> display
(g3 = simplify(expand(f3(m => M+p/3)))) |> display
c3 = sympy.Poly(g3, M).all_coeffs()

# w, p, q, r = symbols("w p q r")
# factor.(solve(w^4 - 2p*w^2 + 8q*w + (p^2-4r), w))

A = [
    w x y z
    x w z y
    y z w x
    z y x w
]

d4 = det(A)

d4 == f4

sympy.Poly(d4, w).all_coeffs()

# ベルヌイ多項式のリスト k = 0,1,2,…,12

B(k,x) = sympy.bernoulli(k,x)
SS(k,x) = (B(k+1,x+1) - B(k+1,Sym(1)))/(k+1)
x = symbols("x")
for k in 0:12
    latexstring("B_{$k}(x)=", sympy.latex(B(k,x))) |> display
end

# べき乗和の公式のリスト k = 0,1,2,…,12

x = symbols("x")
for k in 0:12
    latexstring("S_{$k}(x)=", sympy.latex(factor(SS(k,x)))) |> display
end

# a(a-1)…(a-r+1) の和

a = symbols("a", integer=true)
n = symbols("n", integer=true, positive=true)
ff(a,r) = iszero(r) ? typeof(a)(1) : prod(a-i for i in 0:r-1)

for r in 0:4
    s = sympy.Sum(ff(a,r), (a,0,n))
    latexstring(sympy.latex(s), "=", sympy.latex(s.doit().factor())) |> display
end

# Σ_{a=0}^n a^k を a(a-1)…(a-r+1) の和に帰着する方法で計算

stirlingsecond(n,k) = sympy.functions.combinatorial.numbers.stirling(n, k, kind=2)
n = symbols("n", integer=true, positive=true)
SSS(k,n) = expand(sum(stirlingsecond(k,r)*ff(n+1, r+1)/(r+1) for r in 0:k))

for k in 0:12
    latexstring(raw"\sum_{j=0}^n", "j^{$k} =", sympy.latex(factor(SSS(k,n)))) |> display
end

struct AkiyamaTanigawa{T,S<:Integer}
    a::Array{T,2}
    N::S
end

Base.length(AT::AkiyamaTanigawa) = AT.N

function AkiyamaTanigawa(a0::AbstractArray{T,1}) where T
    N = size(a0,1)
    a = zeros(eltype(a0), N, N)
    a[0+1,:] = a0
    for n in 0:N-2
        for m in 0:N-n-2
            a[(n+1)+1,m+1] = (m+1) * (a[n+1, m+1] - a[n+1, (m+1)+1])
        end
    end
    AkiyamaTanigawa(a, N)
end

function displayAT(AT::AkiyamaTanigawa)
    N = length(AT)
    s = raw"\begin{matrix}" * "\n"
    for n in 0:N-1
        s *= prod(x * raw"&" for x in @.(sympy.latex(Sym(AT.a[n+1,1:N-n]))))
        s = replace(s, r"&$"=>"")
        s *= "& "^n * raw"\\\\" * "\n"
    end
    s *= raw"\end{matrix}" * "\n"
    l = latexstring(raw"\displaystyle", s)
    display(l)
end

L = 10
a0 = 1 .// collect(1:L+1)
AT = AkiyamaTanigawa(a0)
displayAT(AT)

[B(n,1) for n in 0:L]

a, b, c = 0.8, 0.5, 1.0
f(x,y) = a*x + b*y + c
x = -10:0.1:5
y = -10:0.1:15
surface(x, y, f.(x',y), colorbar=false, size=(400,300))
plot!(x, @.(-(a*x + c)/b), zero(x); label="", c=:black, lw=0.5)
plot!(; ylim=extrema(y))
title!("\$z = $a x + $b y + $c\$")

a, b, c = 0.8, 0.5, 1.0
f(x,y) = a*x + b*y + c
x = -10:0.1:5
y = -10:0.1:15
surface(x, y, abs.(f.(x',y)), colorbar=false, size=(400,300))
plot!(title="\$z = |$a x + $b y + $c|\$")

# log x は上に凸な函数

x = 0:0.01:2.0
a, b = 0.3, 1.5
f(x) = log(x)
t = 0:0.01:1.0
g(a,b,t) = (1-t)*f(a) + t*f(b)
h(a,b,t) = (1-t)*a + t*b
plot(size=(400,250), legend=:topleft, xlims=(0,2.0), ylims=(-2.0, 0.8))
plot!(x, f.(x), label=L"y = \log\,x")
plot!(h.(a,b,t), g.(a,b,t), label="")
plot!([a,a], [-10.0, f(a)], label=L"x = a", ls=:dash)
plot!([b,b], [-10.0, f(b)], label=L"x = b", ls=:dashdot)

# exp(x) は下に凸な函数

x = -1:0.01:2
a, b = -0.3, 1.5
f(x) = exp(x)
t = 0:0.01:1.0
g(a,b,t) = (1-t)*f(a) + t*f(b)
h(a,b,t) = (1-t)*a + t*b
plot(size=(400,250), legend=:topleft, xlims=(-1,2), ylims=(0,8))
plot!(x, f.(x), label=L"y = e^x")
plot!(h.(a,b,t), g.(a,b,t), label="")
plot!([a,a], [-0.0, f(a)], label=L"x = a", ls=:dash)
plot!([b,b], [-0.0, f(b)], label=L"x = b", ls=:dashdot)

# 上に凸な函数 f(x) = log x の接線

x = 0:0.01:2.0
μ = 0.7
f(x) = log(x)
g(μ,x) = (1/μ)*(x-μ) + f(μ)
plot(size=(500,350), legend=:topleft, xlims=(0,1.7), ylims=(-2.2, 0.8))
plot!(x, f.(x), label=L"y = \log\,x")
plot!(x, g.(μ,x), label=L"y = a(x-\mu)+f(\mu)", ls=:dashdot)
plot!([μ, μ], [-3.0, f(μ)], label=L"x=\mu", ls=:dash)
plot!(size=(400,250))

# M_p(x,y) = 1 のプロット

f(p,x) = iszero(p) ? 1/x : (2 - x^p)^(1/p)
P = plot(size=(500,500))
ps = [-10, -2, -1, 0, 1, 2, 10]
for p in ps
    a = iszero(p) ? 1/3 : 2^(1/p)
    if p > 0
        Δx = a/1000
        x = Δx:Δx:a
    else
        Δx = 3/1000
        x = a+eps():Δx:3
    end
    plot!(x, f.(p,x), label="p = $p", ls=:auto)
end
plot(P, xlim=(0,3), ylim=(0,3), size=(300,300))

# α_i > π ならば周長と面積をより大きくできること

t = range(0, 2π, length=400)
theta = 2π * [0, 3/20, 14/20, 16/20, 18/20, 1]
phi = copy(theta)
phi[2] = theta[3] - π
phi[3] = theta[2] + π
plot(size=(300,300), aspect_ratio=1, legend=false)
plot!(cos.(t), sin.(t), lw=0.5, color=:black)
plot!(cos.(theta), sin.(theta), color=:blue)
plot!(cos.(phi), sin.(phi), ls=:dash, color=:red)

showimg("image/jpeg", "images/Jikkyo20140125limitsinc.jpg", scale="60%")

showimg("image/jpeg", "images/Jikkyo20140125ArcLength1.jpg", scale="60%")

showimg("image/jpeg", "images/Jikkyo20140125ArcLength2.jpg", scale="60%")

t = symbols("t")
x = √(1-t^2)
y = t
equ = diff(x, t)^2 + diff(y, t)^2
sol = simplify(equ)
latexstring("x=", sympy.latex(x), raw",\quad y=", sympy.latex(y)) |> display
latexstring(raw"\ds\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2=", sympy.latex(sol))

t, x, y = symbols("t x y")
equ = [y-t*x, x^2+y^2-1]
s = solve(equ, [x,y])
latexstring(raw"\text{euqation:}\; ", sympy.latex(equ[1]), raw"=0,\ ", sympy.latex(equ[2]), "=0") |> display
latexstring(raw"\text{solution:}\; x = ", sympy.latex(s[2][1]), ", y = ", sympy.latex(s[2][2])) |> display

X, Y= s[2][1], s[2][2]
sol = simplify(diff(X,t)^2 + diff(Y,t)^2)
latexstring(raw"\ds X=", sympy.latex(X), raw",\quad Y=", sympy.latex(Y)) |> display
latexstring(raw"\ds\left(\frac{dX}{dt}\right)^2+\left(\frac{dY}{dt}\right)^2=", sympy.latex(sol))

showimg("image/jpeg", "images/sin-sinh.jpg", scale="80%")

showimg("image/jpeg", "images/hyperbolicsine.jpg", scale="80%")

# Edwards曲線の公式の確認

k, y = symbols("k y")
x² = (1-y^2)/(1-k^2*y^2)
y² = y^2
simplify(x² + y² - 1 - k^2*x²*y²)

# Edwards曲線の二通りのプロット

k² = -20

y = -1:0.01:1
x = @. √((1-y^2)/(1-k²*y^2))
P1 = plot(title="\$x^2+y^2=1+k^2x^2y^2\$ for \$k^2=$k²\$", titlefontsize=10)
plot!(aspectratio=1, legend=false)
plot!(x, y, color=:red)
plot!(-x, y, color=:red)

u = 0:0.01:3
P2 = plot(title="(cd u, sn u) for \$k^2=$k²\$", titlefontsize=10)
plot!(aspectratio=1, legend=false)
plot!(cd.(u,k²), sn.(u,k²))

plot(P1, P2, size=(500, 260))

# (cd u, sn u)のプロット

k² = -20
u = -5:0.01:5
plot(title="Jacobi's elliptic functions for \$k^2=$k²\$", titlefontsize=10)
plot!(u, cd.(u,k²), label="cd u", ls=:dash)
plot!(u, sn.(u,k²), label="sn u")
plot!(size=(500, 200), legend=:bottomleft)

showimg("image/jpeg", "images/Jikkyo20140125GaussZeta.jpg", scale="80%")

f(x,y) = exp(-(x^2+y^2))
g(x,y,r) = x^2+y^2 ≤ r^2 ? exp(-(x^2+y^2)) : zero(x)
h(x,y,r) = (abs(x) ≤ r && abs(y) ≤ r) ? exp(-(x^2+y^2)) : zero(x)

x = range(-2, 2, length=201)
y = range(-2, 2, length=201)
r = 1/√2

# z = e^{-(x^2+y^2)} のグラフ
surface(x, y, f.(x', y), size=(400,250), colorbar=false)

# A(r) のグラフ (r=1/√2)
surface(x, y, g.(x', y, r), size=(400,250), colorbar=false)

# B(r) のグラフ (r=1/√2)
surface(x, y, h.(x', y, r), size=(400,250), colorbar=false)

# A(√2 r) のグラフ (r=1/√2)
surface(x, y, g.(x', y, √2*r), size=(400,250), colorbar=false)

showimg("image/jpeg", "images/Jikkyo20140125GammaLaplace.jpg", scale="80%")

stirling_approx(n) = n^n * exp(-n) * √(2π*n)
error(x, x₀) = x - x₀            # 誤差
relative_error(x, x₀) = x/x₀ - 1 # 相対誤差

@printf("%2s %10s %13s %13s %13s\n", "n", "n!", "Stirling", "Error", "Rel.Err.")
println("-"^60)
for n in 1:10
    ft = factorial(n)
    s = stirling_approx(n)
    err = error(s, ft)
    relerr = relative_error(s, ft)
    @printf("%2d %10d %13.4f %13.4f %13.4f\n", n, ft, s, err, relerr)
end

stirling_approx1(n) = n^n * exp(-n) * √(2π*n) * (1 + 1/(12n))

@printf("%2s %10s %20s %13s %13s\n", "n", "n!", "Improved Stirling", "Error", "Rel.Err.")
println("-"^65)
for n in 1:10
    ft = factorial(n)
    s₁ = stirling_approx1(n)
    err = error(s₁, ft)
    relerr = relative_error(s₁, ft)
    @printf("%2d %10d %20.5f %13.5f %13.5f\n", n, ft, s₁, err, relerr)
end

@show stirling_approx(1)
@show stirling_approx1(1);

@show √(2π)
@show √(2π) * 13/12
@show float(ℯ);

showimg("image/jpeg", "images/Jikkyo20140125Stirling.jpg", scale="80%")

# 前者の問題の SymPy による解

n = symbols("n", positive=true)
binom(n,k) = gamma(n+1)/(gamma(k+1)*gamma(n-k+1))
sol = limit((binom(3n,n)/binom(2n,n))^(1/n), n=>oo)
latexstring(raw"\ds\lim_{n\to\infty}\left(\frac{\binom{3n}{n}}{\binom{2n}{n}}\right)^{1/n}=", sympy.latex(sol))

# 後者の問題の SymPy による解

n = symbols("n", positive=true)
sol = limit((1/n)*(gamma(2n+1)/gamma(n+1))^(1/n), n=>oo)
latexstring(raw"\ds\lim_{n\to\infty}\frac{1}{n}\left(\frac{(2n)!}{n!}\right)^{1/n}=", sympy.latex(sol))

showimg("image/jpeg", "images/Jikkyo20140125Beta.jpg", scale="80%")

x = symbols("x", real=true)
sol = 2integrate(√(x*(1-x)^2), (x,0,1))
latexstring(raw"\ds 2\int_0^1\sqrt{x(1-x^2)}\,dx=", sympy.latex(sol))

showimg("image/jpeg", "images/Gamma-Beta-01.jpg", scale="60%")

showimg("image/jpeg", "images/Gamma-Beta-02.jpg", scale="60%")

showimg("image/jpeg", "images/Jikkyo20140125Wallis.jpg", scale="50%")

display("text/html", "<h4>ベータ函数の現れ方(1)</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-01.jpg", scale="50%")

display("text/html", "<h4>ベータ函数の現れ方(2)</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-02.jpg", scale="50%")

display("text/html", "<h4>ガンマ函数の基礎</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-03.jpg", scale="80%")

display("text/html", "<h4>Gaussの超幾何函数の現れ方(1)</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-04.jpg", scale="60%")

display("text/html", "<h4>Gaussの超幾何函数の現れ方(2)</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-05.jpg", scale="70%")

display("text/html", "<h4>Kummerの超幾何函数の現れ方</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-06.jpg", scale="70%")

display("text/html", "<h4>ガンマ函数と正弦函数の関係</h4>")
showimg("image/jpeg", "images/Beta-Gamma-Gauss-Kummer-07.jpg", scale="70%")

display("text/html", "<h4>Hurwitzのゼータ函数とガンマ函数の関係</h4>")
showimg("image/jpeg", "images/Hurwitz-Gamma.jpg", scale="70%")

@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;

@show log(2);
@show sum((-1)^(k-1)/k for k in 1:10^8);

@show π/4;
@show sum((-1)^(k-1)/(2k-1) for k in 1:10^8);