算額（その2010）

算額（その2010）

(22) 兵庫県養父町左近山 地蔵堂 明治21年(1888)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円1個，正三角形，正方形

正三角形の頂点を通る円と接する正方形の一辺の長さを求めよ。

正三角形の高さを h とする。h は円と正方形がはみ出ないように，十分に大きいこと（正方形の一辺の長さは h に依存しない）。
正三角形の一辺の長さは a = h/√3
半径 r の円の中心は，正三角形の頂点からおろした垂線上にある。
正方形の一辺の長さを 2b とすると，b = 2r/√3 である。
したがって，正方形の一辺の長さは 4r/√3 = 4√3r/3 である。

include("julia-source.txt");

function draw(r, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   b = 2r/√3
   h = 1.2*(2r + 2b)  # 適当
   a = h/√3
   @printf("円の直径が %g のとき，正方形の一辺の長さは %g である。\n", 2r, 2b)
   plot([a, 0, -a , a], [0, h, 0, 0], color=:green, lw=0.5)
   circle(0, h - r, r)
   rect(-b, h - 2r - 2b, b, h - 2r, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, h - r, "h-r", :red, :center, delta=-delta/2)
       point(0, h - 2r, "h-2r", :red, :center, delta=-delta/2)
       point(0, h - 2r - 2b, "h-2r-2b", :red, :center, delta=-delta/2)
       point(0, h, "h", :red, :center, :bottom, delta=delta/2)
       point(-b, h - 2r - 2b, "(-b,h-2r-2b) ", :blue, :right, :vcenter)
       point(b, h - 2r, " (b,h-2r)", :blue, :left, :vcenter)
   end
end;

draw(2, true)

算額（その2009）

算額（その2009）

(22) 兵庫県養父町左近山 地蔵堂 明治21年(1888)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，長方形，斜線2本

長方形の中に 1 本の対角線ともう 1 本の斜線を引き，等円 4 個を容れる。長方形の長辺と短辺が与えられたとき，等円の直径を求めよ。

長方形の長辺と短辺を a, b
斜線と長方形の長辺の交点座標を (c, 0), (a - c, b)
等円の半径と中心座標を r, (x1, r), (x2, y2), (r, y3), (x3, b - r)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive,
     r::positive, x1::positive, x2::positive,
     y2::positive, y3::positive, x3::positive
eq1 = dist2(0, b, a, 0, x1, r, r)
eq2 = dist2(0, b, a, 0, x2, y2, r)
eq3 = dist2(0, b, a, 0, r, y3, r)
eq4 = dist2(0, b, a, 0, x3, b - r, r)
eq5 = dist2(c, 0, a - c, b, x2, y2, r)
eq6 = dist2(c, 0, a - c, b, x1, r, r)
eq7 = (x2 - r)^2 + (y2 - y3)^2 - 4r^2

using NLsolve

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

function H(u)
   (c, r, x1, x2, y2, y3, x3) = u
   return [
       b*(a^2*b - 2*a^2*r - 2*a*b*x1 + 2*a*r*x1 - b*r^2 + b*x1^2),  # eq1
       a^2*b^2 - 2*a^2*b*y2 - a^2*r^2 + a^2*y2^2 - 2*a*b^2*x2 + 2*a*b*x2*y2 - b^2*r^2 + b^2*x2^2,  # eq2
       a*(a*b^2 - 2*a*b*y3 - a*r^2 + a*y3^2 - 2*b^2*r + 2*b*r*y3),  # eq3
       b*(-2*a*r*x3 - b*r^2 + b*x3^2),  # eq4
       -a^2*r^2 + a^2*y2^2 + 2*a*b*c*y2 - 2*a*b*x2*y2 + 4*a*c*r^2 - 4*a*c*y2^2 + b^2*c^2 - 2*b^2*c*x2 - b^2*r^2 + b^2*x2^2 - 4*b*c^2*y2 + 4*b*c*x2*y2 - 4*c^2*r^2 + 4*c^2*y2^2,  # eq5
       b*(2*a*c*r - 2*a*r*x1 + b*c^2 - 2*b*c*x1 - b*r^2 + b*x1^2 - 4*c^2*r + 4*c*r*x1),  # eq6
       -4*r^2 + (-r + x2)^2 + (y2 - y3)^2,  # eq7
   ]
end;

(a, b) = (10, 5)
iniv = BigFloat[2.6, 1.1, 5.2, 3.1, 2.2, 3.2, 4.8]
res = nls(H, ini=iniv)

   ([2.60785001480277, 1.1239440727077343, 5.238896505102071, 3.1345163522937027, 2.1761341491119186, 3.181420288904903, 4.761103494897929], true)

例えば，長方形の長辺と短辺が 10 寸，5 寸のとき，等円の直径は 2*1.1239440727077343 = 2.2478881454154687 である。

function draw(a, b, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (c, r, x1, x2, y2, y3, x3) = res[1]
   plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5)
   circle(x1, r, r)
   circle(x2, y2, r)
   circle(r, y3, r)
   circle(x3, b - r, r)
   segment(0, b, a, 0, :blue)
   segment(c, 0, a - c, b, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(x1, r, "(x1,r)", :red, :center, delta=-delta/2)
       point(x2, y2, "(x2,y2)", :red, :center, delta=-delta/2)
       point(r, y3, "(r,y3)", :red, :center, delta=-delta/2)
       point(x3, b - r, "(x3,b-r)", :red, :center, delta=-delta/2)
       point(c, 0, "c", :blue, :center, delta=-delta/2)
       point(a, 0, "a", :blue, :center, delta=-delta/2)
       point(0, b, "b ", :blue, :right, :vcenter)
       point(a - c, b, "(a-c,b)", :blue, :center, :bottom, delta=delta/2)
       point(a, b, "(a,b)", :blue, :right, :bottom, delta=delta/2)
       plot!(xlims=(-7delta, a + 5delta), ylims=(-7delta, b + 5delta))
   end
end;

draw(a, b, true)

算額（その2008）

算額（その2008）

(23) 兵庫県姫路市広畑区西蒲田 天満神社 明治21年(1888)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，外円

外円の中に等円 3 個を容れる。外円の直径が 13 寸のとき，等円の直径および（図に示す）矢の長さはいかほどか。

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

include("julia-source.txt");

using SymPy

@syms R::positive, r::positive, y::negative;
eq1 = r^2 + (R - r - y)^2 - 4r^2
eq2 = r^2 + y^2 - (R - r)^2
res = solve([eq1, eq2], (r, y))[1]

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

等円の半径の式は簡約化でき，R*(2√3 - 3) となる。

r = res[1] |> simplify
r |> println

   R*(-3 + 2*sqrt(3))

等円の半径 r は，外円の半径 R の 2√3 - 3 倍である。
外円の直径が 13 寸のとき，等円の直径は 13*(2√3 - 3) = 6.033320996790807 寸である。

等円径を得る「術」は，「外円径/(1.1547 + 1)」である。分母をなぜ 2.1547 としなかったかというのは，分母は 2√3/3 + 1 で，2√3/3 = 1.1547005383792515 を 1.1547 に丸めたことを残したかったためであろう。それにしても，小数の割り算より掛け算の方が計算しやすくはないだろうか。

外円径 = 13
等円径 = 外円径/(1.1547 + 1)

   6.033322504292941

矢は (y + r) - (R - 2r) = R*(7√3 - 12) である。
外円の直径が 13 寸のとき，0.808311744383917 寸である。

y = res[2]
矢 = ((y + r) - (R - 2r)) |> simplify
矢 |> println
矢(R => 13/2).evalf() |> println

   R*(-12 + 7*sqrt(3))
   0.808311744383917

矢を求める「術」は，「等円径 - √3*等円径/2」である。等円径の計算精度が確保できれば完全に正しい式である。

function draw(R, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, y) = (-(-R + R*(-2 + sqrt(3)))*(R*(-2 + sqrt(3)) + R)/(2*R), R*(-2 + sqrt(3)))
   println((y + r) - (R - 2r))
   plot()
   circle(0, 0, R, :blue)
   circle(0, R - r, r)
   circle(r, y, r)
   circle(-r, y, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R, "R", :blue, :center, :bottom, delta=delta/2)
       point(0, R - r, "等円:r,(0,R-r)", :red, :center, delta=-delta/2)
       point(r, y, "等円:r,(r,y)", :red, :center, delta=-delta/2)
       segment(-1.5r, y + r, 1.5r, y + r, :green)
       dimension_line(0, y+ r, 0, R - 2r, " 矢", :black, :left, lw=1)
   end
end;

draw(13/2, true)

算額（その2007）

算額（その2007）

(23) 兵庫県姫路市広畑区西蒲田 天満神社 明治21年(1888)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，直線上

直線上に甲円 2 個，乙円 1 個，丙円 1 個が載っている。甲円，乙円の直径が 10 寸，4 寸のとき，丙円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy

@syms r1::positive, r2::positive, r3::positive
eq1 = (r1 - r2)^2 + (r1 - r3)^2 - (r1 + r3)^2
res = solve(eq1, r3)[1];
res |> println

   (r1 - r2)^2/(4*r1)

甲円，乙円の直径が 10 寸，4 寸のとき，丙円の直径は 0.9 寸（9 分）である。

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

draw(225/2, 81/2, true)

算額（その2006）

算額（その2006）

(25) 兵庫県太子堂町鵤 鵤太子堂 明治26年(1893)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円3個，円弧

円弧（外円）の中に大円　1　個，中円　2　個，小円　2　個を容れる。外円の直径が 225 寸，大円の直径が 81 寸のとき，小円の直径はいかほどか。

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

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive,
     r2::positive, x2::positive,
     r3::positive, x3::positive
eq1 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq2 = x3^2 + (r1 - r3)^2 - (r1 + r3)^2
eq3 = (x2 - x3)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq4 = x2^2 + (R - 2r1 + r2)^2 - (R - r2)^2
res = solve([eq1, eq2, eq3, eq4], (r2, x2, r3, x3))[1]

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

res[3](R => 225/2, r1 => 81/2).evalf() |> println

   8.00000000000000

外円の直径が 225 寸，大円の直径が 81 寸のとき，小円の直径は 16 寸である。

「術」では，小円の直径を求める式は以下のようになっている。

外径 = 225; 大径 = 81
甲 = 外径 - 大径
(sqrt(甲*外径) - 甲)^2/大径

   16.0

甲を代入すると以下のようになる。
SymPy では，自動ではこのような式まで簡約化できない。

@syms 外径, 大径, R, r1
小径 = (sqrt((外径 - 大径)*外径) - (外径 - 大径))^2/大径;
小径(外径 => 225, 大径 => 81) |> println

   16

function draw(R, r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, x2, r3, x3) = (-r1*(-R + r1)/R, 2*r1*sqrt(R - r1)/sqrt(R), -2*R^(3/2)*sqrt(R - r1)/r1 + 2*sqrt(R)*sqrt(R - r1) + 2*R^2/r1 - 3*R + r1, 2*sqrt(R)*sqrt(R - r1) - 2*R + 2*r1)
   @printf("外円，大円の直径が %g，%g のとき，小円の直径は %g である（中円の直径は %g)。\n", 2R, 2r1, 2r3, 2r2)
   @printf("R = %g;  r1 = %g;  r2 = %g;  x2 = %g  r3 = %g;  x3 = %g\n", R, r1, r2, x2, r3, x3)
   y = R - 2r1
   x = sqrt(R^2 - y^2)
   θ = atand(y/x)
   plot()
   circle(0, 0, R, :green, beginangle=θ, endangle=180−θ)
   circle(0, R - r1, r1)
   circle(x2, y + r2, r2, :blue)
   circle(x3, y + r3, r3, :magenta)
   segment(-x, y, x, y, :green)
   segment(0, 0, x, y, :gray90)
   segment(0, 0, -x, y, :gray90)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, R - r1, "大円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
       point(x2, y + r2, "中円:r2\n(x2,y+r2)", :blue, :center, :bottom, delta=delta)
       point(x3, y + r3, " 小円:r3,(x3,y+r3)", :magenta, :left, :vcenter)
       point(0, R, "R", :green, :center, :bottom, delta=delta)
       point(0, y, " y=R-2r1", :green, :center, delta=-delta)
       point(√(R^2-y^2), y, "(√(R^2-y^2),y)", :green, :right, delta=-delta)
   end
end;

draw(225/2, 81/2, true)

算額（その2005）

算額（その2005）

(25) 兵庫県太子堂町鵤 鵤太子堂 明治26年(1893)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：等円4個，正方形，斜線3本

正方形の中に斜線を 3 本引き，区画された領域に等円を 4 個容れる。正方形の一辺の長さが与えられたとき，等円の直径を求めよ。

図形を180度回転させて解く。
正方形の一辺の長さを a
斜線と正方形の一辺の交点座標を (0, b)
等円の半径と中心座標を r, (r, r), (r, y)
とおき，以下の連立方程式を解く。
注：b は正方形の一辺の中点ではない。

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive,
     r::positive, y::positive
eq1 = dist2(a, 0, 0, b, r, r, r)
eq2 = dist2(a, 0, 0, b, r, y, r)
eq3 = dist2(a, 0, 0, a, r, y, r)
res = solve([eq1, eq2, eq3], (r, y, b))[4];

等円の半径 r の式は，以下のように簡約化される。

@syms d
apart(res[1]/a, d)*a |> simplify |> println

   a*(1 - sqrt(-1 + sqrt(2)))/2

等円の半径 r は，正方形の一辺の長さ a の (1 - sqrt(√2 - 1))/2 倍である。
正方形の一辺の長さが 1 のとき，等円の直径は (1 - sqrt(√2 - 1))/2 = 0.17820287354720865 である。

y, b は以下のようになる。

apart(res[2]/a, d)*a |> simplify |> println

   a*(-sqrt(2) + sqrt(-1 + sqrt(2)) + sqrt(-2 + 2*sqrt(2)) + 1)/2

apart(res[3]/a, d)*a |> simplify |> println

   a*sqrt(-2 + 2*sqrt(2))/2

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r = a*(1 - sqrt(√2 - 1))/2
   y = a*(sqrt(√2 - 1) + sqrt(2√2 - 2) + 1 - √2)/2
   b = a*sqrt(2√2 - 2)/2
   @printf("a = %g;  r = %g;  y = %g;  b = %g\n", a, r, y, b)
   plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
   circle(r, r, r)
   circle(a - r, a - r, r)
   circle(r, y, r)
   circle(a - y, a - r, r)
   segment(0, b, a, 0)
   segment(0, a, a, 0)
   segment(a - b, a, a, 0)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, a, "(a,a)", :green, :right, :bottom, delta=delta/2)
       point(r, r, "(r,r)", :red, :center, delta=-delta/2)
       point(r, y, "(r,y)", :red, :center, delta=-delta/2)
       point(0, a, "a ", :green, :right, :vcenter)
       point(0, b, "b ", :green, :right, :vcenter)
       point(a, 0, " a", :green, :center, delta=-delta/2)
       plot!(xlims=(-5delta, a + 5delta), ylims=(-5delta, a + 5delta))
   end
end;

draw(1, true)

 

算額（その2004）

算額（その2004）

(復元) 兵庫県伊丹市伊丹桜崎 猪名野神社 嘉永6年(1853)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円3個，直角三角形，中鈎

直角三角形の中に大円を容れ，（鈎と股が挟む）直角の頂点から対辺（弦）へ垂線（中鈎）をおろす。できた小さな直角三角形の中に中円と小円を容れる。鈎，股，弦が 21 寸，28 寸，35 寸，大円の直径が 14 寸のとき，中円，小円の直径はいかほどか。

類題は多いが，素直な「問」である。

元の直角三角形と，中鈎で区分けされた2つの直角三角形は相似であり，その中にある円もまた相似である。相似比は 1 : 28/35 : 21/35 である。
大円の直径が 14 寸なので，中円，小円の直径は 14*28/35 = 11.2 寸，14*21/35 = 8.4 寸である。

include("julia-source.txt");

using SymPy

@syms 鈎::positive, 股::positive, 弦::positive,
     中鈎::positive, 短弦::positive, 長弦::positive,
     r1::positive, r2::positive, x2::positive,
     r3::positive, y3::positive, x0::positive, y0::positive
eq1 = 鈎 + 股 - 弦 - 2r1
eq2 = 弦 - sqrt(鈎^2 + 股^2)
res1 = solve([eq1, eq2], (弦, r1))
弦 = res1[弦]
r1 = res1[r1]
弦 |> println
弦(鈎 => 21, 股 => 28) |> println
r1 |> println
r1(鈎 => 21, 股 => 28) |> println

   sqrt(股^2 + 鈎^2)
   35

   股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2
   7

eq4 = y0/(股 - x0) - 鈎/股
eq5 = (鈎 - y0)/x0 - 鈎/股
eq6 = 中鈎 - sqrt(x0^2 + y0^2)
eq7 = 弦*中鈎 - 股*鈎
(中鈎, x0, y0) = solve([eq4, eq6, eq7], (中鈎, x0, y0))[1];

中鈎 |> println
中鈎(鈎 => 21, 股 => 28) |> println
x0 |> println
x0(鈎 => 21, 股 => 28) |> println
y0 |> println
y0(鈎 => 21, 股 => 28) |> println

   股*鈎/sqrt(股^2 + 鈎^2)
   84/5

   股*鈎^2/(股^2 + 鈎^2)
   252/25

   股^2*鈎/(股^2 + 鈎^2)
   336/25

eq8 = 短弦 - sqrt(鈎^2 - 中鈎^2)
eq9 = 長弦 + 短弦 - 弦
res3 = solve([eq8, eq9], (短弦, 長弦))
短弦 = res3[短弦]
短弦 |> println
短弦(鈎 => 21, 股 => 28) |> println
長弦 = res3[長弦]
長弦 |> println
長弦(鈎 => 21, 股 => 28) |> println

   鈎^2/sqrt(股^2 + 鈎^2)
   63/5

   股^2/sqrt(股^2 + 鈎^2)
   112/5

eq10 = dist2(股, 0, 0, 鈎, r3, y3, r3)
eq11 = dist2(股, 0, 0, 鈎, x2, r2, r2)
eq12 = 中鈎 + 短弦 - 鈎 - 2r3
eq13 = 中鈎 + 長弦 - 股 - 2r2
solve([eq10, eq11, eq12, eq13], (r2, x2, r3, y3))[1];

(r2, x2, r3, y3) = solve([eq10, eq11, eq12, eq13], (r2, x2, r3, y3))[1]
r2 |> println
r2(鈎 => 21, 股 => 28).evalf() |> println
x2 |> println
x2(鈎 => 21, 股 => 28).evalf() |> println
r3 |> println
r3(鈎 => 21, 股 => 28).evalf() |> println
y3 |> println
y3(鈎 => 21, 股 => 28).evalf() |> println

   股*(股 + 鈎 - sqrt(股^2 + 鈎^2))/(2*sqrt(股^2 + 鈎^2))
   5.60000000000000

   -股*sqrt(-2*鈎*(股^2 + 鈎^2)*(2*股^3 + 股^2*鈎 - 2*股^2*sqrt(股^2 + 鈎^2) + 股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 + 鈎^2*sqrt(股^2 + 鈎^2)) + (股^3 + 2*股^2*鈎 - 股^2*sqrt(股^2 + 鈎^2) + 股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 2*鈎^3)^2)/(2*鈎*(股^2 + 鈎^2)) + 股*(股^3 + 2*股^2*鈎 - 股^2*sqrt(股^2 + 鈎^2) + 股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 2*鈎^3)/(2*鈎*(股^2 + 鈎^2))
   11.2000000000000

   鈎*(股 + 鈎 - sqrt(股^2 + 鈎^2))/(2*sqrt(股^2 + 鈎^2))
   4.20000000000000

   -鈎*sqrt(-2*股*(股^2 + 鈎^2)*(股^3 + 股^2*鈎 + 股^2*sqrt(股^2 + 鈎^2) + 股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 2*鈎^3 - 2*鈎^2*sqrt(股^2 + 鈎^2)) + (2*股^3 + 股^2*鈎 + 2*股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 - 鈎^2*sqrt(股^2 + 鈎^2))^2)/(2*股*(股^2 + 鈎^2)) + 鈎*(2*股^3 + 股^2*鈎 + 2*股*鈎^2 - 股*鈎*sqrt(股^2 + 鈎^2) + 鈎^3 - 鈎^2*sqrt(股^2 + 鈎^2))/(2*股*(股^2 + 鈎^2))
   12.6000000000000

function draw(鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   t = 鈎^2 + 股^2
   弦 = sqrt(t)
   r1 = 股/2 + 鈎/2 - 弦/2
   中鈎 = 股*鈎/弦
   x0 = 股*鈎^2/t
   y0 = 股^2*鈎/t
   短弦 = 鈎^2/弦
   長弦 = 股^2/弦
   r2 = 股*(股 + 鈎 - 弦)/(2弦)
   x2 = -股*sqrt(-2*鈎*t*(2股^3 + 股^2*鈎 - 2股^2*弦 + 股*鈎^2 - 股*鈎*弦 + 鈎^3 + 鈎^2*弦) + (股^3 + 2股^2*鈎 - 股^2*弦 + 股*鈎^2 - 股*鈎*弦 + 2鈎^3)^2)/(2鈎*t) + 股*(股^3 + 2股^2*鈎 - 股^2*弦 + 股*鈎^2 - 股*鈎*弦 + 2鈎^3)/(2鈎*t)
   r3 = 鈎*(股 + 鈎 - 弦)/(2弦)
   y3 = -鈎*sqrt(-2股*t*(股^3 + 股^2*鈎 + 股^2*弦 + 股*鈎^2 - 股*鈎*弦 + 2鈎^3 - 2鈎^2*弦) + (2股^3 + 股^2*鈎 + 2股*鈎^2 - 股*鈎*弦 + 鈎^3 - 鈎^2*弦)^2)/(2股*t) + 鈎*(2股^3 + 股^2*鈎 + 2股*鈎^2 - 股*鈎*弦 + 鈎^3 - 鈎^2*弦)/(2股*t)
   plot([0, 股, 0, 0, x0], [0, 0, 鈎, 0, y0], color=:green, lw=0.5)
   circle(r1, r1, r1, :magenta)
   circle(x2, r2, r2)
   circle(r3, y3, r3, :blue)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(股, 0, " 股", :green, :left, :bottom, delta=delta/2)
       point(0, 鈎, " 鈎", :green, :left, :bottom, delta=delta/2)
       point(x0, y0, " (x0,y0)", :green, :left, :bottom, delta=delta/2)
       point(r1, r1, " 大円:r1\n (r1,r1)", :magenta, :left, :vcenter)
       point(x2, r2, "中円:r2\n(x2,r2)", :red, :center, delta=-delta/2)
       point(r3, y3, "小円:r3\n(r3,y3)", :blue, :center, delta=-delta/2)
   end
end;

draw(21, 28, true)

算額（その2003）

算額（その2003）

(補) 滋賀県栗東町字御園 春日神社 文政7年(1824)
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：円4個，長方形，斜線2本

長方形の中に，対角線 1 本と斜線 1 本を引き，区分された領域に等円を 4 個容れる。長方形の長辺，短辺がそれぞれ 72寸，65 寸のとき等円の直径はいかほどか。

長方形の長辺，短辺を 2a, 2b
斜線と長方形の短辺の交点座標を (a, c), (-a, -c)
等円の半径と中心座標を r, (x0, b - r), (-x0, r - b), (a - r, y0), (-y0, r - a)
とおき，以下の連立方程式を解く。
注：図形は，x 軸，y 軸で対称ではない。（等円の中心は x 軸，y 軸上にはない）

include("julia-source.txt");

using SymPy

@syms a::positive, b::positive, c::positive, r::positive, x0, y0
eq1 = dist2(-a, -c, a, c, x0, b - r, r)
eq2 = dist2(-a, -c, a, c, a - r, y0, r)/a
eq3 = dist2(-a, b, a, -b, x0, b - r, r)/b
eq4 = dist2(-a, b, a, -b, a - r, y0, r)/a;

using NLsolve

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

function H(u)
   (r, c, x0, y0) = u
   return [
       a^2*b^2 - 2*a^2*b*r - 2*a*b*c*x0 + 2*a*c*r*x0 - c^2*r^2 + c^2*x0^2,  # eq1
       a*c^2 - 2*a*c*y0 - a*r^2 + a*y0^2 - 2*c^2*r + 2*c*r*y0,  # eq2
       a^2*b - 2*a^2*r + 2*a*b*x0 - 2*a*r*x0 - b*r^2 + b*x0^2,  # eq3
       a*b^2 + 2*a*b*y0 - a*r^2 + a*y0^2 - 2*b^2*r - 2*b*r*y0,  # eq4
   ]
end;

a = 72//2
b = 65//2

iniv = BigFloat[14, 30, 0.1, 0.1]
res = nls(H, ini=iniv)

   ([14.098029747876524, 28.125, 0.6548773444789615, -0.7794330672778218], true)

等円の直径は 2*14.098029747876524 = 28.196059495753047 である。

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r, c, x0, y0) = res[1]
   @printf("a = %g;  b = %g;  r = %g;  c = %g;  x0 = %g;  y0 = %g\n", a, b, r, c, x0, y0)
   plot([a, a, -a, -a, a], [-b, b, b, -b, -b], color=:green, lw=0.5)
   circle(a - r, y0, r, :blue)
   circle(r - a, -y0, r, :blue)
   circle(x0, b - r, r, :blue)
   circle(-x0, r - b, r, :blue)
   segment(-a, -c, a, c, :red)
   segment(-a, b, a, -b, :red)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(a, 0, " a", :green, :left, :bottom, delta=delta/2)
       point(0, b, "b", :green, :center, :bottom, delta=delta/2)
       point(a - r, y0, "(a-r,y0)", :blue, :center, delta=-delta/2)
       point(r - a, -y0, "(r-a,-y0)", :blue, :center, delta=-delta/2)
       point(x0, b - r, "(x0,b-r)", :blue, :center, delta=-delta/2)
       point(-x0, r - b, "(-x0,r-b)", :blue, :center, delta=-delta/2)
   end
end;

draw(72/2, 65/2, true)