算額（その508）

算額（その508）

秋田県仙北市角館町西長野　熊野神社 嘉永2年(1849)
http://www.wasan.jp/akita/kakunodatekumano1.html
小寺裕: 算額あれこれ
秋月めぐる，遠藤寛子: 算法少女(1)，リイド社，2013.

鈎股弦（直角三角形内）に二等辺三角形と 2 個の正方形が入っている。鈎の長さが 1 のとき，股の長さはいかほどか。
原点と (x,y) を結ぶ線分は，弦と直交する。

図のように記号を定め，連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms 鈎::positive, 股::positive, 弦::positive, x::positive, y::positive, x2::positive, y2::positive

eq1 = y/(股 - x) - 鈎/股
eq2 = 股*鈎 - 弦*sqrt(x^2 + y^2)
eq3 = 2(x - x2) - y2
eq4 = y2/(股 - 2x -y2) - 鈎/股
eq5 = 鈎^2 + 股^2 - 弦^2
eq6 = y/x - y2/x2;
res = solve([eq1, eq2, eq3, eq4, eq5, eq7], (股, 弦, x, y, x2, y2))

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)
   (股, 弦, x, y, x2, y2) = u
   return [
       y/(-x + 股) - 鈎/股,  # eq1
       -弦*sqrt(x^2 + y^2) + 股*鈎,  # eq2
       2*x - 2*x2 - y2,  # eq3
       y2/(-2*x - y2 + 股) - 鈎/股,  # eq4
       -弦^2 + 股^2 + 鈎^2,  # eq5
       -y2/x2 + y/x,  # eq6
   ]
end;

鈎 = 1
iniv = BigFloat[2.0, 2.2, 0.2, 0.8, 0.1, 0.5]
res = nls(H, ini=iniv)

   (BigFloat[1.999999999999999999985926409326636355492543560353991554111136639959569147152117, 2.236067977499789696396585866558016976252953289080026164208793498254395925337324, 0.3999999999999999999976356367668749077227473181394705810906711610880915130403816, 0.7999999999999999999997748225492261816878806969656638648657781909543976392427267, 0.1999999999999999999982548747565029080810754014838949527097809876602829734671478, 0.3999999999999999999987615240207439992833438333111512567617803468556170791464633], true)

   股 = 2;  弦 = 2.23607;  x = 0.4;  y = 0.8;  x2 = 0.2;  y2 = 0.4

鈎の長さが 1 のとき，股の長さは 2 である。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (股, 弦, x, y, x2, y2) = res[1]
   @printf("股 = %g;  弦 = %g;  x = %g;  y = %g;  x2 = %g;  y2 = %g\n", 股, 弦, x, y, x2, y2)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   rect(x2, 0, x2 + y2, y2, :red)
   rect(2x, 0, 2x + y2, y2, :red)
   plot!([0, x, 2x], [0, y, 0], color=:blue, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta)
       point(2x, 0, " 2x", :blue, :left, :bottom, delta=delta)
       point(x2, y2, " (x2,y2)", :red, :left, :bottom, delta=delta)
       point(2x + y2, y2, " (2x+y2,y2)", :red, :left, :bottom, delta=delta)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

 

算額（その507）

算額（その507）

長岡高専和算倶楽部オリジナル算額
山田 章，谷口 悠，加藤 祐樹，長谷川 柊太，平田 大成，中山 雅友美，冨樫 瑠美，涌田 和芳: 長岡高専和算倶楽部オリジナル算額について，長岡工業高等専門学校研究紀要，57 巻，p. 41-45，2021.
https://kinpoku.nagaoka-ct.ac.jp/lib/kiyo/vol_57/57_41yamada.pdf

正三角形に内接する正方形（正方形の3つの頂点が正三角形の 3 つの辺上にある）の一辺の長さが最小になるのはどのようなときか。

正三角形の一辺の長さを a，正方形の一辺の長さを b とする。
正方形の 4 つの頂点座標を (x1, 0 ), (x2, y2), (x3, y3), (x4, y4) とする。
(x1, 0), (x2, y2) を結ぶ線分が正三角形の底辺となす角を θ° として，
以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, θ, x1, x2, y2, x3, y3, x4, y4

x2 = x1 + b*cosd(θ)
y2 =      b*sind(θ)
x3 = x2 - b*sind(θ)
y3 = y2 + b*cosd(θ)
x4 = x1 - b*sind(θ)
y4 =      b*cosd(θ)

eq1 = y4/x4 - sqrt(Sym(3))
eq2 = y3/(a - x3) - sqrt(Sym(3))
res = solve([eq1, eq2], (x1, b))

   Dict{Any, Any} with 2 entries:
     b  => sqrt(3)*a/(2*(sin(pi*(θ/180 + 1/3)) + cos(pi*θ/180)))
     x1 => a*sin(pi*(θ/180 + 1/6))/(sin(pi*(θ/180 + 1/3)) + cos(pi*θ/180))

正方形の一辺の長さは sqrt(3)/(2*(sin(pi*(θ/180 + 1/3)) + cos(pi*θ/180))) である。
0 ≦ θ ≦ 30 で図を描くと以下のようになる。

using Plots
pyplot(size=(400, 200), grid=false, aspectratio=:none, label="", fontfamily="IPAMincho")
plot(res[b]/a, xlims=(0, 30), xlabel="θ", ylabel="b")

b の式の導関数を求め，導関数 = 0 となる θ を求め，そのときの b の値を求める。
正三角形の一辺の長さを a とすれば，θ が 15° のとき，正方形の一辺の長さは最小値 (3√2 - √6)a/4 = 0.448287736084027a になる。

solve(diff(res[b]/a))[1] |> println

   15

eq = (res[b])(θ => 15) |> simplify
eq |> println
eq(θ => 15).evalf() |> println

   -sqrt(6)*a/4 + 3*sqrt(2)*a/4
   0.448287736084027*a

   a = 1;  θ = 15;  x1 = 0.366025;  b = 0.448288

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1
   θ = 15
   b  = sqrt(3)*a/(2*(sin(pi*(θ/180 + 1/3)) + cos(pi*θ/180)))
   x1 = a*sin(pi*(θ/180 + 1/6))/(sin(pi*(θ/180 + 1/3)) + cos(pi*θ/180))
   @printf("a = %g;  θ = %g;  x1 = %g;  b = %g\n", a, θ, x1, b)
   plot([0, a, a/2, 0], [0, 0, √3*a/2, 0], color=:black, lw=0.5)
   x2 = x1 + b*cosd(θ)
   y2 =      b*sind(θ)
   x3 = x2 - b*sind(θ)
   y3 = y2 + b*cosd(θ)
   x4 = x1 - b*sind(θ)
   y4 =      b*cosd(θ)
   plot!([x1, x2, x3, x4, x1], [0, y2, y3, y4, 0], color=:blue, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, 0, "(x1,0) ", :blue, :right, :bottom, delta=delta/2)
       point(x2, y2, "(x2,y2)", :blue, :left, delta=-delta/2)
       point(x3, y3, " (x3,y3)", :blue, :left, :vcenter)
       point(x4, y4, " (x4,y4)", :blue, :left, delta=-delta/2)
       point(a/2, √3a/2, " (a/2, √3a/2)", :black, :left, :vcenter)
       point(a, 0, "a", :black, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その506）

算額（その506）

兵庫県伊丹市 猪名野神社 嘉永6年(1853)
http://www.wasan.jp/hyogo/inano.html
宮崎興二: ある異国人の見た算額，和算，第20号，昭和53年1月1日発行
http://www.wasan.jp/kinkikaisi/wasan020.pdf

正方形内に 5 個の等円が入っている。正方形の一辺の長さが与えられたとき，等円の直径はいかほどか。

正方形の一辺の長さを a，等円の半径を r として，以下の方程式を解く。

include("julia-source.txt");

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

   5*a/26

等円の半径は正方形の一辺の長さの 5/26 倍である。
正方形の一辺の長さが 26 のとき，等円の直径は 10 である。

using Plots

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

Windowsで「RとRStudioをGUIで簡単にインストールする」方法

ChatGPT

WindowsでRとRStudioをGUIで簡単にインストールする方法は以下の通りです。

1. Rのインストール

まず、Rをインストールします。

• R 公式ウェブサイト: https://cran.r-project.org/mirrors.html
  • ウェブサイトにアクセスし、ページ中央部の「Download and Install R」セクションから最新のRをダウンロードします。
  • ダウンロードが完了したら、ダウンロードした .exe ファイルを実行してRをインストールします。

2. RStudioのインストール

次に、RStudioをインストールします。

• RStudio 公式ウェブサイト: https://www.rstudio.com/products/rstudio/download/
  • ウェブサイトにアクセスし、ページ中央部の「RStudio Desktop」セクションから最新のRStudioをダウンロードします。
  • ダウンロードが完了したら、.exeファイルを実行してRStudioをインストールします。

### 注意事項

• インストールの際には、インストーラーが表示する指示に従って進めてください。
• インストールが完了したら、RとRStudioを起動して正しく動作するか確認してください。

以上で、GUIを使用して簡単にRとRStudioをWindowsにインストールすることができます。

macOS で「R と RStudio を超簡単に GUI でインストールする」方法

ChatGPT

macOSでRとRStudioをGUIで簡単にインストールする方法はいくつかあります。以下はその一例です。なお、手順やダウンロード先の情報は変更される可能性があるため、最新の情報を確認することをお勧めします。

1.  Rのインストール

まず、Rをインストールします。

• R 公式ウェブサイト: https://cran.r-project.org/mirrors.html

  • ウェブサイトにアクセスし、ページ中央部の「Download and Install R」セクションから最新のRをダウンロードします。
  • ダウンロードが完了したら、.pkgファイルを実行してRをインストールします。

2.  RStudioのインストール

次に、RStudioをインストールします。

• RStudio 公式ウェブサイト: https://www.rstudio.com/products/rstudio/download/

  • ウェブサイトにアクセスし、ページ中央部の「RStudio Desktop」セクションから最新のRStudioをダウンロードします。
  • ダウンロードが完了したら、.dmgファイルを実行してRStudioをインストールします。

注意事項

• インストールの際には、システムの管理者パスワードが必要な場合がありますので、用意しておいてください。
• インストールが完了したら、RとRStudioを起動して正しく動作するか確認してください。

これで、GUIを使用して簡単にRとRStudioをmacOSにインストールすることができます。

算額（その505）

算額（その505）

愛媛県大洲市新谷字山口 金刀比羅宮（現山口神社） 天明8年(1788)
http://www.wasan.jp/ehime/kotohira1.html

鈎股弦（直角三角形）内に，天地人の 3 個の等円が入っている。天円と地円の中心を結ぶ線分の長さが既知 a，鈎は股と弦の差より b だけ大きい。

角が直角である頂点の座標を (0, 0)
鈎，股，弦の長さを鈎，股，弦
天円，地円，人円の半径と中心座標を r, (r, y1), (x2, r), (x2 - 2r, r)
として，以下の連立方程式の数値解を求める。

一般性を失わずに，a = 1 として解く。

include("julia-source.txt");

using SymPy
@syms 鈎, 股, 弦, r, y1, x2, x3, a, b

x2 = x3 + 2r
弦 = sqrt(鈎^2 + 股^2)
eq1 = (x2 - r)^2 + (y1 - r)^2 - a^2
eq2 = (鈎 - y1) + a + (股 - x2) - 弦
eq3 = (x3 - r)^2 + (y1 - r)^2 - 4r^2
eq4 = (y1 - r)/(x2 - r) - 鈎/股
eq5 = (弦 - 股 + b) - 鈎;
# res = solve([eq1, eq3, eq4, eq5, eq7], (鈎, 股, r, y1, x3))

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, y1, x3) = u
   return [
       -a^2 + (-r + y1)^2 + (r + x3)^2,  # eq1
       a - 2*r - x3 - y1 + 股 + 鈎 - sqrt(股^2 + 鈎^2),  # eq2
       -4*r^2 + (-r + x3)^2 + (-r + y1)^2,  # eq3
       (-r + y1)/(r + x3) - 鈎/股,  # eq4
       b - 股 - 鈎 + sqrt(股^2 + 鈎^2),  # eq5
   ]
end;

a = 1
b = 0.9

iniv = iniv = BigFloat[1.4, 2.9, 0.33, 0.87, 0.77]
res = nls(H, ini=iniv)

   (BigFloat[1.177379233363580301745329510921967893577066939399416499239250156741158075407193, 2.36009488557904520374835794798416053632581553124513131479814373696319717422653, 0.2793822120308020755416636618053848905748184134306487573156428076934921573204049, 0.7257863018027224197024724603072078876900065179854659272910142447388136755685596, 0.6154492741356734514186607085851531396329900169696428080777001398742020097906305], true)

a = 1, b = 0.9 のとき，等円の直径はほぼ 0.5587644240616042 である。

鈎 = 1.17738;  股 = 2.36009;  弦 = 2.63747;  r = 0.279382;  y1 = 0.725786;  x2 = 1.17421;  x3 = 0.615449
天円と地円の中心を結ぶ線分の長さ = sqrt((x2 - r)^2 + (y1 - r)^2)) = a = 1
鈎(1.17738) は 股(2.36009) と 弦(2.63747) の 差(0.277379) より b = 0.9 だけ大きい
等円の直径 = 0.558764　(0.5587644240616042)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (鈎, 股, r, y1, x3) = res[1]
   x2 = x3 + 2r
   弦 = sqrt(鈎^2 + 股^2)
   @printf("鈎 = %g;  股 = %g;  弦 = %g;  r = %g;  y1 = %g;  x2 = %g;  x3 = %g\n", 鈎, 股, 弦, r, y1, x2, x3)
   @printf("天円と地円の中心を結ぶ線分の長さ = sqrt((x2 - r)^2 + (y1 - r)^2)) = a = %g\n", sqrt((x2 - r)^2 + (y1 - r)^2))
   @printf("鈎(%g) は 股(%g) と 弦(%g) の 差(%g) より b = %g だけ大きい\n", 鈎, 股, 弦, 弦 - 股, 鈎 - (弦 - 股))
   @printf("等円の直径 = %g　(%.16g)\n", 2r, 2r)
   plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(r, y1, r)
   circle(x2, r, r)
   circle(x3, r, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, 鈎, " 鈎", :black, :left, :bottom, delta=delta)
       point(股, 0, "股", :black, :left, :bottom, delta=delta)
       point(r, y1, "天円(r,y1)", :red, :center, :bottom, delta=delta)
       point(x2, r, "地円(x2,r)", :red, :center, :top, delta=-delta)
       point(x2 - 2r, r, "人円(x2-2r,r)", :red, :center, :top, delta=-delta)
       segment(r, y1, x2, r, :blue)
   else
       plot!(showaxis=false)
   end
end;

算額（その504）

算額（その504）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8

算額（その164)
https://blog.goo.ne.jp/r-de-r/e/ad36c91a7122981d5a236e8addbe3cff
と基本的に同じであるが，乙円，丙円の半径を未知数として一般解を求める。

二等辺三角形内に甲円，乙円，丙円を入れる。乙円と丙円の直径が与えられたとき，丙円の直径を求めよ。

一般性を失わずに，丙円の直径を 1 とする。乙円と丙円の直径が与えられたとき，乙円の直径は「乙円の直径/丙円の直径」，丙円の直径は 1 として与える。結果は「丙円の直径/乙円の直径」倍になる。

二等辺三角形の底辺の長さを 2a，高さを b とする
甲円の半径と中心座標を r1, (0, r1)
乙円の半径と中心座標を r2, (0, 2r1 + r2)
丙円の半径と中心座標を r3, (x3, r3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, r1, r2, r3, x3
r3 = 1
eq1 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq2 = r3/(a - x3) - r1/a
eq3 = r2/(b - 2r1 - r2) - r1/(b - r1)
eq4 = r1/(b - r1) - a/sqrt(a^2 + b^2)
res = solve([eq1, eq2, eq3, eq4], (a, b, r1, x3))

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

4 組の解が得られるが，4番目の解が適解である。

甲円の直径は， 4r3*(sqrt(r2/r3) + 1/2) である。

(r2, r3) = (1, 2)
t = r2/r3
u = 2sqrt(t) + 1

a = r3*(u^(3/2)/sqrt(t))
b = r3*(2u^2/(u - t))
r1 = r3*u
x3 = r3*(2sqrt(u))
@printf("r2 = %g;  r3 = %g;  a = %g;  b = %g;  r1 = %g;  x3 = %g\n", r2, r3, a, b, r1, x3)
@printf("甲円の直径 = %g\n", 4r3*(sqrt(r2/r3) + 1/2))

   r2 = 1;  r3 = 2;  a = 10.6098;  b = 12.1793;  r1 = 4.82843;  x3 = 6.2151
   甲円の直径 = 9.65685

乙円，丙円の半径を 1, 2 としたとき，甲円の直径は 9.65685 である。

using Plots

function draw(r2, r3, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = r2/r3
   u = 2sqrt(t) + 1

   a = r3*(u^(3/2)/sqrt(t))
   b = r3*(2u^2/(u - t))
   r1 = r3*u
   x3 = r3*(2sqrt(u))
   @printf("r2 = %g;  r3 = %g;  a = %g;  b = %g;  r1 = %g;  x3 = %g\n", r2, r3, a, b, r1, x3)
   @printf("甲円の直径 = %g\n", 4r3*(sqrt(r2/r3) + 1/2))
   plot([a, 0, -a, 0], [0, b, 0, 0,], color=:black, lw=0.5)
   circle(0, r1, r1)
   circle(0, 2r1 + r2, r2, :blue)
   circle(x3, r3, r3, :green)
   circle(-x3, r3, r3, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1, " 甲円:r1,(0,r1)", :magenta, :left, :vcenter)
       point(0, 2r1 + r2, " 乙円:r2,(0,2r1+r2)", :magenta, :left, :vcenter)
       point(x3, r3, " 丙円:r3,(x3,r3)", :magenta, :center, :top, delta=-delta)
   else
       plot!(showaxis=false)
   end
end;

   r2 = 1;  r3 = 4;  a = 22.6274;  b = 18.2857;  r1 = 8;  x3 = 11.3137
   甲円の直径 = 16

 

 

 

算額（その503）

算額（その503）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8
https://dl.ndl.go.jp/pid/3508431/1/9

長方形内に2本の斜線を引き，分割された領域に甲円，乙円を入れる。
長方形の長辺と短辺の長さが与えられたとき，乙円の直径を求めよ。

長方形の長辺と短辺の長さを 2a, 2r1 とする。
甲円の半径と中心座標を r1, (a - r1, 0)
乙円の半径と中心座標を r2, (0, r1 - r2)
斜線と長方形の上辺の交点を (b, r1)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a, b, r1, r2
eq1 = b/sqrt(b^2 + r1^2) - r2/(r1 - r2)
eq2 = r1/(a - r1) - r1/sqrt(b^2 + r1^2);
eq2 = (a - r1) - sqrt(b^2 + r1^2);
res = solve([eq1, eq2], (r2, b))

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

2 組の解が得られるが，2 番目のものが適解である。
SymPy では式を簡約化できないが，容易に以下のようにできる。

b = sqrt(a^2 - 2a*r1)
r2 = (r1*b)/(b + a - r1)

長方形の長辺と短辺が 150, 72 のとき，乙円の半径は 10，直径は 20 である。

長 = 150;  平 = 72;  a = 75;  r1 = 36;  b = 15;  r2 = 10
乙円の直径 = 20

using Plots

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

算額（その502）

算額（その502）

高山忠直編: 算法評論
国立国会図書館 デジタルコレクション
https://dl.ndl.go.jp/pid/3508431/1/8

甲円，乙円，丙円が図のように外接しあっている。甲円，乙円の直径が与えられたとき，丙円の直径を求めよ。

乙円の半径と中心座標を r2, (0, 0 )
甲円の半径と中心座標を r1, (r1, y1)
丙円の半径と中心座標を r3, (0, r2 + r3)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms r1, y1, r2, r3

eq1 = r1^2 + (y1 - r2 - r3)^2 - (r1 + r3)^2
eq2 = r1^2 + y1^2 - (r1 + r2)^2;
res = solve([eq1, eq2], (y1, r3));

2 組の解が得られるが，2 番目のものが適解である。
簡約化すると以下のようになる。

# y1
res[1][1] |> simplify |> println

   sqrt(r2^3*(2*r1 + r2))/r2

# r3: 乙円の半径
res[1][2] |> simplify |> println

   (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2)

   r1 = 5;  r2 = 7;  y1 = 10.9087;  r3 = 0.857477
   甲円の直径が 10，乙円の直径が 14 のとき，丙円の直径は 1.71495 である。

using Plots

function draw(r1, r2, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (y1, r3) = (sqrt(r2^3*(2*r1 + r2))/r2, (r2*(r1 + 2*r2) - 2*sqrt(r2^3*(2*r1 + r2)))/(r1 - 4*r2))
   @printf("r1 = %g;  r2 = %g;  y1 = %g;  r3 = %g\n", r1, r2, y1, r3)
   @printf("甲円の直径が %g，乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r1, 2r2, 2r3)
   plot()
   circle(r1, y1, r1)
   circle(-r1, y1, r1)
   circle(0, 0, r2, :blue)
   circle(0, r2 + r3, r3, :green) 
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(r1, y1, " 甲円:r1\n (r1,y1)", :red, :left, :vcenter)
       point(0, 0.4r2, " 乙円:r2,(0,0)", :blue, :left, mark=false)
       point(0, r2 + r3, " 丙円:r3,(0,r2+r3)", :green, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

算額（その501）

算額（その501）

山形県鶴岡市茨新田 鶴岡神明宮 文久2年(1862)
http://www.wasan.jp/yamagata/sinmei.html

外円内に 4 本の斜線を引き，区画された領域に大円，中円，小円を入れる。
外円，大円，小円の直径がそれぞれ 125寸，99寸，21寸のとき，中円の直径はいかほどか。

外円の半径と中心座標を r0, (0, 0 )
大円の半径と中心座標を r1, (0, r1 - r0)
中円の半径と中心座標を r2, (x2, y2)
小円の半径と中心座標を r3, (0, r0 - r3)
斜線と外円の交点座標を (a, sqrt(r0^2 - a^2), (b, sqrt(r0^2 - b^2)
とおき，以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms a, b, r0, r1, r2, x2, y2, r3

eq1 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), 0, r1 - r0) - r1^2
eq2 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), x2, y2) - r2^2
eq3 = distance(a, sqrt(r0^2 - a^2), -b, sqrt(r0^2 - b^2), 0, r0 -r3) - r3^2
eq4 = distance(-a, sqrt(r0^2 - a^2), b, sqrt(r0^2 - b^2), x2, y2) - r2^2
eq5 = distance(a, sqrt(r0^2 - a^2), b, sqrt(r0^2 - b^2), x2, y2) - r2^2;

using NLsolve

function nls(func, params...; ini = [0.0])
   if typeof(ini) <: Number
       r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=big"1e-40")
       v = r.zero[1]
   else
       r = nlsolve((vout, vin)->vout .= func(vin, params...), ini, ftol=big"1e-40")
       v = r.zero
   end
   return v, r.f_converged
end;

function H(u)
   (a, b, r2, x2, y2) = u
   t1 = sqrt(-a^2 + r0^2)
   t2 = sqrt(-b^2 + r0^2)
   return [
       -r1^2 + (t1 - t2)^2*(a^2*b - a*b^2 + a*r0^2 - a*r0*t1 + a*r0*t2 + a*r1*t1 - a*r1*t2 - a*t1*t2 - b*r0^2 - b*r0*t1 + b*r0*t2 + b*r1*t1 - b*r1*t2 + b*t1*t2)^2/(4*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)^2) + (-r0 + r1 - (-a^3*b - a^2*r0^2 + a^2*r0*t1 - 3*a^2*r0*t2 - a^2*r1*t1 + 3*a^2*r1*t2 + a^2*t1*t2 + a*b^3 + b^2*r0^2 + 3*b^2*r0*t1 - b^2*r0*t2 - 3*b^2*r1*t1 + b^2*r1*t2 - b^2*t1*t2 - 4*r0^3*t1 + 4*r0^3*t2 + 4*r0^2*r1*t1 - 4*r0^2*r1*t2)/(2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)))^2,  # eq1
       -r2^2 + (x2 - (x2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (t1 - t2)*(a^2*b + a^2*x2 - a*b^2 + a*r0^2 + a*y2*t1 - a*y2*t2 - a*t1*t2 + b^2*x2 - b*r0^2 + b*y2*t1 - b*y2*t2 + b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2 + (y2 - (y2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) - (a + b)*(a^2*b + a^2*x2 - a*b^2 + a*r0^2 + a*y2*t1 - a*y2*t2 - a*t1*t2 + b^2*x2 - b*r0^2 + b*y2*t1 - b*y2*t2 + b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2,  # eq2
       -r3^2 + (t1 - t2)^2*(a^2*b - a*b^2 + a*r0^2 + a*r0*t1 - a*r0*t2 - a*r3*t1 + a*r3*t2 - a*t1*t2 - b*r0^2 + b*r0*t1 - b*r0*t2 - b*r3*t1 + b*r3*t2 + b*t1*t2)^2/(4*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)^2) + (r0 - r3 - (-a^3*b - a^2*r0^2 - a^2*r0*t1 + 3*a^2*r0*t2 + a^2*r3*t1 - 3*a^2*r3*t2 + a^2*t1*t2 + a*b^3 + b^2*r0^2 - 3*b^2*r0*t1 + b^2*r0*t2 + 3*b^2*r3*t1 - b^2*r3*t2 - b^2*t1*t2 + 4*r0^3*t1 - 4*r0^3*t2 - 4*r0^2*r3*t1 + 4*r0^2*r3*t2)/(2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2)))^2,  # eq3
       -r2^2 + (x2 - (x2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (t1 - t2)*(-a^2*b + a^2*x2 + a*b^2 - a*r0^2 - a*y2*t1 + a*y2*t2 + a*t1*t2 + b^2*x2 + b*r0^2 - b*y2*t1 + b*y2*t2 - b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2 + (y2 - (y2*(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2) + (a + b)*(-a^2*b + a^2*x2 + a*b^2 - a*r0^2 - a*y2*t1 + a*y2*t2 + a*t1*t2 + b^2*x2 + b*r0^2 - b*y2*t1 + b*y2*t2 - b*t1*t2 - 2*r0^2*x2 + 2*x2*t1*t2)/2)/(a^2*t2 + a*b*t1 - a*b*t2 - b^2*t1 + 2*r0^2*t1 - 2*r0^2*t2))^2,  # eq4
       -r2^2 + (x2 - (a*(a*b - r0^2 + t1*t2) - (a - b)*(a*b + a*x2 - b*x2 - r0^2 + y2*t1 - y2*t2 + t1*t2)/2)/(a*b - r0^2 + t1*t2))^2 + (y2 - (2*t1*(a*b - r0^2 + t1*t2) - (t1 - t2)*(a*b + a*x2 - b*x2 - r0^2 + y2*t1 - y2*t2 + t1*t2))/(2*(a*b - r0^2 + t1*t2)))^2,  # eq5
   ]
end;
(r0, r1, r3) = (125, 99, 21) .// 2
iniv = BigFloat[60.0, 33, 10, 29, 40]
res = nls(H, ini=iniv)

  (BigFloat[60.5769230769230769230769230769230769230769230769230769230769230769230769357718, 31.73076923076923076923076923076923076923076923076923076923076923076923419539377, 10.81730769230769230769230769230769230769230769230769230769230769230769343490884, 28.12500000000000000000000000000000000000000000000000000000000000000000245224407, 40.62499999999999999999999999999999999999999999999999999999999999999999867465516], true)

   r0 = 62.5;  r1 = 49.5;  r3 = 10.5
   a = 60.5769;  b = 31.7308;  r2 = 10.8173;  x2 = 28.125;  y2 = 40.625
   中円の直径 = 2r2 = 21.6346

中円の直径は 21.6346 ほどである。術では「二十一寸五十二分...」とあるので，一致しない。

実際の図は，算額に描かれたものとかなり異なる。
このようなことは算額においてはよくあることであるが，原因の一つは算額の図における各円の直径の割合が問に書かれたものと違うこと。
さらに，この算額の場合には外円の直径が「一百二十五寸」と書かれているのだが，「一」の部分が不鮮明であること。「二」でも「三」でもなさそうではあるのだが。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1, r3) = (125, 99, 21) .// 2
   (a, b, r2, x2, y2) = res[1]
   @printf("r0 = %g;  r1 = %g;  r3 = %g\n", r0, r1, r3)
   @printf("a = %g;  b = %g;  r2 = %g;  x2 = %g;  y2 = %g\n", a, b, r2, x2, y2)
   @printf("中円の直径 = 2r2 = %g\n", 2r2)
   plot()
   circle(0, 0, r0, :black) 
   circle(0, r1 - r0, r1, :red) 
   circle(x2, y2, r2, :blue)
   circle(-x2, y2, r2, :blue)
   circle(0, r0 - r3, r3, :magenta)
   plot!([a, b, -a], [sqrt(r0^2 - a^2), sqrt(r0^2 - b^2), sqrt(r0^2 - a^2)], color=:brown, lw=0.5)
   plot!([-a, -b, a], [sqrt(r0^2 - a^2), sqrt(r0^2 - b^2), sqrt(r0^2 - a^2)], color=:brown, lw=0.5)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1 - r0, " 大円:r1,(0,r1-r0)", :red, :left, :vcenter)
       point(x2, y2, " 中円:r2,(x2,y2)", :blue, :center, :top, delta=-delta)
       point(0, r0 - r3, " 小円:r3\n (0,r0-r3)", :magenta, :left, :vcenter)
       point(r0, 0, "r0 ", :black, :right, :bottom)
       point(a, sqrt(r0^2 - a^2), "(a,√(r0^2-a^2)", :brown, :right)
       point(b, sqrt(r0^2 - b^2), " (b,√(r0^2-b^2)", :brown, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

算額（その500）

算額（その500）

群馬県高崎市 清水寺 文政11年(1828)
http://takasakiwasan.web.fc2.com/gunnsann/pdf/046_2monndai.pdf

長径 7 寸，短径 5 寸の菱形の中に楕円と小円 4 個が入っている。小円の直径はいかほどか。

菱形の中心を原点とし，菱形の長径，短径を 2x，2y とする。
楕円の長径，短径を 2a，2b とする。
等円の半径と中心座標を r1, (a + r, 0), (0, b + r)
として以下の連立方程式の数値解を求める。

include("julia-source.txt");

using SymPy
@syms x, y, a, b, r, x1, y1

z = sqrt(x^2 + y^2)
eq1 = r/(x - a - r) - y/z  # 半径と短径の関係
eq2 = r/(y - b - r) - x/z  # 半径と長径の関係
eq3 = (b^2*x1)/(a^2*y1) - y/x  # 接線の傾きに関して
eq4 = b^2/y1 - y  # 接線の切片に関して
eq5 = x1^2/a^2 + y1^2/b^2 - 1;  # 接点が楕円状にあることに関して
# res = solve([eq1, eq2, eq3, eq4, eq5], (a, b, r, x1, y1))

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)
   (a, b, r, x1, y1) = u
   return [
       r/(-a - r + x) - y/sqrt(x^2 + y^2),  # eq1
       r/(-b - r + y) - x/sqrt(x^2 + y^2),  # eq2
       -y/x + b^2*x1/(a^2*y1),  # eq3
       b^2/y1 - y,  # eq4
       -1 + y1^2/b^2 + x1^2/a^2,  # eq5
   ]
end;
(x, y) = (7, 5) .// 2
iniv = BigFloat[2.55, 1.7, 0.35, 2, 1] .* (x/3.5)
res = nls(H, ini=iniv)

   (BigFloat[2.54414717679825272253216901118052540925755526167623360938083320066487243514035, 1.716860438920218883208445479192299273530396028742001559200850725830125354563721, 0.3513564057748656326259665107244486935073587532847506507495566147571622830042101, 1.849338530631577084618933179122660657165347568695981336767010069217909175976173, 1.179043906691730653843619157769528102024751736645727616594992765005302482977008], true)

x = 3.5;  y = 2.5;  a = 2.54415;  b = 1.71686;  r = 0.351356
円の直径 = 2r = 0.702713

x = 3.5, y = 2.5 のとき，円の直径は 0.702713 である。

算額の図を再現している（かどうかは怪しいが），元の記事に掲載されている図は，実際のものとはかなり異なっている。

一番の問題は，図では分かりづらいが実測してみると左右の円と上下の円の直径が異なっていることである。比は 1.125 程度も違う。
また，菱形の長径と短径の比は，問では 7/5=1.4 であるが，図では 1.64 程度である。
もし菱形の長径と短径が条件通り 7, 5 で，円の直径が 1 とすれば，楕円は菱形に内接しない。

いずれにしろ，得られた数値で図形を描くと条件を満たす図になっているので，得られた結果として「円の直径は 0.702713 である」は正しいと思われる。

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (x, y) = (7, 5) ./ 2
   (a, b, r, x1, y1) = res[1]
   @printf("x = %g;  y = %g;  a = %g;  b = %g;  r = %g\n", x, y, a, b, r)
   @printf("円の直径 = 2r = %g\n", 2r)
   plot([x, 0, -x, 0, x], [0, y, 0, -y, 0], color=:black, lw=0.5)
   circle(0, b + r, r) 
   circle(0, -b - r, r) 
   circle(a + r, 0, r)
   circle(-a - r, 0, r)
   ellipse(0, 0, a, b, color=:blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(x1, y1, " (x1,y1)", :blue, :left, :vcenter)
       point(x, 0, " x", :black, :left, :bottom, delta=delta/2)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta/2)
       point(a + r, 0, "a+r", :red, :center, :top, delta=-delta/2)
       point(0, y, " y", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :blue, :left, :top, delta=-delta/2)
       point(0, b + r, " b+r", :red, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;

算額（その499）

算額（その499）

埼玉県北本市本宿 天神社 明治24年(1891)

北本デジタルアーカイブズ 北本市史 通史編 近世
https://kdas.jp/detail_display.php?t_cd=1&acc_cd=1&aclc_cd=40&chap=3&hp_page=47&bc_cd=1

正方形内に斜線4本で領域を区切り，甲円 1 個，乙円 4 個を入れる。甲円の直径が 1 寸 8 分のとき，乙円の直径はいかほどか。

この算額は，算額（その204）
https://blog.goo.ne.jp/r-de-r/e/698f42011195b818050f5e5bc17ed20e
から最もも小さい円を除いたものである。
算額（その204）では円の中心から斜線までの距離を使って解いたが，もう少し簡単にできる。

正方形の中心を原点とし，正方形の一辺の長さを 2a とする。
甲円の半径と中心座標を r1, (0, 0)
乙円の半径と中心座標を r2, (a/2, a - r2)
として，以下の連立方程式をとき，r2, a を求める（それぞれは r1 を含む式になる）。

「三角函数を援用しなければ解けないかなり高度のもの」という解説があるが，三角関数を使う必要はない。
各円を含む直角三角形と各円の半径を一辺とする直角三角形の相似関係を使えばよい。

include("julia-source.txt");

using SymPy
@syms a::positive, r1::positive, r2::positive;

eq1 = r1/a - a/sqrt(a^2 + (2a)^2)  # 甲円を含む直角三角形との相似関係
eq2 = r2/(a - r2) - (a/2)/sqrt((a/2)^2 + a^2)  # 乙円を含む直角三角形との相似関係
res = solve([eq1, eq2], (r2, a))

   Dict{Any, Any} with 2 entries:
     a  => sqrt(5)*r1
     r2 => r1*(5 - sqrt(5))/4

甲円の半径を 1.8/2 = 0.9 寸とすれば，乙円の半径は 0.6218847050625473, 直径は 1.2437694101250947 である。
ちなみに，正方形の一辺の長さは 4.024922359499622 である。

2*0.9*(5 - √5)/4, 2(√5*0.9)

   (1.2437694101250947, 4.024922359499622)

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1.8/2 # 半径が 0.9 寸のとき
   (r2, a) = (r1*(5 - sqrt(5))/4, sqrt(5)*r1)
   @printf("r1 = %g;  r2 = %g;  a = %g\n", r1, r2, a)
   @printf("乙円の直径 = 2r2 = %g\n", 2r2)
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   plot!([-a, 0, a],[  a, -a, a], color=:green, lw=0.5)
   plot!([-a, 0, a],[ -a,  a, -a], color=:green, lw=0.5)
   circle(0, 0, r1)
   circle4(a/2, a - r2, r2, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r1, " r1", :red, :left, :bottom, delta=delta/2)
       point(0.5r1, 0.2r1, "甲円:r1,(0,0)", :red, :center, mark=false)
       point(a/2, a - r2, "乙円:r2,(a/2,a-r2)", :blue, :center, delta=-delta)
       point(a/2, 0, " a/2", :black, :left, :bottom, delta=delta/2)
       point(a, 0, " a", :black, :left, :bottom, delta=delta/2)
       point(0, a, " a", :black, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その498）

算額（その498）

埼玉県北本市本宿 天神社 明治24年(1891)

北本デジタルアーカイブズ 北本市史 通史編 近世
https://kdas.jp/detail_display.php?t_cd=1&acc_cd=1&aclc_cd=40&chap=3&hp_page=47&bc_cd=1

この算額は，算額（その269）
https://blog.goo.ne.jp/r-de-r/e/1fd7e974eb24c0c5118e308a9e790003
のものと類似している。
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
県内の算額（260）
長野県下高井郡木島平村 天満宮 推定明治21年(1888)

両者が異なる点は，
天満宮（その269）のものは，小円が中円 2 個と大円 1 個に外接するが
天神社のものは，小円が中円 3 個と外接する点である。

天神社のものは天満宮（その269）のものより3年後のものなので，問の写し間違いがあったのかもしれない。

大円 2 個，中円 5 個，小円 4 個が配置されている。小円は 3 個の中円に外接している。
大円の直径が 2 寸 8 分のとき，小円の直径を求めよ。

大円，中円の半径を r1, r2 とおく（r2 = r1/2）。小円の直径を r3, 中心座標を (x3, x3) とおき，以下の方程式を解く。条件式が 2 本で，r3, x3 の二変数について解を求める。その解の式には r1 が変数として含まれる。

なお，x3 = r1/2 は図が正しければ自明だが，一応，未知数として解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive, x3::positive;

r2 = r1//2
eq1 = 2x3^2 - (r2 + r3)^2
eq2 = (2r2 - x3)^2 + x3^2 - (r2 + r3)^2;

solve([eq1, eq2], (r3, x3))

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

ここで面白いことに気づいた。
天満宮の算額（その269）では，SymPy で計算した小円の半径が 「3*r1/14」 になったのに，術では「r1*(sqrt(1/2) - 1/2)」としており不一致であった。
そして，天神社の算額（その498）の術は SymPy で計算した通りの「r1*(sqrt(1/2) - 1/2)」である。
天満宮の算額は，答えも術も正しいが図が間違えていたということで，本来同じ問であった（天神社のものは問，図，答，術すべてが正しい）ということである。

直径が 2.8 寸ならば，小円の直径は (sqrt(1/2) - 1/2)*2.8 = 0.5798989873223331 である。

ここで，今度は天神社の答えが間違っている。(sqrt(1/2) - 1/2)*2.8 を計算した結果が 0.574 に切り捨てられている。

(sqrt(1/2) - 1/2)*2.8

   0.5798989873223331

   r3 = 0.289949;  x3 = 0.700000
   小円の直径 = 2r3 = 0.579899

using Plots

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 28//20  # 半径が 1.4 寸のとき
   r2 = r1//2
   (r3, x3) = (r1*(-1/2 + sqrt(2)/2), r1/2)
   @printf("r3 = %.6f;  x3 = %.6f\n", r3, x3)
   @printf("小円の直径 = 2r3 = %.6f\n", 2r3)
   plot()
   circle(0, r2, 2r2)
   circle(0, -r2, 2r2)
   circle(0, 0, r2, :blue)
   circle42(0, 2r2, r2, :blue)
   circle4(x3, x3, r3, :green)
   if more
       point(r2, 0, " r2")
       point(0, r2, " r2")
       point(2r2, 0, " 2r2")
       point(0, 2r2, " 中円:r2\n 2r2", :blue)
       point(3r2, 0, "3r2 ", :green, :right)
       point(x3, x3, " 小円:r3\n (x3,x3)", :green, :center, :bottom)
       point(1.2r2, 2r2, "大円", :red, mark=false)
       vline!([0], color=:black, lw=0.5)
       hline!([0], color=:black, lw=0.5)
   else
       plot!(showaxis=false)
   end
end;

算額（その497）

算額（その497）

埼玉県さいたま市中央区円阿弥 日枝神社 慶応2年(1866) 正野友三郎一門奉納の算額

与野郷土資料館 展示Web解説（その21）
https://www.city.saitama.jp/004/005/004/005/yono/yonokyodo_tenjikaisetsu/p078011.html

外円の中に2個の三角形と甲円が 2 個，乙円が 6 個入っている。乙円の直径が 1 寸のとき，甲円の直径はいかほどか。

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

ただし，一度に解くと solve() の能力的に解けないことと，数値解を求める場合にも収束が不安定になるので，まずは外円と甲円の r0, r1 を先に求める。

include("julia-source.txt");

using SymPy

@syms r0::positive, r1, r2, x2, y2;

x = sqrt(r0^2 - r2^2)
eq1 = r0*x - (r0 - 2r2)sqrt(x^2 + (r0 - r2)^2)
eq2 = (sqrt(r0^2 - r2^2) + sqrt(r0^2 - r2^2 + (r0 - r2)^2))r1 - sqrt(r0^2 - r2^2)*(r0 - r2);
res1 = solve([eq1, eq2], (r0, r1))

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

解は 2 通り求まるが 2 番目のものが適解である。
すなわち，外円の半径は乙円の半径の 8 倍，甲円の半径は乙円の半径の 3 倍である。
乙円の直径が 1 寸のとき，外円の直径は 8 寸，甲円の直径は 3 寸である。

算額の問の答えとしては以上で十分であるが，図形を描くために右上にある乙円の中心座標を求める。

@syms r0, r1, r2, x2, y2;

(r0, r1) = (8r2, 3r2)
eq3 = y2*r0/(x2*sqrt(r0^2 - r2^2)) - 1
eq4 = x2^2 + y2^2 - (r0 - r2)^2
res2 = solve([eq3, eq4], (x2, y2))

   2-element Vector{Tuple{Sym, Sym}}:
    (-56*sqrt(127)*r2^2/(127*sqrt(r2^2)), -21*sqrt(889)*r2/127)
    (56*sqrt(127)*r2^2/(127*sqrt(r2^2)), 21*sqrt(889)*r2/127)

解は 2 通り求まるが 2 番目のものが適解である。

乙円の直径が 1 寸のとき，
甲円の直径 = 3;  外円の直径 = 8
r0 = 4;  r1 = 1.5;  r2 = 0.5;  x2 = 2.4846;  y2 = 2.46511 である。

using Plots

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r0, r1) = (8*r2, 3*r2)
   (x2, y2) = (56*sqrt(127)*r2^2/(127*sqrt(r2^2)), 21*sqrt(889)*r2/127)
   @printf("甲円の直径 = %g;  外円の直径 = %g;  r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  y2 = %g\n",
       2r1, 2r0, r0, r1, r2, x2, y2)
   x = sqrt(r0^2 - r2^2)
   plot([x, 0, -x, x], [r2, r0, r2, r2], color=:red, lw=0.5)
   plot!([x, 0, -x, x], -[r2, r0, r2, r2], color=:red, lw=0.5)
   circle(0, 0, r0, :gray)
   circle(0, r2 + r1, r1, :blue)
   circle(0, -r2 - r1, r1, :blue)
   circle4(x2, y2, r2, :green)
   circle(r0 - r2, 0, r2, :green)
   circle(r2 - r0, 0, r2, :green)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r2 + r1, " 甲円:r1\n (0,r2+r1)", :blue, :left, :vcenter)
       point(x2, y2, "乙円:r2,(x2,y2)", :green, :center, :bottom, delta=delta/2)
       point(r0 - r2, 0, "r0-r2", :green, :center, :bottom, delta=delta/2)
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta/2)
       point(r0, 0, " r0", :black, :left, :bottom, delta=delta/2)
       point(x, r2, "(√(r0^2-r2^2),r2)  ", :red, :right, :bottom, delta=delta/2)
       point(0, r2, " r2", :green, :left, :bottom, delta=delta/2)
   else
       plot!(showaxis=false)
   end
end;

算額（その496）

算額（その496）

岐阜県郡上市八幡町 八幡神社 嘉永3年(1850)

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

外円内に半径の等しい 6 個の大円が入っている。小円は外円の同心円で，この 6 個の大円に外接する。2個の大円に外接し外円に内接する円を甲円，2個の大円と小円に外接する円を乙円とする。
甲円と乙円の直径が与えられたとき，外円の直径を求めよ。

外円の半径と中心座標を r0, (0, 0)
大円の半径と中心座標を r1, (0, r0 - r1), (x1, y1)
小円の半径と中心座標を r2, (0, 0); r2 = r0 - 2r1
甲円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r4, (x4, y4)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

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

r2 = r0 - 2r1
x1 = (r0 - r1)cosd(Sym(30))
y1 = (r0 - r1)sind(Sym(30))
x3 = (r0 - r3)cosd(Sym(60))
y3 = (r0 - r3)sind(Sym(60))
x4 = (r2 + r4)cosd(Sym(60))
y4 = (r2 + r4)sind(Sym(60))

eq1 = (x1 - x3)^2 + (y1 - y3)^2 - (r1 + r3)^2
eq2 = (x1 - x4)^2 + (y1 - y4)^2 - (r1 + r4)^2

res = solve([eq1, eq2], (r0, r1))

   2-element Vector{Tuple{Sym, Sym}}:
    ((-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), -sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))
    ((sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))

2組の解が得られるが，最初のものが適解である。

甲円，乙円の直径が 8, 2 のとき，外円の直径は約 86.9559 である。
外円の直径 = 86.9559;  r0 = 43.4779;  r1 = 17.3042;  r2 = 8.86948;  r3 = 4,  r4 = 1

using Plots

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r3, r4) = (8, 2) .// 2
   (r0, r1) = ((-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-56*sqrt(3)*r3^2 + 97*r3^2 - 1448*sqrt(3)*r3*r4 + 2508*r3*r4 - 1616*sqrt(3)*r4^2 + 2799*r4^2 + (-sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))*(-2702*r3 + 1560*sqrt(3)*r3 - 2702*r4 + 1560*sqrt(3)*r4))/(-1351*r3^2 + 780*sqrt(3)*r3^2 - 780*sqrt(3)*r4^2 + 1351*r4^2), -sqrt(-16*r3*r4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)^2 + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)^2)/(4*(-26*r3 + 15*sqrt(3)*r3 - 26*r4 + 15*sqrt(3)*r4)) + (-31*sqrt(3)*r3^2 + 54*r3^2 - 58*sqrt(3)*r3*r4 + 100*r3*r4 - 31*sqrt(3)*r4^2 + 54*r4^2)/(4*(-15*sqrt(3)*r3 + 26*r3 - 15*sqrt(3)*r4 + 26*r4)))
   r2 = r0 - 2r1
   x1 = (r0 - r1)cosd(30)
   y1 = (r0 - r1)sind(30)
   x3 = (r0 - r3)cosd(60)
   y3 = (r0 - r3)sind(60)
   x4 = (r2 + r4)cosd(60)
   y4 = (r2 + r4)sind(60)
   @printf("外円の直径 = %g;  r0 = %g;  r1 = %g;  r2 = %g;  r3 = %g,  r4 = %g\n", 2r0, r0, r1, r2, r3, r4)    
   plot()
   circle(0, 0, r0, :black)
   circle(0, 0, r2, :blue)
   rotate(0, r0 - r1, r1, :green, angle=60)
   rotate(r0 - r3,0, r3, :magenta, angle=60)
   rotate(r2 + r4, 0, r4, :red, angle=60)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:black, lw=0.5)
       vline!([0], color=:black, lw=0.5)
       point(0, r0, " r0", :black, :left, :bottom, delta=delta/2)
       point(0, r2, " r2", :blue, :left, :bottom, delta=delta/2)
       point(0, r0 - r1, " r0-r1", :green, :left, :vcenter)
       point(x1, y1, "    大円:r1,(x1, y1)", :green, :center, :bottom, delta=delta/2)
       point(x3, y3, "甲円:r3 \n(x3, y3) ", :black, :right, :vcenter)
       point(x4, y4, " 乙円:r4,(x4, y4) ", :black, :left, :vcenter)
   else
       plot!(showaxis=false)
   end
end;