算額（その866）
岩手県平泉町 中尊寺阿弥陀堂（中尊寺地蔵院にて保管） 安政6年(1859)
http://www.wasan.jp/iwate/chusonji3.html
牧下英世：数学史を取り入れた授業実践―算額の教材化と総合的な学習―，2000筑波大学附属駒場論集第40集
https://tsukuba.repo.nii.ac.jp/record/6486/files/10.pdf
正方形内に対角線（角斜）と 1 本の斜線（甲斜）を引き，区分された領域に四分円，乙円，丙円，丁円，戊円，己円を入れる。正方形の一辺の長さが与えられたとき，甲斜の長さと乙円，丙円，丁円，戊円，己円の直径の和を求めよ。
正方形の一辺の長さを a
甲斜と正方形の辺の交点座標を (0, b)
半円の半径と中心座標を r0, (r0, a)
乙円の半径と中心座標を r1, (r0, a - r1)
丙円の半径と中心座標を r2, (a - r2, y2)
丁円の半径と中心座標を r3, (x3, r3)
戊円の半径と中心座標を r4, (r4, y4)
己円の半径と中心座標を r5, (a - r5, y5)
とおき，以下の連立方程式を解く。
SymPy の能力的に，一度に自動的に解くことができないので，手動で逐次解いていく。
include("julia-source.txt");
using SymPy
@syms a::positivev, b::positice, r0::positive,
r1::positive, r2::positive, y2::positive,
r3::positive, x3::positice, r4::positive,
y4::positive, r5::positice, y5::positice, d
s2 = √Sym(2)
r1 = r0/2
eq1 = numerator(apart(dist(0, 0, a, a, r0, a) - r0^2, d))
eq2 = numerator(apart(dist(0, 0, a, a, a - r2, y2) - r2^2, d))
eq3 = numerator(apart(dist(0, 0, a, a, x3, r3) - r3^2, d))
eq4 = numerator(apart(dist(0, 0, a, a, r4, y4) - r4^2, d))
eq5 = numerator(apart(dist(0, b, a, 0, r0, a) - r0^2, d))
eq6 = numerator(apart(dist(0, b, a, 0, a - r2, y2) - r2^2, d))
eq7 = numerator(apart(dist(0, b, a, 0, x3, r3) - r3^2, d))
eq8 = numerator(apart(dist(0, b, a, 0, r4, y4) - r4^2, d))
eq9 = numerator(apart(dist(0, b, a, 0, a - r5, y5) - r5^2, d))
eq10 = (r2 - r5)^2 + (y2 - y5)^2 - (r2 + r5)^2 |> expand;
1. eq1, eq5 から r0, b を求める。
res_r0 = solve(eq1, r0)[1] # r0 確定
res_r0 |> println
res_r0(a => 1.5) |> N # 0.6213203435596426
a*(-1 + sqrt(2))
0.6213203435596426
res_b = solve(eq5(r0 => res_r0), b)[1] # b 確定
res_b |> println
res_b(a => 1.5) |> N # 1.2016314489305129
-3*a*sqrt(20 - 14*sqrt(2)) - 2*sqrt(2)*a*sqrt(20 - 14*sqrt(2)) + sqrt(2)*a + 2*a
1.201631448930513
2. eq2, eq6 から y2, r2 を求める。
res_r2 = solve(eq2, r2)[2]
res_r2 |> println
-a + y2 + sqrt(2)*(a - y2)
eq16 = eq6(b => res_b, r2 => res_r2) |> expand |> simplify
eq16 |> println
a^2*(-3*a^2 + 2*sqrt(2)*a^2 - 6*sqrt(2)*a*y2 + 2*a*y2*sqrt(20 - 14*sqrt(2)) + 4*a*y2*sqrt(10 - 7*sqrt(2)) + 6*a*y2 - 2*y2^2 - 4*y2^2*sqrt(10 - 7*sqrt(2)) - 2*y2^2*sqrt(20 - 14*sqrt(2)) + 4*sqrt(2)*y2^2)
res_y2 = solve(eq16, y2)[2]
res_y2 |> println
a*(-3*sqrt(2) - sqrt(-4*sqrt(2) - 4*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + 9) + sqrt(20 - 14*sqrt(2)) + 2*sqrt(10 - 7*sqrt(2)) + 3)/(2*(-2*sqrt(2) + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(10 - 7*sqrt(2)) + 1))
res_y2 = apart(res_y2/a, d)*a |> factor # y2 確定
res_y2 |> println
res_y2(a => 1.5) |> N # 0.694654694074625
a*(-9*sqrt(2) + 9*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 14*sqrt(10 - 7*sqrt(2)) + 20)/34
0.694654694074625
res_r2 = res_r2(y2 => res_y2) |> simplify # r2 確定
res_r2 |> println
res_r2(a => 1.5) |> N # 0.33358494810779965
a*(-5*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + 4 + 5*sqrt(2))/34
0.33358494810779965
3. eq3, eq7 から x3，r3 を求める。
res_r3 = solve(eq3, r3)[1]
res_r3 |> println
x3*(-1 + sqrt(2))
eq17 = eq7(b => res_b, r3 => res_r3) |> expand |> simplify
eq17 |> println
a^2*(-20*a^2*sqrt(20 - 14*sqrt(2)) - 28*a^2*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2)*a^2 + 10*a^2 - 20*a*x3 - 14*sqrt(2)*a*x3 + 42*a*x3*sqrt(20 - 14*sqrt(2)) + 60*a*x3*sqrt(10 - 7*sqrt(2)) - 28*x3^2*sqrt(10 - 7*sqrt(2)) - 18*x3^2*sqrt(20 - 14*sqrt(2)) + 4*x3^2 + 10*sqrt(2)*x3^2)
res_x3 = apart(solve(eq17, x3)[1]/a, d)*a |> simplify # x3 確定
res_x3 |> println
res_x3(a => 1.5) |> N # 0.6882058497807045
a*(-2*sqrt(10 - 7*sqrt(2)) - sqrt(20 - 14*sqrt(2)) + 2)/2
0.6882058497807045
res_r3 = res_r3(x3 => res_x3) |> simplify # r3 確定
res_r3 |> println
res_r3(a => 1.5) |> N # 0.2850641966836687
a*(-2 - sqrt(20 - 14*sqrt(2)) + 2*sqrt(2))/2
0.2850641966836687
4. eq4, eq8 から y4，r4 を求める。
res_r4 = solve(eq4, r4)[1]
res_r4 |> println
y4*(-1 + sqrt(2))
eq18 = eq8(b => res_b, r4 => res_r4) |> expand |> simplify
eq18 |> println
a^2*(-20*a^2*sqrt(20 - 14*sqrt(2)) - 28*a^2*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2)*a^2 + 10*a^2 - 10*sqrt(2)*a*y4 - 8*a*y4 + 22*a*y4*sqrt(20 - 14*sqrt(2)) + 32*a*y4*sqrt(10 - 7*sqrt(2)) - 2*y4^2 - 4*y4^2*sqrt(10 - 7*sqrt(2)) - 2*y4^2*sqrt(20 - 14*sqrt(2)) + 4*sqrt(2)*y4^2)
res_y4 = solve(eq18, y4)[2]
res_y4 |> println
a*(-5 - 3*sqrt(2) + 14*sqrt(10 - 7*sqrt(2)) + 10*sqrt(20 - 14*sqrt(2)))/(-2*sqrt(2) + sqrt(2)*sqrt(10 - 7*sqrt(2)) + 2*sqrt(10 - 7*sqrt(2)) + 1)
res_y4 = apart(solve(eq18, y4)[2]/a, d)*a |> simplify # y4 確定
res_y4 |> println
res_y4(a => 1.5) |> N # 0.6451521645656638
a*(-78*sqrt(10 - 7*sqrt(2)) - 55*sqrt(20 - 14*sqrt(2)) + 27 + 21*sqrt(2))/17
0.6451521645656638
res_r4 = res_r4(y4 => res_y4) |> simplify # r4 確定
res_r4 |> println
res_r4(a => 1.5) |> N # 0.26723077635745685
a*(-23*sqrt(20 - 14*sqrt(2)) - 32*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 15)/17
0.2672307763574568
5. eq9, eq10 から r5，y5 を求める。
eq19 = eq9(b => res_b) |> simplify
eq19 |> println
a^2*(-r5^2 - 2*r5*y5*(-3*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + sqrt(2) + 2) + y5^2)
res_y5 = solve(eq19, y5)[1]
res_y5 |> println
r5*(-3*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + sqrt(2) + 2)
eqA = eq10(r2 => res_r2, y2 => res_y2, y5 => res_y5) |> simplify
eqA |> println
a^2*(-9*sqrt(2) + 9*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 14*sqrt(10 - 7*sqrt(2)) + 20)^2/1156 - a*r5*(-9*sqrt(2) + 9*sqrt(20 - 14*sqrt(2)) + 14*sqrt(10 - 7*sqrt(2)) + 20)*(-3*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + sqrt(2) + 2)/17 - 2*a*r5*(-5*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + 4 + 5*sqrt(2))/17 + r5^2*(-3*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + sqrt(2) + 2)^2
res_r5 = solve(eqA, r5)[2] # r5 決定
res_r5 |> println
res_r5(a => 1.5) |> N # 0.13202643690561106
0.13202643690561106
res_y5 = res_y5(r5 => res_r5) |> simplify # y5 決定
res_y5 |> println
res_y5(a => 1.5) |> N # 0.27493082244464034
0.27493082244464034
6. 甲斜の長さを求める。
甲斜 = res_b^2 + a^2 |> expand |> simplify |> sqrt
甲斜 |> println
甲斜(a => 1.5) |> N # 1.281304567672052
sqrt(a^2*(-20*sqrt(20 - 14*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11))
1.921956851508078
7. 6和を求める
(甲斜 + res_r0 + 2(res_r2 + res_r3 + res_r4 + res_r5))(a => 1.5) |> N # 4.579089911176793
4.579089911176793
和 = 甲斜 + res_r0 + 2(res_r2 + res_r3 + res_r4 + res_r5)
和 |> println
a*(-40*sqrt(20 - 14*sqrt(2)) - 56*sqrt(10 - 7*sqrt(2)) - 6*sqrt(-34*sqrt(2) - 8*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 52) - 8*sqrt(-17*sqrt(2) - 4*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 26) + 2*sqrt(-40*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 56*sqrt(10 - 7*sqrt(2)) + 12*sqrt(2) + 22) + 4*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 12*sqrt(2) + 21)*(-2*sqrt(-184*sqrt(2) + 14*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 52*sqrt(10 - 7*sqrt(2)) + 14*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 42*sqrt(2)*sqrt(10 - 7*sqrt(2))*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 26*sqrt(2)*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 156*sqrt(10 - 7*sqrt(2))*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 390) - 9*sqrt(-40*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 56*sqrt(10 - 7*sqrt(2)) + 12*sqrt(2) + 22) - 2*sqrt(20 - 14*sqrt(2)) + 2*sqrt(2) + 12*sqrt(10 - 7*sqrt(2)) + 9*sqrt(-34*sqrt(2) - 8*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 52) + 14*sqrt(-17*sqrt(2) - 4*sqrt(2)*sqrt(10 - 7*sqrt(2)) + 26) + 22 + 20*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11))/(17*(-40*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 56*sqrt(10 - 7*sqrt(2)) - 6*sqrt(2)*sqrt(10 - 7*sqrt(2))*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) - 8*sqrt(10 - 7*sqrt(2))*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 2*sqrt(2)*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 4*sqrt(-20*sqrt(2)*sqrt(10 - 7*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11) + 12*sqrt(2) + 21)^2) + 2*a*(-23*sqrt(20 - 14*sqrt(2)) - 32*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 15)/17 + a*(-2 - sqrt(20 - 14*sqrt(2)) + 2*sqrt(2)) + a*(-1 + sqrt(2)) + a*(-5*sqrt(20 - 14*sqrt(2)) - 4*sqrt(10 - 7*sqrt(2)) + 4 + 5*sqrt(2))/17 + sqrt(a^2*(-20*sqrt(20 - 14*sqrt(2)) - 28*sqrt(10 - 7*sqrt(2)) + 6*sqrt(2) + 11))
6和を求める式は長いので，精一杯簡約化してもなお長い。
a = 1.5
v = 5√2 + 7
s = sqrt(10 - 7√2)
t = sqrt(11 + 6√2 - 4s*v)
u = sqrt(26 - √2(17 + 4s))
a*(2(2 + √2)t + 21 + 12√2 - 8v*s - 2(4 + 3√2)u)*
(-2√2sqrt((26 + 7√2)s + (7 + 13√2)t + (78 + 21√2)s*t + 195 - 92√2) + 2(6 - √2)s + (20 - 9√2)t + (14 + 9√2)u + 2(11 + √2))/
(17(2(2 + √2)t - 8v*s - 2(4 + 3√2)s*t + 21 + 12√2)^2) +
a*(√2(4 - s) - 1 - (68 + 51√2)s/17 + sqrt(11 + 6√2 - 4v*s))
4.579089911176797
正方形の一辺の長さが 1.5 のとき，6和は 4.579089911176797 である。
術では，「置二個開平方以減二個開平法加三個乗方面得六和」とあるが，そんなに短いわけがなく，適切な術でもない。
8. 数値解を求めるのは簡単である。
#=
using NLsolve
function nls(func, params...; ini = [0.0])
if typeof(ini) <: Number
r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
v = r.zero[1]
else
r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
v = r.zero
end
return Float64.(v), r.f_converged
end;
function H(u)
(b, r0, r2, y2, r3, x3, r4, y4, r5, y5) = u
return [
a^2/2 - a*r0 - r0^2/2, # eq1
a^2/2 - a*r2 - a*y2 - r2^2/2 + r2*y2 + y2^2/2, # eq2
-r3^2/2 - r3*x3 + x3^2/2, # eq3
-r4^2/2 - r4*y4 + y4^2/2, # eq4
a^4 - 2*a^3*b + a^2*b^2 + 2*a^2*b*r0 - a^2*r0^2 - 2*a*b^2*r0, # eq5
-a^2*r2^2 + a^2*y2^2 - 2*a*b*r2*y2, # eq6
a^2*b^2 - 2*a^2*b*r3 - 2*a*b^2*x3 + 2*a*b*r3*x3 - b^2*r3^2 + b^2*x3^2, # eq7
a^2*b^2 - 2*a^2*b*y4 - a^2*r4^2 + a^2*y4^2 - 2*a*b^2*r4 + 2*a*b*r4*y4, # eq8
-a^2*r5^2 + a^2*y5^2 - 2*a*b*r5*y5, # eq9
-4*r2*r5 + y2^2 - 2*y2*y5 + y5^2, # eq10
]
end;
a = 1.5
iniv = a .* BigFloat[0.801, 0.414, 0.222, 0.463, 0.19, 0.459, 0.178, 0.43, 0.088, 0.183]
res = nls(H, ini=iniv)
([1.2016314489305129, 0.6213203435596426, 0.33358494810779965, 0.694654694074625, 0.2850641966836687, 0.6882058497807045, 0.26723077635745685, 0.6451521645656638, 0.13202643690561106, 0.27493082244464034], true)
=#
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(b, r0, r2, y2, r3, x3, r4, y4, r5, y5) = [1.2016314489305129, 0.6213203435596426, 0.33358494810779965, 0.694654694074625, 0.2850641966836687, 0.6882058497807045, 0.26723077635745685, 0.6451521645656638, 0.13202643690561106, 0.27493082244464034]
a = 1.5
r1 = r0/2
和 = sqrt(a^2 + b^2) + 2(r1 + r2 + r3 + r4 + r5)
println("和 = $和")
#(b, r0, r2, y2, r3, x3, r4, y4, r5) = [85.716, 44.321, 23.796, 49.552, 20.335, 49.092, 19.062, 46.021, 10]
@printf("b = %g; r0 = %g; r2 = %g; y2 = %g; r3 = %g; x3 = %g; r4 = %g; y4 = %g; r5 = %g; y5 = %g\n", b, r0, r2, y2, r3, x3, r4, y4, r5, y5)
plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:blue, lw=0.5)
circle(r0, a, r0, beginangle=180, endangle=360)
circle(r0, a - r1, r1, :green)
circle(a - r2, y2, r2, :magenta)
circle(x3, r3, r3, :blue)
circle(r4, y4, r4, :brown)
circle(a - r5, y5, r5, :orange)
segment(0, 0, a, a, :tomato)
segment(0, b, a, 0, :tomato)
if more
delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3 # size[2] * fontsize * 2
hline!([0], color=:gray80, lw=0.5)
vline!([0], color=:gray80, lw=0.5)
point(a, 0, " a", :blue, :left, :bottom, delta=delta/2)
point(0, a, "a ", :blue, :right, :vcenter)
point(0, b, "b ", :tomato, :right, :vcenter)
point(r0, a, "半円:r0,(r0,a)", :red, :center, :bottom, delta=delta/2)
point(r0, a - r1, "乙円:r1,(r0,a-r1)", :green, :center, delta=-delta/2)
point(a - r2, y2, "丙円:r2,(a-r2,y2)", :magenta, :center, delta=-delta/2)
point(x3, r3, "丁円:r3,(x3,r3)", :blue, :center, delta=-delta/2)
point(r4, y4, "戊円:r4,(r4,y4)", :brown, :center, delta=-delta/2)
point(a - r5, y5, "己円:r5\n(a-r5,y5)", :black, :center, delta=-delta/2)
plot!(xlims=a.*(-0.05, 1.05))
end
end;
