裏 RjpWiki

算額（その1583）

福島県田村郡三春町御木沢厳島神社
三春まちなか寺子屋：厳島神社新算額
https://miharukoma.com/wp-content/uploads/2018/11/12.22厳島神社-新算額.pdf
https://miharukoma.com/wp-content/uploads/2018/12/寺子屋12.22　問題解説.pdf

キーワード：円5個，半円2個，長方形
#Julia, #SymPy, #算額, #和算, #数学

長方形の中に，半円 2 個，甲円 1 個，乙円 2 個，丙円 2 個を容れる。半円の半径が 48 寸のとき，乙円，丙円の半径はいかほどか。

半円の半径と中心座標を R, (0, R/2)
甲円の半径と中心座標を r1, (0, 0)
丙円の半径と中心座標を r3, (R - r3, 0)
乙円の半径と中心座標を r2, (r2, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r1, r2, r3

eq1 = r2^2 + r1^2 - (R - r2)^2
eq2 = (R - r3)^2 + r1^2 - (R + r3)^2
eq3 = 2r1^2 - R^2
res = solve([eq1, eq2, eq3], (r1, r2, r3))[2] # 2 of 2

(sqrt(2)*R/2, R/4, R/8)

# r1
res[1] |> println

sqrt(2)*R/2

# r2
res[2] |> println

R/4

# r3
res[3] |> println

R/8

甲円，乙円，丙円の半径は，半円の直径の 1/√2, 1/4, 1/8 倍である。

半円の半径が 48 寸のとき，甲円，乙円，丙円の半径は 33.94112549695428 寸， 12 寸，6 寸である。

using Colors
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, r2, r3) = R .* (1/√2, 1/4, 1/8)
@printf("半円の半径が %g のとき，甲円，乙円，丙円の半径は %g, %g, %g である。\n", R, r1, r2, r3)
plot()
circlef(0, r1, R, parse(Colorant, "#FFF565"), beginangle=180, endangle=360)
circlef(0, -r1, R, parse(Colorant, "#FFF565"), beginangle=0, endangle=180)
circlef(0, 0, r1, parse(Colorant, "#FF8B3E"))
x = [R*cosd(z) for z in 45:0.01:135]
y = [R*sind(z) - r1 for z in 45:0.01:135]
append!(x, reverse(x))
append!(y, -reverse(y))
plot!(x, y, seriestype=:shape, color=parse(Colorant, "#FFFCC4"), fcolor=parse(Colorant, "#FFFCC4"), lw=0.5)
circle2f(R - r3, 0, r3, parse(Colorant, "#F02F22"))
circle2f(r2, 0, r2, parse(Colorant, "#757CF7"))
circle(0, r1, R, :black, beginangle=180, endangle=360)
circle(0, -r1, R, :black, beginangle=0, endangle=180)
circle(0, 0, r1, :black)
circle2(R - r3, 0, r3, :black)
circle2(r2, 0, r2, :black)
rect(-R, -r1, R, r1, :black)
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(0, r1, "半円:R,(0,R/2)", :black, :center, :bottom, delta=delta)
point(R - r3, 0, "丙円:r2,(R-r3,0)", :black, :right, delta=-2delta, deltax=4delta)
point(r2, 0, "乙円:r2,(r2,0)", :black, :center, :bottom, delta=delta)
ylims!(-r1 - 2delta, r1 + 2delta)
end
end;

draw(48, true)

頑張ってます。クリックお願いします。

算額（その1582）

算額（その1582）

キーワード：円4個，半円2個，長方形
#Julia, #SymPy, #算額, #和算, #数学

長方形の中に，半円 2 個，乙円 2 個，丙円 2 個を容れる。半円の半径が 48 寸のとき，乙円，丙円の半径はいかほどか。

半円の半径と中心座標を R, (0, R/2)
丙円の半径と中心座標を r3, (R - r3, 0)
乙円の半径と中心座標を r2, (r2, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r2, r3
eq1 = r2^2 + (R/2)^2 - (R - r2)^2
eq2 = (R - r3)^2 + (R/2)^2 - (R + r3)^2;

res = solve([eq1, eq2], (r2, r3));

# r2
res[r2] |> println

3*R/8

# r3
res[r3] |> println

R/16

乙円，丙円の半径は，半円の直径の 3/8, 1/16 倍である。

半円の半径が 48 寸のとき，乙円，丙円の半径は 18 寸，3 寸である。

using Colors
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (3R/8, R/16)
@printf("半円の半径が %g のとき，乙円，丙円の半径は %g, %g である。\n", R, r2, r3)
plot()
circlef(0, R/2, R, parse(Colorant, "#FFF565"), beginangle=180, endangle=360)
circlef(0, -R/2, R, parse(Colorant, "#FFF565"), beginangle=0, endangle=180)
circle2f(R - r3, 0, r3, parse(Colorant, "#F02F22"))
circle2f(r2, 0, r2, parse(Colorant, "#757CF7"))
circle(0, R/2, R, :black, beginangle=180, endangle=360)
circle(0, -R/2, R, :black, beginangle=0, endangle=180)
circle2(R - r3, 0, r3, :black)
circle2(r2, 0, r2, :black)
rect(-R, -R/2, R, R/2, :black)
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(0, R/2, "半円:R,(0,R/2)", :black, :center, :bottom, delta=delta/2)
point(R - r3, 0, "丙円:r2,(R-r3,0)", :black, :right, delta=-2delta, deltax=-2delta)
point(r2, 0, "乙円:r2,(r2,0)", :black, :center, :bottom, delta=delta/2)
ylims!(-R/2 - 2delta, R/2 + 2delta)
end
end;

draw(48, true)

頑張ってます。クリックお願いします。

算額（その1581）

算額（その1581）

キーワード：円4個，半円2個，長方形
#Julia, #SymPy, #算額, #和算, #数学

長方形の中に，半円 2 個，正三角形 4 個，乙円 2 個，丙円 2 個を容れる。半円の半径が 48 寸のとき，乙円，丙円の半径はいかほどか。

半円の半径と中心座標を R, (0, √3R/2)
丙円の半径と中心座標を r3, (R - r3, 0)
乙円の半径と中心座標を r2, (r2, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, r2, r3
eq1 = (R - r3)^2 + (√Sym(3)*R/2)^2 - (R + r3)^2
eq2 = r2^2 + (√Sym(3)*R/2)^2 - (R - r2)^2;

# r2
ans_r2 = solve(eq2, r2)[1]
ans_r2 |> println

R/8

# r3
ans_r3 = solve(eq1, r3)[1]
ans_r3 |> println

3*R/16

乙円，丙円の半径は，半円の直径の 1/8, 3/16 倍である。

半円の半径が 48 寸のとき，乙円，丙円の半径は 6 寸，9 寸である。

using Colors
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (R/8, 3R/16)
@printf("半円の半径が %g のとき，乙円，丙円の半径は %g, %g である。\n", R, r2, r3)
plot()
circlef(0, √3R/2, R, parse(Colorant, "#FF8B3E"), beginangle=180, endangle=360)
circlef(0, -√3R/2, R, parse(Colorant, "#FF8B3E"), beginangle=0, endangle=180)
polygon(R/2, R/√3, R/√3, 3, φ=30, color=parse(Colorant, "#FFF565"), fcolor=parse(Colorant, "#FFF565"))
polygon(-R/2, R/√3, R/√3, 3, φ=30, color=parse(Colorant, "#FFF565"), fcolor=parse(Colorant, "#FFF565"))
polygon(R/2, -R/√3, R/√3, 3, φ=210, color=parse(Colorant, "#FFF565"), fcolor=parse(Colorant, "#FFF565"))
polygon(-R/2, -R/√3, R/√3, 3, φ=210, color=parse(Colorant, "#FFF565"), fcolor=parse(Colorant, "#FFF565"))
x = [R*cosd(z) for z in 60:0.01:120]
y = [R*sind(z) - √3R/2 for z in 60:0.01:120]
append!(x, reverse(x))
append!(y, -reverse(y))
circle2f(R - r3, 0, r3, parse(Colorant, "#F02F22"))
plot!(x, y, seriestype=:shape, color=parse(Colorant, "#FFFCC4"), fcolor=parse(Colorant, "#FFFCC4"), lw=0.5)
circle2f(r2, 0, r2, parse(Colorant, "#757CF7"))
circle(0, √3R/2, R, :black, beginangle=180, endangle=360)
circle(0, -√3R/2, R, :black, beginangle=0, endangle=180)
circle2(R - r3, 0, r3, :black)
circle2(r2, 0, r2, :black)
rect(-R, -√3R/2, R, √3R/2, :black)
plot!([-R, 0, R], √3R/2 .* [-1, 1, -1], color=:black, lw=0.5)
plot!([-R, 0, R], √3R/2 .* [1, -1, 1], color=:black, lw=0.5)
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(0, √3R/2, "半円:R,(0,√3R/2)", :black, :center, :bottom, delta=delta/2)
point(R - r3, 0, "丙円:r2\n(R-r3,0)", :black, :center, delta=-delta/2)
point(r2, 0, "乙円:r2,(r2,0)", :black, :left, delta=-5delta)
end
end;

draw(48, true)

頑張ってます。クリックお願いします。

算額（その1580）

算額（その1580）

キーワード：円4個，半円2個，正方形
#Julia, #SymPy, #算額, #和算, #数学

正方形の中に半円 2 個，乙円 2 個，丙円 2 個を容れる。半円の半径が 48 寸のとき，乙円，丙円の半径はいかほどか。

半円の半径と中心座標を R, (0, R)
丙円の半径と中心座標を r3, (R - r3, 0)
乙円の半径と中心座標を r2, (R - 2r3 - r2, 0)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r2::positive, r3::positive
eq1 = (R - r3)^2 + R^2 - (R + r3)^2
eq2 = (R - 2r3 - r2)^2 + R^2 - (R + r2)^2
res = solve([eq1, eq2], (r2, r3))[1]

(R/12, R/4)

乙円，丙円の半径は，半円の半径の 1/12, 1/4 倍である。

半円の半径が 48 寸のとき，乙円，丙円の半径は 4 寸，12 寸である。

using Colors
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r2, r3) = (R/12, R/4)
@printf("半円の半径が %g のとき，乙円，丙円の半径は %g, %g である。\n", R, r2, r3)
plot()
circlef(0, R, R, parse(Colorant, "#FFF565"), beginangle=180, endangle=360)
circlef(0, -R, R, parse(Colorant, "#FFF565"), beginangle=0, endangle=180)
circle2f(R - r3, 0, r3, parse(Colorant, "#F02F22"))
circle2f(R - 2r3 - r2, 0, r2, parse(Colorant, "#757CF7"))
circle(0, R, R, :black, beginangle=180, endangle=360)
circle(0, -R, R, :black, beginangle=0, endangle=180)
circle2(R - r3, 0, r3, :black)
circle2(R - 2r3 - r2, 0, r2, :black)
rect(-R, -R, R, R, :black)
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(0, R, "半円:R,(0,R)", :black, :center, :bottom, delta=delta/2)
point(R - r3, 0, "丙円:r2\n(R-r3,0)", :black, :center, delta=-delta/2)
point(R - 2r3 - r2, 0, "乙円:r2,(R-2r3-r2,0)", :black, :right, delta=-3delta)
end
end;

draw(48, true)

頑張ってます。クリックお願いします。

算額（その1579）

2025年02月02日 | Julia

算額（その1579）

福島県田村市日渡神社明治21年(1888)
https://www.city.tamura.lg.jp/soshiki/30/bunkazai.html
キーワード：円9個，外円
#Julia, #SymPy, #算額, #和算, #数学

外円の中に，大円 2 個と，小円 6 個を容れる。大円の直径が与えられたときに小円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, y21), (r2, y22)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive,
r2::positive, y21::positive, y22::positive
r1 = R/2
eq1 = r1^2 + y21^2 - (r1 + r2)^2
eq2 = r2^2 + (y22 - y21)^2 - 4r2^2
eq3 = r2^2 + y22^2 - (R - r2)^2
res = solve([eq1, eq2, eq3], (r2, y21, y22))[1]

(R*(-5 + sqrt(33))/4, R*(-5*sqrt(2) + sqrt(66))/(2*sqrt(7 - sqrt(33))), R*sqrt(7/2 - sqrt(33)/2))

小円の半径 r2 は，大円の半径 R の (√33 - 5)/4 倍である。
大円の直径が 3954 のとき，736.0001761028412 である。

3954*(√33 - 5)/4

736.0001761028412

using Colors
function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r1 = R/2
(r2, y21, y22) = (R*(-5 + sqrt(33))/4, R*(-5*sqrt(2) + sqrt(66))/(2*sqrt(7 - sqrt(33))), R*sqrt(7/2 - sqrt(33)/2))
@printf("外円の直径が %g のとき，小円の直径は %g である。\n", 2R, 2r2)
plot()
circlef(0, 0, R, RGB(194/255, 214/255, 236/255))
circle2f(r1, 0, r1, RGB(234/255, 51/255, 35/255))
circle22f(0, y21, r2, RGB(229/255, 74/255, 233/255))
circle4f(r2, y22, r2, RGB(229/255, 74/255, 233/255))
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(0, R, "R", :green, :center, :bottom, delta=delta/2)
point(r1, 0, "大円:r1,(r1,0)", :white, :center, delta=-delta/2)
point(0, y21, "小円:r2,(0,y21)", :black, :left, :bottom, delta=delta)
point(r2, y22, " 小円:r2,(r2,y22)", :black, :left, :bottom, delta=delta)
end
end;

draw(1977, true)

頑張ってます。クリックお願いします。

算額（その1578）

算額（その1578）

『算法續淺問答巻之七』
terasima3.pdf
http://www.wasan.jp/terasima/terasima.html
http://www.wasan.jp/terasima/terasima3.pdf

キーワード：円環，外円
#Julia, #SymPy, #算額, #和算

図のように 8 個の等円が互いに外接し，外円に内接している。外円の直径が 10 寸のとき，等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (R - r, 0)
等円の個数を n
とおき，以下の方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r::positive, n::positive
eq = (R - r)*sind(180/n) - r
res = solve(eq, r)[1]
res |> println

R*sin(pi/n)/(sin(pi/n) + 1)

n = 8, R = 10/2 を代入し，得られた半径を 2 倍すれば直径は 2.76768653914155 である。

2*res(n => 8, R => Sym(10)/2) |> println
2*res(n => 8, R => Sym(10)/2).evalf() |> println

10*sqrt(1/2 - sqrt(2)/4)/(sqrt(1/2 - sqrt(2)/4) + 1)
2.76768653914155

function draw(R, n, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = R*sin(pi/n)/(sin(pi/n) + 1)
p = plot()
circle(0, 0, R, :green)
circle(0, 0, R - r, :gray80)
rotate(0, R - r, r, angle=360/n)
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(0, R, "R", :green, :center, :bottom, delta=delta/2)
point(0, R - r, "R-r", :gray, :center, :bottom, delta=delta/2)
end
return p
end;

draw(10, 8, true)

頑張ってます。クリックお願いします。

算額（その1577）

算額（その1577）

愛知県名古屋市中央区大須観音（北野山真福寺寶生院）天保8年(1837)
藤安順：所掲大須観音堂解義
山田潤：動的幾何ソフトウエアを利用した平面図形の作図についての一考察 - 和算所にある平面図形問題の利用 -
2023018-Yamada.pdf
https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2273-19.pdf

キーワード：円10個，外円
#Julia, #SymPy, #算額, #和算

外円の中に大円 4 個，甲円 3 個，乙円を 2 個容れる。甲円の直径が 1 寸のとき，乙円および外円の直径はいかほどか。

算額（その1575）は，深川が意図的に大円が交差する部分にある円を省略し，甲円と乙円の関係を述べさせるものであったようだ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, R - r1)
甲円の半径と中心座標を r2, (x2, y2)
乙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, x1::positive,
r2::positive, r3::positive
x1 = r1 - r2
eq1 = R - (2r1 + r2)
eq2 = x1^2 + r1^2 - (r1 + r2)^2 # 大円と甲円が接する
eq3 = x1^2 + (R - r3 - r1)^2 - (r1 + r3)^2; # 大円と乙円が接する

res = solve([eq1, eq2, eq3], (R, r1, r3))[1]

(9*r2, 4*r2, r2)

外円の直径は甲円の直径の 9 倍，大円の直径は甲円の直径の 4 倍，乙円の直径は甲円の直径に等しい。

考察：大円の交差するところに甲円が入っているのが，算額（その1575）での y1 の制約条件になっている。

function draw(r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(R, r1, r3) = (9*r2, 4*r2, r2)
plot()
circle(0, 0, R, :green)
circle4(r1 - r2, r1, r1, :blue)
circle(0, 0, r2)
circle22(0, r1, r3)
circle22(0, R - r3, r3, :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(0, R, "R", :green, :center, :bottom, delta=delta/2)
point(r1 - r2, r1, "大円:r1,(r1-r2,r1)", :blue, :center, delta=-delta, deltax=5delta)
point(0, R - r3, "乙円:r3,(0,R-r3)", :magenta, :left, :bottom, delta=delta/2, deltax=4delta)
point(0, r1, "甲円:r2,(0,r1)", :red, :right, :vcenter, deltax=-4delta)
end
end;

draw(1/2, true)

頑張ってます。クリックお願いします。

算額（その1576）

算額（その1576）

佐久間纉：算法起源集,1877.
https://kokusho.nijl.ac.jp/biblio/100234582/69?ln=ja

山田潤：動的幾何ソフトウエアを利用した平面図形の作図についての一考察 - 和算所にある平面図形問題の利用 -
2023018-Yamada.pdf
https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2273-19.pdf

キーワード：円7個，外円
#Julia, #SymPy, #算額, #和算

外円の中に甲円，乙円，丙円を 2 個ずつ容れる。甲円，乙円，丙円の直径が与えられたとき，外円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive,
r2::positive, x2::positive, y2::negative,
r3::positive, x3::positive, y3::positive
eq1 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq2 = x3^2 + y3^2 - (R - r3)^2 |> expand
eq3 = x2^2 + (y2 - r1 + R)^2 - (r1 + r2)^2 |> expand
eq4 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2 |> expand
eq5 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2 |> expand;

まず，eq1, eq3 の連立方程式から x2, y2 を求める。

res2 = solve([eq1, eq3], (x2, y2))[1]

(-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2)/(-R + r1), -(-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1))

次に，eq2, eq4 の連立方程式から x3, y3 を求める。

res3 = solve([eq2, eq4], (x3, y3))[1]

(-2*sqrt(R)*sqrt(r1)*sqrt(r3)*sqrt(R - r1 - r3)/(-R + r1), (-R^2 + R*r1 + R*r3 + r1*r3)/(-R + r1))

eq5 の x2, y2, x3, y3 に求められた解を代入し，整理する。

eq15 = eq5(x2 => res2[1], y2 => res2[2], x3 => res3[1], y3 => res3[2])/R |> simplify |> numerator;

しかし，このままでは SymPy では解けない。

eq15 |> println

4*R^3 - 8*R^2*r1 - 4*R^2*r2 - 4*R^2*r3 + 4*R*r1^2 + 4*R*r1*r2 + 4*R*r1*r3 - 8*r1*sqrt(r2)*sqrt(r3)*sqrt(R - r1 - r2)*sqrt(R - r1 - r3) + 8*r1*r2*r3

8*r1*sqrt(r2)*sqrt(r3)*sqrt(R - r1 - r2)*sqrt(R - r1 - r3) の項とそれ以外の項をそれぞれ二乗し，方程式を再構成する。

eq25 = (4*R^3 - 8*R^2*r1 - 4*R^2*r2 - 4*R^2*r3 + 4*R*r1^2 + 4*R*r1*r2 + 4*R*r1*r3 + 8*r1*r2*r3)^2 - ( - 8*r1*sqrt(r2)*sqrt(r3)*sqrt(R - r1 - r2)*sqrt(R - r1 - r3))^2;

ここまで来ると SymPy でも解ける。

res = solve(eq25)[2]; # 2 of 5
res[R] |> println

r1 + r2 + r3

外円の半径は，甲円，乙円，丙円の半径の合計である。

なお，x2 = x3 であることがわかる。

# x2
res2[1](R => r1 + r2 + r3) |> println
# x3
res3[1](R => r1 + r2 + r3) |> println

-2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(-r2 - r3)
-2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(-r2 - r3)

function draw(r1, r2, r3, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
R = r1 + r2 + r3
(x2, y2) = (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2)/(-R + r1), -(-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1))
(x3, y3) = (-2*sqrt(R)*sqrt(r1)*sqrt(r3)*sqrt(R - r1 - r3)/(-R + r1), (-R^2 + R*r1 + R*r3 + r1*r3)/(-R + r1))
@printf("r1 = %g; r2 = %g; r3 = %g; R = %g; x2 = %g; y2 = %g; x3 = %g; y3 = %g\n",
r1, r2, r3, R, x2, y2, x3, y3)
plot()
circle(0, 0, R, :green)
rotate(0, R - r1, r1, angle=180)
rotate(x2, y2, r2, :blue, angle=180)
rotate(x3, y3, r3, :magenta, angle=180)
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(0, R, "R", :green, :center, :bottom, delta=delta/2)
point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta/2)
point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, delta=-delta/2)
end
end;

draw(0.5, 0.3, 0.2, true)

頑張ってます。クリックお願いします。

算額（その1575）

算額（その1575）

愛知県名古屋市中央区大須観音（北野山真福寺寶生院）天保8年(1837)
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 102，森北出版株式会社，2010年４月22日.

キーワード：円8個，外円
#Julia, #SymPy, #算額, #和算

外円の中に，大円 4 個，小円 3 個を容れる。外円と小円の直径が与えられたとき，大円の直径はいかほどか。

深川は，「中央の小円1と左右の小円2の直径の関係を示せ」という問題としている。
図を見ただけで直径は等しいとわかるが，ちゃんと式で示せということではある。
大円の中心の y 座標はどんな非負の値も取れるが，3個の小円の直径は同じであるということで，算額問題としては，「大円の直径を外円と小円の直径であらわせ」という問題にした。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R::positive, r1::positive, y1::positive, r2::positive, r3::positive
eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r1^2 + y1^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, y1))[1]

(-(-R + r2)/2, sqrt(R)*sqrt(r2))

大円の直径は，外円の直径から小円の直径を差し引いたものの半分である。

外円の直径が同じでも，小円の直径により印象の異なる図が得られる。

function draw(R, r2, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
(r1, y1) = ((R - r2)/2, sqrt(R)*sqrt(r2))
r3 = r2
p = plot()
circle(0, 0, R, :green)
circle4(r1, y1, r1, :magenta)
circle2(R - r2, 0, r2, :blue)
circle(0, 0, r3)
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(R, 0, "R", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
point(R - r2, 0, "R-r2", :blue, :center, delta=-delta)
point(r3, 0, "r3", :red, :right, :vcenter, deltax=-delta/2)
point(r1, y1, "大円:r1,(r1,y1)", :magenta, :center, :bottom, delta=delta/2)
pos = max(R, r2)
str = @sprintf(" R=%.3g, r2=%.3g", R, r2)
point(0, pos, str, :black, :center, :bottom, delta=delta, mark=false)
end
return p
end;

draw(1, 0.1, true)

頑張ってます。クリックお願いします。

算額（その1574）

算額（その1574）

三重県鳥羽市観世音寛政12年(1800)
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 97，森北出版株式会社，2010年４月22日.

キーワード：円1個，二等辺三角形2個
#Julia, #SymPy, #算額, #和算

外円の中に大小の正三角形を容れる。外円の直径が与えられたとき，小正三角形の一辺の長さはいかほどか。

外円の半径と中心座標を R, (0, 0)
小正三角形の一辺の長さを q
とおき，以下の方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms R, q
eq = (q/2)^2 + (-R/2 - √Sym(3)*q/2)^2 - R^2
res = solve(eq, q)[1]
res |> println

R*(-sqrt(3) + sqrt(15))/4

小正三角形の一辺の長さ q は，外円の半径 R の (√15 - √3)/4 倍である。
外円の直径が 100 のとき，小正三角形の一辺の長さは 26.76165673298175 である。

100/2*(√15 - √3)/4

26.76165673298175

function draw(R, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
q = R*(√15 - √3)/4
@printf("外円の直径が = %g のとき，小正三角形の一辺の長さは %g である。\n", 2R, q)
plot()
circle(0, 0, R)
polygon(0, 0, R, 3, color=:blue)
polygon(0, -R/2 - q/√3, q/√3, 3, color=: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(0, R, "R", :red, :center, :bottom, delta=delta/2)
point(0, -R/2, "-R/2", :blue, :center, :bottom, delta=delta/2)
point(√3R/2, -R/2, "(√3R/2,-R/2) ", :blue, :right, :bottom, delta=delta/2)
point(q/2, -R/2 - √3q/2, "(q/2,-R/2-√3q/2)", :blue, :left, delta=-delta/2)
end
end;

draw(100/2, true)

頑張ってます。クリックお願いします。

算額（その1573）

算額（その1573）

岡山県瀬戸内市長船町土師宮森片山日子神社明治6年(1873)
http://www.wasan.jp/okayama/katayamahiko.html
深川英俊，トニー・ロスマン：聖なる数学：算額，p. 96，森北出版株式会社，2010年４月22日.

キーワード：円1個，二等辺三角形
#Julia, #SymPy, #算額, #和算

二等辺三角形に円を容れる。底辺と斜辺が与えられたとき，円の直径はいかほどか。

底辺と斜辺を a, b
円の半径と中心座標を r, (0, r)
とおき，以下の連立方程式を解く。

include("julia-source.txt")
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms a, b, r, h
h = sqrt(b^2 - (a/2)^2)
eq1 = ((a + 2b)*r/2) - (h*a/2)
res = solve(eq1, r)[1]
res |> println

a*sqrt(-a^2 + 4*b^2)/(2*(a + 2*b))

円の半径は a*sqrt(4b^2 - a^2)/(2a + 4b) である。
底辺が 12，斜辺が 10 のとき，内接円の直径は 6 である。

2*res(a => 12, b => 10) |> println

以下のほうが簡単。

@syms a, b, r, h
h = sqrt(b^2 - (a/2)^2)
eq2 = (a/2)/b - r/(h - r)
res2 = solve(eq2, r)[1]
res2 |> println

a*sqrt(-a^2 + 4*b^2)/(2*(a + 2*b))

function draw(a, b, more)
pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
r = a*sqrt(4b^2 - a^2)/(2a + 4b)
@printf("a = %g; b = %g; r = %g\n", a, b, r)
plot([a/2, 0, -a/2, a/2], [0, sqrt(b^2 - (a/2)^2), 0, 0], color=:blue, lw=0.5)
circle(0, r, r)
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(0, sqrt(b^2 - a^2/4), "sqrt(b^2-a^2/4)", :blue, :center, :bottom, delta=delta/2)
point(a/2, 0, "a/2", :blue, :left, :bottom, delta=delta/2)
point(0, r, " r", :red, :left, :vcenter)
end
end;

draw(12, 10, true)

頑張ってます。クリックお願いします。

Julia の歪度と尖度

1.　Julia の歪度と尖度は古典的定義に基づく

1.1.　skewness, kurtosis は StatsBase に含まれる。

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
"""
skewness(v, [wv::AbstractWeights], m=mean(v))

Compute the standardized skewness of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
"""
kurtosis(v, [wv::AbstractWeights], m=mean(v))

Compute the excess kurtosis of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

1.2.　R は 3 通りの定義に基づく関数値を返す

Joanes and Gill (1998) [^1]discuss three methods for estimating kurtosis and skewness.

[1]: D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189.

1.2.1.　Kurtosis

Type 1: This is the typical definition used in many older textbooks.

g2 = m4/m2^2 - 3

Type 2: Used in SAS and SPSS.

G2 = ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))

Type 3: Used in MINITAB and BMDP.

b2 = m4/s^4 − 3 = (g2 + 3)*(1 − 1/n)^2 − 3

Only G2 (corresponding to type = 2) is unbiased under normality.

デフォルトは type=3 である。

1.2.2.　Skewness

Type 1: This is the typical definition used in many older textbooks.

g1 = m3/m2^(3/2)

Type 2: Used in SAS and SPSS.

G1 = g1*sqrt(n*(n - 1))/(n - 2)

Type 3: Used in MINITAB and BMDP.

b1 = m3/s^3 = g1*((n - 1)/n)^(3/2)

All three skewness measures are unbiased under normality.

デフォルトは type=3 である。

1.3.　Julia でも type を選べる関数を定義する

1.3.1.　尖度

using StatsBase

kurt1(x) = kurt(x, type=1)
kurt2(x) = kurt(x, type=2)
kurt3(x) = kurt(x)
function kurt(x; type=3)
doc = """
kurt(x; type=n)
type: an integer between 1 and 3 selecting one of the algorithms for computing kurtosis detailed below.

Type 1: g2. This is the typical definition used in many older textbooks.
Type 2: G2. Used in SAS and SPSS.
Type 3: b2. default. Used in MINITAB and BMDP.
"""
type in 1:3 || error(doc)
n = length(x)
g2 = kurtosis(x)
if type == 2
return ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))
elseif type== 3
return (g2 + 3)*(1 − 1/n)^2 − 3
else
return g2
end
end;

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1];

kurt(x), kurt3(x)

(-1.6335639958376689, -1.6335639958376689)

kurt(x, type=1), kurt1(x)

(-1.3130419701699616, -1.3130419701699616)

kurt(x, type=2), kurt2(x)

(-1.356984911550468, -1.356984911550468)

kurt(x, type=3), kurt3(x)

(-1.6335639958376689, -1.6335639958376689)

1.3.2.　歪度

using StatsBase

skew1(x) = skew(x, type=1)
skew2(x) = skew(x, type=2)
skew3(x) = skew(x)
function skew(x; type=3)
doc = """
skew(x; type=n)
type: an integer between 1 and 3 selecting one of the algorithms for computing skewness detailed below.

Type 1: g1. This is the typical definition used in many older textbooks.
Type 2: G1. Used in SAS and SPSS.
Type 3: b1. default. Used in MINITAB and BMDP.
"""
type in 1:3 || error(doc)
n = length(x)
g1 = skewness(x)
if type == 2
return g1*sqrt(n*(n - 1))/(n - 2)
elseif type== 3
return g1*((n - 1)/n)^(3/2)
else
return g1
end
end;

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1]
skew(x), skew3(x)

(0.10905343709973031, 0.10905343709973031)

skew(x, type=1), skew1(x)

(0.12772490663150174, 0.12772490663150174)

skew(x, type=2), skew2(x)

(0.15146310708295876, 0.15146310708295876)

skew(x, type=3), skew3(x)

(0.10905343709973031, 0.10905343709973031)

頑張ってます。クリックお願いします。

算額（その1572）