算額（その772）

算額（その772）

宮城県白石市小原 小原温泉薬師堂 大正5年(1916)
徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf

外円の中に大円，中円，小円，甲円，乙円が入っている。中円，小円，甲円の直径がそれぞれ 9 寸，6 寸，3 寸のとき，乙円の直径はいかほどか。

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

算額にある図とは大幅に異なる結果になる。そもそも「大円」は「中円」はおろか，「小円」よりも小さい。

include("julia-source.txt");

using SymPy
@syms R, r1, r2, x2, y2, r3, x3, y3, r4, x4, y4, r5, x5, y5
eq1 = x2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
eq2 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2
eq3 = x4^2 + (R - r1 - y4)^2 - (r1 + r4)^2
eq4 = x5^2 + (R - r1 - y5)^2 - (r1 + r5)^2
eq5 = x2^2 + y2^2 - (R - r2)^2
eq6 = x3^2 + y3^2 - (R - r3)^2
eq7 = x4^2 + y4^2 - (R - r4)^2
eq8 = x5^2 + y5^2 - (R - r5)^2
eq9 = (x2 - x3)^2 + (y2 - y3)^2 - (r2 + r3)^2
eq10 = (x2 - x4)^2 + (y2 - y4)^2 - (r2 + r4)^2
eq11 = (x3 - x5)^2 + (y3 - y5)^2 - (r3 + r5)^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)
   (R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = u
   return [
       x2^2 - (r1 + r2)^2 + (R - r1 - y2)^2,  # eq1
       x3^2 - (r1 + r3)^2 + (R - r1 - y3)^2,  # eq2
       x4^2 - (r1 + r4)^2 + (R - r1 - y4)^2,  # eq3
       x5^2 - (r1 + r5)^2 + (R - r1 - y5)^2,  # eq4
       x2^2 + y2^2 - (R - r2)^2,  # eq5
       x3^2 + y3^2 - (R - r3)^2,  # eq6
       x4^2 + y4^2 - (R - r4)^2,  # eq7
       x5^2 + y5^2 - (R - r5)^2,  # eq8
       -(r2 + r3)^2 + (x2 - x3)^2 + (y2 - y3)^2,  # eq9
       -(r2 + r4)^2 + (x2 - x4)^2 + (y2 - y4)^2,  # eq10
       -(r3 + r5)^2 + (x3 - x5)^2 + (y3 - y5)^2,  # eq11
   ]
end;

(r2, r3, r4) = (9, 6, 3) .// 2
iniv = BigFloat[7.5, 2.43, 2.7, -1.31, -4.2, 1.62, 3.9, 4.56, 1.0, -3.4, 5.54]
res = nls(H, ini=iniv)

   ([7.5031904642363045, 2.4327679290250375, 2.7, -1.3149726097831362, -4.2, 1.6244150815566774, 3.9, 4.563802772896491, 1.0, -3.4, 5.543598670009763], true)

乙円の半径は 1 寸である（直径は 2 寸）。

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

   R = 7.50319;  r1 = 2.43277;  x2 = 2.7;  y2 = -1.31497;  x3 = -4.2;  y3 = 1.62442;  x4 = 3.9;  y4 = 4.5638;  r5 = 1;  x5 = -3.4;  y5 =5.5436

function draw(more=false)
   pyplot(size=(400, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r2, r3, r4) = (9, 6, 3) .// 2
   (R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5) = res[1]
   @printf("r2 = %g;  r3 = %g;  r4 = %g のとき，乙円の直径は %g\n", r2, r3, r4, 2r5)
   @printf("R = %g;  r1 = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g;  x4 = %g;  y4 = %g;  r5 = %g;  x5 = %g;  y5 =%g\n", R, r1, x2, y2, x3, y3, x4, y4, r5, x5, y5)
   plot()
   circle(0, 0, R)
   circle(0, R - r1, r1, :blue)
   circle(x2, y2, r2, :green)
   circle(x3, y3, r3, :magenta)
   circle(x4, y4, r4, :purple)
   circle(x5, y5, r5, :orange)
   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\n(0,R-r1)", :blue, :center, :bottom, delta=delta)
       point(x2, y2, "中円:r2\n(x2,y2)", :green, :center, :bottom, delta=delta)
       point(x3, y3, "小円:r3\n(x3,y3)", :magenta, :center, :bottom, delta=delta)
       point(x4, y4, "甲円:r4,(x4,y4)", :purple, :left, delta=-delta)
       point(x5, y5, "乙円:r5,(x5,y5)", :black, :right, delta=-delta)
       point(R, 0, " R", :red, :left, :bottom, delta=delta)
   end
end;

算額（その771）

算額（その771）

神壁算法 天明5年乙巳 關流四傅藤田貞資六弟子 筑後州久留米 佐田甚助方之
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

和算問題あれこれ2 令和６年４月の問題-No.３
https://gunmawasan.web.fc2.com/k-n-mondai.html

円弧の弦の上に，底辺が弦の上で頂点が円弧上にある甲，乙，丙の正三角形が載っている。一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さはいかほどか。

正三角形の底辺の 2 つの x 座標，頂点の x 座標を，左から（丙，乙，甲の順で），a, b, c, d, e, f, g とする。
弦と y 軸の交点の y 座標を h とする。
頂点の y 座標 は h + √3(b - a)，h + √3(d - c)，h + √3(f - e) とする。
A = g - e，B = e - c，C = c - a とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     c::positive, d::positive,
     e::positive, f::positive,
     g::negative, h::positive,
     R::positive
@syms a, b, c, d, e, f, g, h, R, A, B, C
s3 = sqrt(Sym(3))
eq1 = g - e - A
eq2 = e - c - B
eq3 = c - a - C
eq4 = b - (a + c)/2
eq5 = d - (c + e)/2
eq6 = f - (e + g)/2
eq7 = (h + (b - a)s3)^2 + b^2 - R^2
eq8 = (h + (d - c)s3)^2 + d^2 - R^2
eq9 = (h + (f - e)s3)^2 + f^2 - R^2
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9],
   (R, h, a, b, c, d, e, f, g))

   2-element Vector{NTuple{9, Sym{PyCall.PyObject}}}:
    (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
    (sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))

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

弦の長さは 2sqrt(R^2 - h^2) である。
正三角形の一辺の長さがそれぞれ 206 寸 1 分，229 寸，183 寸 2 分のとき，弦の長さは 1029.80004005632 寸である。

2sqrt(res[2][1]^2 - res[2][2]^2)(A => 206.1, B => 229.0, C => 183.2) |> println

   1029.80004005632

その他のパラメータは以下のとおりである。
R = 764.582;  h = 565.211;  a = -337.775;  b = -246.175;  c = -154.575;  d = -40.075;  e = 74.425;  f = 177.475;  g = 280.525

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (A, B, C) = (206.1, 229.0, 183.2)
   (R, h, a, b, c, d, e, f, g) = (-sqrt(3)*sqrt((A^2 - A*B + B^2)*(B^2 - B*C + C^2)*(A^2 + A*B - A*C + B^2 + B*C + C^2))/(3*(A*C - B^2)), sqrt(3)*(-A^2*B - A^2*C - A*B*C - A*C^2 + B^3 - B*C^2)/(6*(A*C - B^2)), (A^2*B - A^2*C - A*B*C - A*C^2 + B^3 + 2*B^2*C - B*C^2)/(2*(A*C - B^2)), (A^2*B - A^2*C - A*B*C + B^3 + B^2*C - B*C^2)/(2*(A*C - B^2)), (B - C)*(A^2 - A*C + B^2 + B*C)/(2*(A*C - B^2)), (A - C)*(A*B - A*C + B*C)/(2*(A*C - B^2)), (A - B)*(A*B - A*C + B^2 + C^2)/(2*(A*C - B^2)), (A^2*B - A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)), (A^2*B + A^2*C - 2*A*B^2 + A*B*C + A*C^2 - B^3 - B*C^2)/(2*(A*C - B^2)))
   @printf("R = %g;  h = %g;  a = %g;  b = %g;  c = %g;  d = %g;  e = %g;  f = %g;  g = %g\n", R, h, a, b, c, d, e, f, g)
   s3 = √3
   plot([a, b, c, d, e, f, g], h .+ [0, (b - a)s3, 0, (d - c)s3, 0, (f - e)s3, 0], color=:blue, lw=0.5)
   θ = asind(h/R)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ, n=1000)
   x = R*cosd(θ)
   segment(-x, h, x, h, :red)
   plot!([-x, 0, x], [h, 0, h], color=:gray50, lw=0.5)
   @printf("弦の長さ = %g\n", 2R*cosd(θ))
   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, b, c, d, e, f, g], h, Char["abcdefg"...], :blue, :center, delta=-delta)
       point(x, h, "(x,h) ", :red, :right, delta=-delta)
       point(x/2+4delta, h/2, "R", :gray50, :center, :vcenter, mark=false)
       point.([b, d, f], h + 75, ["丙三角", "乙三角", "甲三角"], :blue, :center, delta=-delta, mark=false)
   end
end;

算額（その770）

算額（その770）

神壁算法 天明8年戊申 二月 關流藤田貞資門人 細川能登守家士 渡邉金八源耀
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

正五角形に5本の斜線を入れ，内部に小さな正五角形を作り，できた 6 つの区画にそれぞれ同じ直径を持つ円を内接させる。
正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径はいかほどか。

正五角形の一辺の長さを L，正五角形が内接する円の半径を R とする。
内側の正五角形が内接する円の半径を c，その一辺の長さを l とする。
等円の半径を r1 とする。
斜線と正五角形の一辺で作られる三角形 BAR において第二余弦定理を使う。
AR, AB, BR の長さを AR, AB, BR，また，A, B, R の座標を (x, y), (x0, y0), (0, R) とする。
ポイントは，斜線は外側の五角形の一辺と平行ではないということである。

include("julia-source.txt");

using SymPy
@syms L::positive, R::positive,
     x0::positive, y0::positive,
     x1::positive, y1::positive,
     x::negative, y::positive,
     AR::positive, l::positive, r1::positive
d18 = Sym(18)

R = L/2sind(2d18)
x0 = R*cosd(d18)
y0 = R*sind(d18)
c = sqrt(x^2 + y^2)

eq1 = (x0 - x)^2 + (y0 - y)^2 - (AR + l)^2
eq2 = x^2 + (R - y)^2 - AR^2
eq3 = AR^2 + (AR + l)^2 - 2AR*(AR + l)*cosd(4d18) - L^2
eq4 = dist(x0, y0, R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2
eq5 = dist(x0, y0, x, y, x1, y1) - r1^2
eq6 = dist(x0, y0, x, y, 0, 0) - r1^2
eq7 = dist(c*cosd(atand(y, x) - 4d18), c*sind(atand(y, x) - 4d18), R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^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(par)
    (AR, l, x, y, x1, y1, r1) = par
    d18 = 18
    R = L/2sind(2d18)
    x0 = R*cosd(d18)
    y0 = R*sind(d18)
    c = sqrt(x^2 + y^2)
    return [
        (x0 - x)^2 + (y0 - y)^2 - (AR + l)^2,
        x^2 + (R - y)^2 - AR^2,
        AR^2 + (AR + l)^2 - 2AR*(AR + l)*cosd(4d18) - L^2,
        dist(x0, y0, R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2,
        dist(x0, y0, x, y, x1, y1) - r1^2,
        dist(x0, y0, x, y, 0, 0) - r1^2,
        dist(c*cosd(atand(y, x) - 4d18), c*sind(atand(y, x) - 4d18), R*cosd(-3d18), R*sind(-3d18), x1, y1) - r1^2
    ]
end;

L = 5336//10
iniv = BigFloat[347.1, 179.5, -94.4, 119.9, 259.5, 10.0, 123.5]
res = nls(H, ini=iniv)

   ([347.0551765561487, 179.4560303257959, -94.42748155686841, 119.94507841468496, 259.51771893633804, 10.023199995813941, 123.50001782997607], true)

正五角形の一辺の長さが 533 寸 6 分のとき，等円の直径は 247 寸である。

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

AR = 347.055;  x = -94.4275;  y = 119.945;  x1 = 259.518;  y1 = 10.0232;  r1 = 123.5

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   L = 5336//10
   R = L/2sind(2*18)
   x0 = R*cosd(18)
   y0 = R*sind(18)
   (AR, l, x, y, x1, y1, r1) = res[1]
   @printf("等円の直径 = %g\n", 2r1)
   @printf("AR = %g;  x = %g;  y = %g;  x1 = %g;  y1 = %g;  r1 = %g\n", AR, x, y, x1, y1, r1)
   r = sqrt(x^2 + y^2)
   θ = atand(y, x)
   ix = [r*cosd(72i + θ) for i in 0:5]
   iy = [r*sind(72i + θ) for i in 0:5]
   ox = [R*cosd(72i + 18) for i in 0:5]
   oy = [R*sind(72i + 18) for i in 0:5]
   plot()
   circle(0, 0, R, :gray80)
   circle(0, 0, r, :gray90)
   circle(0, 0, r1)
   rotate(x1, y1, r1, angle=72)
   for i in 1:5
       segment(ox[i], oy[i], ix[i], iy[i], :blue)
       segment(ox[i], oy[i], ox[i + 1], oy[i + 1], :blue)
   end
   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(x, y, "A:(x,y) ", :blue, :right, :vcenter)
       point(x0, y0, " B:(x0,y0)", :blue, :left, :vcenter)
       point(x1, y1, "等円:r1\n(x1,y1)", :red, :center, :bottom, delta=delta)
       point(0, R, " R:(0,R)", :black, :left, :bottom, delta=delta)
       point(r, 0, " r", :gray50, :left, :bottom, delta=delta)
   end
end;

算額（その769）

算額（その769）

神壁算法 天明5年乙巳 秋九月 第二術 關流四傅藤田権平定資六弟子 奧州一關 梶山主水次俊
藤田貞資(1789)：神壁算法巻上
http://www.wasan.jp/jinpeki/jinpekisanpo1.pdf

大円，小円，甲円，乙円がある。上下の二線（上線，下線）と斜線を引く。
記述が重複するが正確を期すため，これらは以下のように接する。
大円は二線，斜線，小円，甲円，乙円に接する。
小円は二線，大円に接する。
甲円は上線，斜線，大円に接する。
乙円は下線，斜線，大円に接する。
大円，小円，甲円の直径がそれぞれ 245寸，196寸，25寸のとき，乙円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, r1)
小円の半径と中心座標を r2, (0, r2)
甲円の半径と中心座標を r3, (x3, y3)
乙円の半径と中心座標を r4, (x4, r4)
上線は二点 (x02, y02), (x03, y03) を通る
下線は二点 (x02, y02), (x01, 0) を通る
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, x1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     r4::positive, x4::positive,
     x01::positive, x02::positive, y02::positive,
     x03::positive, y03::positive
eq1 = dist(x02, y02, x03, y03, x1, r1) - r1^2
eq2 = dist(x02, y02, x03, y03, 0, r2) - r2^2
eq3 = dist(x02, y02, x03, y03, x3, y3) - r3^2
eq4 = dist(x01, 0, x02, y02, x1, r1) - r1^2
eq5 = dist(x01, 0, x02, y02, x3, y3) - r3^2
eq6 = dist(x01, 0, x02, y02, x4, r4) - r4^2
eq7 = sqrt((x02 - x03)^2 + (y02 - y03)^2) - x1
eq8 = x1^2 + (r1 - r2)^2 - (r1 + r2)^2
eq9 = (x1 - x3)^2 + (y3 - r1)^2 - (r1 + r3)^2
eq10 = (x1 - x4)^2 + (r1 - r4)^2 - (r1 + r4)^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)
   (x1, x3, y3, r4, x4, x01, x02, y02, x03, y03) = u
   return [
       -r1^2 + (r1 - y02 - (-y02 + y03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-x02 + x1 - (-x02 + x03)*((r1 - y02)*(-y02 + y03) + (-x02 + x03)*(-x02 + x1))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq1
       -r2^2 + (-x02 - (-x02 + x03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (r2 - y02 - (-y02 + y03)*(-x02*(-x02 + x03) + (r2 - y02)*(-y02 + y03))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq2
       -r3^2 + (-x02 + x3 - (-x02 + x03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2 + (-y02 + y3 - (-y02 + y03)*((-x02 + x03)*(-x02 + x3) + (-y02 + y03)*(-y02 + y3))/((-x02 + x03)^2 + (-y02 + y03)^2))^2,  # eq3
       -r1^2 + (r1 - y02*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x1 - (-x01 + x02)*(r1*y02 + (-x01 + x02)*(-x01 + x1))/(y02^2 + (-x01 + x02)^2))^2,  # eq4
       -r3^2 + (-y02*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2) + y3)^2 + (-x01 + x3 - (-x01 + x02)*(y02*y3 + (-x01 + x02)*(-x01 + x3))/(y02^2 + (-x01 + x02)^2))^2,  # eq5
       -r4^2 + (r4 - y02*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2 + (-x01 + x4 - (-x01 + x02)*(r4*y02 + (-x01 + x02)*(-x01 + x4))/(y02^2 + (-x01 + x02)^2))^2,  # eq6
       -x1 + sqrt((x02 - x03)^2 + (y02 - y03)^2),  # eq7
       x1^2 + (r1 - r2)^2 - (r1 + r2)^2,  # eq8
       (-r1 + y3)^2 - (r1 + r3)^2 + (x1 - x3)^2,  # eq9
       (r1 - r4)^2 - (r1 + r4)^2 + (x1 - x4)^2,  # eq10
   ]
end;

(r1, r2, r3) = (245, 196, 25) .// 2
iniv = BigFloat[219.13, 118.51, 212.5, 24.5, 109.57, 82.18, 107.08, 222.73, -106.65, 174.33]
res = nls(H, ini=iniv)

   ([219.1346617949794, 118.51160280748886, 212.5, 24.5, 109.5673308974897, 82.17549817311728, 107.0771642861831, 222.72727272727272, -106.6467651187968, 174.33221099887766], true)

乙円の半径は 24.5寸 である（直径 49寸）。

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

r1 = 122.5;  r2 = 98;  r3 = 12.5;  x1 = 219.135;  x3 = 118.512;  y3 = 212.5;  r4 = 24.5;  x4 = 109.567;  x01 = 82.1755;  x02 = 107.077;  y02 = 222.727;  x03 = -106.647;  y03 = 174.332

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2, r3) = (245, 196, 25) .// 2
   (x1, x3, y3, r4, x4, x01, x02, y02, x03, y03) = res[1]
   @printf("大円，小円，甲円の直径がそれぞれ %g寸, %g寸, %g寸 のとき，乙円の直径は %g寸\n", 2r1, 2r2, 2r3, 2r4)
   @printf("r1 = %g;  r2 = %g;  r3 = %g;  x1 = %g;  x3 = %g;  y3 = %g;  r4 = %g;  x4 = %g;  x01 = %g;  x02 = %g;  y02 = %g;  x03 = %g;  y03 = %g\n", r1, r2, r3, x1, x3, y3, r4, x4, x01, x02, y02, x03, y03)
   plot()
   circle(x1, r1, r1)
   circle(0, r2, r2, :blue)
   circle(x3, y3, r3, :green)
   circle(x4, r4, r4, :magenta)
   segment(x01, 0, x02, y02)
   y04 = (y02 - y03)/(x02 - x03)*(1.25x1 - x03) + y03
   segment(x03, y03, 1.25x1, y04)
   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(x01, 0, "x01 ", :black, :right, :bottom, delta=delta)
       point(x02, y02, "(x02,y02)", :black, :right, :bottom, delta=delta)
       point(x03, y03, "(x03,y03)", :black, :left, :bottom, delta=4delta)
       point(x1, r1, " 大円:r1,(x1,r1)", :red, :center, delta=-delta)
       point(0, r2, " 小円:r2,(0,r2)", :blue, :center, delta=-delta)
       point(x3, y3, "   甲円:r3,(x3,y3)", :green, :left, :vcenter)
       point(x4, r4, "      乙円:r4,(x4,r4)", :magenta, :left, :vcenter)
       segment(x03, 0, 1.4x1, 0)
   end
end;

算額（その768）

算額（その768）

宮城県白石市小原字馬頭山 三瀧神社奉納算額（萬蔵稲荷神社所蔵） 明治8年
徳竹亜紀子，谷垣美保：宮城県白石市小原地区の算額調査，仙台高等専門学校名取キャンパス研究紀要，第57号，2021.
https://www.sendai-nct.ac.jp/natori-library/wp/wp-content/uploads/2021/04/No57_2.pdf

大球の頂きに甲球が載り，大球と甲球の両方に接し，かつ互いにも接し合う乙球を複数個添える。大球，甲球の直径がそれぞれ 39 寸，13 寸のとき，乙球が6個ある場合に，乙球の直径はいかほどか。

大球の半径と中心座標を r1, (0, 0, 0)
甲球の半径と中心座標を r2, (0, 0, r1 + r2)
x-z 平面にある乙球の個数を n，半径と中心座標を r3, (x3, 0, z3); x3 > 0
とおおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, r3::positive,
     x3::positive, z3::positive, n::integer
(r1, r2) = (39, 13) .// 2
eq1 = x3^2 + z3^2 - (r1 + r3)^2
eq2 = x3^2 + (r1 + r2 - z3)^2 - (r2 + r3)^2
eq3 = x3*sind(Sym(180)/n) - r3
res = solve([eq1, eq2, eq3], (r3, x3, z3))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (-78*sin(pi/n)^2/(3*sin(pi/n)^2 - 4), -78*sin(pi/n)/(3*sin(pi/n)^2 - 4), (39*sin(pi/n)^2 - 156)/(6*sin(pi/n)^2 - 8))

t = sin(pi/n)，u = 4 - 3t^2 とすれば，乙球の半径と中心座標は (78t^2/u, 78t/u, (312 - 78t^2)/u) である。

問のように，大球，甲球の直径がそれぞれ 39 寸，13 寸，乙球が 6 個のときには，乙球の直径は 12 寸である。

function draw(which=1, more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   n = 9
   (r1, r2) = (39, 13) .// 2
   t = sin(pi/n)
   u = 4 - 3t^2
   (r3, x3, z3) = (78t^2/u, 78t/u, (312 - 78t^2)/4u)
   plot()
   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)
   end
   if which == 1
       circle(0, 0, r1)
       circle(0, r1+r2, r2, :blue)
       circle(x3, z3, r3, :green)
       point(0, r1, "r1", :red, :center, delta=-0.5)
       point(0, 0, "大球:r1,(0,0,0)", :red, :center, delta=-0.5)
       point(0, r1 + r2, "甲球:r2\n(0,0,r1+r2)", :blue, :center, :bottom, delta=0.5)
       point(x3, z3, "    乙球:r3\n    (x3,0,z3)", :green, :left, :vcenter)
       point(r1, 0, " x軸", :black, :left, delta=-0.25, mark=false) 
       point(0, r1 + 2r2, "z軸", :black, :center, :bottom, delta=0.2, mark=false) 
   else
       circle(0, 0, x3, :gray80)
       rotate(x3, 0, x3*sind(180/n), angle=360/n, :green)
       point(x3, 0, "乙球:r3\n(x3,0,z3)", :green, :center, :bottom, delta=0.25)
       point(x3 + r3, 0, "x軸 ", :black, :right, delta=-0.25, mark=false) 
       point(0, x3 + r3, "y軸", :black, :center, :bottom, mark=false) 
   end
end;

 

算額（その767）

算額（その767）

宮城県大和町落合舞野 舞野正観音堂 慶応4年(1868)
http://www.wasan.jp/miyagi/mainokannondo.html

正方形の中に交差する楕円が入っており，区画された領域に4個の等円が入っている。正方形の一辺の長さが与えられたとき，等円の直径を求めよ。

正方形の一辺の長さを　2a とする。
楕円の長半径と短半径，中心座標を a, b, (0, 0)
等円の半径と中心座標を r, (a - r, a - r), (0, 0)
とする。楕円の長半径は正方形の一辺の長さの半分，等円の半径と楕円の短半径は等しい。
以下の連立方程式を解く。
しかし，たかが3元連立方程式なのに SymPy では答えを得られない。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r::positive,
     x0::positive, y0::positive
@syms a, b, r, x0, y0
a = 1//2
r = b
eq1 = x0^2/a^2 + y0^2/b^2 - 1
eq2 = -b^2*x0/(a^2*y0) + (a - r - x0)/(a - r - y0)
eq3 = (a - r - x0)^2 + (a - r - y0)^2 - r^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)
   (b, x0, y0) = u
   return [
       4*x0^2 - 1 + y0^2/b^2,  # eq1
       -4*b^2*x0/y0 + (-b - x0 + 1/2)/(-b - y0 + 1/2),  # eq2
       -b^2 + (-b - x0 + 1/2)^2 + (-b - y0 + 1/2)^2,  # eq3
   ]
end;
a = 1/2
iniv = BigFloat[0.356, 0.57, 0.295]./2
res = nls(H, ini=iniv)

   ([0.1785016026524881, 0.27969014973090045, 0.14796196722191884], true)

等円の直径は 正方形の一辺の長さの 0.3570032053049762 倍である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   a = 1/2
   (b, x0, y0) = [0.1785016026524881, 0.27969014973090045, 0.14796196722191884]
   r = b
   println("等円の直径は 正方形の一辺の長さの $(2r) 倍である。")
   plot([a, a, -a, -a, a], [-a, a, a, -a, -a], color=:black, lw=0.5)
   ellipse(0, 0, a, b, color=:blue)
   ellipse(0, 0, b, a, color=:blue)
   circle(0, 0, r)
   circle4(a - r, a - r, r)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, b, "r = b", :blue, :center, :bottom, delta=delta)
       point(0, 0, "等円:r,(0,0)", :red, :center, delta=-delta)
       point(a - r, a - r, "等円:r,(a-r,a-r)", :red, :center, delta=-delta)
       point(a, 0, "a ", :blue, :right, :bottom, delta=delta)
   end
end;

算額（その766）

算額（その766）

山形県七日町 長源寺観音堂 大正4年(1915)
http://www.wasan.jp/yamagata/chogenji.html

3 個の甲円が交わり，上の甲円には上の甲円に内接し下の2個の甲円に外接する正方形を入れ，下の甲円には上の甲円に外接し，下の甲円に内接する乙円と，乙円と上の甲円に外接し下の甲円に内接する丙円を入れる。それぞれの大きさを求めよ。

甲円の半径と中心座標を r1, (0, r1), (r1, 0)
乙円の半径と中心座標を r2, (x2, y2)
丙円の半径と中心座標を r3, (x3, y3)
正方形の一辺の長さを a
として，以下の連立方程式を r1 を未知数のままとして解く。

一度に解いてもよいが，乙円・丙円をいれるのと正方形をいれるのは独立なので，まずは乙円と丙円を甲円に入れる。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, x2::positive, y2::negative,
     r3::positive, x3::positive, y3::positive
y2 = r1 - x2
eq1 = (x2 - r1)^2 + y2^2 - (r1 - r2)^2
eq2 = x2^2 + (r1 - y2)^2 - (r1 + r2)^2
eq3 = (x3 - r1)^2 + y3^2 - (r1 - r3)^2
eq4 = (x2 - x3)^2 + (y3 - y2)^2 - (r2 + r3)^2
eq5 = x3^2 + (r1 - y3)^2 - (r1 + r3)^2
solve([eq1, eq2, eq3, eq4, eq5], (r2, x2, r3, x3, y3))

   1-element Vector{NTuple{5, Sym{PyCall.PyObject}}}:
    (sqrt(2)*r1/2, r1*(1 + sqrt(2))/2, sqrt(2)*r1/6, r1*(1 + sqrt(2))/2, r1*(sqrt(2)/6 + 1/2))

右下の甲円の円周上に (sx, sy) を取り，sx = a/2 かつ (sx, sy + a ) が上の甲円の円周上にあるような sx, sy を求める。

@syms a::positive, sx::positive, sy::positive
a = 2sx
eq6 = (sx - r1)^2 + sy^2 - r1^2
eq7 = sx^2 + (sy + a - r1)^2 - r1^2
solve([eq6, eq7], (sx, sy))

   1-element Vector{Tuple{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}}:
    (r1/2, sqrt(3)*r1/2)

まとめると，乙円，丙円の直径はそれぞれ甲円の直径の √2/2，√2/6 倍であり，正方形の一辺の長さは甲円の半径に等しい。

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

甲円の直径 = 1;  乙円の直径 = 0.707107;  丙円の直径 = 0.235702;  正方形の一辺の長さ = 0.5
   r1 = 0.5;  a = 0.5;  r2 = 0.353553;  x2 = 0.603553;  y2 = -0.103553;  r3 = 0.117851;  x3 = 0.603553;  y3 = 0.367851

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 1//2
   (r2, x2, r3, x3, y3) = (sqrt(2)/4, 1/4 + sqrt(2)/4, sqrt(2)/12, 1/4 + sqrt(2)/4, sqrt(2)/12 + 1/4)
   y2 = r1 - x2
   (sx, sy) = (1/4, sqrt(3)/4)
   a = 2sx
   @printf("甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g;  正方形の一辺の長さ = %g\n", 2r1, 2r2, 2r3, a)
   @printf("r1 = %g;  a = %g;  r2 = %g;  x2 = %g;  y2 = %g;  r3 = %g;  x3 = %g;  y3 = %g\n", r1, a, r2, x2, y2, r3, x3, y3)
   plot([sx, sx, -sx, -sx, sx], [sy, sy + a, sy + a, sy, sy], color=:green, lw=0.5)
   circle(0, r1, r1, :blue)
   circle2(r1, 0, r1, :blue)
   circle2(x2, y2, r2)
   circle2(x3, y3, r3, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1, "甲円:r1,(0,r1)", :blue, :center, :bottom, delta=delta)
       point(r1, 0, "甲円:r1,(r1,0)", :blue, :center, :bottom, delta=delta)
       point(x2, y2, "乙円:r2,(x2,y2)", :red, :center, :bottom, delta=delta)
       point(x3, y3, " 丙円:r3,(x3,y3)", :black, :center, delta=-delta/2)
       point(sx, sy, "(sx,sy)", :green, :left, delta=-delta/2)
       point(sx, sy + a, " (sx,sy+a)", :green, :left, :vcenter)
   end
end;

算額（その765）

算額（その765）

山形県七日町 長源寺観音堂 大正4年(1915)
http://www.wasan.jp/yamagata/chogenji.html

外円の中に2本の斜線と大円，中円，小円を入れる。中円の直径が与えられたとき，小円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (0, r1 - R)
中円の半径と中心座標を r2, (0, r2 - R)
小円の半径と中心座標を r3, (x3, y3), (0, 2r1 - R - r3)
斜線と外円の交点座標を (x0, y0); y0 =-sqrt(R^2 - x0^2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, r2::positive,
     r3::positive, x3::positive, y3::positive,
     x0::positive, y0::positive
y0 = -sqrt(R^2 - x0^2)
eq1 = x3^2 + (y3 - r1 + R)^2 - (r1 - r3)^2
eq2 = dist(0, R, x0, y0, 0, 2r1 - R - r3) - r3^2
eq3 = dist(0, R, x0, y0, x3, y3) - r3^2
eq4 = dist(0, R, x0, y0, 0, r2 - R) - r2^2
eq5 = dist(0, R, x0, y0, 0, r1 - R) - (r1 - 2r3)^2
eq6 = r2/(2R - r2) - x0/sqrt((R - y0)^2 + x0^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)
   (R, r1, r3, x3, y3, x0) = u
   t = (-R - sqrt(R^2 - x0^2))
   u = (-2*R + 2*r1 - r3)
   v = (x0^2 + t^2)
   w = -2*R + r1
   x = -R + y3
   y = -2*R + r2
   return [
       x3^2 - (r1 - r3)^2 + (R - r1 + y3)^2,  # eq1
       -r3^2 + x0^2*t^2*u^2/v^2 + (-2*R + 2*r1 - r3 - t^2*u/v)^2,  # eq2
       -r3^2 + (-x0*(x0*x3 + x*t)/v + x3)^2 + (x - t*(x0*x3 + x*t)/v)^2,  # eq3
       -r2^2 + x0^2*y^2*t^2/v^2 + (y - y*t^2/v)^2,  # eq4
       x0^2*w^2*t^2/v^2 - (r1 - 2r3)^2 + (w - w*t^2/v)^2,  # eq5
       r2/(2R - r2) - x0/sqrt(x0^2 + (R + sqrt(R^2 - x0^2))^2),  # eq6
   ]
end;

r2 = 25
iniv = BigFloat[50, 40, 12, 30, 10, 32]
res = nls(H, ini=iniv)

   ([50.39121365997017, 40.15575464114408, 10.077756774346838, 28.39419413598576, -0.3129747007225315, 31.386069651197772], true)

中円の半径が 25 のとき，小円の半径は 10.077756774346838 である。

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

外円の直径 = 100.782;  大円の直径 = 80.3115;  中円の直径 = 50;  小円の直径 = 20.1555
   R = 50.3912;  r1 = 40.1558;  r2 = 25;  r3 = 10.0778;  x3 = 28.3942;  y3 = -0.312975;  x0 = 31.3861;  y0 = -39.4232

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 25
   (R, r1, r3, x3, y3, x0) = res[1]
   y0 = -sqrt(R^2 - x0^2)
   @printf("外円の直径 = %g;  大円の直径 = %g;  中円の直径 = %g;  小円の直径 = %g\n", 2R, 2r1, 2r2, 2r3)
   @printf("R = %g;  r1 = %g;  r2 = %g;  r3 = %g;  x3 = %g;  y3 = %g;  x0 = %g;  y0 = %g\n",
       R, r1, r2, r3, x3, y3, x0, y0)
   plot()
   circle(0, 0, R, :magenta)
   circle(0, r1 - R, r1)
   circle(0, r2 - R, r2, :blue)
   circle2(x3, y3, r3, :green)
   circle(0, 2r1 - R - r3, r3, :green)
   segment(0, R, x0, y0, :orange)
   segment(0, R, -x0, y0, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(0, r1 - R, " 大円:r1\n (0,r1-R)", :red, :left, :vcenter)
       point(0, r2 - R, " 中円:r2\n (0,r2-R)", :blue, :left, :vcenter)
       point(0, 2r1 - R - r3, "小円:r3\n(0,2r1-R-r3)", :green, :center, :bottom, delta=delta/2)
       point(x3, y3, "小円:r3\n(x3,y3)", :green, :center, :bottom, delta=delta/2)
       point(R, 0, " R", :magenta, :left, :bottom, delta=delta/2)
       point(x0, y0, " (x0,y0)")
   end
end;

算額（その764）

算額（その764）

山形県七日町 長源寺観音堂 大正4年(1915)
http://www.wasan.jp/yamagata/chogenji.html

楕円の中に大円（半円）と，中円，小円をいれる。中円，小円の直径が与えられたとき，大円の直径を求めよ。

楕円の長半径と短半径および中心座標を a, b, (0, 0)
大円の半径と中心座標を r0, (0, 2r2 ‐ b)
中円の半径と中心座標を r1, (a - r1, 0)
小円の半径と中心座標を r2, (0, r2 - b)
とおき，以下の連立方程式を解く。

簡単な連立方程式に見えるが，solve() では解けない。

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive,
     r0::positive, r1::positive, r2::positive
@syms a, b, r0, r1, r2
eq1 = r0^2/a^2 + (2r2 - b)^2/b^2 - 1
eq2 = (a - r1)^2 + (2r2 - b)^2 - (r0 + r1)^2
eq3 = r0 + 2r2 - 2b;

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)
   (a, b, r0) = u
   return [
       -1 + (-b + 2*r2)^2/b^2 + r0^2/a^2,  # eq1
       (a - r1)^2 + (-b + 2*r2)^2 - (r0 + r1)^2,  # eq2
       -2*b + r0 + 2*r2,  # eq3
   ]
end;

(r1, r2) = (3, 2)
iniv = BigFloat[22, 11, 17]
res = nls(H, ini=iniv)

   ([22.21485637843628, 10.669764855602098, 17.339529711204197], true)

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

術では「置小圓径三之加中圓径倍」とあるが，正しいものではない。そもそも，中円，小円の直径にかかわらず「三」を加えるのはあきらかにおかしい。

2*17.339529711204197

   34.67905942240839

その他のパラメータは以下のとおりである。
大円の直径 = 34.6791;  中円の直径 = 6;  小円の直径 = 4
r1 = 3;  r2 = 2;  a = 22.2149; b = 10.6698; r0 = 17.3395

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (3, 2)
   (a, b, r0) = res[1]
   @printf("大円の直径 = %g;  中円の直径 = %g;  小円の直径 = %g\n", 2r0, 2r1, 2r2)
   @printf("r1 = %g;  r2 = %g;  a = %g; b = %g; r0 = %g\n", r1, r2, a, b, r0)
   plot()
   ellipse(0, 0, a, b, color=:blue)
   circle(0, 2r2 - b, r0, beginangle=0, endangle=180)
   circle2(a - r1, 0, r1, :magenta)
   circle(0, r2 - b, r2, :green)
   segment(-r0, 2r2 - b, r0, 2r2 - 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", :blue, :left, :bottom, delta=delta/2)
       point(a - r1, 0, "a-r1", :blue, :center, :bottom, delta=delta/2)
       point(0, -b, " -b", :blue, :left, :bottom, delta=delta/2)
       point(0, r2 - b, " r2-b", :blue, :left, :bottom, delta=delta/2)
       point(0, 2r2 - b, " 2r2-b", :red, :left, :bottom, delta=delta/2)
   end
end;

算額（その763）

算額（その763）

山形県七日町 長源寺観音堂 大正4年(1915)
http://www.wasan.jp/yamagata/chogenji.html

外円内に，甲円，乙円，丙円が入っている。乙円の直径が与えられたとき，丙円の直径を求めよ。

引用元の画像では乙円と甲円の関係が不明瞭であるが，図のように，乙円は下部の 2 個の甲円に外接し，外円に内接しているものと思われる。
外円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (0, 3r1 - R), (r1, y1); y1 = 4r1 - R
乙円の半径と中心座標を r2, (x2, 2r1 - R)
丙円の半径と中心座標を r3, (0, R - r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

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

   1-element Vector{NTuple{4, Sym{PyCall.PyObject}}}:
    (44*r2/15, 11*r2/10, 4*sqrt(5)*r2/5, 44*r2/105)

丙円の直径は乙円の直径の 44/105 倍である。

術は一部文字が不鮮明であるが「置四拾四個一百▢除乗乙圓径得丙▢径」と読めるので，「乙円は下部の 2 個の甲円に外接し，外円に内接している」という解釈で間違いなさそうである。

その他，外円の直径は 44/15 倍，甲円の直径は 4√5/5 倍である。

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

   乙円の直径 = 105;  丙円の直径 = 44
   R = 154;  r1 = 57.75;  x2 = 93.9149;  r3 = 22

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r2 = 105//2
   (R, r1, x2, r3) = (44*r2/15, 11*r2/10, 4*sqrt(5)*r2/5, 44*r2/105)
   y1 = 4r1 - R
   @printf("乙円の直径 = %g;  丙円の直径 = %g\n", 2r2, 2r3)
   @printf("R = %g;  r1 = %g;  x2 = %g;  r3 = %g\n", R, r1, x2, r3)
   plot()
   circle(0, 0, R, :blue)
   circle2(r1, y1, r1)
   circle(0, 3r1 - R, r1)
   circle(0, r1 - R, r1)
   circle2(x2, 2r1 - R, r2, :green)
   circle(0, R - r3, r3, :orange)
   if more
       delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
       hline!([0], color=:gray80, lw=0.5)
       vline!([0], color=:gray80, lw=0.5)
       point(R, 0, " R", :blue, :left, :bottom, delta=delta/2)
       point(0, r1 - R, " r1-R", :red, :left, :vcenter)
       point(0, 3r1 - R, "甲円:r1,(0,3r1-R)", :red, :center, delta=-delta)
       point(r1, 4r1 - R, "甲円:r1,(r1,4r1-R)", :red, :center, :bottom, delta=delta)
       point(x2, 2r1 - R, "乙円:r2,(x2,2r1-R)", :green, :center, delta=-delta)
       point(0, R - r3, "丙円:r3,(0,R-r3)", :black, :center, :bottom, delta=delta)
   end
end;

 

算額（その762）

算額（その762）

山形県新庄市 戸沢神社 文政元年(1818)
http://www.wasan.jp/yamagata/tozawa.html

円弧内に斜線2本を入れ，区分された領域に甲円，乙円 2 個ずつを入れる。矢の長さが与えられたときに甲円の径を求めよ。

矢の長さを「矢」
円弧を構成する円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, y1); y1 = 矢/2 + (R - 矢)
乙円の半径と中心座標を r2, (0, R - r2), (0, R - 矢 + r2)
斜線と円弧の交点座標を (sqrt(R^2 - y0^2), y0), (sqrt(R^2 - (R - 矢)^2), R - 矢)
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms R::positive, r1::positive, x1::positive, y1::positive,
     r2::positive, x0::positive, y0::positive, 矢::positive
y1 = 矢/Sym(2) + (R - 矢)
eq1 = x1^2 + y1^2 - (R - r1)^2
eq2 = dist(sqrt(R^2 - (R - 矢)^2), R - 矢, -x0, y0, x1, y1) - r1^2
eq3 = dist(-sqrt(R^2 - (R - 矢)^2), R - 矢, x0, y0, x1, y1) - r1^2
eq4 = dist(-sqrt(R^2 - (R - 矢)^2), R - 矢, x0, y0, 0, R - 矢 + r2) - r2^2
eq5 = dist(-sqrt(R^2 - (R - 矢)^2), R - 矢, x0, y0, 0, R - r2) - r2^2
eq6 = x0 - sqrt(R^2 - y0^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)
   (R, r1, x1, r2, x0, y0) = u
   return [
x1^2 - (R - r1)^2 + (R - 矢/2)^2,  # eq1
-r1^2 + (矢/2 - (矢*(-R + y0 + 矢)/2 + (-x0 - sqrt(R^2 - (R - 矢)^2))*(x1 - sqrt(R^2 - (R - 矢)^2)))*(-R + y0 + 矢)/((-x0 - sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2 + (x1 - sqrt(R^2 - (R - 矢)^2) - (-x0 - sqrt(R^2 - (R - 矢)^2))*(矢*(-R + y0 + 矢)/2 + (-x0 - sqrt(R^2 - (R - 矢)^2))*(x1 - sqrt(R^2 - (R - 矢)^2)))/((-x0 - sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2,  # eq2
-r1^2 + (矢/2 - (矢*(-R + y0 + 矢)/2 + (x0 + sqrt(R^2 - (R - 矢)^2))*(x1 + sqrt(R^2 - (R - 矢)^2)))*(-R + y0 + 矢)/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2 + (x1 + sqrt(R^2 - (R - 矢)^2) - (x0 + sqrt(R^2 - (R - 矢)^2))*(矢*(-R + y0 + 矢)/2 + (x0 + sqrt(R^2 - (R - 矢)^2))*(x1 + sqrt(R^2 - (R - 矢)^2)))/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2,  # eq3
-r2^2 + (r2 - (r2*(-R + y0 + 矢) + sqrt(R^2 - (R - 矢)^2)*(x0 + sqrt(R^2 - (R - 矢)^2)))*(-R + y0 + 矢)/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2 + (sqrt(R^2 - (R - 矢)^2) - (x0 + sqrt(R^2 - (R - 矢)^2))*(r2*(-R + y0 + 矢) + sqrt(R^2 - (R - 矢)^2)*(x0 + sqrt(R^2 - (R - 矢)^2)))/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2,  # eq4
-r2^2 + (sqrt(R^2 - (R - 矢)^2) - (x0 + sqrt(R^2 - (R - 矢)^2))*(sqrt(R^2 - (R - 矢)^2)*(x0 + sqrt(R^2 - (R - 矢)^2)) + (-r2 + 矢)*(-R + y0 + 矢))/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2 + (-r2 + 矢 - (sqrt(R^2 - (R - 矢)^2)*(x0 + sqrt(R^2 - (R - 矢)^2)) + (-r2 + 矢)*(-R + y0 + 矢))*(-R + y0 + 矢)/((x0 + sqrt(R^2 - (R - 矢)^2))^2 + (-R + y0 + 矢)^2))^2,  # eq5
x0 - sqrt(R^2 - y0^2),  # eq6
   ]
end;

矢 = 40
iniv = BigFloat[28.5, 7, 15.9, 6.4, 17, 22.9] .* 矢/28
res = nls(H, ini=iniv)

   ([44.26217677053205, 10.0, 24.191807196045545, 9.531899446242074, 25.79218891003677, 35.97086715238879], true)

甲円の直径は矢の半分である。

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

   甲円の直径 = 20;  乙円の直径 = 19.0638;  矢 = 40
   R = 44.2622;  r1 = 10;  x1 = 24.1918;  y1 = 24.2622;  r2 = 9.5319;  x0 = 25.7922;  y0 = 35.9709

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   矢 = 40
   (R, r1, x1, r2, x0, y0) = res[1]
   y1 =  矢/2 + (R - 矢)
   @printf("甲円の直径 = %g;  乙円の直径 = %g;  矢 = %g\n", 2r1, 2r2, 矢)
   @printf("R = %g;  r1 = %g;  x1 = %g;  y1 = %g;  r2 = %g;  x0 = %g;  y0 = %g\n", R, r1, x1, y1, r2, x0, y0)
   θ = atand(R - 矢, sqrt(R^2 - (R - 矢)))
   plot()
   circle(0, 0, R, :blue, beginangle=θ, endangle=180-θ, n=10000)
   circle2(x1, y1, r1)
   circle(0, R - r2, r2, :green)
   circle(0, R - 矢 + r2, r2, :green)
   segment(-sqrt(R^2 - (R - 矢)^2), R - 矢, sqrt(R^2 - (R - 矢)^2), R - 矢)
   segment(-sqrt(R^2 - (R - 矢)^2), R - 矢, sqrt(R^2 - y0^2), y0)
   segment(sqrt(R^2 - (R - 矢)^2), R - 矢, -sqrt(R^2 - y0^2), y0)
   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-矢", :black, :center, delta=-delta/2)
       point(x1, y1, " 甲円:r1\n (x1,y1)", :red, :center, delta=-delta/2)
       point(0, R - 矢 + r2, " 乙円:r2\n (0,R-矢+r2)", :green, :center, delta=-delta)
       point(0, R - r2, " 乙円:r2\n (0,R-r2)", :green, :center, :bottom, delta=delta)
       point(-sqrt(R^2 - (R - 矢)^2), R - 矢)
       point(R, R - 矢, "(R,R-矢)", :blue, :right, delta=-delta)
   end
end;

算額（その761）

算額（その761）

山形県鶴岡市羽黒町 出羽三山神社 文政6年(1823)
http://www.wasan.jp/yamagata/haguro.html

教材アーカイブス/算数・数学/文献/算額
http://www.s-soarer.jp/?plugin=attach&refer=%E6%95%99%E6%9D%90%E3%82%A2%E3%83%BC%E3%82%AB%E3%82%A4%E3%83%96%E3%82%B9%2F%E7%AE%97%E6%95%B0%E3%83%BB%E6%95%B0%E5%AD%A6%2F%E6%96%87%E7%8C%AE%2F%E7%AE%97%E9%A1%8D&openfile=DSCN1420.JPG

楕円内に2本の斜線と大円（半円）と甲円，乙円，丙円を入れる。乙円と丙円の直径がそれぞれ 3 寸と 2 寸（注）のとき，甲円の直径はいかほどか。

注：前者の URL では数字が判然とせず，検索の結果後者がヒットし，数字が確定できた。

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

include("julia-source.txt");

using SymPy
@syms a::positive, b::positive, r0::positive,
     r1::positive, r2::positive, r3::positive,
     x1::positive, x2::positive,
     y1::positive, y2::positive
eq1 = eq1 = r0^2/a^2 + (b - r0)^2/b^2 -1
eq2 = (a - r3)^2 + (b - r0)^2 - (r0 + r3)^2
eq3 = distance(x1, y1, -x2, -y2, 0, b - r1) - r1^2
eq4 = distance(x1, y1, -x2, -y2, 0, r2 - b) - r2^2
eq5 = distance(x1, y1, -x2, -y2, a - r3, 0) - r3^2
eq6 = x1^2/a^2 + y1^2/b^2 - 1
eq7 = x2^2/a^2 + y2^2/b^2 - 1
eq8 = 2r2 + r0 - 2b

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)
   (a, b, r0, r1, x1, x2, y1, y2) = u
   return [
-1 + (b - r0)^2/b^2 + r0^2/a^2,  # eq1
(a - r3)^2 + (b - r0)^2 - (r0 + r3)^2,  # eq2
-r1^2 + (x1*(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2) - (x1 + x2)*(-b*y1 - b*y2 + r1*y1 + r1*y2 + x1^2 + x1*x2 + y1^2 + y1*y2))^2/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2)^2 + (b - r1 - (b*y1^2 + 2*b*y1*y2 + b*y2^2 - r1*y1^2 - 2*r1*y1*y2 - r1*y2^2 - x1^2*y2 + x1*x2*y1 - x1*x2*y2 + x2^2*y1)/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2))^2,  # eq3
-r2^2 + (x1*(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2) - (x1 + x2)*(b*y1 + b*y2 - r2*y1 - r2*y2 + x1^2 + x1*x2 + y1^2 + y1*y2))^2/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2)^2 + (-b + r2 - (-b*y1^2 - 2*b*y1*y2 - b*y2^2 + r2*y1^2 + 2*r2*y1*y2 + r2*y2^2 - x1^2*y2 + x1*x2*y1 - x1*x2*y2 + x2^2*y1)/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2))^2,  # eq4
-r3^2 + (x1 + x2)^2*(a*y1 + a*y2 - r3*y1 - r3*y2 - x1*y2 + x2*y1)^2/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2)^2 + (a - r3 - ((a - r3)*(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2) - (y1 + y2)*(a*y1 + a*y2 - r3*y1 - r3*y2 - x1*y2 + x2*y1))/(x1^2 + 2*x1*x2 + x2^2 + y1^2 + 2*y1*y2 + y2^2))^2,  # eq5
-1 + y1^2/b^2 + x1^2/a^2,  # eq6
-1 + y2^2/b^2 + x2^2/a^2,  # eq7
-2*b + r0 + 2*r2,  # eq8
   ]
end;

(r2, r3) = (3, 2) .// 2
iniv = BigFloat[14, 7, 12, 6, 15, 6, 2, 7]
res = nls(H, ini=iniv)

   ([10.772733743669523, 6.13102420004887, 9.26204840009774, 4.383884823343673, 10.493949340307568, 5.4542418519895035, 1.3857702087399792, 5.287131918198976], true)

乙円と丙円の直径がそれぞれ 3 寸と 2 寸のとき，甲円の直径は 8.76777 寸となった。

「答」では「甲円の直径は十一寸零三分八厘六毛四一九有奇」となっており，計算結果と合わない。連立方程式に使った条件式のどれかが不適切なのかもしれない。

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

   甲円の直径 = 8.76777;  乙円の直径 = 3;  丙円の直径 = 2
   a = 10.7727;  b = 6.13102;  r0 = 9.26205;  r1 = 4.38388;  x1 = 10.4939;  x2 = 5.45424;  y1 = 1.38577;  y2 = 5.28713

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   #(r2, r3) = (0.17, 0.12)
   (a, b, r0, r1, x1, x2, y1, y2) = res[1]
   @printf("甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r1, 2r2, 2r3)
   @printf("a = %g;  b = %g;  r0 = %g;  r1 = %g;  x1 = %g;  x2 = %g;  y1 = %g;  y2 = %g\n", a, b, r0, r1, x1, x2, y1, y2)
   plot()
   ellipse(0, 0, a, b, color=:red)
   circle(0, b - r1, r1, :blue)
   circle(0, r2 - b, r2, :green)
   circle2(a - r3, 0, r3, :magenta)
   circle(0, b - r0, r0, :orange, beginangle=0, endangle=180)
   segment(-r0, 2r2 - b, r0, 2r2 - b, :orange)
   segment(x1, y1, -x2, -y2)
   segment(-x1, y1, x2, -y2)
   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, b - r1, " 甲円:r1\n (0,b-r1)", :blue, :left, :vcenter)
       point(0, r2 - b, " 乙円:r2\n (0,r2-b)", :green, :left, :vcenter)
       point(r0, b - r0, "(r0,b-r0)", :orange, :right, :bottom, delta=delta/2)
       point(a - r3, 0, "丙円:r3\n(a-r3,0)", :black, :right, delta=-delta/2)
       point(a, 0, " a", :red, :left, :vcenter)
       point(0, b, " b", :red, :left, :bottom, delta=delta)
       point(x1, y1, "(x1,y1) ", :black, :right, :bottom, delta=delta)
       point(x2, -y2, " (x2,-y2)", :black, :left, :vcenter)
   end
end;

算額（その760）

算額（その760）

秋田県仙北市角館町西長野 熊野神社 弘化2年(1846)
http://www.wasan.jp/akita/kakunodatekumano2.html

直角三角形の中に正方形と，直径がそれぞれ 8 寸，6寸の甲円と乙円が入っている。正方形の一辺の長さはいかほどか。

直角三角形の直角を挟むに辺のうち，短い方を「鈎」，長い方を「股」
甲円の半径と中心座標を r1, (股 - a - r1, r1)
乙円の半径と中心座標を r2, (a - r2, a + r2)
正方形の一辺の長さを a
とおき，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive, a::positive,
     鈎::positive, 股::positive
eq1 = r1/a - r2/(鈎 - a)
eq2 = a/(股 - a) - 鈎/股
eq3 = dist(0, 0, 股, 鈎, 股 - a - r1, r1) - r1^2;
res = solve([eq1, eq2, eq3], (a, 鈎, 股))
res[2]

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

2 組の解が得られるが，2 番目のものが適解である。
正方形の一辺の長さは 甲円と乙円の直径の二乗和の平方根を取り，甲円と乙円の直径を加える。
甲円と乙円の直径がそれぞれ 8 寸，6寸のとき，正方形の一辺の長さは 12 寸（1 尺 2 寸）である。

(8 + 6 + sqrt(8^2 + 6^2))/2

   12.0

ちなみに，鈎，股はそれぞれ 21, 28 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (r1, r2) = (8, 6) .// 2
   (a, 鈎, 股) = (r1 + r2 + sqrt(r1^2 + r2^2), (r1 + r2)*(r1 + r2 + sqrt(r1^2 + r2^2))/r1, (r1 + r2)*(r1 + r2 + sqrt(r1^2 + r2^2))/r2)
   plot([0, 股, 股, 0], [0, 0, 鈎, 0], color=:black, lw=0.5)
   circle(股 - a - r1, r1, r1)
   circle(股 - r2, a + r2, r2, :blue)
   segment(股 - a, 0, 股 - a, a, :green)
   segment(股 - a, a, 股, a, :green)
   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 - r1, r1, "甲円:r1\n(股-a-r1,r1)", :red, :center, delta=-delta/2)
       point(股 - r2, a + r2, "乙円:r2\n(股-r2,a+r2)", :blue, :center, :bottom, delta=delta/2)
       point(股, 0, " 股", :black, :left, :bottom, delta=delta/2)
       point(股, 鈎, "(股,鈎)", :black, :right, :bottom, delta=delta)
       point(股, a, "(股,a)", :green, :right, delta=-delta)
       point(股 - a, 0, " 股-a", :green, :left, :bottom, delta=delta)
   end
end;

算額（その759）

算額（その759）

秋田県仙北市角館町西長野 熊野神社 弘化2年(1846)
http://www.wasan.jp/akita/kakunodatekumano2.html

直径が 4 寸の大円の中に，大円と同じ直径を持つ 2 個の弧と乙円 2 個が入っている。乙円の直径はいかほどか。

大円の半径と中心座標を r1, (0, 0), (0, r1), (0, -r1)
乙円の半径と中心座標を r2, (r1 - r2, 0)
として，以下の方程式を解く。

include("julia-source.txt");

using SymPy
@syms r1::positive, r2::positive
eq1 = (r1 - r2)^2 + r1^2 - (r1 + r2)^2
res = solve(eq1, r2)

   1-element Vector{Sym{PyCall.PyObject}}:
    r1/4

乙円の直径は，大円の直径の 1/4 である。
大円の直径が 4 寸のとき，乙円の直径は 1 寸である。

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 4//2
   r2 = r1/4
   @printf("乙円の直径 = %g;  大円の直径 = %g\n", 2r1, 2r2)
   plot()
   circle(0, 0, r1)
   circle(0, r1, r1, beginangle=210, endangle=330)
   circle(0, -r1, r1, beginangle=30, endangle=150)
   circle2(r1 - r2, 0, r2, :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)
   end
end;

 

算額（その758）

算額（その758）

秋田県仙北市角館町西長野 熊野神社 弘化2年(1846)
http://www.wasan.jp/akita/kakunodatekumano2.html

半梯（半分の等脚台形の横倒し）内に，直角三角形と，天円，地円，人円が入っている。
鈎が 63 寸（注），地円と人円の直径がそれぞれ 28 寸，21 寸のとき，天円の直径はいかほどか。

注：算額の写真が不鮮明で，61 寸なのか 62 寸 なのか，はたまた 63 寸なのかが判然としない。63 寸としたときに，「答」と一致した。そうと分かれば，63 寸のようにも見える。

半梯の高さ，下底，上底がそれぞれ a, b, c
直角三角形の直角である頂点が x 軸と接する x 座標を d
天円の半径と中心座標を r1, (x1, y1)
地円の半径と中心座標を r2, (r2, r2)
人円の半径と中心座標を r3, (a - r3, r3)
として，以下の連立方程式を解く。

include("julia-source.txt");

using SymPy
@syms a::positive, c::positive, cd::positive, bd::positive,
     bc::positive, b::positive, d::positive, r1::positive,
     r2::positive, r3::positive, x1::positive, y1::positive
eq1 = c + (a - d) - cd - 2r3
eq2 = cd + bd - bc - 2r1
eq3 = d + b - bd - 2r2
eq4 = a^2 + (b - c)^2 - bc^2
eq5 = c^2 + (a - d)^2 - cd^2
eq6 = b^2 + d^2 - bd^2
eq7 = (bc + cd + bd)*r1 + d*b + (a - d)*c - (c + b)*a
res = solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7], (r1, a, b, c, d, bc, bd));
res[3]

   (cd*(r2 + r3 - sqrt(r2^2 + r3^2))/(2*r3), cd*r2/(2*r3) + cd/2 + r2 - r2*sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/(2*r3) + r3 + sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/2, r2*(cd + 2*r3 + sqrt(cd^2 - 4*cd*r3 - 4*r3^2))/(2*r3), cd/2 + r3 - sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/2, r2*(cd + 2*r3)/(2*r3) - r2*sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/(2*r3), cd*sqrt(r2^2 + r3^2)/r3, cd*r2/r3)

4 組の解が得られるが，3 番目のものが適解である。
なお，4 番目の解は 3 番目の解と線対称ではないが，上底と下底の位置関係が算額のものと逆である。

「術」は「地円と人円それぞれ直径の二乗の和の平方根を，地円と人円の直径の和から引き，鈎を掛け人円の直径で割る」となっており，cd*(r2 + r3 - sqrt(r2^2 + r3^2))/(2*r3) と一致する。

次に，図を描くために天円の中心座標を求める。
4 組の解が得られるが，最初のものが適解である。

天円の直径は 42 寸である。

using SymPy
@syms a::positive, c::positive, cd::positive, bd::positive,
     bc::positive, b::positive, d::positive, r1::positive,
     r2::positive, r3::positive, x1::positive, y1::positive
(r1, a, b, c, d, bc, bd) = (21, 98 - 7*sqrt(Sym(2))/2, 14*sqrt(Sym(2)) + 56, 42 - 21*sqrt(Sym(2))/2, 56 - 14*sqrt(Sym(2)), 105, 84)
eq8 = distance(0, b, a, c, x1, y1) - r1^2
eq9 = distance(0, b, d, 0, x1, y1) - r1^2
res2 = solve([eq8, eq9], (x1, y1));
res2[1]

   (56 - 7*sqrt(2), 28)

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

   天円の直径 = 42;  地円の直径 = 28;  人円の直径 = 21;  鈎 = 63
   r1 = 21;  a = 93.0503;  b = 75.799;  c = 27.1508;  d = 36.201;  bc = 105;  bd = 84;  cd = 63

function draw(more=false)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   cd = 63
   r2 = 28//2
   r3 = 21//2
   (r1, a, b, c, d, bc, bd) = (cd*(r2 + r3 - sqrt(r2^2 + r3^2))/(2*r3), cd*r2/(2*r3) + cd/2 + r2 - r2*sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/(2*r3) + r3 + sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/2, r2*(cd + 2*r3 + sqrt(cd^2 - 4*cd*r3 - 4*r3^2))/(2*r3), cd/2 + r3 - sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/2, r2*(cd + 2*r3)/(2*r3) - r2*sqrt(cd^2 - 4*cd*r3 - 4*r3^2)/(2*r3), cd*sqrt(r2^2 + r3^2)/r3, cd*r2/r3)
   (x1, y1) = (56 - 7*sqrt(2), 28)
   @printf("天円の直径 = %g;  地円の直径 = %g;  人円の直径 = %g;  鈎 = %g\n", 2r1, 2r2, 2r3, cd)
   @printf("r1 = %g;  a = %g;  b = %g;  c = %g;  d = %g;  bc = %g;  bd = %g;  cd = %g\n",
       r1, a, b, c, d, bc, bd, cd)
   plot([0, a, a, 0, 0], [0, 0, c, b, 0], color=:black, lw=0.5)
   circle(x1, y1, r1, :green)
   circle(r2, r2, r2)
   circle(a - r3, r3, r3, :blue)
   segment(0, b, d, 0)
   segment(d, 0, a, c)
   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", :black, :left, :bottom, delta=delta/2)
       point(d, 0, "d  ", :black, :right, :bottom, delta=delta/2)
       point(a, c, " (a,c)", :black, :left, :bottom, delta=delta/2)
       point(0, b, " b", :black, :left, :bottom, delta=delta/2)
       point(x1, y1, "天円:r1,(r1,r1)", :green, :center, delta=-delta/2)
       point(r2, r2, "地円:r2,(r2,r2)", :red, :center, delta=-delta/2)
       point(a - r3, r3, "人円:r3\n(a-r3,r3)", :blue, :center, delta=-delta/2)
   end
end;