算額（その866）

算額（その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;