https://sangaku0418.hatenablog.com/entry/2024/12/30/222817

90 岩手県陸前高田市小友町上の坊常膳寺観音堂天保13年(1842)

安富有恒：和算—岩手の現存算額のすべて，青磁社，東京都，1987.
http://www.wasan.jp/iwatenosangaku_yasutomi.pdf
キーワード：円2個，直角三角形，正方形，斜線
#Julia #SymPy #算額 #和算 #数学

直角三角形の中に甲円，乙円，正方形，斜線を容れる。鈎，股が与えられたとき乙円の直径を求める術を述べよ。

鈎，股をそのまま，「鈎」，「股」
正方形の一辺の長さを \(a\)
斜線と正方形の一辺の交点座標を \( (b,\ 0)\)
甲円の半径と中心座標を \(r_1,\ (a + r_1,\ r_1)\)
乙円の半径と中心座標を \(r_2,\ (a - r_2,\ a - r_2)\)
とおき，以下の方程式を解く。

一度に解けないので，逐次求めていく。

include("julia-source.txt"); # julia-source.txt ソース

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, r1::positive, r2::positive, a::positive, b::positive
eq1 = (股 - a) + a - sqrt( (股 - a)^2 + a^2) - 2r1
eq2 = dist2(0, 鈎, b, 0, a + r1, r1, r1)
eq3 = dist2(0, 鈎, b, 0, a - r2, a - r2, r2)
eq4 = a/(股 - a) - 鈎/股;
# res = solve([eq1, eq4, eq5, eq6], (r2, r1, a, b))

# a
ans_a = solve(eq4, a)[1]
ans_a |> println

股*鈎/(股 + 鈎)

# r1
eq11 = eq1(a => ans_a)
ans_r1 = solve(eq11, r1)[1]#(sqrt(鈎^2 + 股^2) => 弦)
ans_r1 |> println

股*(股 + 鈎 - sqrt(股^2 + 鈎^2))/(2*(股 + 鈎))

# b
eq12 = eq2(a => ans_a, r1 => ans_r1) |> simplify
ans_b = solve(eq12, b)[2](sqrt(鈎^2 + 股^2) => 弦) |> simplify
ans_b |> println

股*鈎^2*(-弦 + 股 + 2*鈎)/(弦*股^2 + 弦*股*鈎 - 股^3 - 股^2*鈎 + 股*鈎^2 + 鈎^3)

# r2
eq13 = eq3(a => ans_a, r1 => ans_r1, b => ans_b)
ans_r2 = solve(eq13, r2)[2]; # 1 of 2
ans_r2 |> println

( (股 + 鈎)^3*(-弦 + 股 + 2*鈎)*(弦*股 - 股^2 + 鈎^2)*sqrt(弦^2*股^4 + 2*弦^2*股^3*鈎 + 2*弦^2*股^2*鈎^2 - 2*弦*股^5 - 4*弦*股^4*鈎 - 2*弦*股^3*鈎^2 + 2*弦*股*鈎^4 + 股^6 + 2*股^5*鈎 + 3*股^2*鈎^4 + 2*股*鈎^5 + 鈎^6) + (弦*股^2 - 股^3 + 3*股*鈎^2 + 鈎^3)*(弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6))*(弦*股^2 + 弦*股*鈎 + 弦*鈎^2 - 股^3 - 股^2*鈎 - 鈎^3)/(2*(股 + 鈎)^2*(-弦 + 股 + 2*鈎)*(弦*股 - 股^2 + 鈎^2)*(弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6))

計算式はかなり長く，SymPy でも簡約化しきれない。
弦を前もって計算しておき，以下により小円の半径を求める。

鈎 = 3, 股 = 4, 弦 = 5 のとき，小円の半径は 0.303296703296703 である。

ans_r2(弦 => sqrt(鈎^2 + 股^2))(鈎 => 3, 股 => 4) |> println
ans_r2(弦 => sqrt(鈎^2 + 股^2))(鈎 => 3, 股 => 4).evalf() |> println

138/455
0.303296703296703

描画関数プログラムのソースを見る

function draw(鈎, 股, more=false)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
弦 = sqrt(鈎^2 + 股^2)
term1 = 弦^2*股^4 + 3*弦^2*股^3*鈎 + 3*弦^2*股^2*鈎^2 + 弦^2*股*鈎^3 - 2*弦*股^5 - 8*弦*股^4*鈎 - 11*弦*股^3*鈎^2 - 5*弦*股^2*鈎^3 + 弦*股*鈎^4 + 弦*鈎^5 + 股^6 + 5*股^5*鈎 + 8*股^4*鈎^2 + 2*股^3*鈎^3 - 7*股^2*鈎^4 - 7*股*鈎^5 - 2*鈎^6
term2 = 弦*股 - 股^2 + 鈎^2
term3 = -弦 + 股 + 2*鈎
term4 = 股 + 鈎
num = (term4^3*term3*term2*sqrt(弦^2*股^4 + 2*弦^2*股^3*鈎 + 2*弦^2*股^2*鈎^2 - 2*弦*股^5 - 4*弦*股^4*鈎 - 2*弦*股^3*鈎^2 + 2*弦*股*鈎^4 + 股^6 + 2*股^5*鈎 + 3*股^2*鈎^4 + 2*股*鈎^5 + 鈎^6) + (弦*股^2 - 股^3 + 3*股*鈎^2 + 鈎^3)*term1)*(弦*股^2 + 弦*股*鈎 + 弦*鈎^2 - 股^3 - 股^2*鈎 - 鈎^3)
den = 2*term4^2*term3*term2*term1
r1 = 股*(-弦 + 股 + 鈎)/(2*(股 + 鈎))
b = 股*鈎^2*(-弦 + 股 + 2*鈎)/(弦*股^2 + 弦*股*鈎 - 股^3 - 股^2*鈎 + 股*鈎^2 + 鈎^3)
a = 股*鈎/(股 + 鈎)
r2 = num/den
println(r2)
plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:magenta, lw=0.5)
plot!([a, a, 0], [0, a, a], color=:magenta, lw=0.5)
segment(0, 鈎, b, 0)
circle(a + r1, r1, r1, :magenta)
circle(a - r2, a - r2, r2, :blue)
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 + r1, r1, "甲円:r1,(a+r1,r1)", :red, :center, delta=-delta)
point(a - r2, a - r2, "乙円:r2,(a-r2,a-r2)", :red, :center, delta=-delta)
point(a, 0, "a", :black, :center, delta=-delta)
point(b, 0, "b", :black, :center, delta=-delta)
point(0, 鈎, "鈎", :magenta, :left, :bottom, delta=delta)
point(股, 0, "股", :magenta, :left, :bottom, delta=delta)
ylims!(-5delta, 鈎 +5delta)
end
end;