算額（その1414）

算額（その1414）

八四 加須市不動岡 総願寺 明治12年(1879)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円4個，正方形，斜線3本
#Julia, #SymPy, #算額, #和算

正方形の中に一つの頂点から対角線を含んで 3 本の斜線を引く。できた 4 個の区画に等円を入れる。

問題は判読不能ということであるが，算額（その1413）にも書いたように，どのような問題であるかは何通りか考えることができる。

一案として，「正方形の一辺の長さが与えられたとき，等円の直径を求める術（すべ）を述べよ。」を吟味しよう。

これは算額（その1205）と同じである。

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

using SymPy
@syms a::positive, b::positive, r::positive, y::positive
eq1 = a + b - sqrt(a^2 + b^2) - 2r
eq2 = dist2(a, 0, 0, b, r, y, r)
eq3 = dist2(a, 0, 0, a, r, y, r)
eq4 = dist2(a, 0, 0, b, r, r, r);

[eq1, eq2, eq3] を使う場合の解と [eq2, eq3, eq4] を使う場合の解では，数値にすると同じであるが，見かけの式は随分違う。

前者では r = a*(1 - sqrt(-1 + sqrt(2)))/2 であるが，
後者では r = a*(-22*sqrt(2) - 3*sqrt(-1 + sqrt(2)) + 33)/(2*(-7*sqrt(2) - 7*sqrt(2)*sqrt(-1 + sqrt(2)) + 8*sqrt(-1 + sqrt(2)) + 11)) であり，SymPy では簡約化できない。

res1 = solve([eq1, eq2, eq3], (r, y, b))[3]  # 3 of 3

    (a*(1 - sqrt(-1 + sqrt(2)))/2, -sqrt(2)*a/2 + sqrt(2)*a*sqrt(-2 + 2*sqrt(2))/4 + a*sqrt(-2 + 2*sqrt(2))/2 + a/2, a*(-20 - 2*sqrt(-2 + 2*sqrt(2)) + 3*sqrt(-1 + sqrt(2)) + 14*sqrt(2))/(-12*sqrt(-1 + sqrt(2)) - 4 + 3*sqrt(2) + 8*sqrt(2)*sqrt(-1 + sqrt(2))))

res1[1](a => 3973).evalf() |> println
res1[2](a => 3973).evalf() |> println
res1[3](a => 3973).evalf() |> println

    708.000016603060
    2263.73675775652
    1808.07201601373

res2 = solve([eq2, eq3, eq4], (r, y, b))[4]  # 4 of 4

    (a*(-22*sqrt(2) - 3*sqrt(-1 + sqrt(2)) + 33)/(2*(-7*sqrt(2) - 7*sqrt(2)*sqrt(-1 + sqrt(2)) + 8*sqrt(-1 + sqrt(2)) + 11)), a*(-3*sqrt(-1 + sqrt(2)) - sqrt(2) + sqrt(-2 + 2*sqrt(2)) + 3)/(2*(-sqrt(2) + sqrt(2)*sqrt(-1 + sqrt(2)) + 1)), a*sqrt(-1/2 + sqrt(2)/2))

res2[1](a => 3973).evalf() |> println
res2[2](a => 3973).evalf() |> println
res2[3](a => 3973).evalf() |> println

    708.000016603060
    2263.73675775652
    1808.07201601373

正方形の一辺の長さが 3973 のとき，等円の直径は 1416.0000 有奇である。

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

  a = 3973;  r = 708;  y = 2263.74;  b = 1808.07

function draw(a, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho") 
    (r, y, b) = (a*(1 - sqrt(-1 + sqrt(2)))/2, -sqrt(2)*a/2 + sqrt(2)*a*sqrt(-2 + 2*sqrt(2))/4 + a*sqrt(-2 + 2*sqrt(2))/2 + a/2, a*(-20 - 2*sqrt(-2 + 2*sqrt(2)) + 3*sqrt(-1 + sqrt(2)) + 14*sqrt(2))/(-12*sqrt(-1 + sqrt(2)) - 4 + 3*sqrt(2) + 8*sqrt(2)*sqrt(-1 + sqrt(2))))
    @printf("正方形の一辺の長さが %.9g のとき，等円の直径は %.9g である。\n", a, 2r)
    @printf("a = %g;  r = %g;  y = %g;  b = %g\n", a, r, y, b)
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    circle(r, r, r)
    circle(r, y, r)
    circle(a - y, a - r, r)
    circle(a - r, a - r, r)
    segment(0, a, a, 0, :magenta)
    segment(0, b, a, 0, :magenta)
    segment(a - b, a, a, 0, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        point(r, r, "等円:r,(r,r)", :red, :left, :center, delta=-delta)
        point(r, y, "等円:r,(r,y)", :red, :left, :center, delta=-delta)
        point(0, b, "b ", :magenta, :right, :vcenter)
        point(a, 0, "a", :green, :center, delta=-delta/2)
        plot!(xlims=(-4delta, a + 5delta), ylims=(-4delta, a + 5delta))
    end
end;

draw(3973, true)

算額（その1413）

算額（その1413）

八四 加須市不動岡 総願寺 明治12年(1879)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：円3個，正方形，斜線2本
#Julia, #SymPy, #算額, #和算

問題は判読不能ということであるが，連立方程式を立ててそれを解くというやり方であれば，何を既知数として何を未知数として解くかさえ決めれば良いので，問題が判読不能であっても問題を設定しその解を求めることは容易である。

算額の図形を文章で記述すると以下のようになる。

正方形の中に一つの頂点から対辺へ 2 本の斜線を引く。区画された領域に大円と小円を入れる。

この図形を描くためにはパラメータが 4 個必要である。正方形の右下隅を原点として描く。

• 正方形の右下の頂点座標 (a, 0)
• 斜線と左側の辺の交点座標 (0, b)
• 大円の直径 r1
• 小円の直径 r2

4 個のパラメータの相互関係は以下の 2 本の連立方程式で記述できる。

1.  大円の中心 (r1, a - r1) と 2 点 (a, 0), (0, b) を結ぶ直線の距離が r1 である。
2.  小円の中心 (r2, r2) と 2 点 (a, 0), (0, b) を結ぶ直線の距離が r2 である。

その他にも何通りかの方程式が考えられることもあるが，見かけは違っても上の式と同じ意味を持つものもある。独立な式は 2 本である。

上の連立方程式には 4 個のパラメータが含まれるが，そのうちの 2 個が既知であれば，残りの2 個のパラメータは連立方程式を解けば得られる。
既知のパラメータに直接数値を代入しても良いが，変数名そのままで方程式を解けば，解の中にも変数名が残るので，一般解が得られる。その時点で変数名にいろいろな値を与えれば，それに応じて数値解が得られる。

4 個のパラメータのうち，どの 2 個のパラメータを既知とするかは 4C2 = 6 通りある。その中には方程式を解くうえで面白いものや難しいものがある。当然「問」として適切なのはいくつかに限られる。

今回の算額図からは，以下のような「問」が考えられる。

1.  正方形の一辺と小円の直径が与えられたとき，大円の直径を求める術（すべ）を述べよ。
2.  正方形の一辺と大円の直径が与えられたとき，小円の直径を求める術（すべ）を述べよ。
3.  大円と小円の直径がが与えられたとき，正方形の一辺を求める術（すべ）を述べよ。

1.　大円の直径を求める

問1. 正方形の一辺と小円の直径が与えられたとき，大円の直径を求める術（すべ）を述べよ。

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

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1)
eq2 = dist2(a, 0, 0, b, r2, r2, r2);
# eq3 = a + b - sqrt(a^2 + b^2) - 2r2;

res1 = solve([eq1, eq2], (r1, b))[2]  # 1 of 2

    ((a^3 - 4*a^2*r2 + 2*a*r2^2)/(2*a^2 - 6*a*r2 + 4*r2^2), 2*r2*(-a + r2)/(-a + 2*r2))

大円の半径は a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2)) で計算できる。

res1[1] |> simplify |> println

    a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2))

正方形の一辺の長さ a = 120，小円の直径 r2 = 48 のとき，大円の直径は 70 である。
具体的な数値が挙げられている場合には，与えられる数，得られる数が整数である場合が多い。

a = 120
r2 = 48/2
a*(a^2 - 4*a*r2 + 2*r2^2)/(2*(a^2 - 3*a*r2 + 2*r2^2))

    35.0

res1[1](a => 120, r2 => 24).evalf() |> println
res1[2](a => 120, r2 => 24).evalf() |> println

    (a^3 - 4.0*a^2*r2 + 2.0*a*r2^2)/(2.0*a^2 - 6.0*a*r2 + 4.0*r2^2)
    2.0*r2*(-a + r2)/(-a + 2.0*r2)

2.　小円の直径を求める

問2. 正方形の一辺と大円の直径が与えられたとき，小円の直径を求める術（すべ）を述べよ。

using SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1)
eq2 = dist2(a, 0, 0, b, r2, r2, r2);
# eq3 = a + b - sqrt(a^2 + b^2) - 2r2;

res2 = solve([eq1, eq2], (r2, b))[1]  # 1 of 4

    ((-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1)), (-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - r1)^2)/(a - 2*r1))

小円の半径はかなり複雑な式になる。

a = 120
r1 = 70/2
(-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1))

    24.0

res2[1](a => 120, r1 => 70/2).evalf() |> println
res2[2](a => 120, r1 => 70/2).evalf() |> println

    0.5*(-2.0*r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 + (a - 2.0*r1)*(a - 0.5*((4.0*a^2*(0.5*a - r1)^2 + 4.0*(r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 - 0.5*(a - r1)^2)^2)/(0.5*a - r1)^2)^0.5) + (a - r1)^2)/(a - 2.0*r1)
    (-2.0*r1*(0.5*a^2 - a*r1 + 0.25*r1^2)^0.5 + (a - r1)^2)/(a - 2.0*r1)

3.　正方形の一辺の長さを求める

問3. 大円と小円の直径がが与えられたとき，正方形の一辺を求める術（すべ）を述べよ。

この場合は，SymPy の能力的に連立方程式を解くのは難しい。逐次的に解いていく。

iusing SymPy
@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = dist2(a, 0, 0, b, r1, a - r1, r1) |> x -> x/a
eq2 = dist2(a, 0, 0, b, r2, r2, r2) |> x -> x/(a*b);

まず，eq2 を解いて b を求める。

ans_b = solve(eq2, b)[1];
ans_b |> println

    2*r2*(a - r2)/(a - 2*r2)

eq1 の b に ans_b を代入し，eq11 とする。eq11 を解いて a を求める。

eq11 = eq1(b => ans_b) |> simplify |> numerator;

ans_a = solve(eq11, a);

5 組の解が得られるが 5 番目のものが形式上は虚数解であるが，虚数部は実質的に 0 なので適解である。a = 120

ans_a[5](r1 => 35, r2 => 24).evalf() |> println

    120.0 - 8.470329472543e-22*I

ans_b に a, r1, r2 を代入すれば b が得られる。

ans_b(r1 => 35, r2 => 24, a => 120) |> println

    64

function draw(type, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho") 
    if type == 1
        (a, r2) = (120, 24)
        (r1, b) = ((a^3 - 4*a^2*r2 + 2*a*r2^2)/(2*a^2 - 6*a*r2 + 4*r2^2), 2*r2*(-a + r2)/(-a + 2*r2))
    elseif type == 2
        (a, r1) = (120, 35)
        (r2, b) = ((-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - 2*r1)*(a - sqrt((a^2*(a - 2*r1)^2 + (r1*sqrt(2*a^2 - 4*a*r1 + r1^2) - (a - r1)^2)^2)/(a - 2*r1)^2)) + (a - r1)^2)/(2*(a - 2*r1)), (-r1*sqrt(2*a^2 - 4*a*r1 + r1^2) + (a - r1)^2)/(a - 2*r1))
    elseif type == 3
        (r1, r2) = (35, 24)
        a = 120
        b = 2*r2*(a - r2)/(a - 2*r2)
    end
    @printf("a = %g;  b = %g;  r1 = %g;  r2 = %g\n", a, b, r1, r2)
    plot([0, a, a, 0, 0], [0, 0, a, a, 0], color=:green, lw=0.5)
    circle(r1, a - r1, r1)
    circle(r2, r2, r2, :blue)
    circle(a - r2, a - r2, r2, :blue)
    segment(0, b, a, 0, :magenta)
    segment(a - b, a, a, 0, :magenta)
    if more
        delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
        hline!([0], color=:gray80, lw=0.5)
        point(r1, a - r1, "大円:(r1,a-r1)", :red, :center, delta=-delta)
        point(r2, r2, "小円:(r2,r2)", :blue, :center, delta=-delta)
        point(0, b, "b ", :magenta, :right, :vcenter)
        point(a, 0, "a", :green, :center, delta=-delta/2)
        plot!(xlims=(-4delta, a + 5delta), ylims=(-4delta, a + 5delta))
    end
end;

draw(3, true)

うどん茶屋 てんてこ舞 TEN TEKO MAI

琴平町 うどん茶屋 てんてこ舞
金毘羅さんの参道（石段の手前）にある

算額（その1412）

算額（その1412）

八二 熊谷市三ケ尻 竜泉寺 明治11年(1878)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：3次元，球3個，円錐
#Julia, #SymPy, #算額, #和算

楕円錐（底面が楕円）の中に上，中，下の 3 個の珠を容れる。上，中，下の珠の直径の和は 120 寸，各直径の差は 10 寸のとき，下珠の直径はいかほどか。

注：本問は 3 個の珠の直径の差が等差数列であることを前提にしている。算額（その1409）では四角錐に内接する球の直径は等比数列であるとした。球を内包する立体は，四角錘でなくて円錐でもよいわけである。もし円錐に等差数列に従う球を入れるとどうなるか。たとえば，図に示すような小球と中球が内接するような円錐では大球はグスグスである（赤い玉と円錐の側面に隙間がある）。しかし，円錐を垂直に立てれば，大球は上にある中球・小球の重さで横方向に押しのけられ，円錐の側面に一点で内接するのである。ということで，円錐と等差数列をなす 3 球ということに矛盾は生じないので，安心して問題を解く。なお，わざわざ楕円錐としたのは長径方向に投影したときに，投影された 3 珠の重心が一直線になるということくらいか？（円錐だとどちら方向に投影するか特定できないから）

蛇足：3つの珠の直径が 30, 40, 50 のとき，初項が 30，公比が 4/3 なので，3番めの珠の直径は 30*(4/3)^2 = 53.33333333333333 になるので，3 番目の珠の直径が 50 なら，グスグスなのである。

中学生レベル（？）の問題であるが，SymPy を使って解くと以下のようになる。
下珠，中珠，上珠の直径を 「上径」，「中径」，「下径」 とおく。直径の和を「径和」，直径の差を「各差」として以下の連立方程式を解く。

using SymPy
@syms 下径, 中径, 上径, 径和, 各差
eq1 = 下径 + 中径 + 上径 - 径和
eq2 = 下径 - 中径 - 各差
eq3 = 中径 - 上径 - 各差
res = solve([eq1, eq2, eq3], (下径, 中径, 上径))
(res[下径], res[中径], res[上径]) |> println

    (各差 + 径和/3, 径和/3, -各差 + 径和/3)

術は「上中下の 3 個の珠の直径の和に直径の差の 3 倍を加え 3 で割ることで下珠の直径を得る」ということで，小学生に説明するなら，このほうがよいかもしれない。でも，式の変形は小学校では扱わないか?

下径 + 中径 + 上径 = 径和
下径 + (下径 - 各差) + (下径 - 2×各差) = 径和
3×下径 - 3×各差 = 径和
3×下径 = 径和 + 3×各差
下径 = (径和 + 3×各差)÷3

算額（その1411）

算額（その1411）

九十二 群馬県富岡市一ノ宮 貫前神社 安政5年(1858)
九十七 群馬県高崎市石原町 清水寺 安政5年(1858)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円12個，外円，正三角形3個
#Julia, #SymPy, #算額, #和算

外円の中に甲円 1 個，乙円 2 個を容れる。甲円，乙円の中には正三角形を容れ，さらに甲円の中には丙円と等円，乙円の中には 3 個ずつの等円を容れる。外円の直径が 9 寸のとき，等円の直径はいかほどか。

注：本問と関連する算額を「算額（その1410）」に掲示した。違いは甲円の中にある丙円の位置と大きさである。

外円の中の甲円と乙円の位置と大きさについて
外円の半径と中心座標を R, (0, )
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (r2, y2)
等円の半径を r
とおき，以下の方程式を解く。

(1) 甲円の中の正三角形と丙円 r3，等円 r について以下の連立方程式を解く。

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

@syms r1, r3, r

eq1 = r/(2r1 - 2r3 - r) - 1//2
eq2 = r3/(2r1 - r3) - 1//2
res = solve([eq1, eq2], (r3, r))
res[r3] |> println
res[r] |> println

    2*r1/3
    2*r1/9

等円の半径は r = 2r1/9 である。

(2) 乙円の中の正三角形と等円の大きさについては，「算額（その918）」に示したように，乙円の中にある等円の半径は乙円の半径の 1/4 である。 r = r2/4

よって，2r1/9 = r2/4 である。外円の中に甲円と乙円が入っているので，以下の連立方程式を解いて甲円と乙円の半径と中心座標を求める。

@syms R, r1, r2, y2
eq1 = 2r1/9 - r2/4
eq2 = r2^2 + y2^2 - (R - r2)^2
eq3 = r2^2 + (R - r1 - y2)^2 - (r1 + r2)^2
res = solve([eq1, eq2, eq3], (r1, r2, y2))[3]  # 3 of 3

    (R/2, 4*R/9, -R/3)

甲円の半径 r1 は外円の半径 R の 1/2 倍，乙円の半径 r2 は外円の半径 R の 4/9 倍である。
甲円の中にある等円の半径は外円の半径 R の (1/2)*(2/9) = 1/9 倍，乙円の中にある等円の半径は外円の半径 R の (1/4)*(4/9) = 1/9 倍で等しくなる。

「問」のように外円の直径が 9 寸のとき，等円の直径は 9*(1/9) = 1 寸となり，「答」，「術」と一致する。

R = 4.5;  r1 = 2.25;  r2 = 2;  y2 = -1.5
甲円の直径 = 4.5，丙円の直径 = 3.0，甲円の中にある等円の直径 = 1.0
乙円の直径 = 4.0，乙円の中にある等円の直径 = 1.0

考察

「九十二 貫前神社」の算額は安政5年戊午3月，「九十七 清水寺」の算額は安政戊午11月に奉納されている。
「百五十一 産泰神社」の算額は年代不明とされている。
多数決ではないが，貫前神社，清水寺のほうが正しいのであろう。
産泰神社のものは，うっかりと図を写し間違えたと思われる。
産泰神社の方が正しいとしても，甲円の中の等円は正三角形の 3 つの隅に 3 個ある方が自然なので，一隅にしかないのはやはり不自然である。貫前神社・清水寺のような場合は甲円の中に等円は 1 個しかないのはもっともである。

function draw_sub(sw, x, y, r, color, color2, color3, color4)
    circle(x, y, r, color)
    (x1, y1) = r .* (cosd(30), -sind(30))
    plot!(x .+ [x1, 0, -x1, x1], y .+ [y1, r, y1, y1], color=color2, lw=0.5)
    if sw == 1
        r1 = 2r/3
        r2 = 2r/9
        circle(x, y - r + r1, r1, color4)
        circle(x, y - r + 2r1 + r2, r2, color3)
        println("甲円の直径 = $(2r)，丙円の直径 = $(2r1)，甲円の中にある等円の直径 = $(2r2)")
    else
        for deg in (30, 150, 270)
            circle(x + 3r/4*cosd(deg), y + 3r/4*sind(deg), r/4, color3)
        end
        sw == 2 && println("乙円の直径 = $(2r)，乙円の中にある等円の直径 = $(2r/4)")
    end
end

function draw(more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 9/2  # 1/2
    (r1, r2, y2) = (R/2, 4*R/9, -R/3)
    @printf("R = %g;  r1 = %g;  r2 = %g;  y2 = %g\n", R, r1, r2, y2)
    plot()
    circle(0, 0, R)
    # circle(0, R - r1, r1, :orange)
    draw_sub(1, 0, R - r1, r1, :green, :chartreuse3, :darkorange, :tomato)
    # circle2(r2, y2, r2, :blue)
    draw_sub(2, r2, y2, r2, :blue, :darkblue, :darkorange, :none)
    draw_sub(3, -r2, y2, r2, :blue, :darkblue, :darkorange, :none)
    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", :red, :center, :bottom, delta=delta)
        point(0, R - r1, "甲円:R/2,(0,R-r1)", :green, :center, delta=-delta/2)
        point(0, R/3, "丙円:R/3,(0,R/3)", :tomato, :center, delta=-delta/2)
        point(0, 7R/9, "等円:R/9,(0,7R/9)", :tomato, :center, delta=-delta/2)
        point(r2, y2, "乙円:4R/9,(r2,y2)", :blue, :center, delta=-delta/2)
        point(r2, -2R/3, "等円:R/9,(r2,y2)", :darkorange, :center, delta=-delta/2)
    end
end;

draw(true)

算額（その1410）

算額（その1410）

百五十一 群馬県前橋市下大屋町 産泰神社 年代不明
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

桜沢門人 武州幡羅郡三箇尾村 権田源之助正賢
山口正義：やまぶき2 第37号,2016年6月20日.
https://yamabukiwasan.sakura.ne.jp/ymbk37.pdf

キーワード：円12個，外円，正三角形3個
#Julia, #SymPy, #算額, #和算

外円の中に甲円 1 個，乙円 2 個を容れる。甲円，乙円の中には正三角形を容れ，さらに甲円の中には丙円と等円，乙円の中には 3 個ずつの等円を容れる。外円の直径が 9 寸のとき，等円の直径はいかほどか。

注：本問と関連する算額を「算額（その1411）」に掲示する。違いは甲円の中にある丙円の位置と大きさである。

外円の中の甲円と乙円の位置と大きさについて
外円の半径と中心座標を R, (0, )
甲円の半径と中心座標を r1, (0, R - r1)
乙円の半径と中心座標を r2, (r2, y2)
等円の半径を r
とおき，以下の方程式を解く。

(1) 甲円の中の正三角形と丙円，等円の大きさについては，「算額（その593）」に示したように，甲円の中にある等円の直径は甲円の直径の 1/6 である。r = r1/6

(2) 乙円の中の正三角形と等円の大きさについては，「算額（その918）」に示したように，乙円の中にある等円の直径は乙円の直径の 1/4である。 r = r2/4

よって，r2 = 2r1/3 である。外円の中に甲円と乙円が入っているので，以下の連立方程式を解いて甲円と乙円の半径と中心座標を求める。

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

    (R*(-4 + sqrt(21)), 2*R*(-4 + sqrt(21))/3, R*(-2 + sqrt(21)/3))

甲円の半径 r1 は外円の半径 R の (√21 - 4) 倍，乙円の半径 r2 は外円の半径 R の (√21 - 4)*2/3 倍である。
甲円の中にある等円の半径は外円の半径 R の (√21 - 4)/6 倍，乙円の中にある等円の半径は外円の半径 R の (√21 - 4)*2/3/4 = (√21 - 4)/6 倍でで等しくなる。

問のように外円の直径が 9 寸のとき，等円の直径は 9*(√21 - 4)/6 = 0.8738635424337597 寸となる。

R = 4.5;  r1 = 2.62159;  r2 = 1.74773;  y2 = -2.12614
甲円の直径 = 5.2431812546025585，丙円の直径 = 1.3107953136506396，甲円の中にある等円の直径 = 0.8738635424337597
乙円の直径 = 3.495454169735039，乙円の中にある等円の直径 = 0.8738635424337597

おかしいですね。
答，術は「外円径を 9 で割れば，等円の直径が得られる。」すなわち，「外円径が 9寸ならば，等円径は 1 寸」のはず。

実は同じような算額が他に 2 件あるのは承知していた（その 2 件の算額の内容は同じ）。（時系列とは逆に）たまたまこちらを先に解いた。
似ている算額について，次の記事「算額（その1411）」を書く。

ここまで書いてきて，だいぶ前に読んだ記事にあった図（山口）を思い出した。再確認した所，「図の写し間違いだろう」という指摘も同じであった。

function draw_sub(sw, x, y, r, color, color2, color3)
    circle(x, y, r, color)
    (x1, y1) = r .* (cosd(30), -sind(30))
    plot!(x .+ [x1, 0, -x1, x1], y .+ [y1, r, y1, y1], color=color2, lw=0.5)
    if sw == 1
        circle(x, y, r/2, color2)
        circle(x, y + (r/2 + r/6), r/6, color3)
        println("甲円の直径 = $(2r)，丙円の直径 = $(r/2)，甲円の中にある等円の直径 = $(2r/6)")
    else
        for deg in (30, 150, 270)
            circle(x + 3r/4*cosd(deg), y + 3r/4*sind(deg), r/4, color3)
        end
        sw == 2 && println("乙円の直径 = $(2r)，乙円の中にある等円の直径 = $(2r/4)")
    end
end

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, r2, y2) = (R*(-4 + sqrt(21)), 2*R*(-4 + sqrt(21))/3, R*(-2 + sqrt(21)/3))
    @printf("R = %g;  r1 = %g;  r2 = %g;  y2 = %g\n", R, r1, r2, y2)
    plot()
    circle(0, 0, R)
    # circle(0, R - r1, r1, :orange)
    draw_sub(1, 0, R - r1, r1, :green, :chartreuse3, :darkorange)
    # circle2(r2, y2, r2, :blue)
    draw_sub(2, r2, y2, r2, :blue, :darkblue, :darkorange)
    draw_sub(3, -r2, y2, r2, :blue, :darkblue, :darkorange)
    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", :red, :center, :bottom, delta=delta)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :green, :center,:bottom, delta=delta/2)
        point(0, R - r1, "丙円:r1/2,(0,R-r1)", :tomato, :center, delta=-delta/2, mark=false)
        point(0, R - r1/3, "等円:r1/6,(0,R-r1/3)", :tomato, :center, delta=-delta/2)
        point(r2, y2, "乙円:r2,(r2,y2)", :blue, :center, delta=-delta/2)
        point(r2, y2 - 3r2/4, "等円:r2/4\n(r2,y2-3r2/4)", :darkorange, :left, delta=-8delta/2, deltax=4delta)
        point(2.1, 0, "r1 = R*(√21 - 4)\nr2 = 2r1/3", :black, :left, :vcenter, deltax=-0.5delta, mark=false)

    end
end;

draw(9/2, true)

 

算額（その1409）

算額（その1409）

百五十 群馬県多野郡新町 稲荷神社 文政3年(1820)
百五十三 多野郡新町 諏訪神社 年代不明（上と同一人物）
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：球3個，四角台
#Julia, #SymPy, #算額, #和算

四角台の中に，上球，中球，下球が互いに接し，側面，上面，下面にも接している。
上球，下球の直径が与えられたとき，中球の直径はいかほどか。

以下の図は四角台の側面を除いて，横倒しにした状態である。

算額（その1405）の拡張なので，おなじように解いてもよい。

しかし，単に球の直径を求めるだけならば，算額（その1393）などに述べた問題の3次元版なので，もっと簡単になる。

2直線に挟まれた互いに外接する円の直径は等比数列をなす。
同じく，底辺が正方形の四角錐の 4 平面に挟まれた互いに外接する球の直径は等比数列をなす。
後者は，1 つの側面と平行な丙面に射影すれば，2直線に挟まれた円が見えるであろうことから同じ問題であることがわかる。

よって，上球，中球，下球 の直径を d1, d2, d3，公比を p とすれば，d2 = d1*p，d3 = d2*p = d1*p^2 である。
d1 * d3 = d1^2 * p^2 = d2^2 なので，d2 = sqrt(d1 * d3) である。

中球の直径は，上球の直径と下球の直径を掛け合わせ，平方根をとれば求まる。

 

 

算額（その1408）

算額（その1408）

百四十九 群馬県多野郡新町 稲荷神社 文政3年(1820)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円4個，半円，直角三角形
#Julia, #SymPy, #算額, #和算

半円の中に三角形（直角三角形）を作り，甲円，乙円，丙円，丁円を容れる。大斜（半円の直径），小斜（三角形の一番短い辺）が与えられたとき，丁円の直径はいかほどか。

半円の半径と中心座標を R, (0, 0)
甲円の半径と中心座標を r1, (x1, r1)
乙円の半径と中心座標を r2, (x2, r2)
丙円の半径と中心座標を r3, (x3, r3)
丁円の半径と中心座標を r4, (x4, r4)
円周上にある直角三角形の頂点の座標を (x0, sqrt(R^2 - x0))
とおき，以下の連立方程式を解く。

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

@syms R, x0, r1, x1, r2, x2, r3, x3, r4, x4,
      大斜, 小斜, 中斜
R = 大斜/2
eq1 = sqrt((R - x0)^2 + (R^2 - x0^2)) - 小斜
eq2 = r2/(x2 + R) - r1/(x1 + R)
eq3 = r3/(R - x3) - r1/(R - x1)
eq4 = x3 - x1 - 2sqrt(r1*r3)
eq5 = x1 - x4 - 2sqrt(r1*r4)
eq6 = x4 - x2 - 2sqrt(r2*r4)
eq7 = x1 - x2 - 2sqrt(r1*r2)
eq8 = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1)
eq9 = 小斜 + sqrt(大斜^2 - 小斜^2) - 大斜 - 2r1;
# eqxx = (大斜 + sqrt(大斜^2 - 小斜^2) + 小斜)*r1 - 小斜*中斜
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x1, r1, r1);
# eqxx = dist2(-R, 0, x0, sqrt(R^2 - x0^2), x2, r2, r2)
# eqxx = dist2(R, 0, x0, sqrt(R^2 - x0^2), x3, r3, r3);
# solve([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq10], (x0, r1, x1, r2, x2, r3, x3, r4, x4))

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)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = u
    return [
        -小斜 + sqrt(-x0^2 + 大斜^2/4 + (-x0 + 大斜/2)^2),  # eq1
        -r1/(x1 + 大斜/2) + r2/(x2 + 大斜/2),  # eq2
        -r1/(-x1 + 大斜/2) + r3/(-x3 + 大斜/2),  # eq3
        -x1 + x3 - 2*sqrt(r1*r3),  # eq4
        x1 - x4 - 2*sqrt(r1*r4),  # eq5
        -x2 + x4 - 2*sqrt(r2*r4),  # eq6
        x1 - x2 - 2*sqrt(r1*r2),  # eq7
        (r1^2*x0 - r1*x1*sqrt(-4*x0^2 + 大斜^2) - x0*x1^2 + 大斜*(-4*r1^2 - 4*r1*sqrt(-4*x0^2 + 大斜^2) - 8*x0*x1 - 2*x0*大斜 + 4*x1^2 + 4*x1*大斜 + 大斜^2)/8)/大斜,  # eq8
        -2*r1 - 大斜 + 小斜 + sqrt(大斜^2 - 小斜^2),  # eq9
    ]
end;

大斜 = 10
R = 大斜/2
小斜 = 4
iniv = BigFloat[3.45, 1.58, 2.58, 1.05, 0.01, 0.46, 4.29, 0.32, 1.16]
res = nls(H, ini=iniv)

    ([3.4, 1.5825756949558398, 2.582575694955839, 1.0456116707250647, 0.009826491126086901, 0.46246227038447935, 4.293577213300662, 0.31816569849213683, 1.163390997458993], true)

たとえば，大斜，小斜，中斜が 10, 4, 9.16515 のとき，丁円の直径は 0.6363313969842737 である。

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

   x0 = 3.4;  r1 = 1.58258;  x1 = 2.58258;  r2 = 1.04561;  x2 = 0.00982649;  r3 = 0.462462;  x3 = 4.29358;  r4 = 0.318166;  x4 =1.16339

function draw(大斜, 小斜, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 大斜/2
    中斜 = sqrt(大斜^2 - 小斜^2)
    (x0, r1, x1, r2, x2, r3, x3, r4, x4) = res[1]
    @printf("大斜，小斜，中斜が %g, %g, %g のとき，丁円の直径は %g である。\n", 大斜, 小斜, 中斜, 2r4)
    @printf("x0 = %g;  r1 = %g;  x1 = %g;  r2 = %g;  x2 = %g;  r3 = %g;  x3 = %g;  r4 = %g;  x4 =%g\n", x0, r1, x1, r2, x2, r3, x3, r4, x4)
    y0 = sqrt(R^2 - x0^2)
    plot([-R, x0, R], [0, y0, 0], color=:blue, lw=0.5)
    circle(0, 0, R, beginangle=0, endangle=180)
    circle(x1, r1, r1)
    circle(x2, r2, r2, :green)
    circle(x3, r3, r3, :magenta)
    circle(x4, r4, r4, :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(x0, sqrt(R^2 - x0^2), "(x0,sqrt(R^2-x0^2))", :blue, :center, :bottom, delta=delta)
        point(x1, r1, "甲円:r1,(x1,r1)", :red, :center, delta=-delta)
        point(x2, r2, "乙円:r2\n(x2,r2)", :green, :center, delta=-delta)
        point(x3, r3, "丙円:r3,(x3,r3)", :magenta, :right, delta=-7delta, deltax=13delta)
        point(x4, r4, "丁円:r4,(x4,r4)", :orange, :right, delta=-5delta, deltax=13delta)
        point(R, 0, "大斜/2", :red, :left, :bottom, delta=delta)
        point(0, 2.6r2, "中斜", :blue, :center, mark=false)
        point((R + x0)/2, r1, "小斜", :blue, :center, mark=false, deltax=7delta)
        plot!(xlims=(-R - 3delta, R + 15delta), ylims!(-10delta, R + 3delta))
    end
end;

draw(10, 4, true)

算額（その1407）

算額（その1407）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：楕円，直角三角形
#Julia, #SymPy, #算額, #和算

直角三角形の中に斜線を 2 本描き，区画された領域に楕円を容れる。鈎が 3 寸，股が 4 寸のとき，楕円の長径，短径はいかほどか。

注：問には明記されておらず，図も正確ではないが，斜線と鈎，股の交点は鈎，股の中点である。

鈎，股，楕円の長径，短径をそのまま変数「鈎」，「股」，「長径」，「短径」とする。

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

三角形に内接する楕円の長径・短径に関する「算法助術」の公式97を援用する。
三角形の底辺, 斜辺1, 斜辺2, 高さを与えて，長径，短径の関係式を返す関数 equation を定義する。

@syms 底辺::positive, 斜辺1::positive, 斜辺2::positive, 高さ::positive, 長径::positive, 短径::positive
equation(底辺, 斜辺1, 斜辺2, 高さ) = -(斜辺1^2 - 斜辺2^2)^2*高さ*短径^2 + (斜辺1^2 - 斜辺2^2)^2*短径^3 + 底辺^4*高さ*(2§高さ - 短径)^2 - 底辺^4*短径*(2高さ - 短径)^2 - 底辺^2*高さ*長径^2*(2高さ - 短径)^2;

同じ楕円が 2 つの直角三角形に内接しているので，以下の連立方程式を解く。

@syms 鈎, 股
eq1 = equation(股, sqrt(股^2 + (鈎/2)^2), 鈎/2, 鈎/2)
eq2 = equation(股/2, sqrt((股/2)^2 + 鈎^2), 鈎, 鈎)
res = solve([eq1, eq2], (長径, 短径))[4];

# 長径
p = res[1] |> simplify |> sympy.sqrtdenest |> simplify
p |> println
p(鈎 => 3, 股 => 4).evalf() |> println

    股*(3 - sqrt(5))/2
    1.52786404500042

# 短径
q = res[2] |> simplify
q |> println
q(鈎 => 3, 股 => 4).evalf() |> println

    鈎*(3 - sqrt(5))/2
    1.14589803375032

長径は股の (3 - √5)/2 倍，短径は鈎の (3 - √5)/2 倍である。
鈎，股が 3 寸，4 寸のとき，長径は 4(3 - √5)/2 = 1.5278640450004204，短径は 3(3 - √5)/2 = 1.1458980337503153 である。

function draw(鈎, 股, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (a, b) = (股, 鈎) .* ((3 - √5)/ 2)  ./ 2
    @printf("鈎，股が %g，%g のとき，楕円の長径，短径は %g，%g である。\n", 鈎, 股, 2a, 2b)
    plot([0, 股, 0, 0], [0, 0, 鈎, 0], color=:blue, lw=0.5)
    ellipse(a, b, a, b, color=:red)
    segment(0, 鈎/2, 股, 0, :green)
    segment(0, 鈎, 股/2, 0, :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, 鈎, "鈎 ", :blue, :right, :vcenter)
        point(0, 鈎/2, "鈎/2 ", :blue, :right, :vcenter)
        point(股, 0, "股", :blue, :center, delta=-delta)
        point(股/2, 0, "股/2", :blue, :center, delta=-delta)
        point(a, b, "楕円:a,b,(a,b)", :red, :center, delta=-delta)
        plot!(xlims=(-10delta, 股 + 5delta), ylims=(-7delta, 鈎 + 5delta))
    end
end;

draw(3, 4, true)

算額（その1406）

算額（その1406）

百四十三 群馬県榛名町榛名山 榛名神社 明治33年(1900)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：円7個，外円，楕円，曲率円
#Julia, #SymPy, #算額, #和算

外円の中に，楕円と 2 個の交わる甲円を容れ，それらによって区切られた領域に乙円，丙円を 2 個ずつ容れる。乙円の直径が 1 寸で，楕円の短径がもっとも大きくなる（乙円が曲率円になる）ときの丙円の直径はいかほどか。

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

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

@syms R, r1, r2, r3
r1 = R - r3 
eq1 = (R - r2)^2 + r3^2 - (r1 + r2)^2  # 乙円は甲円と外接
eq2 = (R - 2r3)^2/R - r2  # 乙円は楕円の曲率円
res = solve([eq1, eq2], (r3, R))[3]  # 3 of 3

    (r2*(sqrt(2) + 2)/2, r2*(2*sqrt(2) + 3))

丙円の半径 r3 は，乙円の半径 r2 の (√2 + 2)/2 倍である。
乙円の直径が 1 寸で，楕円の短径がもっとも大きくなる（乙円が曲率円になる）とき，丙円の直径は (√2 + 2)/2 = 1.7071067811865475 寸である。

術では「天 = √8 + 3; (天 - √天)*乙円径/2」と複雑な計算式を提示しているが，「(√2/2 + 1)*乙円径」と簡約化できる。

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

r2 = 0.5;  r3 = 0.853553;  R = 2.91421;  r1 = 2.06066

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = r2*(2√2 + 3)
    r3 = r2*(√2 + 2)/2
    r1 = R - r3
    @printf("乙円の直径が %g のとき，丙円の直径は %g である。\n", 2r2, 2r3)
    @printf("r2 = %g;  r3 = %g;  R = %g;  r1 = %g\n", r2, r3, R, r1)
    plot()
    circle(0, 0, R)
    circle22(0, R - r3, r3, :blue)
    circle22(0, r3, r1, :green)
    circle(R - r2, 0, r2, :magenta)
    ellipse(0, 0, R, R - 2r3, color=: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, r3, "甲円:r1,(0,r3)", :green, :center ,delta=-delta/2)
        point(R - r2, 0, "乙円:r2\n(R-r2,0)", :magenta, :center, :bottom, delta=delta/2)
        point(0, R - r3, "丙円:r3,(0,R-r3)", :blue, :center, delta=-delta/2)
        point(0, R, "R", :red, :center, :bottom, delta=delta/2)
        point(R, 0, " R", :red, :left, :vcenter)
        point(0, R - 2r3, "R-2r3", :orange, :center, :bottom, delta=delta/2)
    end
end;

draw(1/2, true)

算額（その1405）

算額（その1405）

八七 群馬県碓氷郡松井田町峠 熊野神社 安政4年(1857)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：3次元，球2個，四角台
#Julia, #SymPy, #算額, #和算

上面，下面の一辺が 125 寸，2000 寸の正方形の四角台の中に小球と大球を容れる。小球と大球は互いに接するほかに，小球は側面と上面に 5 点で，大球は側面と下面に 5 点で接している。小球の直径はいかほどか。

上面の正方形，下面の正方形の一辺の長さをそれぞれ 2j, 2k
四角台の高さを i
四角台を含む四角錐の高さを h
大球の半径と中心座標を r1, (0, 0, r1)
小球の半径と中心座標を r2, (0, 0, 2r1 + r2)
とおく。
3次元の図形をx軸の正の方向から原点方向を見た y-z 平面へ投影された図を考える。
上面と下面の正方形の一辺の中点 i, k を結ぶ直線は小球と大球の共通接線になる。o, i, j, k は y-z 平面に平行な平面上にある。これを y-z 平面に投影した図は以下のようになる。

点と線の位置関係は以下の 3 本の方程式で表される。図形を描くために必要な変数は r1, r2, h の 3 変数なので，連立方程式を解けばよい。

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

@syms r1::positive, r2::positive, h::positive, j::positive, k::positive
eq1 = r2/(h - 2r1 - r2) - r1/(h - r1)
eq2 = r1/(h - r1) - k/sqrt(h^2 + k^2)
eq3 = j/(h - 2r1 - 2r2) - k/h
res = solve([eq1, eq2, eq3], (r1, r2, h))[2];  # 2 of 4

r1, r2, h ともに，ルートの中は j と k の高次式になっている。
小球の半径は，以下のようである。

res[2] |> display

res[2] |> println

    sqrt((2*j^(21/2)*k^4 - 12*j^(19/2)*k^5 + 30*j^(17/2)*k^6 - 40*j^(15/2)*k^7 + 30*j^(13/2)*k^8 - 12*j^(11/2)*k^9 + 2*j^(9/2)*k^10 + j^11*k^(7/2) - 5*j^10*k^(9/2) + 9*j^9*k^(11/2) - 5*j^8*k^(13/2) - 5*j^7*k^(15/2) + 9*j^6*k^(17/2) - 5*j^5*k^(19/2) + j^4*k^(21/2))/(j^(19/2)*k^3 - 5*j^(17/2)*k^4 + 9*j^(15/2)*k^5 - 5*j^(13/2)*k^6 - 5*j^(11/2)*k^7 + 9*j^(9/2)*k^8 - 5*j^(7/2)*k^9 + j^(5/2)*k^10 + 2*j^9*k^(7/2) - 12*j^8*k^(9/2) + 30*j^7*k^(11/2) - 40*j^6*k^(13/2) + 30*j^5*k^(15/2) - 12*j^4*k^(17/2) + 2*j^3*k^(19/2)))

SymPy では簡約化できないので，手動で簡約化する。まず根号の中身の分子と分母をそれぞれ変数変換して因子分解し，再度分数式に戻して平方根を取る。

# 分子を変数変換
@syms t, u
num = res[2].args[1](j^(1//2) => t, k^(1//2) => u) |> numerator |> factor
num |>  println

    t^8*u^7*(t - u)^6*(t + u)^8

# 分母を変数変換
den = res[2].args[1](j^(1//2) => t, k^(1//2) => u) |> denominator |> factor
den |> println

    t^5*u^6*(t - u)^6*(t + u)^8

# 分数式に戻す
r2 = num/den
r2 |> println

    t^3*u

# 変数の逆変換
r2 = sqrt(r2(t => j^(1//2), u => k^(1//2))) |> factor
r2 |> println

    j^(3/4)*k^(1/4)

小球の直径は上面の正方形の一辺の長さの 3 乗と下面の正方形の一辺の長さを掛けて，4 乗根を求めることで得られる。

# 実際の数値を代入して答えを求める
r2 = fourthroot(125^3 * 2000)

    250.0

大球の直径は小球の直径の 4 倍，四角台の高さは大球と小球の直径の和（小球の直径の 5 倍）である。

# 大球の直径は小球の直径の何倍か？
res[1]/res[2] |> println
(res[1]/res[2])(j => 125//2, k => 2000//2) |> println

    sqrt(k)/sqrt(j)
    4

# 四角錐の高さは小球の直径の何倍か？
res[3]/res[2] |> println
(res[3]/res[2])(j => 125//2, k => 2000//2) |> println

    2*(sqrt(j)*k + k^(3/2))/(sqrt(j)*(-j + k))
    32/3

1.　連立方程式の解き方（別法）

割り算を避けるなど，なるべく式を簡潔にして，逐次的に解を求めていく。

@syms r1::positive, r2::positive, h::positive, j::positive, k::positive
eq1 = r2*(h - r1) - r1*(h - 2r1 - r2) |> expand
eq2 = r1*sqrt(h^2 + k^2) - k*(h - r1) |> expand
eq3 = j*h - k*(h - 2r1 - 2r2) |> expand;

eq2 は r2 を含まないので，そのぶん簡潔なので，まず eq1, eq2 から r1, h を求める（式に r2 が含まれる）。

res = solve([eq1, eq2], (r1, h))[1]

    (k^(2/3)*r2^(1/3), 2*k^(4/3)*r2^(2/3)/(k^(2/3)*r2^(1/3) - r2))

eq3 の r1, h に代入して r2 を求める。

eq13 = eq3(r1 => res[1], h => res[2]) |> simplify |> numerator;

ans_r2 = solve(eq13, r2)[1]
ans_r2 |> println

    j^(3/4)*k^(1/4)

求められた r2 を，先に求めた r1, h の r2 に代入する。

@syms j, k
j = 125/2
k = 2000/2
r2 = j^(3/4)*k^(1/4)
r1 = k^(2/3)*r2^(1/3)
h  = 2*k^(4/3)*r2^(2/3)/(k^(2/3)*r2^(1/3) - r2)
println("r2 = $r2\nr1 = $r1\nh = $h")

    r2 = 125.0
    r1 = 499.9999999999999
    h = 1333.3333333333328

function draw(j, k, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r2 = fourthroot(j^3 * k)
    r1 = 4r2
    i = 2r1 + 2r2
    h = r2*32/3
    @printf("上面と下面の正方形の一辺の長さが %g, %g のとき，小球，大球の直径は %g, %g，四角台の高さは %g である。\n",2j, 2k, 2r2, 2r1, 2r1 + 2r2)
    plot([k, 0, -k, k], [0, h, 0, 0], color=:blue, lw=0.5)
    circle(0, r1, r1)
    circle(0, 2r1 + r2, r2, :green)
    segment(j, i, -j, i, :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(0, r1, "大球:r1,(0,0,r1)", :red, :center, delta=-delta)
        point(0, 2r1 + r2, "小球:r2,(0,0,2r1+r2)", :green, :left, delta=-delta, deltax=-5delta)
        point(j, i, "(j,0,i)", :blue, :left, :bottom, delta=delta)
        point(k, 0, "(k,0,k)", :blue, :left, :bottom, delta=delta)
        point(0, h, "(0,0,h)", :blue, :left, :bottom, delta=delta)
        xlims!(-k - 3delta, k + 15delta)
    end
end;

draw(125/2, 2000/2, true)

さぬきのおもてなし 麺匠 くすがみ

高松市楠上町 さぬきのおもてなし 麺匠 くすがみ
やや細麺，ぶっかけうどん（冷）にしたせいもあろうが，しっかりした腰があった
あまりに美味しそうで，一口二口食べた後で写真を取り忘れたことに気づいた😋

算額（その1404）

算額（その1404）

依拠した図がでたらめなものであったので，改訂版を書いた。

「算額（その1404）改訂版」を参照のこと。

三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市.
http://www.wasan.jp/yamamura/yamamura.html
キーワード：円3個，円弧，斜線2本
#Julia, #SymPy, #算額, #和算

問の原文は以下のとおりである。
「今有如図半円内設等円二個従其親所洩二斜容等円一個其等円径若干問得半円径術如何」

最初にいっておく。本問は条件不足であり，解けない。

原文には曖昧な点もあるので，補足的説明も加えると以下のようになろう。
「円弧（注1）」の中に等円 2 個を容れ，等円と円弧の接点と円弧の中心を結ぶ半径を 2 本引く（注2）。2 本の半径と円弧の弦が作る三角形の中に，円弧内の等円と同じ大きさの等円が内接する。等円の直径が与えられたとき，円弧の半径を求める方法を述べよ。
注1：原文は「半円」であるが，それでは図が成り立たない。
注2：この半径は，（当然ではあるが）円弧内の等円の中心を通る。

「条件不足」という意味は，以下に示すように，「等円径若干」で与えられた等円の直径が同じでも，円弧の直径が一意に定まらないということである。

円弧の半径と中心座標を R, (0, 0)
弦と y 軸の交点座標を (0, y)
等円の半径と中心座標を r, (x, y + r), (0, y - r)
とおくと以下の 2 方程式が立つ。
特定の r が与えられたとき，図を描くために必要なのは R, x, y の 3 変数であるが条件式が 2 本なので，「条件不足」ということになる。

術は「置等円径三之得半円径」ということで，「等円の直径を 3 倍すれば半円（円弧）の径が得られる」ということであるが，算額の図を見てもそんな比率でないことは明らかである。

r = 1/2 の場合について，連立方程式を解くが，「等円の直径が同じでも，円弧の直径が一意に定まらない」ことを示すために，R を色々変えて，x, y を求める。

R = 3 のときの x, y は以下のようになる。

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

@syms x, y, r, R
r = 1/2
R = 3
eq1 = x^2 + (y + r)^2 - (R - r)^2
eq2 = dist2(0, 0, x, y + r, 0, y - r, r);
res = solve([eq1, eq2], (x, y));

8 組の虚数解が得られるが，虚数部は実質 0 である。また，x < y という条件をつけて，以下の手順で実質的に適切な解を選択・表示する。

for (x, y) in res
   if abs(imag(x)) < 1e-15 && abs(imag(y)) < 1e-15 &&
      real(x) > 0 && real(y) > 0 && real(x) < real(y)
       println("(R, x, y) = ", (R, real(x), real(y)))
   end
end

   (R, x, y) = (3, 0.953431800017409, 1.81105339676857)

等円の直径が 1 のとき，円弧の直径は 6 である。ただし，x = 0.953432，y = 1.81105

しかし，円弧の直径は一意に決まらない。たとえば，他の例として，
等円の直径が 1 のとき，円弧の直径は 5.7 である。ただし，x = 1.08171，y = 1.58624
等円の直径が 1 のとき，円弧の直径は 8 である。ただし，x = 0.721715，y = 2.92478
　　：
当然ではあるが，円弧の直径はどのような値でもよいわけではなく，一定の範囲内に収まっていなければならない。

function draw(r, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   (R, x, y) = (3, 0.953431800017409, 1.81105339676857)
   @printf("等円の直径が %g のとき，円弧の直径は %g である。ただし，x = %g，y = %g である。\n", 2r, 2R, x, y)
   @printf("r = %g;  R = %g;  x = %g;  y = %g\n", r, R, x, y)
   plot()
   x0 = sqrt(R^2 - y^2)
   segment(-x0, y, x0, y)
   θ = atand(y, x0)
   circle(0, 0, R, beginangle=θ, endangle=180 - θ)
   circle2(x, y + r, r, :blue)
   circle(0, y - r, r, :blue)
   segment(0, 0, x*R/(R - r), (y + r)*R/(R - r))
   segment(0, 0, -x*R/(R - r), (y + r)*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(x, y + r, "等円:r,(x,y+r) ", :blue, :right, :vcenter)
       point(0, y - r, "等円:r\n(0,y-r)", :blue, :center, :bottom, delta=2delta)
       point(0, y, "y", :black, :center, :bottom, delta=2delta)
       point(x0, y, "(sqrt(R^2-y^2),y)", :red, :right, delta=-delta)
       point(x - r, y - 1.5r, @sprintf("r = %g\nR = %.15g\nx = %.15g\ny = %.15g", r, R, x, y),
           :black, :left, deltax=5delta, mark=false)
   end
end;

draw(1/2, true)

算額（その1403）

算額（その1403）

依拠した図がでたらめなものであったので，改訂版を書いた。

「算額（その1403）改訂版」を参照のこと。

三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：円2個，直角三角形4個，等脚台形
#Julia, #SymPy, #算額, #和算

合同な直角三角形 4 個を組み合わせて等脚台形を作り，大円と小円を容れる。小円の直径が与えられたときに大円の直径を求める術を述べよ。

最初にいっておく。本問は条件不足であり，解けない。

できる台形の上底，下底，高さを a, 3a, b
大円の半径と中心座標を r1, (a, r1)
小円の半径と中心座標を r2, (2a - r2, b - r2)
とおく。

条件が不足しているので，小円の直径 r2 の他に台形の高さ b も与えられるものとして方程式を解く。

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

@syms a::positive, b::positive, r1::positive, r2::positive
eq1 = a + b - sqrt(a^2 + b^2) - 2r2
eq2 = r1/(b - r1) - a/sqrt(a^2 + b^2)
res = solve([eq1, eq2], (r1, a))[1];

ans_r1 = res[1]
ans_r1 |> println
ans_r1(r2 => 1/2, b => 2) |> println

   2*b*r2*(b - r2)*Abs(b - 2*r2)/(b*sqrt(b^2*(b - 2*r2)^2 + 4*r2^2*(b - r2)^2) + 2*r2*(b - r2)*Abs(b - 2*r2) - 2*r2*sqrt(b^2*(b - 2*r2)^2 + 4*r2^2*(b - r2)^2))
   0.750000000000000

ans_r1 は b > 2r2 を仮定して，以下のように簡約化できる。

ans_r1_ = 2r2*(b - r2)/b
ans_r1_ |> println

   2*r2*(b - r2)/b

ans_a も以下のようになる。

ans_a = res[2]
ans_a |> println

   2*r2*(b - r2)/(b - 2*r2)

術では「41 の平方根から 5 を引いて，小円の直径を掛ければ，大円の直径が得られる」としいているが，小円の直径が同じでも台形の高さが変われば大円の直径も変わる。一律に (√41 - 5) 倍すればよいわけではない。

r1 = 1/2 * (√41 - 5 ) = 0.7015621187164243 であるが，それは b ≒ 1.67539 のときだけである。

ans_r1_(r2 => 1/2, b => 1.67539) |> println

   0.701562024364476

ans_r1_(r2 => 1/2, b => 2) |> println

   0.545454545454546
   0.750000000000000

結論:本問は条件不足で解けない。大円の直径は小円の直径だけでは決まらない。台形の高さの条件が必要である。

function draw(r2, b, more)
   pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
   r1 = 2r2*(b - r2)/b
   a  = 2r2*(b - r2)/(b - 2r2)
   println(b, " ", r1/r2) 
   @printf("r2 = %g;  r1 = %g;  factor = r1/r2 = %g;  a = %g;  b = %g\n", r2, r1, r1/r2, a, b)
   plot([2a, 2a, a, 0, 3a, 2a, a, a],
        [b, 0, b, 0, 0, b, b, 0], color=:green, lw=0.5)
   circle(a, r1, r1)
   circle(2a - r2, b - r2, 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)
       point(a, r1, "大円:r1,(a,r1)", :red, :center, delta=-delta)
       point(2a - r2, b - r2, "小円:r2\n(2a-r2,b-r2)", :blue, :center, delta=-delta)
       point(a, 0, "a", :green, :center, delta=-3delta)
       point(2a, 0, "2a", :green, :center, delta=-3delta)
       point(3a, 0, "3a", :green, :center, delta=-3delta)
       point(a, b, "(a,b)", :green, :left, :bottom, delta=delta/2)
       point(2.1a, 0.8r1, "r2 = $r2\nb = $b\nr1 = $(round(r1, digits=5))", :black, :left, mark=false)
   end
end;

draw(1/2, 1.67539, true)

算額（その1402）

算額（その1402）

三十二 一関市舞川 観福寺内地蔵堂 明治34年(1901)
山村善夫：現存 岩手の算額，昭和52年1月30日，熊谷印刷，盛岡市. http://www.wasan.jp/yamamura/yamamura.html
キーワード：3次元，球11個，高さ
#Julia, #SymPy, #算額, #和算

盤の上に 7 個の等球を互いに外接するように置き，その上に 3 個，更にその上に 1 個の等球を載せる。盤面から一番上の球の頂点までの高さはいかほどか。

 

算額（その1023）の発展版である。
算額（その1023）では，一番上に球が 1 個，それを支える下の層に球が 3 個であった。
本問は更にその下の層にそれを支える 7 個の球がある。
下3個と上1個の球は互いに外接しており，上の層にある球の中心と下の層にある球の中心の z 座標値の差は球の半径 r の 2√6/3 = sqrt(8/3) 倍である。
よって，「高さ」は，一番下の層にある球の下端から一番上の層にある球の上端までの距離である。つまり，r + 2r * 2√6/3 + r = 2.63299316185545 である。

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

@syms r, xa, ya, za, xb, yb, zb, xc, yc, zc
(xa, ya, za) = (0, 0, 0)
(xb, yb, zb) = (2r, 0, 0);

上から x-y 平面に投影した図で，一番下の層の7個の球のうちの 3 つの球の中心座標を求める。
球A の中心座標を (0, 0, 0)
球B の中心座標を (2r, 0, 0)
球C の中心座標を (xc, yc, 0)
とおく。
球C の中心座標は（図を描けば）暗算でもとまるが，あえて方程式を立てて解を求める。
3 つの球が互いに概説しているので，球A と 球C の中心間距離は，三平方の定理により以下のようになり，これを解くと xc = r, yc = √3r　である。

(xc, zc) = (r, 0)
eq1 = (xc - xa)^2 + (yc - ya)^2 + (zc - za)^2- 4r^2
# yc
solve(eq1, yc)[2] |> println

   sqrt(3)*r

x-z 平面に投影した図で，球A, 球B, 球C に載っている 球D の中心座標 (r, yd, zd) を求める。

@syms r, xd, yd, zd
yc = √Sym(3)r
xd = r
yd = r/√Sym(3)
eq2 = (xd - xc)^2 + (yd - yc)^2 + (zd - zc)^2 - 4r^2
# zd
solve(eq2, zd)[2] |> println

   2*sqrt(6)*r/3

zd = 2√6r/3 = sqrt(8/3)r である。
求める「高さ」は r + 2(2√6r/3) + r = 2r(3 + 2√6)/3 で，球の直径の (3 + 2√6)/3 倍である。
球の直径が 1 寸のとき，高さは (3 + 2√6)/3 = 2.632993161855452 寸である。

術は，「球の直径の √(8/3) + 1 倍」としており，同じである。