算額（その814）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
外円内に 5 個の等円を入れる。外円と等円に挟まれる部分の面積（黒積）が 2330.89 歩のとき，等円の直径を求めよ。
黒積は，「扇形 OCD - 直角三角形 OAB - 扇形ACB」の面積の 10 倍である。
等円の半径 OB = R，AB = r とすれば，∠AOB = 36° より，R = r/sind(36)
OC = OB + BC = R + r = Rr とすれば，
扇形 OCD = πRr^2/10
直角三角形 OAB = r*R*cosd(36)/2
扇形 ACB = πr^2*(126/360)
よって，黒積 = 10*(π*Rr^2/10 - r*R*cosd(36)/2 - π*r^2*126/360)
黒積 = 2332.89 を解いて r を求めると，r = 21.5000025771553，等円の直径は 43.0000051543107 寸である。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r::positive, R::positive, Rr::positive,
黒積::positive, BLK::positive
R = solve(R*sind(Sym(36)) - r, R)[1]
Rr = R + r
R |> println
Rr |> println
(R.evalf(), Rr.evalf()) |> println
2*sqrt(2)*r/sqrt(5 - sqrt(5))
r + 2*sqrt(2)*r/sqrt(5 - sqrt(5))
(1.70130161670408*r, 2.70130161670408*r)
黒積 = PI*Rr^2/10 - R*r*cosd(Sym(36))/2 - PI*r^2*126//360
eq = 10黒積 - BLK # 2332.89;
eq |> println
-BLK - 7*pi*r^2/2 - 10*sqrt(2)*r^2*(1/4 + sqrt(5)/4)/sqrt(5 - sqrt(5)) + pi*(r + 2*sqrt(2)*r/sqrt(5 - sqrt(5)))^2
方程式を解いて，等円の直径を求める。
res = 2solve(eq, r)[1]
res |> println
2*sqrt(2)*sqrt(BLK)*(5 - sqrt(5))^(3/4)/sqrt(-8*sqrt(10)*pi - 20*sqrt(10) - 9*pi*sqrt(5 - sqrt(5)) + 5*sqrt(5)*pi*sqrt(5 - sqrt(5)) + 40*sqrt(2)*pi)
黒積が 2332.89 のときの等円の直径を求める。
res(BLK => 2332.89) |> println
96.6*sqrt(2)*(5 - sqrt(5))^(3/4)/sqrt(-8*sqrt(10)*pi - 20*sqrt(10) - 9*pi*sqrt(5 - sqrt(5)) + 5*sqrt(5)*pi*sqrt(5 - sqrt(5)) + 40*sqrt(2)*pi)
等円の直径は 43 寸である。
res(BLK => 2332.89).evalf() |> println
43.0000051543107
function transform(x, y; deg=0, dx=0, dy=0)
return [cosd(deg) -sind(deg); sind(deg) cosd(deg)] * [x, y]
end;
function rotatefigure(x2, y2)
for i in 1:5
plot!(x2, y2, color=:gray80, seriestype=:shape, fillcolor=:gray80)
i == 5 && return
(x2, y2) = transform(x2, y2, deg = 72)
end
end;
function fillblack(r, Rr, x, y)
x2 = []
y2 = []
for i in 180:0.5:306
append!(x2, r*cosd(i) + x[1])
append!(y2, r*sind(i) + y[1])
end
for i in 306:-0.5:270
append!(x2, Rr*cosd(i))
append!(y2, Rr*sind(i))
end
x2 = vcat(x2, -reverse(x2))
y2 = vcat(y2, reverse(y2))
rotatefigure(x2, y2)
end;
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
x = Vector{Float64}(undef, 6)
y = Vector{Float64}(undef, 6)
r = 21.5000025771553
R = 2√2*r/sqrt(5 - √5)
Rr = R + r
blk = 10(π*Rr^2/10 - R*r*cosd(36)/2 - π*r^2*126/360)
@printf("等円の直径 = %g; 黒積 = %g; R = %g; Rr = %g\n", 2r, blk, R, Rr)
for i = 1:6
θ = 306 + (i - 1)*72
x[i] = R*cosd(θ)
y[i] = R*sind(θ)
end
plot()
fillblack(r, Rr, x, y)
for i = 1:5
segment(x[i], y[i], x[i + 1], y[i + 1], :blue)
circle(x[i], y[i], r, :green)
end
circle(0, 0, R)
circle(0, 0, Rr, :magenta)
segment(0, 0, 0, -Rr, :black, lw=1)
segment(0, 0, Rr*cosd(306), Rr*sind(306), :black, lw=1)
if more == true
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(x[1], y[1])
point(R, 0, " R", :red, :left, :bottom, delta=delta)
point(Rr, 0, " Rr", :magenta, :left, :bottom, delta=delta)
point(0, -R*cosd(36), "A ", :black, :right, delta=-delta/2)
point(x[1], y[1], " B", :black, :left, delta=-delta/2)
point(Rr*cosd(306), Rr*sind(306), "C", :black, :center, delta=-delta/2)
point(0, -Rr, "D ", :black, :right, :bottom, delta=delta/2)
point(0, 0, "O", :black, :center, :bottom, delta=delta/2)
end
end;
算額（その813）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
不等辺三角形内に全円，大円，中円，小円を入れる。大円，中円，小円の直径がそれぞれ 256 寸，225 寸，144 寸のとき，全円の直径はいかほどか。
底辺の長さ a，頂点の座標を (b, h)
全円の半径と中心座標を r0, (x0, r0)
大円の半径と中心座標を r1, (x1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, h::positive,
r0::positive, x0::positive,
r1::positive, x1::positive,
r2::positive, x2::positive,
r3::positive, x3::positive, y3::positive,
d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x3 - x2)^2 + (y3 - r2)^2 - (r2 + r3)^2
eq4 = (a + sqrt(b^2 + h^2) + sqrt((a - b)^2 + h^2))r0 - a*h
eq5 = r0/(a - x0) - r1/(a - x1)
eq6 = r0/x0 - r2/x2
eq7 = dist(0, 0, b, h, x2, r2) - r2^2
eq8 = dist(0, 0, b, h, x3, y3) - r3^2
eq9 = dist(a, 0, b, h, x3, y3) - r3^2;
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)
(a, b, h, x0, r0, x1, x2, x3, y3) = u
return [
(r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2, # eq1
(-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2, # eq2
(-r2 + y3)^2 - (r2 + r3)^2 + (-x2 + x3)^2, # eq3
-a*h + r0*(a + sqrt(b^2 + h^2) + sqrt(h^2 + (a - b)^2)), # eq4
r0/(a - x0) - r1/(a - x1), # eq5
r0/x0 - r2/x2, # eq6
-r2^2 + (-b*(b*x2 + h*r2)/(b^2 + h^2) + x2)^2 + (-h*(b*x2 + h*r2)/(b^2 + h^2) + r2)^2, # eq7
-r3^2 + (-b*(b*x3 + h*y3)/(b^2 + h^2) + x3)^2 + (-h*(b*x3 + h*y3)/(b^2 + h^2) + y3)^2, # eq8
-r3^2 + (-h*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2) + y3)^2 + (-a + x3 - (-a + b)*(h*y3 + (-a + b)*(-a + x3))/(h^2 + (-a + b)^2))^2, # eq9
]
end;
(r1, r2, r3) = (256, 225, 144) .// 2
iniv = BigFloat[1000, 355, 358, 385, 162, 513, 269, 368, 269]
res = nls(H, ini=iniv)
res |> println
2res[1][5] |> println
([1014.0, 354.88757396449705, 357.8698224852071, 384.0, 160.0, 510.0, 270.0, 367.9881656804734, 268.8284023668639], true)
320.0
大円，中円，小円の直径がそれぞれ 256 寸，225 寸，144 寸のとき，全円の直径は 320 寸である。
その他のパラメータは以下のとおりである。
a = 1014; b = 354.888; h = 357.87; x0 = 384; r0 = 160; x1 = 510; x2 = 270; x3 = 367.988; y3 = 268.828
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = (256, 225, 144) .// 2
(a, b, h, x0, r0, x1, x2, x3, y3) = res[1]
@printf("全円の直径 = %g\n", 2r0)
@printf("a = %g; b = %g; h = %g; x0 = %g; r0 = %g; x1 = %g; x2 = %g; x3 = %g; y3 = %g\n", a, b, h, x0, r0, x1, x2, x3, y3)
plot([0, a, b, 0], [0, 0, h, 0], color=:black, lw=0.5)
circle(x0, r0, r0, :orange)
circle(x1, r1, r1, :blue)
circle(x2, r2, r2, :green)
circle(x3, y3, r3)
if more == true
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(x0, r0, " 全円:r0,(x0,r0)", :black, :left, :vcenter)
point(x1, r1, "大円:r1\n(x1,r1)", :blue, :center, delta=-delta)
point(x2, r2, "中円:r2\n(x2,r2)", :green, :center, delta=-delta)
point(x3, y3, " 小円:r3,(x3,y3)", :red, :left, :vcenter)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(b, h, "(b,h)", :black, :left, :bottom, delta=delta)
end
end;
算額（その812）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
長方形内に大円，中円，小円が入っている。長方形の長辺，短辺がそれぞれ 86634 寸，77008 寸のとき，大円の直径はいかほどか。
長方形の長辺，短辺を a, b
大円の半径と中心座標を r1, (a - r1, r1)
中円の半径と中心座標を r2, (r2, b - r2)
小円の半径と中心座標を r3, (r3, r3)
とおいて以下の連立方程式を解く。
なお，SymPy の能力的に，a, b を変数のままにして解くことはできない。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, r3::positive,
a::positive, b::positive
(a, b) = (86634, 77008)
eq1 = (a - r1 - r2)^2 + (b - r2 - r1)^2 - (r1 + r2)^2
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (r2 - r3)^2 + (b - r2 - r3)^2 - (r2 + r3)^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))
res[2]
(-57756*sqrt(338 + 240*sqrt(2)) + 24065 + 81821*sqrt(169 + 120*sqrt(2)), -81821*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(338 + 240*sqrt(2)), -4813*sqrt(169 + 120*sqrt(2)) + 24065 + 57756*sqrt(2))
6 通りの解が得られるが，2 番目のものが適解である。
長方形の長辺，短辺がそれぞれ 86634 寸，77008 寸のとき，大円の直径は 53345.0001130072 寸である。
2res[2][1].evalf() |> println
53345.0001130072
その他のパラメータは以下のとおりである。
r1 = 26672.5; r2 = 21457.5; r3 = 17166.1
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (86634, 77008)
t = sqrt(169 + 120√2)
u = 81821 - 57756√2
(r1, r2, r3) = 24065 .+ (t*u, -t*u, 57756√2 - 4813t)
@printf("大円の直径 = %g; 中円の直径 = %g; 小円の直径 = %g\n", 2r1, 2r2, 2r3)
@printf("r1 = %g; r2 = %g; r3 = %g\n", r1, r2, r3)
plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
circle(a - r1, r1, r1)
circle(r2, b - r2, r2, :blue)
circle(r3, r3, r3, :green)
if more == true
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(r1, b - r2, "中円:r2,(r2,b-r2)", :blue, :center, delta=-delta)
point(r3, r3, "小円:r3,(r3,r3)", :green, :center, delta=-delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
point(0, b, " b", :black, :left, :bottom, delta=delta)
end
end;
算額（その811）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
長方形に楕円と円が入っている。長方形の長辺と短辺が 16 寸，9 寸のとき，円の直径はいかほどか。
算法助術の公式90 そのものである。
楕円の長径，短径を a，b，円の直径を d とすれば，
d^2 +ab -2d√ab - 2d(a + b) = 0
が成り立つ。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, d::positive
formula90 = d^2 + a*b - 2d*sqrt(a*b) - 2d*(a + b);
d について解けば，a, b を与えて d を求める公式になる。
d = solve(formula90, d)[1]
d |> println
sqrt(a)*sqrt(b) + a + b - sqrt(2*a^(3/2)*sqrt(b) + 2*sqrt(a)*b^(3/2) + a^2 + 2*a*b + b^2)
楕円の長径，短径を与えれば円の直径が求まる。
d(a => 16, b => 9) |> println
2
関数にしてみる。
func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);
長径，短径が 16, 9 のとき，円の直径は 2 である。
func(16, 9)
2.0
この公式，関数を用いるとき，長径，短径をあたえれば直径が得られ，長半径，短半径を与えれば半径が得られる。
func(a, b) = sqrt(a*b) + a + b - sqrt(2a^(3/2)*√b + 2√a*b^(3/2) + a^2 + 2a*b + b^2);
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b) = (16, 9) .// 2
r = func(a, b)
@printf("円の直径 = %g\n", 2r)
# 図形描画システムと統一を図るため長半径，短半径，半径にする
plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:blue, lw=0.5)
ellipse(0, 0, a, b, color=:red)
circle4(a - r, b - r, r, :green)
if more == true
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 - r, b - r, "(a-r,b-r) ", :black, :right, :bottom, delta=delta)
point(a, 0, " a", :red, :left, :bottom, delta=delta)
point(0, b, " b", :red, :left, :bottom, delta=delta)
end
end;
算額（その810）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
直角三角形内に，大円，中円，小円を入れる。中円，小円の直径がそれぞれ 4635 寸，2060 寸のとき，大円の直径はいかほどか。
直角三角形の中に3個の円を入れる類題は多い。算額（その740），算額（その453），算額（その188）
大円の半径と中心座標を r1, (r1, r2)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (r3, y3)
直角を挟む短い方の辺と長いほうの辺の長さを a, b
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, x2::positive,
r3::positive, y3::positive, a::positive,
b::positive
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (r1 - r3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = r1/(b - r1) - r2/(b - x2)
eq4 = r1/(a - r1) - r3/(a - y3)
eq5 = a + b - sqrt(a^2 + b^2) - 2r1;
#res = solve([eq1, eq2, eq3, eq4, eq5], (r1, x2, y3, a, b))
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)
(r1, x2, y3, a, b) = u
return [
(-r1 + x2)^2 + (r1 - r2)^2 - (r1 + r2)^2, # eq1
(-r1 + y3)^2 + (r1 - r3)^2 - (r1 + r3)^2, # eq2
r1/(b - r1) - r2/(b - x2), # eq3
r1/(a - r1) - r3/(a - y3), # eq4
a + b - 2*r1 - sqrt(a^2 + b^2), # eq5
]
end;
(r2, r3) = (4635, 2060) .// 2
iniv = BigFloat[3500, 9200, 7300, 8900, 20300]
res = nls(H, ini=iniv)
([3515.5003002014273, 9224.150559489277, 7321.26713972666, 8898.390201990349, 20267.383999142992], true)
大円の直径は 7031 寸である。
その他のパラメータは以下のとおりである。
r1 = 3515.5; x2 = 9224.15; y3 = 7321.27; a = 8898.39; b = 20267.4
function draw(more)
pyplot(size=(700, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (4635, 2060) .// 2
(r1, x2, y3, a, b) = res[1]
@printf("大円の直径 = %g\n", 2r1)
@printf("r1 = %g; x2 = %g; y3 = %g; a = %g; b = %g\n", r1, x2, y3, a, b)
plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
circle(r1, r1, r1)
circle(x2, r2, r2, :green)
circle(r3, y3, r3, :magenta)
if more == true
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(r1, r1, "大円:r1\n(r1,r1)", :red, :center, :bottom, delta=3delta)
point(x2, r2, " 中円:r2\n(x2,r2)", :green, :center, delta=-delta)
point(r3, y3, "小円:r3,(r3,y3) ", :blue, :left, delta=-3delta)
point(0, a, " a", :black, :left, :bottom)
point(b, 0, " b", :black, :left, :bottom)
end
end;
算額（その809）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
不等辺三角形内に全円，大円，中円，小円を入れる。大円，中円，小円の直径がそれぞれ 9 寸，4 寸，1 寸のとき，全円の直径はいかほどか。
三角形の底辺の長さを a，頂点の座標を (x0, y0)
全円の半径と中心座標を r1, (x1, r1)
大円の半径と中心座標を r2, (x2, r2)
中円の半径と中心座標を r3, (x3, r3)
小円の半径と中心座標を r4, (x4, y4)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, x1::positive, r2::positive,
x2::positive, r3::positive, x3::positive,
r4::positive, x4::positive, y4::positive,
a::positive, x0::positive, y0::positive,
d
eq1 = (x1 - x2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x1 - x4)^2 + (r1 - y4)^2 - (r1 + r4)^2
# eq4 = r1/(a - x1) - r3/(a - x3)
eq5 = r2/x2 - r1/x1
eq6 = (a + sqrt(x0^2 + y0^2) + sqrt((a - x0)^2 + y0^2))*r1 - y0*a
eq7 = dist(0, 0, x0, y0, x2, r2) - r2^2
eq8 = dist(0, 0, x0, y0, x4, y4) - r4^2
eq9 = dist(a, 0, x0, y0, x3, r3) - r3^2
eq10 = dist(a, 0, x0, y0, x4, y4) - r4^2
eq7 = numerator(apart(eq7, d));
eq8 = numerator(apart(eq8, d));
eq9 = numerator(apart(eq9, d));
eq10 = numerator(apart(eq10, d));
# res = solve([eq1, eq2, eq3, eq5, eq6, eq7, eq8, eq9, eq10],
# (r1, x1, x2, x3, x4, y4, a, x0, y0))
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)
(r1, x1, x2, x3, x4, y4, a, x0, y0) = u
return [
(r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2, # eq1
(r1 - r3)^2 - (r1 + r3)^2 + (x1 - x3)^2, # eq2
-(r1 + r4)^2 + (r1 - y4)^2 + (x1 - x4)^2, # eq3
-r1/x1 + r2/x2, # eq5
-a*y0 + r1*(a + sqrt(x0^2 + y0^2) + sqrt(y0^2 + (a - x0)^2)), # eq6
-r2^2*y0^2 - 2*r2*x0*x2*y0 + x2^2*y0^2, # eq7
-r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2, # eq8
-2*a^2*r3*y0 + a^2*y0^2 + 2*a*r3*x0*y0 + 2*a*r3*x3*y0 - 2*a*x3*y0^2 - r3^2*y0^2 - 2*r3*x0*x3*y0 + x3^2*y0^2, # eq9
-a^2*r4^2 + a^2*y0^2 - 2*a^2*y0*y4 + a^2*y4^2 + 2*a*r4^2*x0 + 2*a*x0*y0*y4 - 2*a*x0*y4^2 - 2*a*x4*y0^2 + 2*a*x4*y0*y4 - r4^2*x0^2 - r4^2*y0^2 + x0^2*y4^2 - 2*x0*x4*y0*y4 + x4^2*y0^2, # eq10
]
end;
(r2, r3, r4) = (9, 4, 1) .// 2
iniv = BigFloat[5.5, 54.7, 44.8, 61.4, 57, 11.1, 65.1, 57.21, 11.6]
res = nls(H, ini=iniv)
([5.5, 54.7243090408641, 44.7744346697979, 61.357558621574896, 56.97961389830577, 11.06, 65.14798695340964, 57.20514438404994, 11.616], true)
全円の直径は 11 寸である。
算額にはよくあることであるが，実際の図と算額に描かれている図には大きな違いがある。
その他のパラメータは以下のとおりである。
r1 = 5.5; x1 = 54.7243; x2 = 44.7744; x3 = 61.3576; x4 = 56.9796; y4 = 11.06; a = 65.148; x0 = 57.2051; y0 = 11.616
function draw(more)
pyplot(size=(800, 170), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3, r4) = (9, 4, 1) .// 2
(r1, x1, x2, x3, x4, y4, a, x0, y0) = res[1]
@printf("全円の直径 = %g\n", 2r1)
@printf("r1 = %g; x1 = %g; x2 = %g; x3 = %g; x4 = %g; y4 = %g; a = %g; x0 = %g; y0 = %g\n", r1, x1, x2, x3, x4, y4, a, x0, y0)
plot([0, a, x0, 0], [0, 0, y0, 0], color=:black, lw=0.5)
circle(x1, r1, r1)
circle(x2, r2, r2, :green)
circle(x3, r3, r3, :magenta)
circle(x4, y4, r4, :blue)
if more == true
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\n(x1,r1)", :red, :center, :bottom, delta=3delta)
point(x2, r2, "大円:r2\n(x2,r2)", :green, :center, delta=-3delta)
point(x3, r3, "中円:r3,(x3,r3) ", :black, :right, :vcenter)
point(x4, y4, "小円:r4,(x4,y4) ", :blue, :right, :bottom, delta=2delta)
point(x0, y0, " (x0,y0)", :black, :left, :vcenter)
point(a, 0, " a", :black, :left, :bottom)
end
end;
算額（その808）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
台形（半梯）の中に甲円と乙円が入っている。大頭，小頭，縦（台形の高さ）がそれぞれ 8.5 寸，4.9 寸，4.8 寸，甲円の直径が 4 寸のとき，乙円の直径はいかほどか。
大頭，小頭，縦（台形の高さ）を b, c, a
甲円の半径と中心座標を r1, (r1, r1)
乙円の半径と中心座標を r2, (a - r2, y2)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b::positive, c::positive,
r1::positive, r2::positive, y2::positive, d
(a, b, c, r1) = (48//10, 85//10, 49//10, 4//2)
eq1 = (a - r2 - r1)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = distance(0, b, a, c, a - r2, y2) - r2^2
eq2 = numerator(apart(eq2, d))
res = solve([eq1, eq2], (r2, y2))
1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
(1, 22/5)
大頭，小頭，縦（台形の高さ）がそれぞれ 8.5 寸，4.9 寸，4.8 寸，甲円の直径が 4 寸のとき，乙円の直径は 2 寸である。
実際の図は，算額の図とは大違いである。
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(a, b, c, r1) = (48//10, 85//10, 49//10, 4//2)
(r2, y2) = (1, 22/5)
@printf("乙円の直径 = %.9g\n", 2r2)
@printf("r2 = %g; y2 = %g\n", r2, y2)
plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
circle(r1, r1, r1)
circle(a - r2, y2, r2, :blue)
if more == true
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(r1, r1, "甲円:r1,(r1,r1)", :red, :center, delta=-delta)
point(a - r2, y2, "乙円:r2\n(a-r2,y2)", :blue, :center, delta=-delta)
point(0, b, " b", :black, :left, :bottom)
point(a, c, "(a,c) ", :black, :right, delta=-delta/2)
point(a, 0, "a ", :black, :right, :bottom, delta=delta/2)
end
end;
算額（その807）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
二等辺三角形内に隔線を入れ，3 個の等円を入れる。底辺が 120 寸，斜辺が 87 寸のとき，等円の直径はいかほどか。
二等辺三角形の底辺の長さを 2a, 高さを h
3 本の隔線の交点座標を (0, b)
等円の半径と中心座標を r, (0, r), (r, y)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, b, h::positive, r::positive, y::positive, d
a = 120//2
eq1 = a^2 + h^2 - 87^2
eq2 = dist(0, h, a, 0, r, y) - r^2
eq3 = (b - r)*a/sqrt(a^2 + b^2)- r
eq4 = dist(0, b, a, 0, r, y) - r^2
res = solve([eq1, eq2, eq3, eq4], (r, y, b, h))
2-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
(12, 33, 25, 63)
(180/7, 513/7, 63, 63)
2 組の解が得られるが，最初のものが適解である。
底辺が 120 寸，斜辺が 87 寸のとき，等円の直径は 24 寸である。
その他のパラメータは以下のとおりである。
r = 12; y = 33; b = 25; h = 63
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
a = 120//2
(r, y, b, h) = (12, 33, 25, 63)
@printf("等円の直径 = %.9g\n", 2r)
@printf("r = %g; y = %g; b = %g; h = %g\n", r, y, b, h)
plot([a, 0, -a, a], [0, h, 0, 0], color=:black, lw=0.5)
circle(0, r, r)
circle2(r, y, r)
segment(-a, 0, 0, b, :blue)
segment(a, 0, 0, b, :blue)
segment(0, b, 0, h, :blue)
if more == true
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(r, y, "等円:r\n(r,y)", :red, :center, delta=-delta)
point(0, r, "等円:r\n(0,r)", :red, :center, delta=-delta)
point(a, 0, " a", :blue, :center, :bottom, delta=delta)
point(0, b, " b", :blue, :center, delta=-2delta)
point(0, h, " h", :black, :left, :bottom, delta=delta/2)
end
end;
算額（その806）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
直角三角形内に，大，中，小の 3 円を入れる。直角を挟む二辺の短い方（鈎），長い方（股）の長さが 2352 寸，3689 寸のとき，小円の直径を求めよ。
鈎，股の長さを a, b
大円の半径と中心座標を r1, (r1, r1)
中円の半径と中心座標を r2, (x2, r2)
小円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, x2::positive, r3::positive, x3::positive,
a::positive, b::positive
(a, b) = (2352, 3689)
eq1 = (x2 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (x3 - r1)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = a + b - sqrt(a^2 + b^2) - 2r1
eq4 = (a + b - 2r1)^2 - (a^2 + b^2)
eq5 = r1/(b - r1) - r2/(b - x2)
eq5 = r1*(b - x2) - r2*(b - r1)
res = solve([eq1, eq2, eq3, eq4, eq5], (r1, r2, x2, r3, x3));
res[1]
(833, 7497/16, 4165/2, 153, 1547)
5 組の解が得られるが，最初のものが適解である。
鈎，股が 2352 寸，3689 寸のとき，小円の直径は 306 寸である。
その他のパラメータは以下のとおりである。
r1 = 833; r2 = 468.562; x2 = 2082.5; r3 = 153; x3 = 1547
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, x2, r3, x3) = (833, 7497/16, 4165/2, 153, 1547)
@printf("小円の直径 = %.9g\n", 2r3)
@printf("r1 = %g; r2 = %g; x2 = %g; r3 = %g; a = %gx3\n", r1, r2, x2, r3, x3)
plot([0, b, 0, 0], [0, 0, a, 0], color=:black, lw=0.5)
circle(r1, r1, r1)
circle(x2, r2, r2, :blue)
circle(x3, r3, r3, :green)
if more == true
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(0, a, " a", :black, :left, :bottom)
point(b, 0, " b", :black, :left, :bottom)
point(r1, r1, "大円:r1,(r1,r1)", :red, :center, delta=-delta)
point(x2, r2, "中円:r2,(x2,r2)", :blue, :center, delta=-delta)
point(x3, r3, " 小円:r3,(x3,r3)", :black, :left, :vcenter)
end
end;
算額（その805）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
台形内に甲円，乙円，丙円を入れる。大頭（下底）が 782 寸，小頭（上底）が 291 寸のとき，甲円の直径はいかほどか。
横倒しになっているが，台形の大頭，小頭，高さを b, c, a とおく。
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, y2)
丙円の半径と中心座標を r3, (r3, r3)
とおいて，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms r1::positive, r2::positive, y2::positive, r3::positive, x3::positive,
a::positive, b::positive, c::positive, d
eq1 = (a - r1 - r2)^2 + (y2 - r1)^2 - (r1 + r2)^2 |> expand
eq2 = (a - r1 - r3)^2 + (r1 - r3)^2 - (r1 + r3)^2 |> expand
eq3 = (r2 - r3)^2 + (y2 - r3)^2 - (r2 + r3)^2 |> expand
eq4 = dist(0, b, a, c, a - r1, r1) - r1^2
eq4 = numerator(apart(eq4, d))
eq5 = dist(0, b, a, c, r2, y2) - r2^2
eq5 = numerator(apart(eq5, d))
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)
(r1, r2, y2, r3, a) = u
return [
a^2 - 2*a*r1 - 2*a*r2 + r1^2 - 2*r1*y2 + y2^2, # eq1
a^2 - 2*a*r1 - 2*a*r3 + r1^2 - 2*r1*r3 + r3^2, # eq2
-4*r2*r3 + r3^2 - 2*r3*y2 + y2^2, # eq3
a^2*c^2 - 2*a^2*c*r1 + 2*a*b*c*r1 - 2*a*b*r1^2 - 2*a*c^2*r1 + 2*a*c*r1^2, # eq4
a^2*b^2 - 2*a^2*b*y2 - a^2*r2^2 + a^2*y2^2 - 2*a*b^2*r2 + 2*a*b*c*r2 + 2*a*b*r2*y2 - 2*a*c*r2*y2, # eq5
]
end;
(b, c) = (782, 391)
iniv = BigFloat[250, 200, 500, 175, 780]
res = nls(H, ini=iniv)
([242.00048108871312, 183.47778926651384, 484.00096217742623, 151.04933058240556, 775.4318754192511], true)
甲円の直径は 242 寸有奇である。
その他のパラメータは以下のとおりである。
r1 = 242; r2 = 183.478; y2 = 484.001; r3 = 151.049; a = 775.432
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, y2, r3, a) = res[1]
@printf("甲円の直径 = %.9g\n", 2r1)
@printf("r1 = %g; r2 = %g; y2 = %g; r3 = %g; a = %g\n", r1, r2, y2, r3, a)
plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
circle(a - r1, r1, r1)
circle(r2, y2, r2, :blue)
circle(r3, r3, r3, :green)
if more == true
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(r2, y2, "乙円:r2,(r2,y2)", :blue, :center, delta=-delta)
point(r3, r3, "丙円:r3,(r3,r3)", :green, :center, delta=-delta)
point(0, b, " b", :black, :left, :bottom)
point(a, 0, " a", :black, :left, :bottom)
point(a, c, " (a,c)", :black, :left, :bottom)
plot!(xlims=(0, 850))
end
end;
算額（その804）
藤田貞資：精要算法（下巻），天明元年(1781)
http://www.wasan.jp/seiyou/seiyou.html
長方形内に甲，乙，丙，丁の 4 円を入れる。甲円の直径が 5140 寸のとき，乙円の直径はいかほどか。
甲円の半径と中心座標を r1, (a - r1, r1)
乙円の半径と中心座標を r2, (r2, 2r1 - r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (r4, r4)
長方形の長辺を a とおく。短辺は 2r1 である。
以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms a::positive, r1::positive, r2::positive, r3::positive, x3::positive, r4::positive
eq1 = (a - r1 - r2)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = (a - r1 - x3)^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x3 - r2)^2 + (2r1 - r2 - r3)^2 - (r2 + r3)^2
eq4 = (r2 - r4)^2 + (2r1 - r2 - r4)^2 - (r2 + r4)^2
eq5 = (x3 - r4)^2 + (r3 - r4)^2 - (r3 + r4)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (a, r2, r3, x3, r4));
res[3]
(r1*(3*sqrt(28*sqrt(2) + 185) + 34*sqrt(2) + 8*sqrt(56*sqrt(2) + 370) + 187)/136, r1*(-43*sqrt(28*sqrt(2) + 185) - 187 + 10*sqrt(56*sqrt(2) + 370) + 272*sqrt(2))/(68*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-19*sqrt(28*sqrt(2) + 185) - 68*sqrt(2) + 6*sqrt(56*sqrt(2) + 370) + 221)/(17*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-58*sqrt(56*sqrt(2) + 370) - 391 + 25*sqrt(28*sqrt(2) + 185) + 612*sqrt(2))/(68*(-sqrt(28*sqrt(2) + 185) + 2*sqrt(2) + 7)), r1*(-sqrt(7*sqrt(2)/16 + 185/64) + sqrt(2)/4 + 15/8))
4 組の解が得られるが，3 番目のものが適解である。
乙円の直径は，t = √(28√2 + 185) とおくと，甲円の直径の ((7 - 4√2)t + 119 - 34√2)/136 倍である。
甲円の直径が 5140 寸のとき，乙円の直径は 3441.0001373070395 寸である。
「答」では「三千四百四十寸◯◯◯有奇」とあるが，誤記であろう。
t = √(28√2 + 185)
5140*((7 - 4√2)t + 119 - 34√2)/136
3441.0001373070395
その他のパラメータは以下のとおりである。
a = 8496.06; r2 = 1720.5; r3 = 959.736; x3 = 2785.03; r4 = 912.939
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = 5140//2
t = sqrt(28*sqrt(2) + 185)
u = 2√2 + 7 - t
(a, r2, r3, x3, r4) = r1 .* (
(3*t + 8*√2t + 187 + 34√2)/136,
(-43*t + 10*√2t - 187 + 272√2)/68u,
(-19*t + 6*√2t + 221 - 68√2)/17u,
(25*t -58*√2t - 391 + 612√2)/68u,
(-sqrt(7√2/16 + 185/64) + √2/4 + 15/8))
@printf("乙円の直径 = %g\n", 2r2)
@printf("a = %g; r2 = %g; r3 = %g; x3 = %g; r4 = %g\n", a, r2, r3, x3, r4)
plot([0, a, a, 0, 0], [0, 0, 2r1, 2r1, 0], color=:black, lw=0.5)
circle(a - r1, r1, r1)
circle(r2, 2r1 - r2, r2, :blue)
circle(x3, r3, r3, :green)
circle(r4, r4, r4, :magenta)
if more == true
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(r2, 2r1 - r2, "乙円:r2,(r2,2r1-r2)", :blue, :center, delta=-delta)
point(x3, r3, "丙円:x3\n(x3,r3)", :green, :center, delta=-delta)
point(r4, r4, "丁円:r4\n(r4,r4)", :magenta, :center, delta=-delta)
point(a, 0, " a", :black, :left, :bottom, delta=delta)
end
end;
神壁算法
- 算額（その771）
https://blog.goo.ne.jp/r-de-r/e/d86f8eafad5635b089840548e2577bfb - 算額（その770）
https://blog.goo.ne.jp/r-de-r/e/28303bccaf422cc2177040e433f30ab2 - 算額（その769）
https://blog.goo.ne.jp/r-de-r/e/1fde4dc2613ad46677b38c699f071794 - 算額（その697）
https://blog.goo.ne.jp/r-de-r/e/ab7929468b8b6664735d652fee031327 - 算額（その696）
https://blog.goo.ne.jp/r-de-r/e/6b46a88a2e9e196cf9e0c035d3794167 - 算額（その695）
https://blog.goo.ne.jp/r-de-r/e/5ce9b6cd24b57fd7d01497b6aa9a1331 - 算額（その694）
https://blog.goo.ne.jp/r-de-r/e/d8227dd71b21273a0effe1e59029b573 - 算額（その693）
https://blog.goo.ne.jp/r-de-r/e/61cdf4ecbf7104c1f85b9f86e3079710 - 算額（その692）
https://blog.goo.ne.jp/r-de-r/e/1b3e7a7043254fc1b5c573ba784967f8 - 算額（その691）
https://blog.goo.ne.jp/r-de-r/e/17c7f558bedea88e621af02d3541043f - 算額（その690）
https://blog.goo.ne.jp/r-de-r/e/30af0ab7a44f52601537b571991c8fac - 算額（その689）
https://blog.goo.ne.jp/r-de-r/e/4ae4441cff6e545596ad2a364157d53d - 算額（その687）
https://blog.goo.ne.jp/r-de-r/e/ebbd04fd41bc039e531d8e67f0eea1ae - 算額（その686）
https://blog.goo.ne.jp/r-de-r/e/20fbf15bdacfff2a94c780f3cf0fa0f3 - 算額（その650）
https://blog.goo.ne.jp/r-de-r/e/141375a967a505289641a98d59bd5f58 - 算額（その649）
https://blog.goo.ne.jp/r-de-r/e/a89c101131a973c5fcd4bd2c2b31feaf - 算額（その575）
https://blog.goo.ne.jp/r-de-r/e/36a7228a5c3c1b7ae0189f08334ed5e3 - 算額（その447）
https://blog.goo.ne.jp/r-de-r/e/e3b4d9e0d5fdc6c578c6925d586b183d
続神壁算法
- 算額（その792）
https://blog.goo.ne.jp/r-de-r/e/5fc1e2f9b4eeee6953c3cce3295a5fc0 - 算額（その791）
https://blog.goo.ne.jp/r-de-r/e/d74e1016f980b8f6556f77af065a2b52 - 算額（その790）
https://blog.goo.ne.jp/r-de-r/e/dbb15ebd39534aeb3b00793f7d864dc3 - 算額（その789）
https://blog.goo.ne.jp/r-de-r/e/c1965824c3ffb56527ca0f472779b132 - 算額（その788）
https://blog.goo.ne.jp/r-de-r/e/19175b3fed4b6ba45257a043104d7a49 - 算額（その787）
https://blog.goo.ne.jp/r-de-r/e/01ca15ba267aea8104461af6a02a4afe - 算額（その446）
https://blog.goo.ne.jp/r-de-r/e/0743744b8d4a5067267df6f06f239e8b - 算額（その445）
https://blog.goo.ne.jp/r-de-r/e/9f35b3c1b37302b50c61355f43956949 - 算額（その444）
https://blog.goo.ne.jp/r-de-r/e/62b830b5a07cbc2f693a80ab9d29d1e0 - 算額（その443）
https://blog.goo.ne.jp/r-de-r/e/08f51fefb42a8dc2a610944e746bc186 - 算額（その442）
https://blog.goo.ne.jp/r-de-r/e/6cad986da6c2eaad06db77453c767746 - 算額（その441）
https://blog.goo.ne.jp/r-de-r/e/7effc64587864aa350523addf0442972 - 算額（その440）
https://blog.goo.ne.jp/r-de-r/e/59a28abd7846fcf2d0cb8e257673a330 - 算額（その439）
https://blog.goo.ne.jp/r-de-r/e/8a2c65a3828a2ca879d7fb5b86444a12 - 算額（その438）
https://blog.goo.ne.jp/r-de-r/e/5411ee85b446d73ca00e82a2e26be622 - 算額（その437）
https://blog.goo.ne.jp/r-de-r/e/7d06276e0130f0c09bb89c49d1fc668d - 算額（その436）
https://blog.goo.ne.jp/r-de-r/e/25ceda8b47f762c54a441ebd059cd781 - 算額（その394）
https://blog.goo.ne.jp/r-de-r/e/10fe593496313a52b495d0080109b0c1
神壁算法追加
- 算額（その515）
https://blog.goo.ne.jp/r-de-r/e/894cc97a5779e7592d5f0a58cabf578e - 算額（その514）
https://blog.goo.ne.jp/r-de-r/e/548cd85b7970aa08c54aa0cfc63cc96d
- 算額（その832）
埼玉県加須市外野 棘脱地蔵堂 第1問
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その771）
神壁算法
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その769）
神壁算法
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その742）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その726）
群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm - 算額（その725）
群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm - 算額（その724）
群馬県桐生市梅田町 日枝神社 明治12年
https://gunmawasan.web.fc2.com/files/sangak-corner/hiejin.htm - 算額（その723）
群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm - 算額（その722）
群馬県桐生市天神町 桐生天満宮 明治11年
https://gunmawasan.web.fc2.com/files/sangak-corner/tenman-g1.htm - 算額（その718）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その717）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その715）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その714）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その713）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その712）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その711）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その710）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その709）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その708）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その707）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その706）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その705）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その704）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その703）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その702）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その701）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その700）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その689）
和算問題あれこれ 2 令和5年4月の問題-No.2（『神壁算法』15問）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その688）
和算問題あれこれ 2 令和5年10月の問題-No.3（『五明算法』34問）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その687）
和算問題あれこれ 2 令和5年7月の問題-No.2（『神壁算法』第53問）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その686）
和算問題あれこれ 2 令和5年11月の問題-No.2（『神壁算法』第34問）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その685）
和算問題あれこれ 2 令和5年2月の問題-No.4（『算法求積通考』第6問）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その684）
和算問題あれこれ 2 令和5年2月の問題-No.2（茨城県 板橋不動願成寺）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その683）
和算問題あれこれ 2 令和5年2月の問題-No.1（茨城県 板橋不動願成寺）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その682）
和算問題あれこれ2 令和6年2月の問題-No.2（東京都 福徳神社）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その635）
和算問題あれこれ2 令和6年3月の問題-No.1（東京都 福徳神社）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その633）
和算問題あれこれ2 令和6年3月の問題-No.2（東京都 福徳神社）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その629）
和算図形問題あれこれ 令和4年5月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その632）
和算図形問題あれこれ 令和3年11月の問題2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その631）
和算図形問題あれこれ
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その630）
和算図形問題あれこれ 令和4年8月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その628）
和算図形問題あれこれ 令和3年12月の問題2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その627）
和算図形問題あれこれ 令和4年9月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その626）
和算図形問題あれこれ 令和4年9月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その625）
和算図形問題あれこれ 令和4年10月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その624）
和算図形問題あれこれ
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その623）
和算図形問題あれこれ
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その622）
和算図形問題あれこれ 令和 2 年 12 月の問題 - No.2?
飯塚文庫『関流叚書見題部円類五十問本書』より
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その621）
和算図形問題あれこれ 令和 2 年 10 月の問題 - No.2?（『絵本工夫之錦』より））
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その620）
和算図形問題あれこれ 令和 4 年 3 月の問題 - No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その619）
和算図形問題あれこれ 令和 4 年 4 月の問題 - No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その604）
和算図形問題あれこれ 令和4年6月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その603）
和算図形問題あれこれ 令和4年11月の問題-No.2 額題輯録より
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その599）
埼玉の算額ほか
https://gunmawasan.web.fc2.com/files/saitama-sangaku-h24.html
愛宕神社の復元算額 明治１３年（部分拡大図）（加須市）
https://gunmawasan.web.fc2.com/files/sangak-corner/atago-3s.jpg - 算額（その540）
和算図形問題あれこれ - 『両式容題問』より
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その539）
和算図形問題あれこれ-加須市騎西での算額展示会（氷川神社の復元算額より）
https://gunmawasan.web.fc2.com/k-n-mondai.html - 算額（その538）
和算図形問題あれこれ - 加須市騎西での算額展示会（氷川神社の復元算額より）
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その537）
和算図形問題あれこれ - 令和4年3月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その536）
和算図形問題あれこれ - 加須市騎西での算額展示会（氷川神社の復元算額より）
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その535）
和算図形問題あれこれ - 令和3年11月の問題？
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その534）
和算図形問題あれこれ - 令和3年11月の問題3
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その533）
和算図形問題あれこれ - 令和3年11月の問題1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その532）
和算図形問題あれこれ - 令和4年1月の問題
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その531）
和算図形問題あれこれ - 令和4年2月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その530）
和算図形問題あれこれ - 令和4年2月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その529）
和算図形問題あれこれ - 令和4年4月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その528）
和算図形問題あれこれ - 令和4年5月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その527）
和算図形問題あれこれ - 令和4年8月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その526）
和算図形問題あれこれ - 令和4年8月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その525）
和算図形問題あれこれ - 令和4年7月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その524）
和算図形問題あれこれ - 令和4年10月の問題-No.2
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その523）
和算図形問題あれこれ - 令和4年11月の問題-No.1
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html - 算額（その98）
東京都渋谷区渋谷 金王八幡神社 安政6年(1859)4月
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf - 算額（その493）
東京都渋谷区 金王神社（金王八幡宮） 元治元年(1864)
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf - 算額（その396）
東京都渋谷区 金王神社（金王八幡宮） 元治元年(1864)
https://gunmawasan.web.fc2.com/files/konnou-HP.pdf - 算額（その384）
埼玉県北埼玉郡根古屋 氷川神社
https://gunmawasan.web.fc2.com/kongetu-no-mondai.html
算額（その803）
藤田貞資：精要算法（下巻），天明元年(1781) http://www.wasan.jp/seiyou/seiyou.html
三辺（大斜，中斜，小斜）が 28 寸，25 寸，17 寸の不等辺三角形の内部の 1 点から 3 つの頂点を結ぶ線分を引き，できる三角形の面積が等しくなるようにするとき，大斜と小斜の交点から件の点までの距離（甲斜）を求めよ。
大斜と小斜の交点を原点とし，大斜，中斜，小斜の長さを L, M, S とおき，三角形内部の 1 点の座標を (x1, y1)，中斜と小斜の交点を (x2, y2) として，以下の連立方程式を解く。
include("julia-source.txt");
# julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
using SymPy
@syms x1, y1, x2, y2, L, M, S, d
eq1 = x2^2 + y2^2 - S^2
eq2 = (L - x2)^2 + y2^2 - M^2
eq3 = S*sqrt(dist(0, 0, x2, y2, x1, y1)) - L*y1
eq3 = dist(0, 0, x2, y2, x1, y1) - (L*y1/S)^2
eq3 = numerator(apart(eq3, d))
eq4 = M*sqrt(dist(L, 0, x2, y2, x1, y1)) - L*y1
eq4 = dist(L, 0, x2, y2, x1, y1) - (L*y1/M)^2
eq4 = numerator(apart(eq4, d))
res = solve([eq1, eq2, eq3, eq4], (x1, y1, x2, y2));
res[2]
((3*L^2 - M^2 + S^2)/(6*L), sqrt(-L^4 + 2*L^2*M^2 + 2*L^2*S^2 - M^4 + 2*M^2*S^2 - S^4)/(6*L), (L^2 - M^2 + S^2)/(2*L), sqrt(-(L - M - S)*(L - M + S)*(L + M - S)*(L + M + S))/(2*L))
2 組の解が得られるが， 2 番目のものが適解である。
甲斜はピタゴラスの定理で計算する。
甲斜は，大斜の 2 乗の 2 倍と，小斜の 2 乗の 2 倍の和から中斜の 2 乗を引き，平方根を取ったものを 3 で割ることで得られる。
甲斜 = sqrt(res[2][1]^2 + res[2][2]^2) |> simplify
甲斜 |> println
sqrt(2*L^2 - M^2 + 2*S^2)/3
大斜，中斜，小斜が 28 寸，25 寸，17 寸のとき，甲斜は 13 寸である。
甲斜(L => 28, M => 25, S => 17).evalf() |> println
13.0000000000000
function draw(more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(L, M, S) = (28, 25, 17)
(x1, y1, x2, y2) = ((3*L^2 - M^2 + S^2)/(6*L), sqrt(-L^4 + 2*L^2*M^2 + 2*L^2*S^2 - M^4 + 2*M^2*S^2 - S^4)/(6*L), (L^2 - M^2 + S^2)/(2*L), sqrt(-(L - M - S)*(L - M + S)*(L + M - S)*(L + M + S))/(2*L))
plot([0, 28, x2, 0], [0, 0, y2, 0], color=:magenta, lw=0.5)
segment(0, 0, x1, y1, :red)
segment(28, 0, x1, y1, :blue)
segment(x2, y2, x1, y1, :blue)
if more == true
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, y1, "(x1,y1)", :blue, :left, :bottom, delta=delta)
point(x2, y2, "(x2,y2)", :magenta, :left, :bottom, delta=delta)
point(L/2, 0, "大斜", :black, :center, :bottom, delta=delta, mark=false)
point((L + x2)/2, y2/2, "中斜", :black, :left, :bottom, delta=delta, mark=false)
point(x2/2, y2/2, "小斜", :black, :right, :bottom, delta=delta, mark=false)
point(x1/2, y1/2, "甲斜", :red, :right, :bottom, delta=delta, mark=false)
end
end;
PVアクセスランキング にほんブログ村
