算額（その1521）

算額（その1521）

二十八 一関市萩荘 赤萩観音寺 天保2年(1831)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html

今有如図 03025
https://w.atwiki.jp/sangaku/pages/143.html

キーワード：円5個，直線上
#Julia, #SymPy, #算額, #和算

問：今有如図直線上下添大円其交容小円二個只謂小円径一寸問大円径幾何

直線の上に互いに外接する大円 2 個が載っており，その大円 2 個に下側から外接する大円が 1 個ある。
左右の大円と直線および下の大円が作る隙間に小円が左右に 2 個ある。
小円の直径が 1 寸のとき，大円の直径はいかほどか。

村山の図は本人自身も問がどのような図形を表現しているのかわからなかったのであろう。図には，あるはずの直線が描かれておらず，そんな描写はないにも関わらず 3 個の大円を内接する四辺形が描かれており（しかも四辺形は正方形であるとして「方」という文字を書き添え，結果として問題をややこしくしている），また，小円の位置は不適切だ。不適切にも程がある。

「今有如図」はちゃんとした図を示している。この図があれば，解を求めるのは簡単である。

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

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

using SymPy
@syms r1::positive, r2::positive, x2::positive
eq1 = x2^2 + ((r2 - r1) + √Sym(3)*r1)^2 - (r1 + r2)^2
# eq1 = x2^2 + ((r2 - r1) + √3*r1)^2 - (r1 + r2)^2  # 記号計算を行わない場合
eq2 = (x2 - r1)^2 + (r2 - r1)^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, x2))[1];

大円の半径は以下のようになる。桁数の大きい整数係数を持つ，二重根号を含む長い式になる。

res[1] |> println

    (-25413135*sqrt(3)*r2 + 44016841*r2 + (-28103451 + 16225535*sqrt(3))*(4*sqrt(3)*r2 + 7*r2 - 4*sqrt(2)*r2*sqrt(1939428460839052729369227845122062114242 - 1119729543939448663603914987253430016841*sqrt(3))/(164355571104683741015 - 94890733220103835857*sqrt(3))))*(4*sqrt(3)*r2 + 7*r2 - 4*sqrt(2)*r2*sqrt(1939428460839052729369227845122062114242 - 1119729543939448663603914987253430016841*sqrt(3))/(164355571104683741015 - 94890733220103835857*sqrt(3)))^2/(2*r2^2*(36060146 - 20819335*sqrt(3)))

SymPy では自動的に簡約化できないので，SymPy の手を借りながら手作業で簡約化する。
まず二重根号を外すと，4 つの項に分解される。

ans_r1 = res[1] |> sympy.sqrtdenest;

それぞれの項を個別に簡約化すると以下のようになる。

t1 = ans_r1.args[1] |> simplify
t1 |> println

    3008680048456320*sqrt(3)*r2/121 + 5211186707645184*r2/121

t2 = ans_r1.args[2] |> simplify
t2 |> println

    -3008680048441074*sqrt(3)*r2/121 - 5211186707618685*r2/121

t3 = ans_r1.args[3] |> simplify
t3 |> println

    52055161703399520*r2/121 + 30054061622167230*sqrt(3)*r2/121

t4 = ans_r1.args[4] |> simplify
t4 |> println

    -52055161703423478*r2/121 - 30054061622181024*sqrt(3)*r2/121

4 項を足し合わせ，更に簡約化すると， r1 = 3*r2*(4*sqrt(3) + 7) = r2*3(4√3 + 7) = r2*(12√3 + 21) = r2*(√432 + 21) となる。

最後の式は，術に述べられているものである。すなわち，「432 の平方根に 21 を加え，小円の直径を掛けると大円の直径が得られる」。

t1 + t2 + t3 + t4 |> simplify |> println

    3*r2*(4*sqrt(3) + 7)

大円の半径 r1 は，小円の半径 r2 の (√432 + 21) 倍である。
小円の直径が 1 寸のとき，大円の直径は √432 + 21 = 41.78460969082653 寸である。

もちろん方程式を解いて得られた解に直接 r2 = 1/2 を代入して 2 倍すれば，上のように迂遠な手順を踏まなくても大円の直径が得られる。

2res[1](r2 => 1/2).evalf() |> println

    41.7846096908265

なお，当初の方程式 eq1 を，eq1 = x2^2 + ((r2 - r1) + √3*r1)^2 - (r1 + r2)^2 とすれば，記号計算は行われず 41.7846096908267*r2 という解が得られる。

function draw(r2, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = 41.7846096908267*r2  # = r2*(√432 + 21)
    x2 = 28.8564064605511*r2
    @printf("小円の直径が %g のとき，大円の直径は %g である。\n", 2r2, 2r1)
    plot()
    circle2(r1, 0, r1, :blue)
    circle(0, -√3r1, r1, :blue)
    circle2(x2, r2 - r1, r2, :red)
    segment(-2r1, -r1, 2r1, -r1, :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(r1, 0, "大円:r1,(r1,0)", :blue, :center, :bottom, delta=delta)
        point(0, -√3r1, "大円:r1,(0,-√3r1)", :blue, :center, delta=-delta)
        point(x2, r2 - r1, "小円:r2,(x2,r2-r1)", :red, :left, delta=-delta, deltax=3delta, mark=false)
    end
end;

draw(1/2, true)