算額（その417）

算額（その417）

北海道函館市住吉町 住吉大明神（住三吉神社） 文化2年(1805)

中村信弥(2001)：幻の算額
http://www.wasan.jp/maborosi/maborosi.html

五芒星内に 5 個の等円を入れる。内側の正五角形の一辺の長さが 1 尺のとき，等円の直径を求めよ。

半径 1/(2*sind(36)) の円に内接する正五角形の一辺の長さが 1 である。
等円の半径と中心座標を r, (x, y) とし，以下の方程式の数値解を求める。

include("julia-source.txt");

using SymPy

@syms x::positive, y::positive, r::positive, a::positive;

(cosd36, sind36, cosd18, sind18) = (cosd(Sym(36)), sind(Sym(36)), cosd(Sym(18)), sind(Sym(18)))
eq1 = distance(0, a, a*cosd18, a*sind18, x, y) - r^2
eq2 = distance(0, a, -a*cosd18, a*sind18, x, y) - r^2
eq3 = distance(a*sind36, -a*cosd36, a*cosd18, a*sind18, x, y) - r^2;
# res = solve([eq1, eq2, eq3], (x, y, r))

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 v, r.f_converged
end;

function H(u)
   (r, x, y) = u
   return [
-r^2 + (y - (-5*sqrt(2)*x - sqrt(10)*x + sqrt(sqrt(5) + 5)*(sqrt(5)*a + 3*a - sqrt(5)*y + 5*y))/(8*sqrt(sqrt(5) + 5)))^2 + (-a*sqrt(2*sqrt(5) + 10)/8 - sqrt(5)*x/8 + 5*x/8 + y*sqrt(2*sqrt(5) + 10)/8)^2,  # eq1
-r^2 + (a*sqrt(2*sqrt(5) + 10)/8 - sqrt(5)*x/8 + 5*x/8 - y*sqrt(2*sqrt(5) + 10)/8)^2 + (-3*a/8 - sqrt(5)*a/8 - sqrt(2)*x*sqrt(sqrt(5) + 5)/8 + sqrt(5)*y/8 + 3*y/8)^2,  # eq2
-r^2 + (x - (-5*sqrt(10)*a/2 + 5*sqrt(2)*a/2 - 5*x*sqrt(sqrt(5) + 5) - 2*x*sqrt(25 - 5*sqrt(5)) + 5*x*sqrt(5 - sqrt(5)) + 2*x*sqrt(5*sqrt(5) + 25) - 5*sqrt(10)*y + 10*sqrt(2)*y)/(-15*sqrt(sqrt(5) + 5) + sqrt(5)*sqrt(5 - sqrt(5)) + sqrt(5)*sqrt(sqrt(5) + 5) + 15*sqrt(5 - sqrt(5))))^2 + (y - (-15*a*sqrt(5 - sqrt(5))/4 - a*sqrt(5*sqrt(5) + 25)/4 - a*sqrt(25 - 5*sqrt(5))/4 + 15*a*sqrt(sqrt(5) + 5)/4 - 5*sqrt(10)*x + 10*sqrt(2)*x - 5*y*sqrt(sqrt(5) + 5) + 5*y*sqrt(5 - sqrt(5)))/(-15*sqrt(sqrt(5) + 5) + sqrt(5)*sqrt(5 - sqrt(5)) + sqrt(5)*sqrt(sqrt(5) + 5) + 15*sqrt(5 - sqrt(5))))^2,  # eq3
   ]
end;

a = 1/2sind(36)
iniv = [big"0.36", 0.62, 0.85]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[0.3632712640026803262857421710703028576320274738964989772105963682762334090546295, 0.6180339887498946353591452335230846232288953659873432429140776623397096012159547, 0.8506508083520398708822131798566121390753734825678071178183687667072162785925422], true)

r = 0.363271;  x = 0.618034;  y = 0.850651;  a = 0.850651
正5角形の一辺の長さを 1 としたとき，等円の直径 = 0.726543

等円の半径は 0.36327126400268034，直径は 0.7265425280053607 である。

なお，0.7265425280053607 = √(5 - 2√5) = tand(36) である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (cosd36, sind36, cosd18, sind18) =
   ((1 + √5)/4, √(10 - 2√5)/4, √(10 + 2√5)/4, (√5 - 1)/4)
   a = 1/2sind(36)
   ax = zeros(5)
   ay = zeros(5)
   r0 = a*cosd36 + a*sind(36)/sind(18)*cosd(18)
   for i = 1:5
       ax[i] = cosd(72i-18)r0
       ay[i] = sind(72i-18)r0
   end
   (r, x, y) = res[1]
   @printf("r = %g;  x = %g;  y = %g;  a = %g\n", r, x, y, a)
   @printf("正5角形の一辺の長さを 1 としたとき，等円の直径 = %g\n", 2r)
   plot(a.*[0, -cosd18, -sind36, sind36, cosd18, 0], a.*[1, sind18, -cosd36, -cosd36, sind18, 1], color=:green, lw=1)
   plot!(ax[[1, 3, 5, 2, 4, 1]], ay[[1, 3, 5, 2, 4, 1]], color=:blue, lw=0.5)
   rotate(x, y, r, angle=72)
   circle(0, 0, r0, :gray85)
   circle(0, 0, a, :gray85)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, " r0", :gray, delta=-delta)
       point(x, y, "r,(x,y)", :red, :center, delta=-delta)
       point(0, a, "  a", :blue, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その416）

算額（その416）

群馬県高崎市榛名山町 榛名神社 明治33年(1900)

早坂平四郎: 算額の一考察，苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

埼玉県加須市　秋葉山本坊(秋葉宮) 文化5年(1808)2月
http://www.wasan.jp/saitama/akibahonbo1.html

 

円内に大円1個，中円2個，小円2個が入っているこれらは互いに内接・外接している他2,弦を共通接線としている。中円，小円の半径がわかっているとき，大円の半径を求めよ。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (0, y1)
中円の半径と中心座標を r2, (r2, y2)
小円の半径と中心座標を r3, (r3,  y3)
弦の端点座標を (xa, ya), (xb, yb) ただし，ya^2 = r0^2 - xa^2, yb^2 = r0^2 - xb^2
solve() では解ききらないので nlsolve() で数値解を求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, y1::positive,
     r2::positive, y2::negative, r3::positive, y3::positive,
     xa::positive, ya::negative, xb::positive, yb::positive;

(r2, r3) = (15, 5)
eq1 = r3^2 + y3^2 - (r0 - r3)^2
eq2 = r2^2 + y2^2 - (r0 - r2)^2
eq3 = r3^2 + (y3 - y1)^2 - (r1 + r3)^2
eq4 = r2^2 + (y1 - y2)^2 - (r1 + r2)^2
eq5 = distance(xa, ya, xb, yb, 0, y1) - r1^2
eq6 = distance(xa, ya, xb, yb, r2, y2) - r2^2
eq7 = distance(xa, ya, xb, yb, r3, y3) - r3^2;
eq8 = xa^2 + ya^2 - r0^2
eq9 = xb^2 + yb^2 - r0^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9], (r0, r1, y1, y2, y3, xa, xb))

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 v, r.f_converged
end;

function H(u)
   (r0, r1, y1, y2, y3, xa, ya, xb, yb) = u
   return [
       y3^2 - (r0 - 5)^2 + 25,  # eq1
       y2^2 - (r0 - 15)^2 + 225,  # eq2
       -(r1 + 5)^2 + (-y1 + y3)^2 + 25,  # eq3
       -(r1 + 15)^2 + (y1 - y2)^2 + 225,  # eq4
       -r1^2 + (y1 - (xa^2*yb - xa*xb*ya - xa*xb*yb + xb^2*ya + y1*ya^2 - 2*y1*ya*yb + y1*yb^2)/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (xa*y1*ya - xa*y1*yb - xa*ya*yb + xa*yb^2 - xb*y1*ya + xb*y1*yb + xb*ya^2 - xb*ya*yb)^2/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2)^2,  # eq5
       (y2 - (ya*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (ya - yb)*(xa^2 - xa*xb - 15*xa + 15*xb - y2*ya + y2*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (-(xa*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (xa - xb)*(xa^2 - xa*xb - 15*xa + 15*xb - y2*ya + y2*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) + 15)^2 - 225,  # eq6
       (y3 - (ya*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (ya - yb)*(xa^2 - xa*xb - 5*xa + 5*xb - y3*ya + y3*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2))^2 + (-(xa*(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) - (xa - xb)*(xa^2 - xa*xb - 5*xa + 5*xb - y3*ya + y3*yb + ya^2 - ya*yb))/(xa^2 - 2*xa*xb + xb^2 + ya^2 - 2*ya*yb + yb^2) + 5)^2 - 25,  # eq7
       -r0^2 + xa^2 + ya^2,  # eq8
       -r0^2 + xb^2 + yb^2,  # eq9
   ]
end;

(r2, r3) = (15, 5)
iniv = [big"40.0", 18, 10, -20.4, 34.1, 32, -24.4, 8.2, 39.1]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[38.7638837486628373465473353663136173085470569162791234706045756107292743498459, 17.99038105676657970145584756129404275207103940357785471041855234588949762681627, 10.95152380775283700472801642392206965400276275976845771998475084075123109744924, -18.4315536735230681094936333192018549711528600474534715446198618905612484280433, 33.39161340506353031902486710976142560545305462282349919389008399018128308922975, 32.33934584554869449875985865589356317819349516144509459022026784764907086157268, -21.37300618915924098536769340774952612678627228671598130988515134391456899286098, 8.622548477793335355655385214702161614702564480769434367924497563401166456827728, 37.79272867931013343463889124083200695874231476480976753169233139115218150527541], true)

   r0 = 38.7639;  r1 = 17.9904;  y1 = 10.9515;  y2 = -18.4316;  y3 = 33.3916;  xa = 32.3393;  ya = -21.373;  xb = 8.62255;  yb = 37.7927

中円，小円の半径をそれぞれ 15, 5 としたとき，大円の半径は 17.99038105676658 である。

Float64(res[1][2])

   17.99038105676658

術によれば，中円，小円の半径を a, b とすると，大円の「直径」は (a + b + 6√(a*b))*(a + b + 2√(a*b))/(4a*b) であるとしている。しかしこれは上述の大円の半径を2倍したものではない。何らかの不適切な式であると思われる。

(a, b) = (15, 5)
(a + b + 6√(a*b))*(a + b + 2√(a*b))/(4a*b)

   8.952135486850342

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3) = (15, 5)
   (r0, r1, y1, y2, y3, xa, ya, xb, yb) = res[1]
   @printf("r0 = %g;  r1 = %g;  y1 = %g;  y2 = %g;  y3 = %g;  xa = %g;  ya = %g;  xb = %g;  yb = %g\n", r0, r1, y1, y2, y3, xa, ya, xb, yb)
   plot()
   circle(0, 0, r0, :black)
   circle(0, y1, r1)
   circle(r2, y2, r2, :orange)
   circle(-r2, y2, r2, :orange)
   circle(r3, y3, r3, :green)
   circle(-r3, y3, r3, :green)
   segment(xa, ya, xb, yb, :blue)
   segment(-xa, ya, -xb, yb, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta)
       point(0, y1, " y1  大円:r1", :red, :left, :vcenter)
       point(r2, y2, " 中円:r2,(r2,y2)", :orange, :center, :top, delta=-delta)
       point(r3, y3, " 小円:r3,(r3,y3)", :green, :left, :vcenter)
       point(xa, ya, "(xa,ya)", :blue, :left, :top, delta=-delta)
       point(xb, yb, "(xb,yb)", :blue, :left, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その415）

算額（その415）

東京都府中市 大國魂神社 明治18年(1885)

早坂平四郎: 算額の一考察，苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

甲円と乙円が交わっており，共通部分に丙円1こと丁円2個を入れる。それぞれの円は互いに内接・外接している。
甲円，乙円，丙円の直径が与えられたとき，丁円の直径を求めよ。

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

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::negative,
     r4::positive, x4::negative, y4::positive;
eq1 = (r1 - r3 - x4)^2 + y4^2 - (r1 - r4)^2
eq2 = x4^2 + y4^2 - (r3 + r4)^2
eq3 = (r3 - r2 - x4)^2 + y4^2 - (r2 - r4)^2
res = solve([eq1, eq2, eq3], (r4, x4, y4))

   2-element Vector{Tuple{Sym, Sym, Sym}}:
    (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), -2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))
    (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), 2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))

最初の組が適解である。

res[1][1] |> simplify

甲円，乙円，丙円の半径 r1, r2, r3 がそれぞれ 55, 45, 17 のとき，丁円の半径は
(r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2) = 8.274473924977126 である。

   r1 = 55;  r2 = 45;  r3 = 17
   r4 = 8.27447;  x4 = -1.32205;  y4 = 25.2399

(r1, r2, r3) = (55, 45, 17)
r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2)

   8.274473924977126

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (55, 45, 17)
   (r4, x4, y4) = (r3*(r1 - r3)*(r2 - r3)/(r1*r2 - r3^2), -r3^2*(r1 - r2)/(r1*r2 - r3^2), 2*sqrt(r1)*sqrt(r2)*r3*sqrt(r1 - r3)*sqrt(r2 - r3)/(r1*r2 - r3^2))
   @printf("r1 = %g;  r2 = %g;  r3 = %g\n", r1, r2, r3)
   @printf("r4 = %g;  x4 = %g;  y4 = %g\n", r4, x4, y4)
   plot()
   circle(r1 - r3, 0, r1, :blue)
   circle(r3 - r2, 0, r2, :red)
   circle(0, 0, r3, :green)
   circle(x4, y4, r4, :orange)
   circle(x4, -y4, r4, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(r1 - r3, 0, "r1-r3  甲円:r1", :blue, delta=-delta)
       point(r3 - r2, 0, "乙円:r2  r3-r2 ", :red, :right, delta=-delta)
       point(0, 0, "丙円:r3", :green, :center, delta=-delta)
       point(x4, y4, "    丁円:r4,(x4,y4)", :orange, :left, :vcenter)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その414）

算額（その414）

愛媛県松山市桜谷町 伊佐爾波神社 嘉永3年(1850)

早坂平四郎: 算額の一考察，苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

直線の上に等円 3 個と大円が互いに接している。さらに，大円の中に4本の弦を引き分割された領域に，等円と同じ半径を持ち弦と接する 2 個の円と，中央の菱形に内接する中円を入れる（菱形の長い方の対角線は，上述の直線と並行である）。

大円の半径と中心座標を a, (0, a)
等円の半径と中心座標を r1, (x1, r1), ( x1 - r1, y1) ただし，y1 = r1(1 + √3)
菱形の中の円の半径と中心座標を r2, (0, a)
弦の1つの端点の座標を (bx, by)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, x1::negative, y1::positive,
     r2::positive, bx::positive;
a = 10
eq1 = x1^2 + (a - r1)^2 - (a + r1)^2
eq2 = (x1 - r1)^2 + (a - r1*(1+sqrt(Sym(3))))^2 - (a + r1)^2
eq3 = distance(-bx, a - sqrt(a^2 - bx^2), a, a, 0, a) - r2^2
eq4 = distance(-bx, a - sqrt(a^2 - bx^2), a, a, 0, r1) - r1^2;
res = solve([eq1, eq2, eq3, eq4], (r1, x1, bx, r2))

   1-element Vector{NTuple{4, Sym}}:
    (210 - 120*sqrt(3), 60 - 40*sqrt(3), -5891176044722613657600*sqrt(-154227 + 89043*sqrt(3))/52283326179841 - 3401272075258036224000*sqrt(3)*sqrt(-154227 + 89043*sqrt(3))/52283326179841 - 18210720*sqrt(3)/7230721 + 47757270/7230721 + 5891187146667514675200*sqrt(3)*sqrt(-51409 + 29681*sqrt(3))/52283326179841 + 10203835454942977476480*sqrt(-51409 + 29681*sqrt(3))/52283326179841, 10*sqrt(-185511737623 - 132968544*sqrt(-51409 + 29681*sqrt(3)) + 76769280*sqrt(3)*sqrt(-51409 + 29681*sqrt(3)) + 107105251656*sqrt(3))/sqrt(3836935845121 - 2215255943040*sqrt(3)))

   r1 = 2.1539;  x1 = -9.28203
   bx = 5.48676;  by = 1.63965;  r2 = 4.75039

菱形の中の円の半径は 4.75038687805787 である。

res[1][4] |> simplify

res[1][4].evalf()

   4.75038687805787

術では以下の式が示されているが，計算される値は不適切なもののようだ。

47/(16*sqrt(57√3 - 73) - 141(7 - 4√3)) * 10

   6.616796128691446

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   (r1, x1, bx, r2) = res[1]
   by = a - sqrt(a^2 - bx^2)
   @printf("r1 = %g;  x1 = %g\n", r1, x1)
   @printf("bx = %g;  by = %g;  r2 = %g\n", bx, by, r2)
   plot()
   circle(0, a, a, :blue)
   circle(x1, r1, r1)
   circle(x1 - 2r1, r1, r1)
   circle(x1 - r1, (1 + √3)r1, r1)
   circle(0, r1, r1)
   circle(0, 2a - r1, r1)
   circle(0, a, r2, :green)
   segment(a, a, -bx, by)
   segment(-a, a, bx, by)
   segment(a, a, -bx, 2a-by)
   segment(-a, a, bx, 2a-by)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(x1, r1, "(x1,r1)", :red, :center, delta=-delta)
       point(x1 - r1, r1*(1 + √3), "(x1-r1,(1+√3)r1)", :red, :center, delta=-delta)
       point(0, a, " a", :green)
       point(float(bx), float(by), " (bx,by)", :black)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その413）

算額（その413）

石川県七尾市中島町宮前（旧熊木村） 久麻加夫都阿良加志比古神社 文政6年(1823)

早坂平四郎: 算額の一考察，苫小牧工業専門学校紀要
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

半径 a の外円の中に，交わる 2 個の大円と，半径の等しい 4 個の小円が入っている。a を知って，大円と小円の半径を求めよ。

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

include("julia-source.txt");

using SymPy

@syms a::positive, r1::positive, r2::positive;
eq1 = a - 2r1 + 2r2
eq2 = (a - r1)^2 + (a - r2)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, r2))

   1-element Vector{Tuple{Sym, Sym}}:
    (a*(-5/4 + sqrt(33)/4) + a/2, a*(-5/4 + sqrt(33)/4))

res[1][1] |> simplify |> println
res[1][2] |> simplify |> println

   a*(-3 + sqrt(33))/4
   a*(-5 + sqrt(33))/4

大円の半径は (√33 - 3)a/4，小円の半径は (√33 - 5)a/4 である。

a = 10 のとき，6.8614066163450715 と 1.8614066163450715。

a = 10
a*(-3 + sqrt(33))/4, a*(-5 + sqrt(33))/4

   (6.8614066163450715, 1.8614066163450715)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 10
   r1 = a*(-3 + sqrt(33))/4
   r2 = a*(-5 + sqrt(33))/4
   @printf("r1 = %g;  r2 = %g\n", r1, r2)
   plot()
   circle(0, 0, a, :black)
   circle(a - r1, 0, r1)
   circle(r1 - a, 0, r1)
   circle(0, a - r2, r2, :blue)
   circle(0, r2 - a, r2, :blue)
   circle(r2, 0, r2, :blue)
   circle(-r2, 0, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, a - r2, " a-r2", :blue)
       point(r2, 0, "r2", :blue, delta=-delta)
       point(a - r1, 0, "a-r1", :red, :right, :bottom, delta=delta)
       point(a, 0, "a ", :black, :right, :bottom, delta=delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その412）

算額（その412）

百十五 群馬県富岡市一ノ宮 貫前神社 明治4年(1871)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

早坂平四郎: 算額の一考察
https://www.tomakomai-ct.ac.jp/wp01/wp-content/uploads/2014/06/kiyou5-8.pdf

キーワード：円6個

図のように，相接する等円 3 個はその中心が直線上にある。中央の等円の上に大円が乗り，大円と等円に外接する小円が 2 個ある。
大円と小円の半径がわかっているとき，等円の半径を求める式を作れ。

大円の半径と中心座標を r2, (0, r1 + r2)
小円の半径と中心座標を r3, (x3, y3) ただし x3 = r1
等円の半径と中心座標を r1, (0, 0), (2r1, 0)
とおき以下の連立方程式を解く。

「大円と小円の半径がわかっているとき，等円の半径を求める式」は，sympy では簡略化できない長い式なので，大円と小円の半径を与えて等円の半径を求める。」

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive,
      r3::positive, x3::positive, y3::positive;

(r2, r3) = (5, 1)
x3 = r1
eq1 = (2r1 - x3)^2 + y3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - y3)^2 - (r2 + r3)^2
res = solve([eq1, eq2], (r1, y3))

   3-element Vector{Tuple{Sym, Sym}}:
    (((-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))*(2 + (-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))/2, 1 + (-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3) - 10/(3*(-1/2 - sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)))
    ((-10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))*(2 - 10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))/2, 1 - 10/(3*(-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3)) + (-1/2 + sqrt(3)*I/2)*(10 + 10*sqrt(111)/9)^(1/3))
    ((-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + (10 + 10*sqrt(111)/9)^(1/3))*(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 2 + (10 + 10*sqrt(111)/9)^(1/3))/2, -10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 1 + (10 + 10*sqrt(111)/9)^(1/3))

3 番目の組が適解である。
大円，小円の半径を 5 寸, 1 寸としたとき，
等円の半径は 2.865876288528763 寸である。

(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + (10 + 10*sqrt(111)/9)^(1/3))*(-10/(3*(10 + 10*sqrt(111)/9)^(1/3)) + 2 + (10 + 10*sqrt(111)/9)^(1/3))/2

   2.865876288528763

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, y3) = res[3]
   r2 = 5
   r3 = 1
   x3 = r1
   @printf("r1 = %g;  x3 = %g;  y3 = %g\n", r1, x3, y3)
   plot()
   circle(0, 0, r1)
   circle(2r1, 0, r1)
   circle(-2r1, 0, r1)
   circle(x3, y3, r3, :blue)
   circle(-x3, y3, r3, :blue)
   circle(0, r1 + r2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1 + r2, " r1+r2  大円")
       point(r1, 0, " r1", :red, delta=-delta)
       point(2r1, 0, " 2r1 等円", :red, delta=-delta)
       point(x3, y3, "小円:r3,(x3,y3)", :blue, :left, :bottom, delta=2delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その411）

算額（その411）

会田安明(1811):「算法天生法指南」
永井信一:「アルベロスに関連した問題」
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/arbe.pdf

長方形内に斜線を引き，2 区分された領域に大円，小円を入れる。大円の直径が 2 寸，長方形の長辺が 3 寸のとき，小円の直径はいかほどか。

長方形の右下を原点とし，長辺と短辺を a, b，斜線と長辺の交点座標を y とする。
大円の半径と中心座標を　r1, (r1, a - r1)
小円の半径と中心座標を　r2, (r2, r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, y::positive,
     r1::positive, r2::positive;

r1 = 1
a = 3
b = 2r1
eq1 = y + b - sqrt(y^2 + b^2) - 2r2
eq2 = distance(0, y, b, 0, r1, a - r1) - r1^2
res = solve([eq1, eq2], (r2, y))

   1-element Vector{Tuple{Sym, Sym}}:
    (1/2, 3/2)

小円の直径は 1 寸である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   a = 3
   b = 2r1
   (r2, y) = (1/2, 3/2)
   @printf("r2 = %g;  y = %g\n", r2, y)
   @printf("小円の直径 = %g\n", 2r2)
   plot([0, b, b, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   circle(r1, a - r1, r1)
   circle(r2, r2, r2, :blue)
   segment(0, y, b, 0, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, a, " a", :black, :left, :bottom)
       point(b, 0, " b", :black, :left, :bottom)
       point(0, y, " y", :green, :left, :bottom)
       point(r1, a - r1, "大円:r1,(r1,a-r1)", :red, :center, delta=-delta)
       point(r2, r2, "小円:r2,(r2,r2)", :blue, :center, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その410）

算額（その410）

会田安明(1811):「算法天生法指南」

永井信一:「アルベロスに関連した問題」
http://www2.ttcn.ne.jp/~nagai/waseda/wasan/arbe.pdf

円内に4個の円を入れる。甲円，乙円の直径がそれぞれ 14 寸，7 寸のとき，丙円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive;

eq1 = r1 + r2 - r0
eq2 = x3^2 + (y3 - r1 + r0)^2 - (r1 + r3)^2
eq3 = x3^2 + (r0 - r2 - y3)^2 - (r2 + r3)^2
eq4 = x3^2 + y3^2 - (r0 -r3)^2
res = solve([eq1, eq2, eq3, eq4], (r0, r3, x3, y3))

   1-element Vector{NTuple{4, Sym}}:
    (r1 + r2, r1*r2*(r1 + r2)/(r1^2 + r1*r2 + r2^2), 2*r1*r2*(r1 + r2)/(r1^2 + r1*r2 + r2^2), (r1 - r2)*(r1 + r2)^2/(r1^2 + r1*r2 + r2^2))

r0 = 10.5;  r3 = 3;  x3 = 6;  y3 = 4.5
丙円の直径 = 6

甲円，乙円の直径がそれぞれ 14 寸，7 寸のとき，丙円の直径は 6 寸である。

(r1, r2) = (14, 7) ./ 2
2r1*r2*(r1 + r2)/(r1^2 + r1*r2 + r2^2)

   6.0

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (14, 7) .// 2
   (r0, r3, x3, y3) = (r1 + r2, r1*r2*(r1 + r2)/(r1^2 + r1*r2 + r2^2), 2*r1*r2*(r1 + r2)/(r1^2 + r1*r2 + r2^2), (r1 - r2)*(r1 + r2)^2/(r1^2 + r1*r2 + r2^2))
   @printf("r0 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r0, r3, x3, y3)
   @printf("丙円の直径 = %g\n", 2r3)
   plot()
   circle(0, 0, r0, :black)
   circle(0, r1 - r0, r1)
   circle(0, r0 - r2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(-x3, y3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
       point(0, r1 - r0, " r1-r0\n 甲円:r1", :red, :left, :vcenter)
       point(0, r0 - r2, " r0-r2\n 乙円:r2", :blue, :left, :vcenter)
       point(x3, y3, " 丙円,r3,(x3,y3)", :green, :center, delta=-delta)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;