算額（その1441）

算額（その1441）

二本松亀谷町坂 観世音 天保10年(1839)
『諸流算額起原一』（佐久間庸軒）

～落書き帳「○△□」～  81. 佐久間文庫
http://streetwasan.web.fc2.com/math15.9.5.html
キーワード：円3個，正三角形，接線
#Julia, #SymPy, #算額, #和算

正三角形内に内接円と甲円，乙円を容れ，3円の共通接線を引く。内接円の直径が 675 寸，甲円の直径が 81 寸のとき，乙円の直径はいかほどか。

最初は，後半に示すような全ての描画パラメータを求める方法で解いた。かなり複雑な7元連立方程式になったので，数値解しか求められなかった。
それでも，乙円は内接円の 1/27 倍の大きさであることは分かった。

その後で，「落書き帳」の著者が「三角関数の加法定理を用いてすんなり？解けます」というヒントを与えていたので，以下のように「簡単に」解けた。図を描くのは別として。

1.　乙円の直径のみを得る

内接円，甲円，乙円の半径を r0, r1, r2 とおく。
甲円の中心 A，内接円の中心 O，内接円と正三角形の斜辺との接点 B を結ぶ直角三角形において，∠AOB = θ1
乙円の中心 C，内接円の中心 O，内接円と正三角形の斜辺との接点 D を結ぶ直角三角形において，∠COD = θ2
とすれば θ1 + θ2 = 60°なので，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms r0, r1, r2, θ1, θ2
θ2 = PI/3 - θ1
eq1 = (r0 + r1)*cos(θ1) - (r0 - r1)
eq2 = (r0 + r2)*cos(θ2) - (r0 - r2)
res = solve([eq1, eq2], (r2, θ1))[2]  # 2 of 2

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

res[1](r0 => 675, r1 => 81).evalf() |> println
res[2](r0 => 675, r1 => 81).evalf() |> println

    25.0000000000000
    0.666946344503664

内接円，甲円の直径がそれぞれ 675 寸，81 寸のとき，乙円の直径は 25 寸である。

2.　全ての描画パラメータを得る

内接円の半径と中心座標を r0, (0, r0)
正三角形の一辺の長さを 2a; a = √3r0
甲円の半径と中心座標を r1, (x1, y1)
乙円の半径と中心座標を r2, (x2, y2)
内接円，甲円，乙円の共通接線と正三角形の斜辺との交点座標を (cx, cy), (dx, dy)
cy = √3(a + cx), dy = √3(a - dx)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, r0::positive,
      r1::positive, x1::positive, y1::positive,
      r2::positive, x2::positive, y2::positive,
      cx::negative, cy::positive, dx::positive, dy::positive
a = r0*√Sym(3)
cy = √Sym(3)*(a + cx)
dy = √Sym(3)*(a - dx)

eq1 = x1^2 + (y1 -  r0)^2 - (r0 + r1)^2
eq2 = x2^2 + (y2 -  r0)^2 - (r0 + r2)^2
eq3 = dist2(cx, cy, dx, dy, x1, y1, r1)
eq4 = dist2(cx, cy, dx, dy, x2, y2, r2)
eq5 = dist2(cx, cy, dx, dy, 0, a/√Sym(3), r0)
eq6 = dist2(a, 0, 0, √Sym(3)a, x1, y1, r1)
eq7 = dist2(-a, 0, 0, √Sym(3)a, x2, y2, r2);

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, y1, r2, x2, y2, cx, dx) = u
    return [
        x1^2 + (-r0 + y1)^2 - (r0 + r1)^2,  # eq1
        x2^2 + (-r0 + y2)^2 - (r0 + r2)^2,  # eq2
        (12*cx^2*dx^2 - 12*sqrt(3)*cx^2*dx*r0 - 12*cx^2*dx*x1 + 4*sqrt(3)*cx^2*dx*y1 + 9*cx^2*r0^2 + 6*sqrt(3)*cx^2*r0*x1 - 6*cx^2*r0*y1 - 4*cx^2*r1^2 + 3*cx^2*x1^2 - 2*sqrt(3)*cx^2*x1*y1 + cx^2*y1^2 + 12*sqrt(3)*cx*dx^2*r0 - 12*cx*dx^2*x1 - 4*sqrt(3)*cx*dx^2*y1 - 18*cx*dx*r0^2 + 12*cx*dx*r0*y1 - 4*cx*dx*r1^2 + 6*cx*dx*x1^2 - 2*cx*dx*y1^2 + 9*dx^2*r0^2 - 6*sqrt(3)*dx^2*r0*x1 - 6*dx^2*r0*y1 - 4*dx^2*r1^2 + 3*dx^2*x1^2 + 2*sqrt(3)*dx^2*x1*y1 + dx^2*y1^2)/(4*(cx^2 + cx*dx + dx^2)),  # eq3
        (12*cx^2*dx^2 - 12*sqrt(3)*cx^2*dx*r0 - 12*cx^2*dx*x2 + 4*sqrt(3)*cx^2*dx*y2 + 9*cx^2*r0^2 + 6*sqrt(3)*cx^2*r0*x2 - 6*cx^2*r0*y2 - 4*cx^2*r2^2 + 3*cx^2*x2^2 - 2*sqrt(3)*cx^2*x2*y2 + cx^2*y2^2 + 12*sqrt(3)*cx*dx^2*r0 - 12*cx*dx^2*x2 - 4*sqrt(3)*cx*dx^2*y2 - 18*cx*dx*r0^2 + 12*cx*dx*r0*y2 - 4*cx*dx*r2^2 + 6*cx*dx*x2^2 - 2*cx*dx*y2^2 + 9*dx^2*r0^2 - 6*sqrt(3)*dx^2*r0*x2 - 6*dx^2*r0*y2 - 4*dx^2*r2^2 + 3*dx^2*x2^2 + 2*sqrt(3)*dx^2*x2*y2 + dx^2*y2^2)/(4*(cx^2 + cx*dx + dx^2)),  # eq4
        cx*dx*(3*cx*dx - 2*sqrt(3)*cx*r0 + 2*sqrt(3)*dx*r0 - 3*r0^2)/(cx^2 + cx*dx + dx^2),  # eq5
        9*r0^2/4 - 3*sqrt(3)*r0*x1/2 - 3*r0*y1/2 - r1^2 + 3*x1^2/4 + sqrt(3)*x1*y1/2 + y1^2/4,  # eq6
        9*r0^2/4 + 3*sqrt(3)*r0*x2/2 - 3*r0*y2/2 - r2^2 + 3*x2^2/4 - sqrt(3)*x2*y2/2 + y2^2/4,  # eq7
    ]
end;

r0 = 675/2
r1 = 81/2
iniv = BigFloat[140, 689, 12, -217, 612, -225, 159]
res = nls(H, ini=iniv)

    ([140.29611541307898, 688.5, 12.499999999999389, -216.5063509461113, 612.499999999998, -224.83351829019105, 159.42740387849864], true)

内接円，甲円の直径が 675 寸，81 寸のとき，乙円の直径は 25 寸である。
2r2 = 25，2r0 = 675 なので，乙円の直径は内接円の直径の 1/27 倍である。
ちなみに，甲円の直径は内接円の直径の 3/25 倍になっている。

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r0 = 675/2
    r1 = 081/2
    (x1, y1, r2, x2, y2, cx, dx) = [140.29611541307898, 688.5, 12.499999999999389, -216.5063509461113, 612.499999999998, -224.83351829019105, 159.42740387849864]
    a =  √3r0
    cy = √3(a + cx)
    dy = √3(a - dx)
    plot([a, 0, -a, a], [0, √3a, 0, 0], color=:brown, lw=0.5)
    circle(0, a/√3, a/√3, :orange)
    circle(x1, y1, r1, :blue)
    circle(x2, y2, r2)
    segment(cx, cy, dx, dy, :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, 0, "a", :brown, :left, :bottom, delta=delta)
        point(0, r0, "O", :red, :center, delta=-delta)
        point(x1, y1, "A", :blue, :center, :vcenter, deltax=5delta)
        point(x2, y2, "C", :red, :right, :vcenter, deltax=-2delta)
        (x01, y01) = foot(a, 0, 0, √3a, 0, r0)
        point(x01, y01, " B", :orange, :left, :vcenter)
        (x02, y02) = foot(-a, 0, 0, √3a, 0, r0)
        point(x02, y02, "D ", :orange, :right, :vcenter)
    end
end;

draw(true)

エビス ウドンファクトリー

高松市香川町浅野 エビス ウドンファクトリー
さぬきうどんの冬季メニュー「しっぽくうどん」（上の写真の中央に，ドンと書かれています）
いわゆる「卓袱料理」から連想されるものとは大違い

算額（その1440）

算額（その1440）

四 岩手県花巻市南笹間 東光寺 慶応2年(1866)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円1個，外円，正方形4個
#Julia, #SymPy, #算額, #和算

円の中に大正方形 1 個，小正方形 3 個を容れる。大正方形の一辺の長さが 1 寸のとき，小正方形の一辺の長さはいかほどか。

円の半径と中心座標を R, (0, 0)
大正方形の一辺の長さを a
小正方形の一辺の長さを b
正方形の頂点座標を (ax, ay), ..., (fx, fy)
とおき，方程式を解く。
座標点のうちいくつかは特定できる。結局，求めるべきパラメータは b, R, cx である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms R::positive, a::positive, b::positive,
      ax::positive, ay::positive, bx::positive, by::positive,
      cx::positive, cy::positive, dx::positive, dy::positive,
      ex::positive, ey::positive, fx::positive, fy::positive
ax = a/2
ay = sqrt(R^2 - ax^2)
bx = ax
by = ay - a
cy = -sqrt(R^2 - cx^2)
dx = a//2 ### 数値解から得られた関係式
dy = -sqrt(R^2 - dx^2)
ex = b/2
ey = -a//2 ### 数値解から得られた関係式
fx = ex
fy = ey - b

eq1 = (cx - dx)^2 + (cy - dy)^2 - b^2 |> expand
eq2 = (ex - bx)^2 + (ey - by)^2 - b^2 |> expand
eq3 = (bx - dx)^2 + (by - dy)^2 - 2b^2 |> expand;

res = solve([eq1, eq2, eq3], (b, R, cx))[4]  # 4 of 4

    (a*sqrt(-3*sqrt(2) - sqrt(11 - 6*sqrt(2)) + 6), a*sqrt(12 - 6*sqrt(2))/2, a*(7*sqrt(34 - 24*sqrt(2)) + sqrt(22 - 12*sqrt(2)) + 14*sqrt(17 - 12*sqrt(2)) + 2*sqrt(11 - 6*sqrt(2)))/12)

# b  二重根号を外す
res[1] |> sympy.sqrtdenest |> println

    a*(-1 + sqrt(2))

# R
res[2] |> println

    a*sqrt(12 - 6*sqrt(2))/2

# cx  二重根号を外す
res[3] |> sympy.sqrtdenest |> simplify |> println

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

小正方形の一辺の長さ b は，大正方形の一辺の長さ a の (√2 - 1) 倍である。
大正方形の一辺の長さが 1 寸のとき，小正方形の一辺の長さは (√2 - 1) = 0.41421356237309515 寸である。

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    b = a*(√2 - 1)
    R = a*sqrt(12 - 6√2)/2
    cx = a*(3 - √2)/2
    ax = a/2
    ay = sqrt(R^2 - ax^2)
    bx = ax
    by = ay - a
    cy = -sqrt(R^2 - cx^2)
    dx = a//2 ### 数値解から得られた関係式
    dy = -sqrt(R^2 - dx^2)
    ex = b/2
    ey = -a//2 ### 数値解から得られた関係式
    fx = ex
    fy = ey - b
    @printf("a = %g;  b = %g;  R = %g;  cx = %g\n", a, b, R, cx)
    plot([bx, ax, -ax, -bx, bx], [by, ay, ay, by, by], color=:blue, lw=0.5)
    circle(0, 0, R)
    plot!([cx, bx, ex, dx, cx], [cy, by, ey, dy, cy], color=:green, lw=0.5)
    plot!(-[cx, bx, ex, dx, cx], [cy, by, ey, dy, cy], color=:green, lw=0.5)
    plot!([fx, ex, -ex, -fx, fx], [fy, ey, ey, fy, fy], color=:green, lw=0.5)
    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(ax,  ay, " (ax,ay)", :green, :left, :vcenter)
        point(bx,  by, " (bx,by)", :green, :left, :vcenter)
        point(cx,  cy, " (cx,cy)", :green, :left, :vcenter)
        point(dx,  dy, " (dx,dy)", :green, :left, :vcenter)
        point(ex,  ey, " (ex,ey)", :green, :left, :vcenter)
        point(fx,  fy, " (fx,fy)", :green, :left)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
    end
end;

draw(1, true)

算額（その1439）

算額（その1439）

郡山市西田町 旧土棚村 鹿嶋神社 明治16年(1883）
～落書き帳「○△□」～ 79. 金木犀の花びら
http://streetwasan.web.fc2.com/math15.9.3.html
キーワード：正方形，正五角形
#Julia, #SymPy, #算額, #和算

正方形と正五角形が図のように重なっている。正方形の一辺の長さが与えられたとき，正五角形の一辺の長さを求めよ。

正方形の一辺の長さを a
正五角形の一辺の長さを b
正五角形の外接円の半径を r
正五角形の高さを h
とおき，以下の方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, b::positive, r::positive, h::positive
eq1 = r*sind(Sym(36)) - b/2
eq2 = (r + r*cosd(Sym(36))) - h
eq3 = (a/2 - b /2)/(h - a) - tand(Sym(18))
res = solve([eq1, eq2, eq3], (b, r, h))

それぞれの式は長いが，かなり簡約化できる。

# b
@syms d
res[b]/a |> x -> apart(x, d) |> factor |> x -> a*x |> println
res[b](a => 1).evalf() |> println

    a*(2*sqrt(5)*sqrt(5 - 2*sqrt(5)) + 5)/10
    0.824919696232906

正五角形の一辺の長さは，正方形の一辺の長さの (2√5sqrt(5 - 2√5) + 5)/10 倍である。
正方形の一辺の長さが 1 のとき，正五角形の一辺の長さは 0.824919696232906 である。

# r
res[r]/a |> x -> apart(x, d) |> factor |> x -> a*x |> println
res[r](a => 1).evalf() |> println

    a*(-2*sqrt(5) + 5*sqrt(5 - 2*sqrt(5)) + 3*sqrt(5)*sqrt(5 - 2*sqrt(5)) + 10)/20
    0.701718606426041

# h
res[h]/a |> x -> apart(x, d) |> factor |> x -> a*x |> println
res[h](a => 1).evalf() |> println

    a*(2*sqrt(5 - 2*sqrt(5)) + sqrt(5)*sqrt(5 - 2*sqrt(5)) + 2)/4
    1.26942088429381

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    t = sqrt(5 - 2√5)
    b = a*(2*sqrt(5)*sqrt(5 - 2*sqrt(5)) + 5)/10
    r = a*(-2*sqrt(5) + 5*sqrt(5 - 2*sqrt(5)) + 3*sqrt(5)*sqrt(5 - 2*sqrt(5)) + 10)/20
    h = a*(2*sqrt(5 - 2*sqrt(5)) + sqrt(5)*sqrt(5 - 2*sqrt(5)) + 2)/4
    x = [r*cosd(18), 0, -r*cosd(18), -r*sind(36), r*sind(36), r*cosd(18)]
    y = [r*cosd(36) + r*sind(18), h, r*cosd(36) + r*sind(18), 0, 0, r*cosd(36) + r*sind(18)]
    @printf("正方形の一辺の長さが %g のとき，正五角形の一辺の長さは %g である。\n", a, b)
    plot(x, y, color=:blue, lw=0.5)
    circle(0, r*cosd(36), r)
    rect(-a/2, h - a, a/2, h, :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(0, h, "h", :red, :center, :bottom, delta=delta/2)
        point(b/2, 0, "b/2 ", :blue, :right, :bottom, delta=delta/2)
        point(a/2, h - a, "(a/2,h-a) ", :green, :right, :bottom, delta=delta/2)
        point(0, r*cosd(36), "外接円の中心", :red, :center, delta=-delta/2)
        point(0.82492/2,0,"***")
    end
end;

draw(1, true)

 

算額（その1438）

算額（その1438）

福島県福島市飯野町青木小手神森 小手神社 明治19年(1886)
～落書き帳「○△□」～ 75．変なお供え
http://streetwasan.web.fc2.com/math15.8.30.html
キーワード：円周，円周率
#Julia, #SymPy, #算額, #和算

直線の上に 1 個の正方形と 3 個の等円が載っている。正方形の一辺の長さが与えられたとき，等円の直径を求めよ。

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

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms r::positive, a::positive
eq = (a/2 + r)^2 + a^2 - (2r)^2
solve(eq, r)[1] |> println

    5*a/6

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

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = a*(5/6)
    plot([a/2, a/2, -a/2, -a/2, a/2], [0, a, a, 0, 0], color=:blue, lw=0.5)
    circle(0, a + r, r)
    circle2(a/2 + r, 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, a + r, "等円:r,(0,a+r)", :red, :center, delta=-delta)
        point(a/2 + r, r, "等円:r,(a/2+r,0)", :red, :center, delta=-delta)
        point(0, a, "a", :blue, :center, :bottom, delta=delta)
        point(a/2, 0, " a/2", :blue, :left, :bottom, delta=delta)
    end
end;

draw(1, true)

 

算額（その1437）

算額（その1437）

田村郡三春町御木沢 厳島神社 明治18年(1885)
http://www.wasan.jp/fukusima/miharuitukusima6.html
キーワード：円周，円周率
#Julia, #SymPy, #算額, #和算

長さが 10 尺のヘビがいる。これを「の」の字に巻いたとき，直径はいかほどか。ただし，円周率としては 3．16 を用いよ。

算額の写真（ヘビが怖い人はクリックしないように）

ちょっと変わった問題で，スキップしていたのであるがやってみた。

『「の」の字』というのは，算額の問を見るとよく分かるが，「円の中に垂直な直径を引いた図形」ということで，10尺のヘビの尻尾の先端を A 点に置き，まっすぐ下に引っ張り，円周上の B 点に届いたら，時計回りに丸く円周に沿って伸ばしていき，最後に頭が B 点に届く」ということである。

AB の長さを 2r とすれば，円周は 3.16*2r で，これに直径 2r を加えればヘビの長さになるということで，方程式を立てて解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms r::positive
eq = (3.16 * 2r + 2r) ⩵ 10
ans_r = solve(eq, r)[1]
2*ans_r |> println

    2.40384615384615

直径は 2.40384615384615 ≒ 二尺四寸零分三厘余

術は，10/(3.16 + 1) = 2.4038461538461537

function draw(more)
    pyplot(size=(150, 150), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 1
    plot(showaxis=false)
    circle(0, 0, R, beginangle=275, endangle=630, lw=1)
    segment(0,R, 0, -R, lw=1, :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(0, R, "A", :red, :center, :bottom, delta=2delta)
        point(0, -R, "B", :red, :center, delta=-2delta)
    end
end;

draw(true)

算額（その1436）

算額（その1436）

福島市立子山 篠葉沢稲荷 明治10年(1877)
街角の数学 ～落書き帳「○△□」～ 68．ウルトラマン・ピーチ
http://streetwasan.web.fc2.com/math15.8.19.html
キーワード：円3個，半円2個，長方形
#Julia, #SymPy, #算額, #和算

外円の中に 2 本の斜線（等斜）を引き，区画された領域に 4 個の等円を容れる。外円の直径が 13 寸のとき，等円の直径はいかほどか。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (x1, y1), (r, y2)
斜線と円の交点座標を (x0, -sqrt(r^2 - x0^2))
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms R::positive, r::positive, x1::positive, y1::positive, y2::positive, x0::positive
t = sqrt(R^2 - x0^2)
eq1 = x1^2 + y1^2 - (R - r)^2
eq2 = y1* (t + R) - x0*x1  # y1/x1 - x0/(t + R)
eq3 = dist2(0, R, x0, -t, x1, y1, r) |> numerator
eq4 = dist2(0, R, x0, -t, r, y2, r) |> numerator
eq5 = sqrt(x0^2 + (R - t)^2) + sqrt(x0^2 + (R + t)^2) - 2R - 2r;
# res = solve([eq1, eq2, eq3, eq4, eq5], (r, x1, y1, y2, x0))

using NLsolve

function nls(func, params...; ini = [0.0])
    if typeof(ini) <: Number
        r = nlsolve((vout, vin) -> vout[1] = func(vin[1], params..., [ini]), ftol=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, x1, y1, y2, x0) = u
    t = sqrt(R^2 - x0^2)
    return [
        x1^2 + y1^2 - (R - r)^2,  # eq1
        -x0*x1 + y1*(R + sqrt(R^2 - x0^2)),  # eq2
        -4*R^3*r^2 + R^3*x0^2 - 4*R^3*x0*x1 + 4*R^3*x1^2 - 4*R^2*r^2*sqrt(R^2 - x0^2) - 2*R^2*x0^2*y1 + R^2*x0^2*sqrt(R^2 - x0^2) + 4*R^2*x0*x1*y1 - 4*R^2*x0*x1*sqrt(R^2 - x0^2) + 4*R^2*x1^2*sqrt(R^2 - x0^2) + 2*R*r^2*x0^2 + 2*R*x0^3*x1 - 3*R*x0^2*x1^2 + R*x0^2*y1^2 - 2*R*x0^2*y1*sqrt(R^2 - x0^2) + 4*R*x0*x1*y1*sqrt(R^2 - x0^2) - 2*x0^3*x1*y1 - x0^2*x1^2*sqrt(R^2 - x0^2) + x0^2*y1^2*sqrt(R^2 - x0^2),  # eq3
        x0*(-4*R^3*r + R^3*x0 + 4*R^2*r*y2 - 4*R^2*r*sqrt(R^2 - x0^2) - 2*R^2*x0*y2 + R^2*x0*sqrt(R^2 - x0^2) - R*r^2*x0 + 2*R*r*x0^2 + 4*R*r*y2*sqrt(R^2 - x0^2) + R*x0*y2^2 - 2*R*x0*y2*sqrt(R^2 - x0^2) - r^2*x0*sqrt(R^2 - x0^2) - 2*r*x0^2*y2 + x0*y2^2*sqrt(R^2 - x0^2)),  # eq4
        -2*R - 2*r + sqrt(x0^2 + (R - sqrt(R^2 - x0^2))^2) + sqrt(x0^2 + (R + sqrt(R^2 - x0^2))^2),  # eq5
    ]
end;

R = 13/2
iniv = BigFloat[1.9, 4.1, 1.7, -3.4, 4.6]
res = nls(H, ini=iniv)

    ([2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615], true)

外円の直径が 13 寸のとき，等円の直径は 4 寸である。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r, x1, y1, y2, x0) = [2.0, 4.153846153846154, 1.7307692307692308, -3.5, 4.615384615384615]
    t = sqrt(R^2 - x0^2)
    @printf("外円の直径が %g のとき，等円の直径は %g である。\n", 2R, 2r)
    @printf("R = %g;  r = %g;  x1 = %g;  y1 = %g;  y2 = %g;  x0 = %g\n", R, r, x1, y1, y2, x0)
    plot()
    circle(0, 0, R)
    circle2(x1, y1, r, :blue)
    circle2(r, y2, r, :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)
        plot!([0, x0, 0, -x0, 0, 0], [R, -t, -R, -t, R, -R], color=:green, lw=0.7)
        point(x1, y1, "等円:r,(x1,y1)", :blue, :center, delta=-delta/2)
        point(r, y2, "等円:r,(r,y2)", :blue, :center, delta=-delta/2)
        point(x0, -t, " (x0,sqrt(R^2-x0^2)", :green, :left, :vcenter)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
    end
end;

draw(13/2, true)

算額（その1435）

算額（その1435）

福島県伊達市梁川町粟野堂内 粟野地蔵堂 明治16年(1883)
街角の数学 ～落書き帳「○△□」～ 43．ルージュの伝言
http://streetwasan.web.fc2.com/math15.6.25.html
キーワード：円3個，半円2個，長方形
#Julia, #SymPy, #算額, #和算

長方形の中に，大半円 2 個，中円 1 個，小円 2 個が入っている。小円の半径が与えられたとき，ピンクの部分の面積を求めよ。

中円の半径と中心座標を r2, (0, 0)
大円の半径と中心座標を r1, (0, r2), (0 -r2)
小円の半径と中心座標を r3, (r1 - r3, 0)
とおき，以下の連立方程式を解いて r1, r2 を得る。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms r1::positive, r2::positive, r3::positive
eq1 = (r1 - r3)^2 + r2^2 - (r1 + r3)^2
eq2 = 2r2^2 - r1^2
res = solve([eq1, eq2], (r1, r2))[1];

# r1 大円の半径
res[1] |> println

    8*r3

# r2 中円の半径
res[2] |> println

    4*sqrt(2)*r3

求める面積の 1/4 は第1象限にある A, B, Cで囲まれる図形である。
その図形は，「四分円 OAB」 - 「八分円 DAC」 + 「二等辺三角形 OAD」である。
OA = r2, DA = r1 より，π*r2^2/4 - πr1^2/8 + r2^2/2 = 16r3^2ゆえ，求める全体の面積はその 4 倍の 64r3^2 である。
r3 が 1 寸のとき，面積は 64 平方寸である。

(r1, r2) = (8*r3, 4*sqrt(Sym(2))*r3)
S = 4 * (PI*r2^2/4 - PI*r1^2/8 + r2^2/2)
S |> println

    64*r3^2

定積分を求めるのが簡単で間違いがない。

@syms x
S = 4 * integrate(sqrt(r2^2 - x^2) - sqrt(r1^2 - x^2) + r2, (x, 0, r2))
S |> println

    64*r3^2

function draw(r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2) = (8*r3, 4*sqrt(2)*r3)
    @printf("r3 = %g;  r1 = %g;  r2 = %g\n", r3, r1, r2)
    plot([r1, r1, -r1, -r1, r1], [-r2, r2, r2, -r2, -r2], color=:green, lw=0.5)
    circlef(0, 0, r2, :pink)
    circle(0, 0, r2)
    circlef(0, r2, r1, :white, beginangle=225, endangle=315)
    circlef(0, -r2, r1, :white, beginangle=45, endangle=135)
    circle(0, r2, r1, beginangle=180, endangle=360)
    circle(0, -r2, r1, beginangle=0, endangle=180)
    circle2(r1 - r3, 0, 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)
        point(0, 0, "O ", :black, :right, :vcenter)
        point(r2, 0, "A  ", :black, :right, :vcenter)
        point(0, r2, "B", :black, :center, :bottom, delta=delta)
        point(0, r1 - r2, "C ", :black, :right, :bottom, delta=delta)
        point(0, -r2, "D ", :black, :right, :bottom, delta=delta)
        point(r1, r2, "(r1,r2)", :green, :right, :bottom, delta=delta)
        point(r1 - r3, 0, "小円\n\nr3", :blue, :center, :vcenter)
        vline!([0], color=:gray70, lw=0.5)
    end
end;

draw(1, true)

算額（その1434）

算額（その1434）

福島県伊達市梁川町二野袋庭渡 庭渡神社 明治27年(1894)
街角の数学 ～落書き帳「○△□」～ 40．正方形と四分円（5）
http://streetwasan.web.fc2.com/math15.6.22.html
キーワード：円2個，四分円2個，長方形，斜線2本
#Julia, #SymPy, #算額, #和算

長方形の中に 2 個の四分円を容れ，長方形の頂点から接戦を引く。区画された領域に乾円，坤円を 1 個ずつ容れる。坤円の半径を与えて，乾円の半径を求めよ。

長方形の長辺，短辺を 2a, b
接線と長方形の長辺の交点座標を (c, 0)
四分円の半径と中心座標を b, (a, b), (-a, b)
乾円の半径と中心座標を r1, (0, b - r1)
坤円の半径と中心座標を r2, (0, r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, b::positive, c::positive, r1::positive, r2::positive
eq1 = a^2 + r1^2 - (r1 + b)^2
eq2 = dist2(-a, b, c, 0, 0, b - r1, r1)
eq3 = dist2(-a, b, c, 0, 0, r2, r2)
eq4 = dist2(-a, b, c, 0, a, b, b);
res = solve([eq1, eq2, eq3, eq4], (a, b, c, r1))[1]  # 1 of 2

    (-r2*sqrt(11 + 8*sqrt(2))*sqrt(56*sqrt(2) + 81)/6 + r2*sqrt(11 + 8*sqrt(2))/3 + sqrt(2)*r2*sqrt(11 + 8*sqrt(2))*sqrt(56*sqrt(2) + 81)/6, 3*sqrt(2)*r2/2 + 5*r2/2, r2*sqrt(11 + 8*sqrt(2)), -r2*sqrt(112*sqrt(2) + 162)/9 - r2*sqrt(56*sqrt(2) + 81)/9 - 5*r2/18 + sqrt(2)*r2/18 + sqrt(2)*r2*sqrt(56*sqrt(2) + 81)/9 + sqrt(2)*r2*sqrt(112*sqrt(2) + 162)/9)

# a
res[1] |> simplify |> println
res[1](r2 => 1).evalf() |> println

    r2*(-sqrt(1787 + 1264*sqrt(2)) + 2*sqrt(11 + 8*sqrt(2)) + sqrt(3574 + 2528*sqrt(2)))/6
    5.70205716976694

# b
res[2] |> simplify |> println
res[2](r2 => 1).evalf() |> println

    r2*(3*sqrt(2) + 5)/2
    4.62132034355964

# c
res[3] |> println
res[3](r2 => 1).evalf() |> println

    r2*sqrt(11 + 8*sqrt(2))
    4.72373882628843

# r1
res[4] |> simplify |> sympy.sqrtdenest |> simplify |> println
res[4](r2 => 1).evalf() |> println

    r2*(1 + sqrt(2))/2
    1.20710678118655

乾円の半径 r1 は，坤円の半径 r2 の (1 + √2)/2 倍である。かなりきれいな倍数だと思う。
坤円の半径 r2 が 1 のとき，乾円の半径 r1 は (1 + √2)/2 = 1.2071067811865475 である。

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

    r2 = 1;  r1 = 1.20711;  a = 5.70206;  b = 4.62132;  c = 4.72374

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    a = r2*(-sqrt(1787 + 1264√2) + 2*sqrt(11 + 8√2) + sqrt(3574 + 2528√2))/6
    b = r2*(3√2 + 5)/2
    c = r2*sqrt(11 + 8√2)
    r1 = r2*(1 + √2)/2
    @printf("r2 = %g;  r1 = %g;  a = %g;  b = %g;  c = %g\n", r2, r1, a, b, c)
    plot([a, a, -a, -a, a], [0, b, b, 0, 0], color=:green, lw=0.5)
    circle(a, b, b, :magenta, beginangle=180, endangle=270)
    circle(-a, b, b, :magenta, beginangle=270, endangle=360)
    circle(0, b - r1, r1)
    circle(0, r2, r2, :blue)
    segment(-a, b, c, 0, :orange)
    segment(a, b, -c, 0, :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(a, 0, "a", :green, :center, delta=-delta)
        point(c, 0, "c", :black, :center, delta=-delta)
        point(a - b, b, "(a-b,b)", :magenta, :center, :bottom, delta=delta)
        point(0, b - r1, "乾円:r1\n(0,b-r1)", :red, :center, delta=-delta)
        point(0, r2, "坤円:r2\n(0,r2)", :blue, :center, delta=-delta)
        ylims!(-8delta, b + 2delta)
    end
end;

draw(1, true)

算額（その1433）

算額（その1433）

（福島県福島市岡部 春日神社 明治28年(1895)）
街角の数学 ～落書き帳「○△□」～ 39．正方形と四分円（4）
http://streetwasan.web.fc2.com/math15.6.20.html
キーワード：円2個，四分円2個，正方形2個
#Julia, #SymPy, #算額, #和算

外正方形（外方）と 2 つの四分円とで囲まれた上下の隙間に，小円と大円が内接している。
外正方形の一辺の長さを 1 寸として大円と小円の半径をそれぞれ求めよ。

外正方形の一辺の長さ（=四分円の半径）を R
内正方形の一辺の長さを a
大円の半径と中心座標を r1, (0, r1)
小円の半径と中心座標を r2, (R - r2)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, R::positive, r1::positive, r2::positive
eq1 = (R/2)^2 + r1^2 - (R - r1)^2
eq2 = (R/2)^2 + (R - r2)^2 - (R + r2)^2
eq3 = (a/2 + R/2)^2 + a^2 - R^2
res = solve([eq1, eq2, eq3], (r1, r2, a))[1];

# r1
res[1] |> println
res[1](R => 1) |> println

    3*R/8
    3/8

# r2
res[2] |> println
res[2](R => 1) |> println

    R/16
    1/16

# a
res[3] |> println
res[3](R => 1) |> println

    3*R/5
    3/5

外方の一辺の長さが R のとき，大円の半径 r1 はその 3/8 倍，小円の半径 r2 は 1/16 倍，内方の一辺の長さ a は 3/5 倍である。
外方の一辺の長さが 1 のとき，大円の直径は 0.75，小円の直径は 0.125，, 内方の一辺の長さは 0.6 である。
外方の一辺の長さが 80 のとき，大円の直径は 60，小円の直径は 10，内方の一辺の長さは 48 である。

四分円と大円の接点は，内方の頂点と一致する。

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, r1, r2) = (3*R/5, 3*R/8, R/16)
    @printf("外方の一辺の長さが %g のとき，大円の直径は %g，小円の直径は %g，, 内方の一辺の長さは %g である。\n", R, 2r1, 2r2, a)
    plot([R/2, R/2, -R/2, -R/2, R/2], [0, R, R, 0, 0], color=:green, lw=0.5)
    circle(-R/2, 0, R, beginangle=0, endangle=90)
    circle(R/2, 0, R, beginangle=90, endangle=180)
    circle(0, r1, r1, :orange)
    circle(0, R - r2, r2, :blue)
    plot!([a, a, -a, -a, a] ./2, [0, a, a, 0, 0], color=:magenta, lw=0.5)
    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/2, 0, " a/2", :magenta, :left, :bottom, delta=delta/2)
        point(0, R - r2, " 小円:r2,(0,R-r2)", :blue, :center, delta=7delta)
        point(0, r1, "大円:r1,(0,r1)", :orange, :center, delta=-delta/2)
        point(a/2, a, " (a/2,a)", :magenta, :left, :vcenter)
        point(R/2, R, "(R/2,R) ", :red, :right, :bottom, delta=delta/2)
    end
end;

draw(1, true)

 

算額（その1432）

算額（その1432）

（福島県田村市船引町石森戸屋 稲荷神社 明治18年(1885)）
街角の数学 ～落書き帳「○△□」～ 36．正方形と四分円（2）
http://streetwasan.web.fc2.com/math15.6.18.html
キーワード：四分円2個，正方形3個
#Julia, #SymPy, #算額, #和算

算額（その1431）の拡張で，正方形をもう 1 個加える。
正方形の中に四分円を 2 個を容れ，交差した部分に正方形を 2 個容れる。外側の正方形の一辺の長さが 5 寸のとき，内側の一番小さな正方形の一辺の長さはいかほどか。

外側の正方形の左下の頂点を原点とし，一辺の長さを a
内側の正方形の右下の頂点の座標を (b, 0)
内側の一番小さな正方形の右下の頂点の座標を (c, 2(b - a/2))
とおき，以下の連立方程式を解く。
内側の 2 つの正方形の右上の頂点いずれもが，原点を中心とする四分円の円周上にあるということを利用する。
内側の正方形の一辺の長さは a - 2(a - b) = 2b - a である。
内側の一番小さい正方形の一辺の長さは a - 2(a - c) = 2c - a である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, b::positive, c::positive
eq1 = b^2 + (2b - a)^2 - a^2
eq2 = c^2 + ((2b - a) + (2c - a))^2 - a^2
res = solve([eq1, eq2], (b, c))[1]
res |> println

    (4*a/5, 3*a/5)

# b
res[1](a => 5) |> println

    4

# c
res[2](a => 5) |> println

    3

a = 5 のとき b = 4, c = 3 なので，内側の一番小さな正方形の一辺の長さは 2c - a = 1 寸である。

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    b = 4a/5
    c = 3a/5
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    plot!([a - b, b, b, a - b], [0, 0, 2(b - a/2), 2(b - a/2), 0], color=:blue, lw=0.5)
    plot!([a - c, c, c, a - c], 2(b - a/2) .+ [0, 0, 2(c - a/2), 2(c - a/2), 0], color=:magenta, lw=0.5)
    circle(0, 0, a, beginangle=0, endangle=90)
    circle(a, 0, a, beginangle=90, endangle=180)
    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(b, 0, " b", :blue, :left, :bottom, delta=delta/2)
        point(b, 2(b - a/2), "(b,2(b-a/2))", :blue, :left, :bottom, delta=delta/2)
        point(c, 2(b - a/2) + 2(c - a/2), "(c,2(b-a/2)+2(c-a/2))", :magenta, :left, :bottom, delta=delta/2)
    end
end;

draw(5, true)

算額（その1431）

算額（その1431）

福島県田村市船引町石森戸屋 稲荷神社 明治18年(1885)
街角の数学 ～落書き帳「○△□」～ 36．正方形と四分円
http://streetwasan.web.fc2.com/math15.6.17.html
キーワード：四分円2個，正方形2個
#Julia, #SymPy, #算額, #和算

正方形の中に四分円を 2 個を容れ，交差した部分に正方形を 1 個容れる。外側の正方形の一辺の長さが 5 寸のとき，内側の正方形の一辺の長さはいかほどか。

外側の正方形の左下の頂点を原点とし，一辺の長さを a
内側の正方形の右下の頂点の座標を (b, 0)
とおき，以下の連立方程式を解く。
内側の正方形の右上の頂点が原点を中心とする四分円の円周上にあるということを利用する。
内側の正方形の一辺の長さは a - 2(a - b) = 2b - a である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a::positive, b::positive
eq = b^2 + (2b - a)^2 - a^2
res = solve(eq, b)[1]
res |> println
res(a => 5) |> println

    4*a/5
    4

a = 5 のとき b = 4 なので，内側の正方形の一辺の長さは 2b - a = 3 寸である。

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    b = 4a/5
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    plot!([a - b, b, b, a - b], [0, 0, 2(b - a/2), 2(b - a/2), 0], color=:blue, lw=0.5)
    circle(0, 0, a, beginangle=0, endangle=90)
    circle(a, 0, a, beginangle=90, endangle=180)
    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(b, 0, " b", :blue, :left, :bottom, delta=delta/2)
        point(b, 2(b - a/2), "(b,2(b-a/2))", :blue, :left, :bottom, delta=delta/2)
    end
end;

draw(5, true)

算額（その1430）

算額（その1430）

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円12個，外円，正三角形3個
#Julia, #SymPy, #算額, #和算

正三角形の中に斜線を 2 本，全円 1 個，甲円 1 個，乙円 3 個，全円 1 個を容れる。甲円の直径が与えられたとき，丙円の直径を得る術を述べよ。

正三角形の一辺の長さを 2a
甲円の半径と中心座標を r1, (0, 2r0 - r1)
乙円の半径と中心座標を r2, (x2, y2), (0, 2r0 + r2)
丙円の半径と中心座標を r3, (x3, y3)
全円の半径と中心座標を r0, (0, r0)
とおき，以下の連立方程式を解く。

なお，r0 = √3a/3, r2 = √3a/9, y3 = y2 である。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cfusing SymPy
@syms a, r0, r1, r2, x2, y2, r3, x3, y3, b
# eq1 = r0 - a/√Sym(3)
# eq2 = 2r2 - (√Sym(3)a - 2r0 - r2)
# res = solve([eq1, eq2], (r0, r2))
# r0 = √Sym(3)a/3
# r2 = √Sym(3)a/9

r0 = √Sym(3)a/3
r2 = √Sym(3)a/9
y3 = y2
# dist2() は使えない
eq3 = dist(-a, 0, b, √Sym(3)*(a - b), x2, y2) - r2^2
eq4 = dist(-a, 0, b, √Sym(3)*(a - b), 0, 2r0 - r1) - r1^2
eq5 = dist(a, 0, -b, √Sym(3)*(a - b), x2, y2) - r2^2
eq6 = dist(a, 0, 0, √Sym(3)a, x2, y2) - r2^2
eq7 = dist(-a, 0, b, √Sym(3)*(a - b), x3, y3) - r3^2
eq8 = (x2 - x3) - (r2 + r3);

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)
    (x2, y2, r1, b, r3, x3) = u
    return [
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x2 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x2))/(3*(a - b)^2 + (a + b)^2))^2,  # eq3
-r1^2 + (a - (a + b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r1))/(3*(a - b)^2 + (a + b)^2))^2 + (2*sqrt(3)*a/3 - r1 - sqrt(3)*(a - b)*(a*(a + b) + sqrt(3)*(a - b)*(2*sqrt(3)*a/3 - r1))/(3*(a - b)^2 + (a + b)^2))^2,  # eq4
-a^2/27 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/((-a - b)^2 + 3*(a - b)^2))^2 + (-a + x2 - (-a - b)*(sqrt(3)*y2*(a - b) + (-a - b)*(-a + x2))/((-a - b)^2 + 3*(a - b)^2))^2,  # eq5
-a^2/27 + (y2 - sqrt(3)*a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2 + (-a + x2 + a*(sqrt(3)*a*y2 - a*(-a + x2))/(4*a^2))^2,  # eq6
-r3^2 + (y2 - sqrt(3)*(a - b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2 + (a + x3 - (a + b)*(sqrt(3)*y2*(a - b) + (a + b)*(a + x3))/(3*(a - b)^2 + (a + b)^2))^2,  # eq7
-sqrt(3)*a/9 - r3 + x2 - x3,  # eq8
    ]
end;

a = 10
r0 = √3a/3  # a*0.5773502691896257
r2 = √3a/9  # a*0.19245008972987523
iniv = BigFloat[0.575411, 0.34912, 0.387466, 0.662547, 0.095137, 0.284131] .* a
res = nls(H, ini=iniv)

    ([5.711594854986474, 3.5787338000056996, 3.8641388667743293, 6.575250125060548, 0.9544492311458473, 2.832644726541875], true)

正三角形の一辺の長さが 2a = 20 のとき，
全円の半径 r0 = 5.773502691896257
甲円の半径 r1 = 3.8641388667743293
乙円の半径 r2 = 1.9245008972987523
丙円の半径 r3 = 0.9544492311458473

丙円の直径/甲円の直径 = 0.9544492311458473/3.8641388667743293 = 0.2470017936862072
甲円の直径/丙円の直径 = 3.8641388667743293/0.9544492311458473 = 4.048553595810754

術は，「15572 の平方根から 121 を引き，甲円の直径を掛ければ全円の直径がえられる」とあるが，√15572 - 121 = 3.787819918451973 であり，題意を満たさない。

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (x2, y2, r1, b, r3, x3) = [5.711594854986474, 3.5787338000056996, 3.8641388667743293, 6.575250125060548, 0.9544492311458473, 2.832644726541875]
    y3 = y2
    r0 = √3a/3
    r2 = √3a/9
    plot([a, 0, -a, a], [0, sqrt(3)a, 0, 0], color=:green, lw=0.5)
    circle(0, a/sqrt(3), r0)
    circle(0, 2r0 + r2, r2, :blue)
    circle2(x2, y2, r2, :blue)
    circle2(x3, y3, r3, :brown)
    circle(0, 2r0 - r1, r1, :orange)
    segment(-a, 0, b, (a - b)*√3)
    segment(a, 0, -b, (a - b)*√3)
    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, √3a, " √3a", :green, :left, :vcenter)
        point(a, 0, "a", :green, :left, :bottom, delta=delta/2)
        point(0, r0, "全円:r0,(0,r0)", :red, :center, delta=-delta)
        point(0, 2r0 - r1, "甲円:r1,(0,2r0-r1)", :orange, :center, delta=-delta)
        point(x2, y2, "乙円:r2\n(x2,y2)", :blue, :left, delta=-delta, deltax=-3delta)
        point(0, 2r0 + r2, " 乙円:r2\n(x2,r2)", :blue, :center, delta=-delta)
        point(x3, y3, "丙円:r3,(x3,y3)", :brown, :right, delta=-4delta, deltax=-delta)
    end
end;

draw(10, true)

橋本製麺所

高松市仏生山町 橋本製麺所
旧二級国道沿いにあるが，店名表示も暖簾も表札すらもない。何百回となく店の前を通っているのに，製麺所とは知らなかった。
いわゆるうどん店ではない。うどん玉1個90円で売っており，丼と箸を持っていけば店内で食べることもできる（店内に数脚のパイプ椅子があるがテーブルはない）。
家に帰り，スーパーで買ってきた天ぷらを乗せて掛けうどんでいただく。

算額（その1429）

算額（その1429）

八十九 陸前高田市小友町 常膳寺観音堂 天保13年(1842)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円6個，半円1個，円弧1個
#Julia, #SymPy, #算額, #和算

半円と円弧に挟まれた領域に，甲円，乙円，丙円を 2 個ずつ容れる。丙円の直径が 4 寸のとき，甲円の直径はいかほどか。

半円の半径は甲円の半径の 4 倍である。
円弧の半径と中心座標を r0, (0, 0)
円弧と半円の共通弦と y 軸の交点座標を (0, y); y = r0 - 2r1
甲円の半径と中心座標を r1, (0, y + r1), (0, (y + 3r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
とおき，以下の連立方程式を解く。

include("julia-source.txt");
# julia-source.txt ソース　https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf

using SymPy
@syms r0::positive, r1::positive, r2::positive, x2::positive,
      r3::positive, x3::positive, y::positive
y = r0 - 2r1
eq1 = x2^2 + (y +  r2)^2 - (r0 - r2)^2
eq2 = x3^2 + (y +  r3)^2 - (r0 - r3)^2
eq3 = x2^2 + (r1 - r2)^2 - (r1 + r2)^2
eq4 = (x3 - x2)^2 + (r2 - r3)^2 - (r2 + r3)^2
eq5 = (4r1)^2 + y^2 - r0^2
res = solve([eq1, eq2, eq3, eq4, eq5], (r0, r1, r2, x2, x3))[1];

# r0
res[1] |> println

    45*r3/4

# r1
res[2] |> println

    9*r3/4

# r2
res[3] |> println

    9*r3/5

# x2
res[4] |> println

    9*sqrt(5)*r3/5

# x3
res[5] |> println

    3*sqrt(5)*r3

甲円の半径 r1 は，丙円の半径 r3 の 9/4 倍である。
丙円の直径が 4 寸のとき，甲円の直径は 4*9/4 = 9 寸である。

答，術ともに「甲円径六寸」としているが，誤りである。

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

    r3 = 2;  r0 = 22.5;  r1 = 4.5;  r2 = 3.6;  x2 = 8.04984;  x3 = 13.4164

function draw(r3, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r0, r1, r2, x2, x3) = r3 .* (45/4, 9/4, 9/5, 9√5/5, 3√5)
    @printf("乙円の直径が %g のとき，甲円の直径は %g である。\n", 2r3, 2r1)
    @printf("r3 = %g;  r0 = %g;  r1 = %g;  r2 = %g;  x2 = %g;  x3 = %g\n", r3, r0, r1, r2, x2, x3)
    y = r0 - 2r1
    θ = atand(y, 4r1)
    plot()
    circle(0, y, 4r1, beginangle=0, endangle=180, n=500)
    circle(0, 0, r0, :blue, beginangle=θ, endangle=180 - θ, n=500)
    circle(0, y + r1, r1, :magenta)
    circle(0, y + 3r1, r1, :magenta)
    circle2(x2, y + r2, r2, :green)
    circle2(x3, y + r3, r3, :orange)
    segment(-4r1, y, 4r1, y)
    segment(0, 0, 4r1, y, :skyblue)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /  3  # size[2] * fontsize * 2
        vline!([0], color=:gray80, lw=0.5)
        hline!([0], color=:gray80, lw=0.5)
        point(0, y + r1, "甲円:r1\n(0,y+r1)", :magenta, :center, delta=-delta)
        point(0, y + 3r1, "甲円:r1\n(0,y+3r1)", :magenta, :center, delta=-delta)
        point(x2, y + r2, "乙円:r2\n(x2,y+r2)", :green, :center, delta=-delta)
        point(x3, y + r3, "丙円:r3,(x3,y+r3)", :orange, :center, delta=-5delta, deltax=-3delta)
        point(4r1, y, "(4r1,y)", :red, :left, :bottom, delta=delta)
        point(0, y + 4r1, "y+4r1", :red, :center, :bottom, delta=delta)
        xlims!(-4r1 - 3delta, 4r1 + 10delta)
    end
end;

draw(4/2, true)