算額（その1041）
八十三 藤沢町保呂羽保呂羽山 保呂羽神社 明治26年(1893)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
直線の上に甲円，丙円，丁円が隣同士互いに接して載っている。さらに，丙円と丁円の上に乙円が載っている。乙円のてっぺんまでの高さはいかほどか。
乙円の半径がわかれば，算額1040 と同じ問題になる。
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (0, r3)
丁円の半径と中心座標を r4, (x4, r4)
とおき，以下の連立方程式を解く。
高さは y2 + r2 である。
r1, r3, r4 が既知なので，x1 と x4 は簡単に計算できる。ちなみに，「算法助術」の公式40ではある。
include("julia-source.txt");
using SymPy
@syms r1::positive, r3::positive, x1::positive,
r3::positive, r4::positive, x4::positive
r2::positive, x2::positive, y2::positive,
eq2 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^2 |> expand
eq5 = x4^2 + (r3 - r4)^2 - (r3 + r4)^2 |> expand;
res1 = solve([eq2, eq5], (x1, x4))[2] # 2 of 2
(2*sqrt(r1)*sqrt(r4) + 2*sqrt(r3)*sqrt(r4), 2*sqrt(r3)*sqrt(r4))
x1, x4 も既知となったので， x2, y2 も求めることができる。x4, x1 を展開するかしないかで，簡約化の程度が変わるが最終的な結果は同じである。
include("julia-source.txt");
using SymPy
@syms r1::positive, r3::positive, x1::positive,
r3::positive, r4::positive, x4::positive,
r2::positive, x2::positive, y2::positive
#x4 = 2sqrt(r3*r4)
#x1 = x4 + 2sqrt(r1*r4)
eq3 = x2^2 + (y2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = (x4 - x2)^2 + (y2 - r4)^2 - (r2 + r4)^2 |> expand
res2 = solve([eq3, eq4], (x2, y2))[2] # 2 of 2
((2*r2*r3 - 2*r2*r4 + x4^2 + (2*r3 - 2*r4)*(x4^2*sqrt(4*r2^2 + 4*r2*r3 + 4*r2*r4 + 4*r3*r4 - x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2)) + (-2*r2*r3^2 + 4*r2*r3*r4 - 2*r2*r4^2 + r3*x4^2 + r4*x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2))))/(2*x4), x4^2*sqrt(4*r2^2 + 4*r2*r3 + 4*r2*r4 + 4*r3*r4 - x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2)) + (-2*r2*r3^2 + 4*r2*r3*r4 - 2*r2*r4^2 + r3*x4^2 + r4*x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2)))
x2, y2 を eq1 に代入し f(r2) = 0 となる方程式を抽出する（分数式の分母のみ，因数分解で0でないことが明らかな項を除く）。
eq1 = (x1 - x2)^2 + (y2 - r1)^2 - (r1 + r2)^2
eq01 = eq1(x2 => res2[1], y2 => res2[2]) |> simplify;
eq02 = numerator(eq01) |> expand |> factor
res3 = solve(eq02.args[4], r2)[1]
res3(x1 => x4 + 2sqrt(r1*r4))(x4 => 2sqrt(r3*r4)) |> simplify |> println
r4^2*(r1^(3/2)*sqrt(r3) + sqrt(r1)*r3^(3/2) + r1^2/4 + 3*r1*r3/2 + r3^2/4)/(2*r1^(3/2)*r3^(3/2) - r1^(3/2)*sqrt(r3)*r4 - sqrt(r1)*r3^(3/2)*r4 + r1^2*r3 + r1*r3^2 - 2*r1*r3*r4)
以上をまとめる。
注：SymPy では式の定義は，ほかの式にも引き継がれるので，最終的な式は定義したものよりも複雑になる。プログラム言語では式は即時に評価されるので，式の長さ（複雑さ）が変わることはない。
using SymPy
@syms r1::positive, r3::positive, x1::positive,
r3::positive, r4::positive, x4::positive,
r2::positive, x2::positive, y2::positive
x4 = 2sqrt(r3*r4)
x1 = x4 + 2sqrt(r1*r4)
r2 = r4^2*(r1^(3/2)*sqrt(r3) + sqrt(r1)*r3^(3/2) + r1^2/4 + 3r1*r3/2 + r3^2/4)/(2r1^(3/2)*r3^(3/2) - r1^(3/2)*sqrt(r3)*r4 - sqrt(r1)*r3^(3/2)*r4 + r1^2*r3 + r1*r3^2 - 2r1*r3*r4)
x2 = (2r2*r3 - 2r2*r4 + x4^2 + (2r3 - 2r4)*(x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2))))/2x4
y2 = x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2*(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2));
r2 はかなり長い式になり，SymPy ではうまく簡約化できないが，手作業でやると以下のようになる。
num = r4^2*(r1^(3/2)*sqrt(r3) + sqrt(r1)*r3^(3/2) + r1^2/4 + 3r1*r3/2 + r3^2/4)
den = (2r1^(3/2)*r3^(3/2) - r1^(3/2)*sqrt(r3)*r4 - sqrt(r1)*r3^(3/2)*r4 + r1^2*r3 + r1*r3^2 - 2r1*r3*r4);
@syms t1, t3
num.args[2](sqrt(r1) => t1)(sqrt(r3) => t3) |> simplify |> println
t1^4/4 + 3*t1^2*t3^2/2 + t1*t3^3.0 + t1^3.0*t3 + t3^4/4
t1^4/4 + 3t1^2*t3^2/2 + t1*t3^3 + t1^3*t3 +t3^4/4 |> factor |> println
(t1 + t3)^4/4
num2 = r4^2*(√r1 + √r3)^4/4
num2 |> println
r4^2*(sqrt(r1) + sqrt(r3))^4/4
num(r1 => 5, r3 => 3, r4 => 2).evalf(), num2(r1 => 5, r3 => 3, r4 => 2).evalf()
(247.935467078637, 247.935467078637)
den(sqrt(r1) => t1)(sqrt(r3) => t3) |> println
-2*r4*t1^2*t3^2 - r4*t1*t3^3.0 - r4*t1^3.0*t3 + t1^4*t3^2 + t1^2*t3^4 + 2*t1^3.0*t3^3.0
den2 = -2*r4*t1^2*t3^2 - r4*t1*t3^3 - r4*t1^3*t3 + t1^4*t3^2 + t1^2*t3^4 + 2*t1^3*t3^3;
den2 |> factor |> println
-t1*t3*(r4 - t1*t3)*(t1 + t3)^2
den2 = -√r1*√r3*(r4 - √r1*√r3)*(√r1 + √r3)^2
den2 |> println
-sqrt(r1)*sqrt(r3)*(sqrt(r1) + sqrt(r3))^2*(-sqrt(r1)*sqrt(r3) + r4)
den(r1 => 5, r3 => 3, r4 => 2).evalf(), den2(r1 => 5, r3 => 3, r4 => 2).evalf()
(114.221766846904, 114.221766846904)
r2 は最終的に以下のように比較的短い式で表すことができる。
r2 = num2/den2
r2 = r2 |> simplify
r2 |> println
r4^2*(2*sqrt(r1)*sqrt(r3) + r1 + r3)/(4*(-sqrt(r1)*sqrt(r3)*r4 + r1*r3))
甲円，丙円，丁円の直径が 8, 5, 4 のとき，乙円のてっぺんまでの高さは 10.883036880224509 である。
(r1, r3, r4) = (8, 5, 4) ./ 2
x4 = 2sqrt(r3*r4)
x1 = x4 + 2sqrt(r1*r4)
r2 = r4^2*(2*sqrt(r1)*sqrt(r3) + r1 + r3)/(4*(-sqrt(r1)*sqrt(r3)*r4 + r1*r3))
x2 = (2r2*r3 - 2r2*r4 + x4^2 + (2r3 - 2r4)*(x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2))))/2x4
y2 = x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2*(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2))
(x4, x1, r2, x2, y2) |> println
height = y2 + r2
height |> println
(4.47213595499958, 10.12899020449196, 3.4892527130926094, 3.452828490790752, 7.393784167131899)
10.883036880224509
術は，以下のように非常に短い式を提示している。
using SymPy
@syms d1, d3, d4 # r1, r3, r4 の 2 倍（直径）
(d1, d3, d4) = (8, 5, 4)
A = sqrt(d1*d3)
B = d4/A
C = 1 - B
高 = d4/C
高 |> println
10.883036880224505
function draw(more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r3, r4) = (8, 5, 4) ./ 2
x4 = 2sqrt(r3*r4)
x1 = x4 + 2sqrt(r1*r4)
r2 = r4^2*(2*sqrt(r1)*sqrt(r3) + r1 + r3)/(4*(-sqrt(r1)*sqrt(r3)*r4 + r1*r3))
x2 = (2r2*r3 - 2r2*r4 + x4^2 + (2r3 - 2r4)*(x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2(r3^2 - 2r3*r4 + r4^2 + x4^2))))/2x4
y2 = x4^2*sqrt(4r2^2 + 4r2*r3 + 4r2*r4 + 4r3*r4 - x4^2)/(2*(r3^2 - 2r3*r4 + r4^2 + x4^2)) + (-2r2*r3^2 + 4r2*r3*r4 - 2r2*r4^2 + r3*x4^2 + r4*x4^2)/(2*(r3^2 - 2*r3*r4 + r4^2 + x4^2))
@printf("甲円，丙円，丁円の直径がそれぞれ %g, %g, %g のとき，乙円のてっぺんまでの高さは %g である。\n", 2r1, 2r3, 2r4, y2 + r2)
@printf("r1 = %g; x1 = %g; r2 = %g; x2 = %g; y2 = %g; r3 = %g; r4 = %g; x4 = %g\n", r1, x1, r2, x2, y2, r3, r4, x4)
plot()
circle(x1, r1, r1)
circle(x2, y2, r2, :green)
circle(0, r3, r3, :blue)
circle(x4, r4, r4, :magenta)
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(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
point(x2, y2, "乙円:r2,(x2,y2)", :green, :center, delta=-delta)
point(0, r3, "丙円:r3,(0,r3)", :blue, :center, delta=-delta)
point(x4, r4, "丁円:r4,(x4,r4)", :magenta, :center, delta=-delta)
point(x2, y2 + r2, "高さ:y2+r2", :green, :center, :bottom, delta=delta)
end
end;
