算額（その1491）

算額（その1491）

五十七 岩手県花泉町 花泉天満宮 文政13年(1830)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円14個，直線上
#Julia, #SymPy, #算額, #和算

直線上に，甲円 7 個を内包する大円，乙円 2 個，丙円 4 個が載っている。甲円の直径が 1 寸のとき，丙円の大きさはいかほどか。

大円の中に甲円を 7 個容れるのは本筋ではない。以下に示す「算法助術の公式53」は，大きさの異なる 3 種類の円の直径の間に成り立つ関係式である。この公式は，2 種の円の直径が分かっているとき，残りの 1 種の円の直径を求めるために使う。

本問と対応付けると，大円が O₂，乙円が O₃ 丙円が O₁ で，甲円を容れた後に大円の直径は決まるが，乙円が決まらない限り丙円は決まらない。
つまり，本問は条件不足で解けないというのが結論である。

試みに 乙円の直径を何通りか変えて図を描いてみる。大円とその中の甲円の大きさは同じでも，乙円の大きさの違いにより丙円の大きさも変わるのが明らかである。

元の算額では読み取れない文字があり，答曰「丙円径六□□□□余□」とあり，山村はこれを「6分9厘9毛」と解釈して「術」を解説しているが，そもそも不毛な解釈である。

最初に示した図は，甲円の直径が 1，乙円の直径が 0.26 のときのものであり，丙円の直径がほぼ 0.7 （山村の解釈に近いもの）となっている。しかし，いつも通り，算額の図とは少し違う。

大円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (r1, y1), (r1, y1 + 2r1), (x1, y2), (0, R1 - R); y2 < 0
乙円の半径と中心座標を r4, (2r3, y4); y4 < 0
丙円の半径と中心座標を r3, (r3, y3), (3r3, y3); y3 < 0
とおき，r1, r4 が既知のとき，R, x1, y1, y2, r3, y3, y4 を未知数として，以下の連立方程式を解く。

甲円，乙円の直径が 1, 0.26 のとき，丙円の直径は 0.700673 である。大円の直径は 3.4305 である。

追記：

https://w.atwiki.jp/sangaku/pages/166.html に，同じ算額図があるが，全く異なったものになっている。（山村の図は分が悪いが）
確かに，山村の問では各円の個数が読み取り不能で □ になっており，山村の解釈で個数を書き入れたものとも取れる。

今有線上如圖設大圓容附甲圓（五ケ）乙圓（六ケ）丙円（二ケ）只云甲圓径一寸答丙圓径幾何
答曰丙圓径六分七厘九毛余
術曰置一個（名天）九七個開平方減天一十六除之（名地）天興地二段和因天開平方（名人）天地興人二段和（名角）天地和乗人二段減天地差因地（名元）乗角加天興地二段和因地冪（名氐）人冪因地冪加充冪（名房）乗角冪内減氐冪開平方加氐以房除之乗甲圓径因地得丙圓径合問

「額文は現物によるが、不明箇所は「和算 岩手の現存算額のすべて」を参考とした。」と記載されており，（術の一部も）山村のものとは異なる。
「和算 岩手の現存算額のすべて」は，「和算-岩手の現存算額のすべて 安富有恒 青磁社 1987」であると思う。

include("julia-source.txt");

using SymPy

@syms R::positive, r1::positive, x1::positive, y1::positive, y2::negative
# r1 = 1
eq1 = x1^2 + y2^2 - (R - r1)^2 |> expand
eq2 = r1^2 + (y1 + 2r1)^2 - (R - r1)^2 |> expand
eq3 = (x1 - r1)^2 + (y1 - y2)^2 - 4r1^2 |> expand
eq4 = x1^2 + (y2 - r1 + R)^2 - 4r1^2 |> expand;
# res = solve([eq1, eq2, eq3, eq4], (R, x1, y1, y2))

@syms r3, y3, r4, y4
eq5 = 4r3^2 + y4^2 - (R + r4)^2 |> expand
eq6 = r3^2 + y3^2 - (R + r3)^2 |> expand
eq7 = r3^2 + (y4 - y3)^2 - (r3 + r4)^2 |> expand
#res2 = solve([eq5, eq6, eq7], (r3, y3, y4))

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 driver(r4)
    function H(u)
        (R, x1, y1, y2, r3, y3, y4) = u
        return [
            -R^2 + 2*R*r1 - r1^2 + x1^2 + y2^2,  # eq1
            -R^2 + 2*R*r1 + 4*r1^2 + 4*r1*y1 + y1^2,  # eq2
            -3*r1^2 - 2*r1*x1 + x1^2 + y1^2 - 2*y1*y2 + y2^2,  # eq3
            R^2 - 2*R*r1 + 2*R*y2 - 3*r1^2 - 2*r1*y2 + x1^2 + y2^2,  # eq4
            -R^2 - 2*R*r4 + 4*r3^2 - r4^2 + y4^2,  # eq5
            -R^2 - 2*R*r3 + y3^2,  # eq6
            -2*r3*r4 - r4^2 + y3^2 - 2*y3*y4 + y4^2,  # eq7
        ]
    end;
    r1 = 1/2
    iniv = BigFloat[3.41, 1.8, 0.24, -1.6, 1.172, -4.476, -3.167]
    res = nls(H, ini=iniv)
    return res[1]
end;

function draw(r1, r4, more=false)
    pyplot(size=(500, 400), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, x1, y1, y2, r3, y3, y4) = driver(r4)
    # r1 = 1 のときは以下が成り立つ
    # (R, x1, y1, y2) = (5/3 + 2*sqrt(7)/3, 1/2 + sqrt(7)/2, -2/3 + sqrt(7)/3, -7/6 - sqrt(7)/6)
    @printf("甲円，乙円の直径が %g, %g のとき，丙円の直径は %g である。大円の直径は %g である。", 2r1, 2r4, 2r3, 2R)
    p = plot()
    circle(0, 0, R, :green)
    circle2(r1, y1, r1)
    circle2(r1, y1 + 2r1, r1)
    circle2(x1, y2, r1)
    circle(0, r1 - R, r1)
    circle2(r3, y3, r3, :blue)
    circle2(3r3, y3, r3, :blue)
    circle2(2r3, y4, r4, :magenta)
    str = replace(@sprintf("r4/r1 = %s, r3 = %g", rationalize(r4/r1), r3), "//" => "/")
    point(0, 2.15, str, :magenta, :center, mark=false)
    hline!([0], color=:gray80, lw=0.5)
    vline!([0], color=:gray80, lw=0.5)
    if more        
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        point(r1, y1 + 2r1, "r1", :red, :center, :vcenter, mark=false)
        point(2r3, y4, "r4", :magenta, :center, :vcenter, mark=false)
        point(r3, y3, "r3", :blue, :center, :vcenter, mark=false)
        point(R, 0, " R", :green, :left, :bottom, delta=delta/2)
    end
    plot!(xlims=(-3.5, 3.5), ylims=(-3.5, 2))  #, titlefontsize=10, title=@sprintf("r4/r1 = %s", rationalize(r4/r1)))
    return p
end;

draw(1/2, 0.13, true)