算額（その379）

算額（その379）

二八 都川村西平 坂東九番観音堂（慈光寺境内観音堂） 文政13年(1830)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県都幾川町西平 慈光寺 文政13年(1830)
やまぶき　和算と歴史随想
https://yamabukiwasan.sakura.ne.jp/page3.html

キーワード：円5個，円弧2個

二本の等弧の間に，天円 2 個，人円 1 個，地円 2 個が互いに接している。天円の径は 68 寸(注），地円の径が 17 寸のとき，人円の径はいかほどか。

注：「六寸八寸」と描かれているが「六十八寸」の誤記である。

その一部が「弧」になる円の半径と中心座標を r1, (r1, 0)
天円の半径と中心座標を r2, (r2, y2)
人円の半径と中心座標を r3, (0, y3)
地円の半径と中心座標を r4, (r4, y4)
とし，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, y2::positive,  r3::positive,
     y3::positive, r4::positive, y4::positive;

(r2, r4) = (68//2, 17//2)
eq1 = (r1 - r2)^2 + y2^2 - (r1 + r2)^2
eq2 = r1^2 + y3^2 - (r1 + r3)^2
eq3 = (r1 - r4)^2 + y4^2 - (r1 + r4)^2
eq4 = r2^2 + (y2 - y3)^2 - (r2 + r3)^2
eq5 = r4^2 + (y3 - y4)^2 - (r3 + r4)^2

res = solve([eq1, eq2, eq3, eq4, eq5], (r1, y2, r3, y3, y4))

   1-element Vector{NTuple{5, Sym}}:
    (272, 136*sqrt(2), 32, 96*sqrt(2), 68*sqrt(2))

   r1 = 272;  y2 = 192.333;  r3 = 32;  y3 = 135.765;  y4 = 96.1665
   人円の直径 = 2r3 = 64 寸

人円の半径は 32 寸である（直径は 64 寸）

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = (68//2, 17//2)
   (r1, y2, r3, y3, y4) = res[1]
   @printf("r1 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  y4 = %g\n", r1, y2, r3, y3, y4)
   @printf("人円の直径 = 2r3 = %g 寸\n", 2r3)
   plot()
   circle(r1, 0, r1)
   circle(-r1, 0, r1)
   circle(r2, y2, r2, :blue)
   circle(-r2, y2, r2, :blue)
   circle(0, y3, r3, :green)
   circle(r4, y4, r4, :magenta)
   circle(-r4, y4, r4, :magenta)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(r1, 0, " r1", :red)
       point(r2, y2, "  天円:r2,(r2,y2)", :blue, :left, :vcenter)
       point(0, y3, "  人円:r3,(0,y2)", :green, :left, :vcenter)
       point(r4, y4, "  地円:r4,(r4,y4)", :magenta, :left, :vcenter)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       plot!(xlims=(-100, 350), ylims=(-100, 275))
   else
       plot!(showaxis=false)
   end
end;

算額（その378）

算額（その378）

三六 川越市古谷本郷 古尾谷八幡神社 天保12年(1841)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県川越市 古尾谷八幡神社 天保12年
http://tamtom.blog44.fc2.com/blog-entry-406.html

算額[古尾谷八幡神社]
http://blog.livedoor.jp/koike631/archives/50527866.html

等脚台形内に，乾円，坤円が入っている。等脚台形の上底の長さが 9 寸，乾円と坤円の直径がそれぞれ 7寸，9 寸のとき下底の長さを求めよ。

上底と下底の長さをそれぞれ 2w1，2w2，高さを h
乾円の半径と中心座標を r1, (x1, h - r1)
坤円の半径と中心座標を r2, (x2, r2)
とする。

r1 = 7/2, r2 = 4/2, w1 = 9/2 として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms h::positive, w1::positive, w2::positive,  r1::positive,
     r2::positive, x1::positive, x2::positive;

eq1 = distance(-w2, 0, -w1, h, x1, r1) - r1^2
eq2 = distance(w2, 0, -w1, h, x2, h - r2) - r2^2
eq3 = distance(w2, 0, -w1, h, x1, r1) - r1^2
eq4 = distance(w2, 0, w1, h, x2, h - r2) - r2^2;

# solve([eq1, eq2, eq3, eq4], (h, w2, x1, x2))

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (h, w2, x1, x2) = u
   return [
       -r1^2 + (x1 - (-h*(h*w2 + h*x1 + r1*w1 - r1*w2) + x1*(h^2 + w1^2 - 2*w1*w2 + w2^2))/(h^2 + w1^2 - 2*w1*w2 + w2^2))^2 + (-h*(h*r1 - w1*w2 - w1*x1 + w2^2 + w2*x1)/(h^2 + w1^2 - 2*w1*w2 + w2^2) + r1)^2,  # eq1
       -r2^2 + (x2 - (-h^2*w1 + h*r2*w1 + h*r2*w2 + w1^2*x2 + 2*w1*w2*x2 + w2^2*x2)/(h^2 + w1^2 + 2*w1*w2 + w2^2))^2 + (h - h*(h^2 - h*r2 + w1*w2 - w1*x2 + w2^2 - w2*x2)/(h^2 + w1^2 + 2*w1*w2 + w2^2) - r2)^2,  # eq2
       -r1^2 + (x1 - (h^2*w2 - h*r1*w1 - h*r1*w2 + w1^2*x1 + 2*w1*w2*x1 + w2^2*x1)/(h^2 + w1^2 + 2*w1*w2 + w2^2))^2 + (-h*(h*r1 + w1*w2 - w1*x1 + w2^2 - w2*x1)/(h^2 + w1^2 + 2*w1*w2 + w2^2) + r1)^2,  # eq3
       -r2^2 + (x2 - (h*(h*w1 - h*x2 - r2*w1 + r2*w2) + x2*(h^2 + w1^2 - 2*w1*w2 + w2^2))/(h^2 + w1^2 - 2*w1*w2 + w2^2))^2 + (h - h*(h^2 - h*r2 - w1*w2 + w1*x2 + w2^2 - w2*x2)/(h^2 + w1^2 - 2*w1*w2 + w2^2) - r2)^2,  # eq4
   ]
   end;

(r1, r2, w1) = [7, 4, 9] .// 2
iniv = [big"8.1", 10, -3, 3]
res = nls(H, ini=iniv);
println(res);

   (BigFloat[8.000000000000000289915981508997928978680584913763671176144672868231933897568114, 10.49999999999999484138720074297806292880734825561173528954255155770257311298397, -3.499999999999999683399841807629998190876460019780784537186673263801523533784222, 3.499999999999996870048545983071833654565390837799776925526657173046905812471916], true)

h = 8;  w2 = 10.5;  x1 = -3.5;  x2 = 3.5
下底の長さ = 2w2 = 20.99999999999999

算額の図は実際の数値に基づいていない。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (h, w2, x1, x2) = [9.0, 7.5, 2, -2.8]
   (h, w2, x1, x2) = res[1]
   @printf("h = %g;  w2 = %g;  x1 = %g;  x2 = %g\n", h, w2, x1, x2)
   println("下底の長さ = $(Float64(2w2))")
   plot([-w2, w2, w1, -w1, -w2], [0, 0, h, h, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :red)
   circle(x2, h - r2, r2, :green)
   segment(-w1, h, w2, 0)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1, "乾円:r1(x1,r1)", :red, :center, :top, delta=-delta)
       point(x2, h - r2, "坤円:r2\n(x2,h-r2)", :green, :center, :bottom, delta=delta)
       point(w1, h, "(w1,h)", :black, :right, :bottom, delta=delta)
       point(w2, 0, "w2", :black, :left, :bottom, delta=delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その377）

算額（その377）

三六 川越市古谷本郷 古尾谷八幡神社 天保12年(1841)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県川越市 古尾谷八幡神社 天保12年
http://tamtom.blog44.fc2.com/blog-entry-406.html

長方形内に，木，火，土，金，水の 5 円が入っている。木円の直径が 17 寸のとき，水円の直径を求めよ。

術に曰く

8 の開平方に 3 を加えたものを「極」と名付ける。
5 から 24 の開平方を引いたものに木円の径(17)を掛け，更に「極」を掛けることで水円の径を得る。

極 = √8 + 3
(5 - √24)*17*極 |> println
(5 - √24)*17*(√8 + 3) |> println

   10.009442010174723
   10.009442010174723

木円の半径と中心座標を r1, (x1, r1)
火円の半径と中心座標を r2, (x2, 2r1 - r2)
土円の半径と中心座標を r3, (x3, y3)
金円の半径と中心座標を r4, (r4, 2r1 - r4)
水円の半径と中心座標を r5, (r5, r5)
とする。

r1 = 17 として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms x1::positive, r1::positive, r2::positive,  x2::positive,
     r3::positive, x3::positive, y3::positive, r4::positive, r5::positive;

eq1 = (x1 - x2)^2 + (2r1 - r2 - r1)^2 - (r1 + r2)^2
eq2 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq3 = (x1 - r5)^2 + (r1 - r5)^2 - (r1 + r5)^2
eq4 = (x2 - x3)^2 + (2r1 - r2 - y3)^2 - (r2 + r3)^2
eq5 = (x2 - r4)^2 + (r4 - r2)^2 - (r2 + r4)^2
eq6 = (x3 - r4)^2 + (2r1 - r4 - y3)^2 - (r3 + r4)^2
eq7 = (x3 - r5)^2 + (y3 - r5)^2 - (r3 + r5)^2
eq8 = (r5 - r4)^2 + (2r1 -r4 - r5)^2 - (r4 + r5)^2;

# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (x1, r2, x2, r3, x3, y3, r4, r5) = u
   return [
       (r1 - r2)^2 - (r1 + r2)^2 + (x1 - x2)^2,  # eq1
       (-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq2
       (r1 - r5)^2 - (r1 + r5)^2 + (-r5 + x1)^2,  # eq3
       -(r2 + r3)^2 + (x2 - x3)^2 + (2*r1 - r2 - y3)^2,  # eq4
       (-r2 + r4)^2 - (r2 + r4)^2 + (-r4 + x2)^2,  # eq5
       -(r3 + r4)^2 + (-r4 + x3)^2 + (2*r1 - r4 - y3)^2,  # eq6
       -(r3 + r5)^2 + (-r5 + x3)^2 + (-r5 + y3)^2,  # eq7
       (-r4 + r5)^2 - (r4 + r5)^2 + (2*r1 - r4 - r5)^2,  # eq8
   ]
end;

r1 = 17/2
iniv = [big"2.1", 0.25, 1, 0.2, 0.96, 1.3, 0.4, 0.6] .* r1
res = nls(H, ini=iniv);
println(res);

   (BigFloat[18.04927980072966600990743424177031286841046841344571292514916108734929792090129, 2.277568135664543005016640483057481026967779743580825713557355507836459659630661, 9.249432267243960090044883726521819586682889562263977714737446982268376356193549, 1.688090125635917522526267428949100510249937740184158307062589102275226843795516, 8.150820151636668050007516015640830020184999772933719726083627782869396869278436, 10.91198610762122951754068288011827079712431459719790662797327282314006745314288, 3.556929041116716190186087773344970300658574809456580227043502433088470316206073, 5.004721005087340180215705712235570136210476406875180438641404697785532644015275], true)

x1 = 18.0493;  r2 = 2.27757;  x2 = 9.24943;
r3 = 1.68809;  x3 = 8.15082;  y3 = 10.912;
r4 = 3.55693;  r5 = 5.00472
木円の直径 = 17.0;  水円の直径 = 10.00944201017468

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 17/2
   (x1, r2, x2, r3, x3, y3, r4, r5) = res[1]
   b = 2r1
   @printf("x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  r5 = %g\n", x1, r2, x2, r3, x3, y3, r4, r5)
   println("木円の直径 = $(Float64(2r1));  水円の直径 = $(Float64(2r5))\n")
   plot([0, x1 + r1, x1 + r1, 0, 0], [0, 0, b, b, 0], color=:black, lw=0.5)
   circle(x1, r1, r1, :green)
   circle(x2, b - r2, r2, :red)
   circle(x3, y3, r3, :brown)
   circle(r4, b - r4, r4, :orange)
   circle(r5, r5, r5, :blue)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x1, r1, "木円:r1,(x1,r1)", :green, :center, delta=-delta/2)
       point(x2, b - r2, "火円:r2\n(x2,b-r2)", :red, :center, :bottom, delta=delta/2)
       point(x3, y3, "土円:r3\n(x3,y3)", :brown, :center, :vcenter)
       point(r4, b - r4, "金円:r4\n(r4,b-r4)", :orange, :center, delta=-delta/2)
       point(r5, r5, "水円:r5,(r5,r5)", :blue, :center, delta=-delta/2)
       point(x1 + r1, 2r1, "(x1+r1,2r1)", :black, :right, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その376）

算額（その376）

東京都 住吉神社 嘉永8年(1851)
http://www.wasan.jp/tokyo/sumiyosi.html

厥円（弓型）内に，互いに外接する 3 個の等円が入っている。また，それぞれは外円に内接している。
矢と弦それぞれの長さが与えられたとき，等円の直径を求めよ。

過去に，算額（その11） で掲載したが，それは「矢と弦それぞれの長さが与えられたとき」という条件を無視したものであった。また，この問題においては，区切られた斜線に内接しているのではなく，互いに外接しているようである。
https://blog.goo.ne.jp/r-de-r/e/0dac8f68b8bbc6af25eb9e60427a9c59
今回新たに解を求める。

厥円のもとである外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r) および (x, r + a)
とする。
直線と外円の交点座標を (x0, y0), y0 = sqrt(R^2 - x0^2)
とおき，以下の方程式を nlsolve() で解き，数値解を求める。

なお，問では「矢」と「弦」を与えるとあるので，方程式を解くに当たって，「矢」と「弦」から外円の径 R と切り取り位置 a（図参照）を求める式を求めておく。

using SymPy
@syms 矢::positive, 弦::positive, R::positive,  a::positive;

eq1 = R - a - 矢
eq2 = 2sqrt(R^2 - a^2) - 弦
solve([eq1, eq2], (R, a))

   1-element Vector{Tuple{Sym, Sym}}:
    ((弦^2 + 4*矢^2)/(8*矢), (弦 - 2*矢)*(弦 + 2*矢)/(8*矢))

include("julia-source.txt");

using SymPy
@syms 矢::positive, 弦::positive,
     R::positive,  a::positive, b::positive,
     x::positive, r::positibe;

eq1 = x^2 + (r + a)^2 - (R - r)^2
eq2 = x^2 + ((R - r) - (r + a))^2 - 4r^2

res = solve([eq1, eq2], (r, x))

   2-element Vector{Tuple{Sym, Sym}}:
    (R*(-R + a)/(-3*R + a), -sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))
    (R*(-R + a)/(-3*R + a), sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))

どちらも実質的に同じである。左右の円のどちらかを与えることになっている。

等円の半径は R*(-R + a)/(-3*R + a) であるが，これを 矢と弦を使って表そう。

eq = R*(-R + a)/(-3*R + a);

eq の R, a に以下を代入する。
R = (弦^2 + 4*矢^2)/(8*矢)
a = (弦 - 2*矢)*(弦 + 2*矢)/(8*矢)

eq(R => (弦^2 + 4*矢^2)/(8*矢), a => (弦 - 2*矢)*(弦 + 2*矢)/(8*矢)) |> simplify |> println

   矢*(弦^2 + 4*矢^2)/(2*(弦^2 + 8*矢^2))

等円の直径はこれを 2 倍して，矢*(弦^2 + 4*矢^2)/(弦^2 + 8*矢^2) である。

(矢, 弦) = (3, 8)
矢*(弦^2 + 4*矢^2)/(2*(弦^2 + 8*矢^2))

   2.2058823529411766

関数定義
等円直径(矢, 弦) = 矢*(弦^2 + 4*矢^2)/(弦^2 + 8*矢^2);

等円直径(3, 8)

   2.2058823529411766

等円直径(2, 10)

   1.7575757575757576

矢 = 3;  弦 = 8 のとき

   (25//6, 7//6)
   矢 = R - a = 3;  弦 = 2b = 8;  R = 4.16667;  a = 1.16667;  b = 4
   r = 1.10294;  x = 2.05798

using Plots

function draw(矢, 弦, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
   println((R, a))
   b = sqrt(R^2 - a^2)
   θ = round(Int64, atand(a, b))
   (r, x) = (R*(-R + a)/(-3*R + a), -sqrt(-1/(-3*R + a))*(-R + a)*sqrt(R + a))
   @printf("矢 = R - a = %g;  弦 = 2b = %g;  R = %g;  a = %g;  b = %g\n", 矢, 弦, R, a, b)
   @printf("r = %g;  x = %g\n", r, x)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180 - θ)
   circle(0, R - r, r, :blue)
   circle(x, r + a, r, :blue)
   circle(-x, r + a, r, :blue)
   segment(-b, a, b, a, :black)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, a, " a", :black, delta=-delta/2)
       point(b, a, "(b,a)", :black, :right, delta=-delta/2)
       point(0, R - r, " R - r", :blue)
       point(0, R, " R", :black, :left, :bottom)
       point(x, r + a, "(x,r+a)", :blue, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

矢 = 2，弦 = 10 のとき

   (29//4, 21//4)
   矢 = R - a = 2;  弦 = 2b = 10;  R = 7.25;  a = 5.25;  b = 5
   r = 0.878788;  x = 1.74078

算額（その375）

算額（その375）

山形県山形市山寺 立石寺根本中堂
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

厥円（弓型）内の 3 個の等円が 2 本の直線で隔てられている。
矢（R - a）と弦（2b）それぞれの長さが与えられたとき，等円の直径を求めよ。

過去に，算額（その11）で掲載したが，それは「矢と弦それぞれの長さが与えられたとき」という条件を無視したものであった。
https://blog.goo.ne.jp/r-de-r/e/0dac8f68b8bbc6af25eb9e60427a9c59
今回新たに解を求める。

厥円のもとである外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (0, R - r) および (x, r + a)
とする。
直線と外円の交点座標を (x0, y0), y0 = sqrt(R^2 - x0^2)
とおき，以下の方程式を nlsolve() で解き，数値解を求める。

なお，問では「矢」と「弦」を与えるとあるので，方程式を解くに当たって，「矢」と「弦」から外円の径 R と切り取り位置 a（図参照）を求める式を求めておく。

using SymPy
@syms 矢::positive, 弦::positive, R::positive,  a::positive;

eq1 = R - a - 矢
eq2 = 2sqrt(R^2 - a^2) - 弦
solve([eq1, eq2], (R, a))

   1-element Vector{Tuple{Sym, Sym}}:
    ((弦^2 + 4*矢^2)/(8*矢), (弦 - 2*矢)*(弦 + 2*矢)/(8*矢))

include("julia-source.txt");

using SymPy
@syms 矢::positive, 弦::positive,
     R::positive,  a::positive, b::positive, y0::positive,
     x::positive, r::positibe, x0::positive;

y0 = sqrt(R^2 - x0^2)
eq1 = x^2 + (r + a)^2 - (R - r)^2
eq2 = distance(0, a, x0, y0, x, r + a) - r^2
eq3 = distance(0, a, x0, y0, 0, R - r) - r^2;

# solve([eq1, eq2, eq3], (r, x, x0))

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r, x, x0) = u
   return [
       x^2 - (R - r)^2 + (a + r)^2,  # eq1
       -r^2 + (x - x0*(-a*r + r*sqrt(R^2 - x0^2) + x*x0)/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2 + (a + r - (-x0*(a*x + r*x0 - x*sqrt(R^2 - x0^2)) + (a + r)*(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2,  # eq2
       -r^2 + x0^2*(a - sqrt(R^2 - x0^2))^2*(-R + a + r)^2/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2))^2 + (R - r - (x0^2*(-R + a + r) + (R - r)*(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))/(R^2 + a^2 - 2*a*sqrt(R^2 - x0^2)))^2,  # eq3        
   ]
end;

(矢, 弦) = (3, 8)
(R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
iniv = [big"1.0", 2, 2]
res = nls(H, ini=iniv)
println(Float64.(res[1]), ", converge = ",res[2])

   [1.0909090909090908, 2.088931871468374, 1.8031774011398194], converge = true

   矢 = R - a = 3;  弦 = 2b = 8;  R = 4.16667;  a = 1.16667;  b = 4
   r = 1.09091;  x = 2.08893;  x0 = 1.80318;  y0 = 3.75628

(矢, 弦) = (3, 8) のとき，等円の半径は 12/11 ≒ 1.090909 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (矢, 弦) = (3, 8)
   (R, a) = (弦^2//(8*矢) + 矢//2, 弦^2//(8*矢) - 矢//2)
   println((R, a))
   b = sqrt(R^2 - a^2)
   θ = round(Int64, atand(a, b))
   (r, x, x0) = res[1]
   y0 = sqrt(R^2 - x0^2)
   @printf("矢 = R - a = %g;  弦 = 2b = %g;  R = %g;  a = %g;  b = %g\n", 矢, 弦, R, a, b)
   @printf("r = %g;  x = %g;  x0 = %g;  y0 = %g\n", r, x, x0, y0)
   plot()
   circle(0, 0, R, :black, beginangle=θ, endangle=180 - θ)
   segment(-b, a, b, a, :black)
   circle(0, R - r, r, :blue)
   circle(x, r + a, r, :blue)
   circle(-x, r + a, r, :blue)
   segment(0, a, x0, y0, :black)
   segment(0, a, -x0, y0, :black)
   if more == true
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(x0, y0, " (x0,y0)", :black, :left, :bottom)
       point(0, a, " a", :black, delta=-delta/2)
       point(b, a, "(b,a)", :black, :right, delta=-delta/2)
       point(0, R - r, " R - r", :blue)
       point(0, R, " R", :black, :left, :bottom)
       point(x, r + a, "(x,r+a)", :blue, :center, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その374）

算額（その374）

長野県埴科郡坂城町 諏訪社 文化2年(1805)
中村信弥「改訂増補長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額1（78）

大小の半円が交わる中に，甲円，乙円を入れる。小円の径が 27 寸，乙円の径が 9寸のとき，甲円の径を求めよ。

算額（その373）に似ているが，この問題ではヒントとなる (0, 2r1)-(2r1, 0) の線分が描かれていないので，条件が足りないなあと思って先送りしていた。

追記: 2024/03/16
乙円の中心座標 (x4, y4) については特に条件を設定しなかったのであるが，暗黙の了解で（?）「乙円の中心は，大円と小円の中心を結ぶ線分上にある」とするようである。そのような条件にすると，中村が指摘した「二百五十分之七十の誤写」のとおりになる。

大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, r2);  r2 = 27/2
甲円の半径と中心座標を r3, (x3, r3)
乙円の半径と中心座標を r4, (x4, y4); r4 = 9/2
大円と小円の交点座標を (x0, y0)
として，以下の連立方程式解く。

include("julia-source.txt");

using SymPy
@syms r1::positive,  r2::positive, r3::positive, r4::positive,
   x3::positive, x4::positive, y4::positibe, x0::positive, y0::positive;

eq1 = (x3 - r1)^2 + r3^2 - (r1 - r3)^2  # 甲円が大円に内接する
eq2 = (r1 - x4)^2 + y4^2 - (r1 - r4)^2  # 乙円が大円に内接する
eq3 = x4^2 + (r2 - y4)^2 - (r2 - r4)^2  # 乙円が小円に内接する
eq4 = x3^2 + (r2 - r3)^2 - (r2 + r3)^2  # 甲円が小円に外接する
eq5 = (2r2 - y0)*r1 - r2*x0  # (x0, y0) は大円と小円の交点
eq5 = y4/(r1 - x4) - r2/r1  # こちらの条件を使う
eq6 = (x0 - r1)^2 + y0^2 - r1^2;
eq7 = x0^2 + (y0 - r2)^2 - r2^2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, r3, x3, x4, y4, x0, y0))

   1-element Vector{NTuple{7, Sym{PyCall.PyObject}}}:
    (2*r4*(r2 - r4)/(r2 - 2*r4), 4*r2*r4^2*(r2 - r4)^2/(-r2^2 + r2*r4 + r4^2)^2, 4*r2*r4*(r2 - r4)/(r2^2 - r2*r4 - r4^2), 2*r4*(r2 - r4)^2/(r2^2 - 2*r2*r4 + 2*r4^2), r2^2*r4/(r2^2 - 2*r2*r4 + 2*r4^2), 4*r2^2*r4*(r2 - 2*r4)*(r2 - r4)/(r2^2 - 2*r2*r4 + 2*r4^2)^2, 8*r2*r4^2*(r2 - r4)^2/(r2^2 - 2*r2*r4 + 2*r4^2)^2)

2res[1][2] |> println
2res[1][2](r2 => 27//2, r4 => 9//2) |> println
2res[1][2](r2 => 27//2, r4 => 9//2).evalf() |> println

   8*r2*r4^2*(r2 - r4)^2/(-r2^2 + r2*r4 + r4^2)^2
   432/25
   17.2800000000000

小円の径が 27 寸，乙円の径が 9 寸のとき，甲円の径は 432/25 =  17.28 寸 = 17 寸 7/25 分である。

その他のパラメータは以下のとおりである。

r1 = 18;  r3 = 8.64;  x3 = 21.6;  x4 = 7.2;  y4 = 8.1;  x0 = 12.96;  y0 = 17.28

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r4) = (27, 9) .// 2
   (r1, r3, x3, x4, y4, x0, y0) = (
       2*r4*(r2 - r4)/(r2 - 2*r4),
       4*r2*r4^2*(r2 - r4)^2/(-r2^2 + r2*r4 + r4^2)^2,
       4*r2*r4*(r2 - r4)/(r2^2 - r2*r4 - r4^2),
       2*r4*(r2 - r4)^2/(r2^2 - 2*r2*r4 + 2*r4^2),
       r2^2*r4/(r2^2 - 2*r2*r4 + 2*r4^2),
       4*r2^2*r4*(r2 - 2*r4)*(r2 - r4)/(r2^2 - 2*r2*r4 + 2*r4^2)^2,
       8*r2*r4^2*(r2 - r4)^2/(r2^2 - 2*r2*r4 + 2*r4^2)^2)
   @printf("小円の径が 27 寸，乙円の径が 9 寸のとき，甲円の直径は %.10g 寸\n", 2r3)
   @printf("r1 = %g;  r3 = %g;  x3 = %g;  x4 = %g;  y4 = %g;  x0 = %g;  y0 = %g\n", r1, r3, x3, x4, y4, x0, y0)
   plot()
   circle(r1, 0, r1, :black, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue, beginangle=270, endangle=450)
   circle(x3, r3, r3)
   circle(x4, y4, r4, :orange)
   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)
       segment(0, 2r2, 2r1, 0, :gray)
       segment(0, r2, r1, 0, :green)
       point(r1, 0, "大円:r1 ", :black, :right, :bottom)
       point(x3, r3, "甲円:r3\n(x3,r3)\n", :red, :center, :top, delta=-delta/2)
       point(x4, y4, "乙円:r4=$r4\n(x4,y4)", :orange, :center, :bottom, delta=delta/2)
       point(0, r2, " 小円\n r2=$r2", :blue, :left, :bottom)
       point(x0, y0, "(x0,y0)\n", :black, :left, :bottom)
   end
end;

算額（その373）

算額（その373）

山形県山形市旅籠町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

直角三角形の 2 つの辺それぞれを直径とする大半円と中半円があり，2 つの等円が入っている。大半円の直径が与えられたとき，中半円の直径を求めよ。

大半円の半径と中心座標を r1, (r1, 0)
中半円の半径と中心座標を r2, (0, r2)
等円の半径と中心座標を r3, (r3, y31) および (x3, y32)
大半円と中半円の交点座標を (x0, y0)
とする。

以下の連立方程式を nlsolve() で解き，数値解を求める。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, y31::positive,
     x3::positive, y32::positive, x0::positive, y0::positive;

eq1 = (r1 - r3)^2 + y31^2 - (r1 + r3)^2
eq2 = (r1 - x3)^2 + y32^2 - (r1 - r3)^2
eq3 = x3^2 + (r2 - y32)^2 - (r2 - r3)^2
eq4 = sqrt(x0^2 + y0^2) * sqrt(4r1^2 + 4r2^2) - 2r1*2r2
eq5 = (x0^2 + y0^2)/((2r1 - x0)^2 + y0^2) - r2^2/r1^2
eq6 = y0/(2r1 - x0) - r2/r1
eq7 = distance(0, 2r2, 2r1, 0, r3, y31) - r3^2;
# res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7])

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=1e-14)
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=1e-14)
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (r2, r3, x3, y31, y32, x0, y0) = u
   return [
       y31^2 + (r1 - r3)^2 - (r1 + r3)^2,  # eq1
       y32^2 - (r1 - r3)^2 + (r1 - x3)^2,  # eq2
       x3^2 - (r2 - r3)^2 + (r2 - y32)^2,  # eq3
       -4*r1*r2 + sqrt(4*r1^2 + 4*r2^2)*sqrt(x0^2 + y0^2),  # eq4
       (x0^2 + y0^2)/(y0^2 + (2*r1 - x0)^2) - r2^2/r1^2,  # eq5
       y0/(2*r1 - x0) - r2/r1,  # eq6
       -r3^2 + (-r1*(r1*r3 + 2*r2^2 - r2*y31)/(r1^2 + r2^2) + r3)^2 + (-r2*(2*r1^2 - r1*r3 + r2*y31)/(r1^2 + r2^2) + y31)^2,  # eq7
   ]
end;

r1 = 10
iniv = [big"0.735", 0.245, 0.430, 0.990, 0.495, 0.700, 0.955] .* r1
res = nls(H, ini=iniv);
println(res);
   (BigFloat[7.364773268277131843543802900503553345300188559250709059836081841354806617265332, 2.445735829117338727593748564700272182455934087389288548663153970404609644379193, 4.27028686542385572015136287764767660467637682250915731571476513356507761504113, 9.890876258688789657450970376693568285214500451690718312230118896022365345155555, 4.923138689794054227478994166323836685186339000292746456092958149436942045604101, 7.033185377446014888320211825529028530741347528474622758676476217579965993488618, 9.549764970689061667860842518627745663621779873444347024592820593189972532890017], true)

大円の半径が 10 のとき，中円の半径は 7.36477 あまりである。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, x3, y31, y32, x0, y0)  = res[1]
   @printf("r2 =　%g;  r3 =　%g;  x3 =　%g;  y31 =　%g;  y32 =　%g;  x0 =　%g;  y0 = %g\n",  r2, r3, x3, y31, y32, x0, y0)
   plot()
   circle(r1, 0, r1, beginangle=0, endangle=180)
   circle(0, r2, r2, :blue, beginangle=-90, endangle=90)
   circle(x3, y32, r3, :green)
   circle(r3, y31, r3, :green)
   segment(0, 2r2, 2r1, 0, :black)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(r1, 0, " r1", :red, :left, :bottom, delta=delta/2)
       point(0, r2, " r2", :blue, :left, :vcenter)
       point(r3, y31, "(r3,y31)", :green, :center, delta=-delta/2)
       point(x3, y32, "(x3,y32)", :green, :center, delta=-delta/2)
       point(x0, y0, "  (x0,y0)", :green, :left, :vcenter, delta=-delta/2)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その372）

算額（その372）

山形県山形市旅籠町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

円を折って弓形を作り，弓形の中に大円と小円4個を入れる。大円の直径が与えられたとき，小円の直径を求めよ。

外円の半径と中心座標を r0, (0, 2r1 - r0)
大円の半径と中心座標を r1, (0, r1)
等円の半径と中心座標を r2, (x21, r2) および (x22, y22)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x21::positive,
     x22::positive, y22::positive;

eq1 = x21^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x22^2 + (r1 - y22)^2 - (r1 + r2)^2
eq3 = x22^2 + (2r1 - r0 - r2)^2 - (r0 -r2)^2
eq4 = x21^2 + (r2 - r0 + 2r1)^2 - (r0 - r2)^2
eq5 = x22^2 + (y22 - r0 + 2r1)^2 - (r0 + r2)^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r2, x21, x22, y22))

   2-element Vector{NTuple{5, Sym}}:
    (3*r1/2, r1/4, r1, sqrt(6)*r1/2, 3*r1/4)
    (r1*(9 - sqrt(33))/2, r1*(7 - sqrt(33))/4, r1*sqrt(7 - sqrt(33)), r1*sqrt(-27 + 5*sqrt(33)), r1*(-13/4 + 3*sqrt(33)/4))

1 組目のものが適解である。

大円の半径 r1 と　等円の半径 r2 の比は r2/r1 = 1/4。したがって，大円の直径の 1/4 が等円の直径である。

res[1][2] / r1 |> println

   1/4

r1 = 1 としたとき，以下のようになる。
r0 = 1.5;  r1 = 1;  r2 = 0.25; , x21 = 1; , x22 = 1.22474; , y22 = 0.75

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1
   (r0, r2, x21, x22, y22) = (3*r1/2, r1/4, r1, sqrt(6)*r1/2, 3*r1/4)
   @printf("r0 = %g;  r1 = %g;  r2 = %g; , x21 = %g; , x22 = %g; , y22 = %g\n", r0, r1, r2, x21, x22, y22)
   x0 = sqrt(r0^2 - (2r1 - r0)^2)
   θ = round(Int64, atand(2r1 - r0, x0))
   plot()
   circle(0, r1, r1)
   #circle(0, 2r1 - r0, r0, :blue, beginangle=330, endangle=570)
   circle(0, 2r1 - r0, r0, :blue, beginangle=-θ, endangle=180 + θ)
   circle(0, 2r1 - r0, r0, :gray85, beginangle=180 + θ, endangle=360 - θ)
   circle(0, r0 - 2r1, r0, :blue, beginangle=θ, endangle=180 - θ)
   circle(0, r0 - 2r1, r0, :gray85, beginangle=180 - θ, endangle=360 + θ)
   circle(x21, r2, r2, :green)
   circle(-x21, r2, r2, :green)
   circle(x22, y22, r2, :green)
   circle(-x22, y22, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, 2r1, " 2r1", :red, :left, :bottom)
       point(0, r1, " r1", :red)
       point(0, -r1, " -r1", :gray)
       point(0, 2r1 - r0, " 2r1-r0", :blue)
       point(0, r0 - 2r1, " r0-2r1", :gray)
       point(x21, r2, "(x21,r2)", :green, :center, delta=-delta/3)
       point(x22, y22, "(x22,y22)", :green, :center, delta=-delta/3)
       point(x0, 0, "  x0", :blue, :left, :bottom)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その371）

算額（その371）

山形県山形市旅籠町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

外円の中に 3 個の円弧を描き，等円 3 個を入れる。外円の直径が与えられたとき，等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
上の等円の半径と中心座標を r, (0, R - r)
右下の等円の半径と中心座標を r, (√3(R-r)/2, (r-R)/2)
右下の円弧の半径と中心座標を R, (√3R/2, -R/2)
とする。

以下の連立方程式を解く。
eq1 と eq2 はr2 を含まないので先に解く。得られた解に基づいて eq2 を解く。

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, a::positive;

l = sqrt(a^2 + 4r1^2)
eq1 = r3/(a - (2r1 + r3)) - 2r1/l
eq2 = r1/(a - r1) - 2r1/l
eq3 = distance(0, 2r1, a, 0, 2r1 - r2, 2r1 - r2) - r2^2
res = solve([eq1, eq2], (r3, a))

   1-element Vector{Tuple{Sym, Sym}}:
    (r1/4, 8*r1/3)

@syms r1::positive, r2::positive, r3::positive, a::positive;

(r3, a) = (r1/4, 8*r1/3)
eq3 = distance(0, 2r1, a, 0, 2r1 - r2, 2r1 - r2) - r2^2
res = solve([eq3], (r2))

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

r2 = r1/2 が適解である。

r2 = r1/2
r3 = r1/4
r3/r2 |> println

   1/2

乙円の半径は甲円の半径の 1/2，丙円の半径は乙円の半径の 1/2

 r1 = 1;  r2 = 0.5;  r3 = 0.25;  a = 2.66667

using Plots

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

算額（その370）

算額（その370）

一〇一 大宮市高鼻町 氷川神社 明治31年(1898)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

山形県山形市旅籠町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

外円の中に 3 個の円弧を描き，等円 3 個を入れる。外円の直径が与えられたとき，等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
上の等円の半径と中心座標を r, (0, R - r)
右下の等円の半径と中心座標を r, (√3(R-r)/2, (r-R)/2)
右下の円弧の半径と中心座標を R, (√3R/2, -R/2)
とする。

以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive;

eq = (sqrt(Sym(3))R/2)^2 + (R-r + R//2)^2 - (R + r)^2
res = solve([eq], (r))

   Dict{Any, Any} with 1 entry:
     r => 2*R/5

等円の半径は外円の半径の 2/5 である。

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   R = 1
   r = (2/5)R
   @printf("R = %g;  r = %g\n", R, r)
   plot()
   circle(0, 0, R, :black)
   rotate(0, R - r, r)
   rotate(0, R, R, :blue, beginangle=210, endangle=330)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, R, " R", :red, :left, :bottom)
       point(0, R - r, " R-r", :red, :left, :vcenter)
       point(√3R/2, -R/2, "(√3R/2,-R/2)  ", :blue, :center, :bottom, delta=delta/1.5)
       point(√3(R-r)/2, (r-R)/2, "(√3(R-r)/2,(r-R)/2)", :red, :center, :bottom, delta=delta/1.5)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

 

算額（その369）

算額（その369）

番外編 2023 鳥取大学工学部等前期日程数学
https://blog.goo.ne.jp/mh0920-yh/e/8075ae29b935640983c8fcc18a38a222

⊿ABC において，∠A = 60°，AB = 8，AC = 6 とする。⊿ABC の垂心を H とするとき，ベクトルAH をベクトルAB，ベクトルAC を用いて表わせ。

include("julia-source.txt");

using SymPy

@syms xa::positive, ya::positive,
     xb::positive, yb::positive,
     xc::positive, yc::positive,
     xh::positive, yh::positive,
     x::positive, y::positive;

図のように記号を定める。
頂点 A, B, C の座標をそれぞれ (xa, ya), (xb, yb), (xc, yc) とする。
頂点 B を原点に，頂点C が x 軸上に来るように三角形を置く。すなわち，xb = yb = yc = 0 とする。
xa, ya, xc を決める。

xb = yb = yc = 0
eq1 = (xa - xb)^2 + ya^2 - 8^2  # ピタゴラスの定理
eq2 = (xc - xa)^2 + ya^2 - 6^2  # ピタゴラスの定理
eq3 = ((8^2 + 6^2 - (xc - xb)^2))/(2*8*6) - cos(PI/3)  # 第2余弦定理
res1 = solve([eq1, eq2, eq3], (xa, ya, xc))
(xa, ya, xc) = res1[1]

   (20*sqrt(13)/13, 12*sqrt(39)/13, 2*sqrt(13))

(xa, ya, xc) = (20*sqrt(13)/13, 12*sqrt(39)/13, 2*sqrt(13)) である。

次に，BH⊥AC，CH⊥AB であることから，垂心 H の座標 (xh, yh) を求める。

eq1 = y/(x - xb) * ya/(xc - xa) - 1
eq2 = y/(xc - x) * ya/(xa - xb) - 1
res2 = solve([eq1, eq2], (x, y))

   Dict{Any, Any} with 2 entries:
     y => 10*sqrt(39)/39
     x => 20*sqrt(13)/13

(x, y) = (20*sqrt(13)/13, 10*sqrt(39)/39) である。

次に，ベクトルAB の a 倍とベクトルAC の b 倍の和がベクトルAH になるので，以下の方程式が成り立つ。

@syms a, b
(x, y) = (res2[x], res2[y])
(xh, yh) = ((xa, ya) .- (xb, 0)).*a .+ ((xa, ya) .- (xc, 0)).*b .- ((xa, ya) .- (x, y))

   (20*sqrt(13)*a/13 - 6*sqrt(13)*b/13, 12*sqrt(39)*a/13 + 12*sqrt(39)*b/13 - 2*sqrt(39)/3)

連立方程式を解き，a, b を求める。

solve([xh, yh], (a, b))

   Dict{Any, Any} with 2 entries:
     b => 5/9
     a => 1/6

a = 1/6, b = 5/9 である。

   xa = 5.547;  ya = 5.76461;  xc = 7.2111;
   AB = 8;  BC = 6;  AC = 7.2111
   xh = 5.547;  yh = 1.60128
   a = 0.166667, b = 0.555556

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   xb = yb = yc = 0
   (xa, ya, xc) = (20*sqrt(13)/13, 12*sqrt(39)/13, 2*sqrt(13))
   (xh, yh) = (20*sqrt(13)/13, 10*sqrt(39)/39)
   (a, b) = (1/6, 5/9)
   @printf("xa = %g;  ya = %g;  xc = %g;\n", xa, ya, xc)
   @printf("AB = %g;  BC = %g;  AC = %g\n", sqrt((xa - xb)^2 + ya^2), sqrt((xc - xa)^2 + ya^2), xc - xb)
   @printf("xh = %g;  yh = %g\n", xh, yh)
   @printf("a = %g, b = %g\n", a, b)
   plot([0, xc, xa, 0], [0, 0, ya, 0], color=:black, lw=0.5)
   #segment(xa - (xa - xb)*a, ya - ya*a, xh, yh)
   plot!([xa, xa - (xa - xb)*a], [ya, ya - ya*a], arrow=0.6, color=:red, lw=0.5)
   plot!([xa - (xa - xb)*a, xh], [ya - ya*a, yh], arrow=0.6, color=:blue, lw=0.5)
   plot!([xa - (xa - xc)*b, xh], [ya - ya*b, yh], arrow=0.5, color=:red, lw=0.5)
   plot!([xa, xa - (xa - xc)*b], [ya, ya - ya*b], arrow=0.5, color=:blue, lw=0.5)
   plot!([xa, xh], [ya, yh], arrow=0.6, color=:green, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(xa, ya, " A")
       point(0, 0, "  B", :green, :left, :bottom)
       point(xc, 0, " C", :green, :left, :bottom)
       point(x, y, " H")
       point(xa - (xa - xb)*a, ya - ya*a, " a*AB")
       point(xa - (xa - xc)*b, ya - ya*b, " b*AC")
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その368）

算額（その368）

山形県山形市旅籠町 湯殿山神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

小円は外円（団扇）の弦と部分円，外円に接する。小円の直径の最大値を求めよ。

外円の半径と中心座標を r0, (0, r0)
部分円の半径と中心座標を r1, (0, 0)
小円の半径と中心座標を r2, (x2, r1 -r0 - r2)
とおき，以下の連立方程式を解く。
eq3, eq4 は外円と部分円の交点座標 (x, y) を求めるためのもので，必須ではない。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, x2::positive,
     x::positive, y::positive;

eq1 = x2^2 + (r0 - (r1 - r2))^2 - (r0 - r2)^2
eq2 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = x^2 + (y - r0)^2 - r0^2
eq4 = x^2 + y^2 - r1^2
res = solve([eq1, eq2, eq3, eq4], (r2, x2, x, y))

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

小円の半径 r2 は部分円の半径 r1 の関数になっている。0 < r1 < r0 において，上に凸の曲線である。

diff(res[1][1], r1) |> simplify |> println

   (-r1*(2*r0 - r1)/2 + (r0 - r1)*(2*r0 + r1))/(2*r0 + r1)^2

微分して，0 になるときの　r1 を求める。2*r0*(-1 + √2) である。

solve(diff(res[1][1], r1), r1) |> println

   Sym[2*r0*(-1 + sqrt(2))]

そのときの r2 を求める。r0*(3 - 2√2) である。

res[1][1](r1 => 2*r0*(-1 + sqrt(Sym(2)))) |> simplify |> println

   r0*(3 - 2*sqrt(2))

まとめると，外円の半径を r0 としたとき，r1 = 2r0*(-1 + √2) = 0.8284271247461903 のとき，r2 は最大値 r0*(3 - 2*√2) をとる。r0 = 1 のとき， r1 = 2*r0*(-1 + 2) のとき r2 = r0*(3 - 22) = 0.1715728752538097

r0 = 1
2r0*(-1 + √2)  # r1

   0.8284271247461903

r0*(3 - 2√2)  # r2

   0.1715728752538097

ついでに， x2, x, y も求めておく。

r0 = 1
r1 = r1 = 2r0*(-1 + √2)
(sqrt(2)*r1*sqrt(2*r0 - r1)/sqrt(2*r0 + r1)) |> println  # x2
r1*sqrt(2*r0 - r1)*sqrt(2*r0 + r1)/(2*r0)  |> println  # x
r1^2/(2*r0) |> println  # y

   0.7540175693734211
   0.7540175693734212
   0.34314575050762

   r0 = 1;  r1 = 0.828427;  r2 = 0.171573;  x2 = 0.754018  x = 0.754018;  y = 0.343146

using Plots
pyplot(size=(400, 250), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(-r1*(r1 - 2)/(2*(r1 + 2)), xlims=(0,2), ylims=(0, 0.18), xlabel="r1", ylabel="r2")
vline!([2√2 - 2])
point(2√2 - 2, 3 - 2√2, " (2√2-2, 3-2√2)", :green, :left, :bottom)

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (1, 2√2 - 2)
   (r2, x2, x, y) = (r1*(2*r0 - r1)/(2*(2*r0 + r1)), sqrt(2)*r1*sqrt(2*r0 - r1)/sqrt(2*r0 + r1), r1*sqrt(2*r0 - r1)*sqrt(2*r0 + r1)/(2*r0), r1^2/(2*r0))
   (r2, x2, x, y) = (r0*(3 - 2*sqrt(2)), sqrt(2)*r1*sqrt(2*r0 - r1)/sqrt(2*r0 + r1), r1*sqrt(2*r0 - r1)*sqrt(2*r0 + r1)/(2*r0), r1^2/(2*r0))
   @printf("r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g  x = %g;  y = %g\n", r0, r1, r2, x2, x, y)
   θ = round(Int64, atand(y, x))
   plot()
   circle(0, r0, r0)
   circle(0, 0, r1, :blue, beginangle=θ, endangle=180 - θ)
   circle(x2, r1 - r2, r2, :magenta)
   x00 = sqrt(r0^2 - (r0 - r1)^2)
   segment(-x00, r1, x00, r1)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       point(0, r0, " r0")
       point(0, r1, " r1")
       point(x2, r1 - r2, "(x2,r1-r2)")
       point(x, y, "(x,y)")
       hline!([0, r1], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その367）

算額（その367）

山形県山形市緑町 三島稲荷神社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

外円の中に等斜 2 個と弦，大円，中円，小円の 4 円がある。中円の直径が 4 寸であるとき，小円の直径を求めよ。

外円の下端を原点とする。
外円の半径と中心座標を r0, (0, 2r1 + r2 + r3)
大円の半径と中心座標を r1, (0, 2r0 - r1) および (0, 2r3 + r1)
中円の半径と中心座標を r2, (0, 2r0 - 2r1 - r2)
小円の半径と中心座標を r3, (r3, r3)
斜線と外円の交点座標を (a, r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1::positive, r2::positive, r3::positive, a::positive;

eq1 = r2/(2r1 + r2) - r1/(3r1 + 2r2)
eq2 = r2/(2r1 + r2) - a/sqrt((2r0 - 2r3)^2 + a^2)
eq3 = 2r1 + r2 + r3 - r0
eq4 = (2r1 + r2 - r3)^2 + a^2 - r0^2

res = solve([eq1, eq2, eq3, eq4], (r0, r1, r3, a))

   1-element Vector{NTuple{4, Sym}}:
    (r2*(4*sqrt(5) + 9)/4, r2*(1 + sqrt(5))/2, r2/4, r2*sqrt(2 + sqrt(5)))

r3 = r2/4。小円の直径は中円の直径の 1/4 である。

中円の直径が 4 寸であるとき，小円の直径は 1 寸である。

 r0 = 8.97214;  r1 = 3.23607;  r2 = 2;  r3 = 0.5;  a = 4.11634

using Plots

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

算額（その366）

算額（その366）

山形県山形市小白川町 天満社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

正方形内に二等辺三角形，大円，中円，小円が入っている。大円の直径を与えられたときに，中円の直径を求めよ。

正方形の左下隅の座標を (0, 0)
大円の半径と中心座標を r1, (0, r1)
中円の半径と中心座標を r2, (0, r2)
小円の半径と中心座標を r3, (r1 - r3, r3)
二等辺三角形の底辺の長さを 2b とする。

算額（その361）の簡単なものである。
https://blog.goo.ne.jp/r-de-r/e/ff37569366bf01f983d6cec055ecc257

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive, b::positive;

eq1 = 2(r1 - r3)^2 - (r1 + r3)^2
eq2 = distance(0, 2r1, b, 0, 0, r2) - r2^2
eq3 = distance(0, 2r1, b, 0, r1 - r3, r3) - r3^2

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

   3-element Vector{Tuple{Sym, Sym, Sym}}:
    (r1/2, r1*(3 - 2*sqrt(2)), sqrt(2)*r1/2)
    (-11*r1/4 - 19*r1*sqrt(203753 - 142008*sqrt(2))/1156 - 3*sqrt(2)*r1*sqrt(203753 - 142008*sqrt(2))/578 + 3*sqrt(2)*r1/2, r1*(3 - 2*sqrt(2)), r1*(-2 + 3*sqrt(2))/2)
    (-11*r1/4 + 3*sqrt(2)*r1*sqrt(203753 - 142008*sqrt(2))/578 + 19*r1*sqrt(203753 - 142008*sqrt(2))/1156 + 3*sqrt(2)*r1/2, r1*(3 - 2*sqrt(2)), r1*(-2 + 3*sqrt(2))/2)

1番目のものが適解である。
r2 = r1/2 で，中円の直径は大円の直径の 1/2 である。

 r2 = 0.5;  r3 = 0.171573;  b = 0.707107;  r2/r1 = 0.5

using Plots

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

算額（その365）

算額（その365）

山形県山形市小白川町 天満社
山形算額勝負-湯殿山神社を目指せ-
https://www.sci.yamagata-u.ac.jp/wasan/pdf/20180711SSEP.pdf

正方形内に大小の正三角形と甲円，乙円がある。甲円の直径が与えられたときに，乙円の直径を求めよ。

正方形の辺の長さを a とする。

甲円の半径と中心座標を r1, (x1, r1)
容円の半径と中心座標を r2, (a - r2, y2)
とする。

以下の連立方程式を解く。a = 1 としても乙円と甲円の比を求める場合には一般性を失わない。

include("julia-source.txt");

using SymPy

@syms r1::positive, x1::positive, r2::positive,
     y2::positive, a::positive;

a = 1
eq1 = distance(0, 0, a/2, sqrt(Sym(3))a/2, x1, r1) - r1^2
eq2 = distance(0, (sqrt(Sym(3)) - 1)a, (sqrt(Sym(3)) - 1)a, 0, x1, r1) - r1^2
eq3 = distance(a, 0, a/2, sqrt(Sym(3))a/2, a - r2, y2) - r2^2
eq4 = distance(a, a, (sqrt(Sym(3)) - 1)a, 0, a - r2, y2) - r2^2

res = solve([eq1, eq2, eq3, eq4], (r1, x1, r2, y2))

   4-element Vector{NTuple{4, Sym}}:
    ((-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3), -1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5), 1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3), -1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5) + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16), 1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + sqrt(3) + (19*sqrt(3)/3 + 11)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3, -1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5), 1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))
    ((-11 - 19*sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + sqrt(3) + (19*sqrt(3)/3 + 11)*(-1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3, -1/2 + sqrt(-10*sqrt(3) + sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2, (-22 - 10*sqrt(3))*(1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5) + (1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16), 1/2 + sqrt(-sqrt(3)/8 + sqrt(2 - sqrt(3))/4 + 3/8))

1 番目のものが適解である。

甲円の半径

res[1][1] |> println

   (-11 - 19*sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^2 + (-1 + sqrt(3)/3)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2) + (19*sqrt(3)/3 + 11)*(-1/2 - sqrt(-10*sqrt(3) - sqrt(6)*sqrt(97 - 56*sqrt(3)) + 35/2) + sqrt(3)/2)^3 + sqrt(3)

乙円の半径

res[1][3] |> println

   (-22 - 10*sqrt(3))*(1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^2 + 1/2 + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))^3*(8*sqrt(3) + 16) + (1/2 - sqrt(-sqrt(3)/8 - sqrt(2 - sqrt(3))/4 + 3/8))*(2*sqrt(3) + 5)

甲円と乙円の式は複雑で長いが，二重根号を外して簡約化すると短い式になる。

r1 = res[1][1] |> sympy.sqrtdenest |> simplify
r2 = res[1][3] |> sympy.sqrtdenest |> simplify
r1 |> println
r2 |> println

   -sqrt(2) - 1/2 + sqrt(3)/2 + sqrt(6)/2
   -sqrt(2)/2 - 1/4 + sqrt(3)/4 + sqrt(6)/4

乙円の半径 / 甲円の半径 を求める。

r2/r1 |> simplify |> println

   1/2

乙円の直径は甲円の直径の 1/2 である。

   a = 1;  r1 = 0.176557;  r2 = 0.0882784;  r2/r1 = 0.5

x1, y2 も簡単な式になる。

res[1][2] |> sympy.sqrtdenest |> simplify |> println
res[1][4] |> sympy.sqrtdenest |> simplify |> println

   -sqrt(6) - sqrt(3)/2 + 3/2 + 3*sqrt(2)/2
   -sqrt(2)/4 + 1/4 + sqrt(3)/4

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   (r1, x1, r2, y2) = res[1]
   @printf("a = %g;  r1 = %g;  r2 = %g;  r2/r1 = %g\n", a, r1, r2, r2/r1)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:black, lw=0.5)
   plot!([a, 0, (√3 - 1)a, a], [a, (√3 - 1)a, 0, a], color=:blue, lw=0.5)
   plot!([0, a, a/2, 0], [0, 0, √3a/2], color=:red, lw=0.5)
   circle(x1, r1, r1, :magenta)
   circle(a - r2, y2, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       r3m1a = (√3 - 1)a
       point(0, r3m1a, "   (√3-1)a", :blue, :left, :vcenter)
       point(r3m1a, 0, " (√3-1)a", :blue, :left, :bottom)
       point(a, 0, " a", :blue, :left, :bottom)
       point(0, a, " a", :blue, :left, :bottom)
       point(a/2, √3a/2, "(a/2,√3a/2)", :red, :right, :bottom)
       point(x1, r1, "甲円:r1,(x1,r1)", :magenta, :center, delta=-delta/2)
       point(a - r2, y2, "乙円:r2\n(a-r2,y2)", :green, :center, delta=2.3delta)
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;