算額（その1543）

算額（その1543）

千葉県富津市寺尾 六所神社 明治4年(1871)
山口正義：やまぶき 第 57 号
https://yamabukiwasan.sakura.ne.jp/ymbk57.pdf
山口正義：やまぶき 第 59 号
https://yamabukiwasan.sakura.ne.jp/ymbk59.pdf

キーワード：球9個，3次元
#Julia, #SymPy, #算額, #和算

外球の中に甲球 1 個，乙球 7 個を容れる。外球の直径が 12 寸，甲球の直径が 8 寸のとき，乙球の直径はいかほどか。

外球の半径と中心座標を R, (0, 0, 0)
甲球の半径と中心座標を r1, (0, 0, R - r1)
乙球の半径と中心座標を r2, (x2, y2, z2)
とおき，以下の連立方程式を解く。

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

@syms R, r1, r2, x2, y2, z2
eq1 = x2^2 + (R - r1 - z2)^2 - (r1 + r2)^2
eq2 = x2^2 + z2^2 - (R - r2)^2
eq3 = x2*sind(Sym(360)/14) - r2
res = solve([eq1, eq2, eq3], (r2, x2, z2))[2]  # 2 of 2

    (2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))

# 2r2
2res[1](R => 12/2, r1 => 8/2).evalf() |> println

    3.27511575021090

外球，甲球の直径がそれぞれ 12 寸，8 寸のとき，乙球の直径は 3.27511575021090 である。
答・術は 4.7071 寸としているが，山口も 3.275 有奇として異を唱えている。

function draw(R, r1, more=false)
    pyplot(size=(500, 250), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, x2, z2) = (2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2), 2*R*r1*(-R*cos(2*pi/7) + R - r1 + r1*cos(2*pi/7))/((R^2 - 2*R*r1*cos(2*pi/7) + r1^2)*sin(pi/7)), R*(R^2 - 2*R*r1 - r1^2 + 2*r1^2*cos(2*pi/7))/(R^2 - 2*R*r1*cos(2*pi/7) + r1^2))
    p1 = plot()
    circle(0, 0, R, :cyan)
    circle(0, 0, x2, :blue)
    rotate(x2, 0, r2, angle=360/7)        
    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, :bottom, delta=delta/2)
        point(x2, 0, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
    end
    p2 = plot()
    circle(0, 0, R, :cyan)
    circle(0, R- r1, r1, :blue)
    circle(x2, z2, r2)
    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, :bottom, delta=delta/2)
        point(0, R - r1, "甲球:(0,0, R-r1)", :blue, :center, delta=-delta/2)
        point(x2, z2, "乙球:r2,(x2,0,z2)", :red, :right, delta=-delta/2, deltax=15delta)
    end
    plot(p1, p2)
end;

draw(12/2, 8/2, true)

算額（その1542）

算額（その1542）

四十 埼玉県熊谷市玉井 玉井大神社 嘉永元年(1848)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

山口正義：やまぶき 第 4 号
https://yamabukiwasan.sakura.ne.jp/ymbk4.pdf

キーワード：正五角形，面積分割
#Julia, #SymPy, #算額, #和算

正五角形の一辺が 1 寸のとき，この面積を三等分するときの截面（一辺の分割した長さ；図参照）はいかほどか。

計算しなくても，図を見ればわかる。
四角形OAEB = 五角形OECDE' ならば，三角形OBE が三角形OCE の二倍，すなわち CE が BC の 1/3 である。

あえて計算すると，以下のようになる。
正五角形が内接する円の半径を r，一辺の長さを「五角面」，截面を「截面」とおき，以下の連立方程式を解く。

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

@syms r, x, y, 五角面
s18 = Sym(18)
s36 = Sym(36)
r = 五角面/2sind(s36)
S = r^2*sind(s36)*cosd(s36)  # 正五角形の面積/5
s1 = area([0 0; r*cosd(s18) r*sind(s18); x y])  # 面積/5 の更に 2/3
s2 = area([0 0; x y; r*sind(s36) -r*cosd(s36)])  # 面積/5 の更に 1/3
eq1 = s1 - 2S/3  # S/5 + s1 = S/5 + s2
eq2 = s2 - S/3;
res = solve([eq1, eq2], (x, y));

@syms d
ans_x = apart(res[x]) |> simplify
ans_x |> println

    五角面*(sqrt(5) + 5)/12

ans_y = apart(res[y]) |> simplify
ans_y|> println

    -五角面*(5*sqrt(2) + 3*sqrt(10))*sqrt(sqrt(5) + 5)/120

截面 = sqrt((r*sind(s36) - ans_x)^2 + (-r*cosd(s36) - ans_y)^2) |> simplify
截面 |> println
截面(五角面 => 1).evalf() |> println

    sqrt(10)*sqrt(-15*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) - 6*sqrt(5)*sqrt(5 - sqrt(5))*sqrt(sqrt(5) + 5) + 20*sqrt(5) + 110)*sqrt(五角面^2)/(30*sqrt(5 - sqrt(5)))
    0.333333333333333

截面は五角面の 1/3 である。

function draw(五角面, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = 五角面/2sind(36)
    x = 五角面*(√5 + 5)/12
    y = -五角面*(5√2 + 3√10)*sqrt(√5 + 5)/120
    (x1, y1) = r.*(sind(36), -cosd(36))
    (x2, y2) = r.*(cosd(18), sind(18))
    length = sqrt((x - x1)^2 + (y - y1)^2)
    @printf("正五角形の一辺の長さが %g のとき，面積を 3 分割するときの\"截面\"は %g である。\n", 五角面, length)
    plot([-x, 0, x], [y, 0, y], color=:blue, lw=0.5, showaxis=false)
    polygon(0, 0, r, 5, color=:red)
    segment(0, 0, 0, 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)
        point(x, y, " (x,y)", :blue, :left, :vcenter)
        point(x1, y1, " (x1,y1)", :red, :left, :bottom, delta=delta/2)
        dimension_line(-x1, y1 - delta, x1, y1 - delta, "五角面", :red, :center, delta=-2delta)
        dimension_line(-x1 - 1.5delta, y1 - delta/2, -x - 1.5delta, y - delta/2, "截面", :magenta, :right)
        segment(0, 0, -x2, y2, :gray80)
        segment(0, 0, -x1, y1, :gray80)
        segment(0, 0, x1, y1, :gray80)
        point(0, 0, " O", :blue, :left, :bottom, delta=delta/2)
        point(0, r, "A", :blue, :center, :bottom, delta=delta/2)
        point(-x2, y2, "B ", :blue, :right, :vcenter)
        point(-x1, y1, " C", :blue, :left, :bottom, delta=delta/2)
        point(x1, y1, "D ", :blue, :right, :bottom, delta=delta/2)
        point(-x, y, "E  ", :blue, :right, :bottom, delta=delta/2)
        point(x, y, "E' ", :blue, :right, delta=-delta/2)
    end
end;

draw(1, true)

算額（その1541）

算額（その1541）

埼玉県本庄市都島 正観寺 戸塚盛政の算額 享保11年(1726)
山口正義：やまぶき 第 26 号
https://yamabukiwasan.sakura.ne.jp/ymbk26.pdf
山口正義：やまぶき 第 40 号
https://yamabukiwasan.sakura.ne.jp/ymbk40.pdf

甲，乙，丙，丁，戊の立方体がある。甲と乙の体積の和が 189，丙と丁と戊の体積の和が 36 である。辺の長さの差は一定である。それぞれの一辺の長さはいかほどか。

立方体:甲の一辺の長さを「甲」，辺の差を「差」とおき，以下の連立方程式を解く。
問題文に出てくる 2 つの定数を記号のままとして解こうとすると，"RecursionError('maximum recursion depth exceeded')" が発生するので解けない。それほど複雑だということであろう。

using SymPy
@syms 甲::positive, 差::positive
eq1 = 甲^3 + (甲 - 差)^3 - 189
eq2 = (甲 -2差)^3 + (甲 -3差)^3 + (甲 -4差)^3 - 36
res = solve([eq1, eq2], (甲, 差))[1]

    (5, 1)

甲の一辺の長さは 5，差が 1 なので，順次 4, 3, 2, 1 となる。

山口の解説は「（まともにやれば）9 次方程式が得られるが，これをとくのはやっかいである」と書いているが，SymPy では何の苦も無く解ける（当たり前だが）。

丙，丁，戊の体積の和が 36 なので，丙^3 + 丁^3 + 戊^3 = 36, 丙 > 丁 > 戊, 各辺は整数なので，3, 2, 1 が解であることはすぐにわかり，5^3 + 4^3 = 189 となるのを確認すれば解を得るのは簡単である。直感で解いていはいけないのかな？

直感で解かれるのを避けるなら，以下のようなテストデータを作ればよい。
甲^3 + 乙^3 = 3699397765
丙^3 + 丁^3 + 戊^3 = 5120681355

using SymPy
@syms 甲::positive, 差::positive
eq1 = 甲^3 + (甲 - 差)^3 - 3699397765
eq2 = (甲 -2差)^3 + (甲 -3差)^3 + (甲 -4差)^3 - 5120681355
res = solve([eq1, eq2], (甲, 差))[1]

    (1234, 13)

甲の一辺の長さは 1234，差は 13 である。

1.　ブルートフォースで解く

甲の一辺の長さと差が未知なので，探索範囲を誤ると答えが出ないが，実行時間を見ながら範囲を広げていけば解が見つかる可能性は高い。
実行時間を短くするために色々策を弄することもできるが，まずは，単純なプログラムでやってみる。
以下のプログラムは 2 秒以内で正解を出す。

@time begin
    K1 = 3699397765
    K2 = 5120681355
    for 差 = 1:100
        for 甲 = 5:999999
            乙 = 甲 - 差
            丙 = 甲 - 2差
            丁 = 甲 - 3差
            戊 = 甲 - 4差
            戊 > 0 && 甲^3 + 乙^3 == K1 && 丙^3 + 丁^3 + 戊^3 == K2 && (println("甲 = $甲, 差 = $差"); break)
        end
    end
end

    甲 = 1234, 差 = 13
      1.396169 seconds (99.00 M allocations: 1.476 GiB, 5.47% gc time, 0.70% compilation time)