算額（その751）第二版

算額（その751）第二版

解析解を求めるようにした

三五 大宮市中釘 秋葉神社 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.

埼玉県さいたま市西区中釘 秋葉神社 天保11年(1840)
山口正義：やまぶき，第20号
https://yamabukiwasan.sakura.ne.jp/ymbk20.pdf

キーワード：円5個，不等辺三角形

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

三角形の頂点を (0, 0), (x, y), (x2, 0)
甲円の半径と中心座標を r1, (x1, r1), (x1 + 2r1, r1), (x1 + 4r1, r1)
乙円の半径と中心座標を r2, (x1 + 3r1, y2)
丙円の半径と中心座標を r3, (x1 + r1, y3)
とおき，以下の連立方程式を解く。

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

using SymPy
@syms a, b, c, d,
      r1::positive, x1::positive,
      r2::positive, y2::positive,
      r3::positive, y3::positive,
      x::positive, y::positive, x2::positive;

真ん中の甲円と乙円の中心の垂直距離を a とする。
甲円と乙円が外接することから，a^2 = (r1 + r2)^2 - r1^2 = r2*(2r1 + r2) である。 

eq1 = a^2 + r1^2 - (r1 + r2)^2
ans_a = solve(eq1, a)[2]  # 2 of 2
ans_a |> println

    sqrt(r2)*sqrt(2*r1 + r2)

乙円と丙円の中心の垂直距離を b とする。
乙円と丙円が外接することから b^2 = (r2 + r3)^2 - (2r1)^2 = -4r1^2 + r2^2 + 2r2*r3 + r3^2 である。

eq2 = b^2 + (2r1)^2 - (r2 + r3)^2
ans_b = solve(eq2, b)[2]  # 2 of 2
ans_b |>  println

    sqrt(-4*r1^2 + r2^2 + 2*r2*r3 + r3^2)

真ん中の甲円と丙円の中心の垂直距離を c とする。
甲円と丙円が外接することから，c^2 = (r1 + r3)^2 - r1^2 = r3*(2r1 + r3) である。 

eq3 = c^2 + r1^2 - (r1 + r3)^2
ans_c = solve(eq3, c)[2]  # 2 of 2
ans_c |> println

    sqrt(r3)*sqrt(2*r1 + r3)

a = b + c なのでこれを方程式としてもよいが，SymPy の能力的には両辺を二乗して a^2 = (b + c)^2 = b^2 + 2b*c + c^2 を方程式とするほうがよい。

eq123 = ans_a^2 - (ans_b + ans_c)^2 |> expand
eq123 |> println

    4*r1^2 + 2*r1*r2 - 2*r1*r3 - 2*r2*r3 - 2*sqrt(r3)*sqrt(2*r1 + r3)*sqrt(-4*r1^2 + r2^2 + 2*r2*r3 + r3^2) - 2*r3^2

次に，左の甲円と乙円の中心間の距離を 2 通りの方法で記述する。
左の甲円と乙円の中心を斜辺　d1 とする直角三角形で，d1^2 = a^2 + (3r)^2 = r2*(2r1 + r2) + 9r1^2 である。

@syms d1, d2
eq4 = d1^2 - (ans_a^2 + (3r1)^2)
ans_d1 = solve(eq4, d1)[2]  # 2 of 2
ans_d1 |> println

    sqrt(9*r1^2 + 2*r1*r2 + r2^2)

乙円と甲円の半径の差 (r2 - r1) と，左の甲円と乙円が三角形の斜辺の接点間の距離 2sqrt(r1*r3) + 2sqrt(r2*r3) を直角を挟む 2 辺とする直角三角形において，d2^2 = (2sqrt(r1*r3) + 2sqrt(r2*r3))^2 + (r2 - r1)^2 = 4*r3*(sqrt(r1) + sqrt(r2))^2 + (r1 - r2)^2 である。

eq5 = d2 - ((r2 - r1)^2 + (2sqrt(r2*r3) + 2sqrt(r1*r3))^2)
ans_d2 = solve(eq5, d2)[1]
ans_d2 |> println

    4*r3*(sqrt(r1) + sqrt(r2))^2 + (r1 - r2)^2

d1 = d2 を 2 つ目の方程式とする。

eq45 = (ans_d1)^2 - ans_d2 |> expand |> simplify
eq45 |> println

    -8*sqrt(r1)*sqrt(r2)*r3 + 8*r1^2 + 4*r1*r2 - 4*r1*r3 - 4*r2*r3

eq123, eq45 の 2 つの方程式を解き r2, r3 を求める。

res = solve([eq123, eq45], (r2, r3))[1]

    (2*r1*(sqrt(7) + 4)/9, 2*r1*(sqrt(7) + 13)/(2*sqrt(7) + 17 + 6*sqrt(2)*sqrt(sqrt(7) + 4)))

# r2
ans_r2 = res[1]
ans_r2 |> println
ans_r2(r1 => 10).evalf() |> println

    2*r1*(sqrt(7) + 4)/9
    14.7683362468102

乙円の直径は，甲円の直径の (2√7 + 8)/9 倍である。
甲円の直径が 20 寸のとき，乙円の直径は 29.5366724936204 寸である。

# r3
@syms d
ans_r3 = apart(res[2], d) |> sympy.sqrtdenest |> simplify |> factor
ans_r3 |> println
ans_r3(r1 => 10).evalf() |>  println

    -2*r1*(-3 + sqrt(7))
    7.08497377870819

丙円の直径は，甲円の直径の (6 - 2√7) 倍である。
甲円の直径が 20 寸のとき，乙円の直径は 14.169947557416371 寸である。

乙円，丙円の中心座標は以下のようになる。

# y2
eq3 = r1^2 + (y2 - r1)^2 - (r1 + r2)^2
ans_y2 = solve(eq3, y2)[2]  # 2 of 2
ans_y2 |> println
ans_y2(r2 => ans_r2)(r1 => 10).evalf() |> println

    r1 + sqrt(r2)*sqrt(2*r1 + r2)
    32.6598870349137

左の甲円の中心座標が (x1, r1) のとき，乙円の中心座標は (x1 + 3r1, r1 + sqrt(r2)*sqrt(2r1 + r2)) である。

# y3
eq4 = r1^2 + (y3 - r1)^2 - (r1 + r3)^2
ans_y3 = solve(eq4, y3)[2]  # 2 of 2
ans_y3 |> println
ans_y3(r3 => ans_r3)(r1 => 10).evalf() |> println

    r1 + sqrt(r3)*sqrt(2*r1 + r3)
    23.8526650511426

左の甲円の中心座標が (x1, r1) のとき，丙円の中心座標は (x1 + r1, r1 + sqrt(r3)*sqrt(2r1 + r3)) である。

三角形の頂点の座標を決定するのは複雑である。
別途数値計算を行い，甲円の直径が 20 寸のとき以下を得る。

x1 = 24.534049509916088
x  = 63.120699623258076
y  = 61.70734973660003
x2 = 77.04124269850618

function driver(r1, r2, r3, y2, y3)
    function H(u)
        (x1, x, y, x2) = u
        return [
            dist(0, 0, x, y, x1, r1) - r1^2,
            dist(0, 0, x, y, x1 + r1, y3) - r3^2,
            dist(x, y, x2, 0, x1 + 3r1, y2) - r2^2,
            dist(x, y, x2, 0, x1 + 4r1, r1) - r1^2,
        ]
    end;

    iniv = BigFloat[24, 62.5, 62.5, 78] .* r1/10
    res = nls(H, ini=iniv)
end;

function draw(r1, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r2 = 2r1*(√7 + 4)/9
    r3 = 2r1*(3 - √7)
    y2 = r1 + sqrt(r2)*sqrt(2*r1 + r2)
    y3 = r1 + sqrt(r3)*sqrt(2*r1 + r3)
    res = driver(r1, r2, r3, y2, y3)
    (x1, x, y, x2) = res[1]
    @printf("甲円の直径 = %g;  乙円の直径 = %g;  丙円の直径 = %g\n", 2r1, 2r2, 2r3)
    @printf("x1 = %g;  r2 = %g;  y2 = %g;  r3 = %g;  y3 = %g;  x = %g;  y = %g;  x2 = %g\n", x1, r2, y2, r3, y3, x, y, x2)
    plot([0, x2, x, 0], [0, 0, y, 0], color=:black, lw=0.5)
    circle(x1, r1, r1, :green)
    circle(x1 + 2r1, r1, r1, :green)
    circle(x1 + 4r1, r1, r1, :green)
    circle(x1 + r1, y3, r3, :blue)
    circle(x1 + 3r1, y2, r2, :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)
        plot!(showaxis=false, xlims=(-5delta, x2+5delta))
        point(x, y, "(x,y)", :black, :right, :bottom, delta=delta)
        point(x1, r1, " 甲円:r1\n (x1,r1)", :green, :center, delta=-delta)
        point(x1 + 2r1, r1, "甲円:r1\n(x1+2r1,r1)", :green, :center, delta=-delta)
        point(x1 + 4r1, r1, "甲円:r1\n(x1+4r1,r1)", :green, :center, delta=-delta)
        point(x1 + 3r1, y2, "乙円:r2,(x1+3r1,y2)", :red, :center, :bottom, delta=delta)
        point(x1 + r1, y3, "丙円:r3,(x1+r1,y3) ", :blue, :right, :vcenter)
        point(x2, 0, "(x2,0)", :black, :center, delta=-delta)
        point(x1 + 3r1, r1)
        point(0, 0, "(0,0)", :black, :center, delta=-delta)
        dimension_line(x1, r1, x1 + 3r1, y2, "d", deltax=-2.5delta, length=0.2r1)
        dimension_line(x1 + 3r1, r1, x1 + 3r1, y2, "a", deltax=1.5delta, length=0.2r1)
    end
end;

draw(20/2, true)

 

算額（その1643）

算額（その1643）

三四 武蔵国埼玉郡下忍村遍照院境内 金毘羅社（神楽堂） 天保11年(1840)
埼玉県立図書館：埼玉県史料集 第二集『埼玉の算額』，昭和44年，誠美堂印刷所，埼玉県与野市.
キーワード：楕円3個，二等辺三角形，外円
#Julia, #Julia, #SymPy, #算額, #和算, #数学

円弧の中に等楕円 3 個を容れる。楕円は短径で円弧に接し，長径が最長のものである。円弧の直径が 5 寸，短径が 1 寸のとき，矢の長さ（弦の中点から円弧への最短線分の長さ）を求めよ。

円弧の直径と中心座標を R, (0, 0)
中央の楕円の長半径，短半径を a, b, (0, R - b)
円弧の中心を通る楕円の接線の楕円との接点を (x0, y0)
接線と円弧の弦がなす角を θ
とおく。

「楕円は短径で円弧に接し，長径が最長のもの」とは，円弧は曲率円の一部であり a^2/b = R が成り立つ。

接点 (x0, y0) を求める。

接線と弦のなす角は θ = atand(y0, x0) = 60° である。

中央の楕円を見込む角も 60° であり，弦と接線の間にも等楕円が入り，y 軸を線対称として反対側にも等楕円が入る。

3 つの楕円を見込む角は合計 180° で，弦は円弧をなす円の直径である。

よって，矢は円の半径に当たるので 2.5 寸である。

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

using SymPy
@syms R::positive, a::positive, b::positive,
      x0::positive, y0::positive
a = sqrt(b*R)
eq1 = x0^2/a^2 + (y0 - (R - b))^2/b^2 - 1
eq2 = -b^2*x0/(a^2*(y0 - (R - b))) - y0/x0
res = solve([eq1, eq2], (x0, y0))[2]

    (R*sqrt(b)*sqrt(R - 2*b)/(R - b), R*(R - 2*b)/(R - b))

tan(θ) = y0/x0 から θ を求めると θ = 60° である。

atand(res[2]/res[1])(R => 5/2, b => 1/2).evalf() |> println

    60.0000000000000

function draw(R, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    a = sqrt(b*R)
    (x0, y0) = (R*sqrt(b)*sqrt(R - 2*b)/(R - b), R*(R - 2*b)/(R - b))
    θ = atand(y0, x0)
    plot(showaxis=false)
    circle(0, 0, R, beginangle=0, endangle=180)
    segment(-R, 0, R, 0, :red)
    circle(0, 0, R, :gray, beginangle=180, endangle=360)
    ellipse(0, R - b, a, b, color=:blue)
    ellipse(0, b - R, a, b, color=:gray80)
    ellipse((R - b)*cosd(θ/2), (R - b)*sind(θ/2), a, b, color=:gray80, φ=2θ)
    ellipse(-(R - b)*cosd(θ/2), (R - b)*sind(θ/2), a, b, color=:gray80, φ=θ)
    ellipse((R - b)*cosd(θ/2), -(R - b)*sind(θ/2), a, b, color=:gray80, φ=-2θ)
    ellipse(-(R - b)*cosd(θ/2), -(R - b)*sind(θ/2), a, b, color=:gray80, φ=-θ)
    segment(0, 0, R*cosd(θ), R*sind(θ), :green)
    segment(0, 0, -R*cosd(θ), R*sind(θ), :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(x0, y0, "(x0,y0)", :blue, :left, delta=-delta/2)
        point(0, R - b, "楕円:a,b\n(0,R-b)", :blue, :right, :vcenter, deltax=-delta)
        circle(0, 0, 0.1R, beginangle=0, endangle=θ)
        circle(0, 0, 0.12R, beginangle=0, endangle=θ)
        point(0.12R, 0.1R, "θ", :red, mark=false)
        dimension_line(0, 0, 0, R, " 矢", :magenta, :left, :vcenter)
    end
end;

draw(5/2, 0.5, true)

2025/03/31

肺炎のため緊急入院となりました

退院までブログもお休みです

算額（その1673）

算額（その1673）

長野県山内町 渋医王殿 文化15年(1818)
中村信弥「幻の算額」（319）
http://www.wasan.jp/maborosi/maborosi3-3.pdf
キーワード：整数係数式、係数推定、最小値を取るｘの値
#Julia, #SymPy, #算額, #和算，#数学

x についての 3 次整数係数式がある。
f(x) = a + b*x + c*x^2 + d*x^3
f（1）= 58,f（3）= 42,f(5) = 10, f(9) = 90
このとき，各項の係数と f(x) が極小値をとるときの x をもとめよ。

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

using SymPy

@syms a, b, c, d, x
f = a + b*x + c*x^2 + d*x^3
eq1 = f(x => 1) - 58
eq2 = f(x => 3) - 42
eq3 = f(x => 5) - 10
eq4 = f(x => 9) - 90;

res1 = solve([eq1, eq2, eq3, eq4], (a, b, c, d))

係数は a = 45，b = 23，c = -11，d = 1 である。

極小値をとる x は，f の導関数をとり，導関数 = 0 を解く。2番目のものが適解である。

res2 = solve(diff(f, x), x)[2]

b, c, d に実値を入れれば x =6.07036751697599 を得る。

res2(b => 23, c => -11, d => 1).evalf() |> println

6.07036751697599

ちなみにそのときの f の値は 2.96464202605517 である。

(x^3 - 11*x^2 + 23*x + 45)(x => 6.07036751697599).evalf() |> println

2.96464202605517

plot(x^3 - 11*x^2 + 23*x + 45, xlims=(0, 10), label="f(x)")
plot!(23 + -22x + 3x^2, xlims=(0, 10), label="g(x)")
vline!([6.07036751697599], label="")

はなまるうどん 発祥の地 木太町店

高松市木田町 セルフ はなまるうどん
香川県内の1号店。最寄り駅は林道だが，はなまるうどん駅になった。

算額（その1672）

算額（その1672）

和算で遊ぼう！！ 「三春まちなか寺子屋」2017レポート
９月 狐田稲荷神社 
https://miharukoma.com/wp-content/uploads/2018/01/%E4%B8%89%E6%98%A5%E3%81%BE%E3%81%A1%E3%81%AA%E3%81%8B%E5%AF%BA%E5%AD%90%E5%B1%8B2017%E3%83%AC%E3%83%9D%E3%83%BC%E3%83%88.pdf
キーワード：中国剰余定理
#Julia, #SymPy, #算額, #和算，#数学

狐が田植えをする。苗の束数はわからないが，5 把ずつ植えると 1 把余り，7
 把ずつ植えると 2 把余るという。総束数をもとめよ。

using SymPy

#from sympy.ntheory.modular import crt

# 剰余の条件
moduli = [5, 7]   # 除数のセット
remainders = [1, 2]  # 剰余のセット

# 中国剰余定理を使用して解を求める
solution = sympy.ntheory.modular.crt(moduli, remainders)  # 最小解と周期が返される

    (16, 35)

総束数は 16 把である。

算額（その1671）

算額（その1671）

長野県錯視臼田町広河原 禅昌寺 文政9年(1826)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
キーワード：球7個，減弧円錐，3次元
#Julia, #SymPy, #算額, #和算，#数学

減弧円錐の中に甲球 3 個，乙球 1 個，丙球 3 個を容れる。甲球は互いに外接し，底面と側面に接し，乙球とも外接する。丙球は乙球と外接し，側面とも接している。乙球は甲球，丙球と外接し，側面にも接している。
丙球の直径が甲球の直径の 1/5(2 分)のとき，乙球の直径の最小値を求めよ。

減弧円錐とは，左側の図において，赤で描いた 2 個の円を z 軸を中心として回転させたときにできる先の尖った円錐状の立体である。

大球の半径と中心座標を R, (R, 0, R)
甲球の半径と中心座標を r1, (2r1/√3, 0, r1)
乙球の半径と中心座標を r2, (0, 0, z2)
丙球の半径と中心座標を r3, ((2r3/√3, 0, z3)
とおき，以下の連立方程式を解く。

r1, r2, z2, r3, z3 は R の関数として得られる。

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

using SymPy

@syms R, r1, r2, z2, r3, z3
eq1 = (R -2r1/√Sym(3))^2 + (R - r1)^2 - (R + r1)^2
eq2 = R^2 + (R - z2)^2 - (R + r2)^2
eq3 = (2r1/√Sym(3))^2 + (z2 - r1)^2 - (r1 + r2)^2
eq4 = (R - 2r3/√Sym(3))^2 + (R - z3)^2 - (R + r3)^2
eq5 = (2r3/√Sym(3))^2 + (z3 - z2)^2 - (r2 + r3)^2;

ans_r1 = solve(eq1, r1)[2]  # 2 of 2
ans_r1 |> println

    R*(-sqrt(9 + 6*sqrt(3)) + sqrt(3) + 3)/2

eq12 = eq2(r1 => ans_r1)
eq13 = eq3(r1 => ans_r1);
res = solve([eq12, eq13], (r2, z2))[1]

    (R*(-1940*sqrt(3 + 2*sqrt(3)) - 1120*sqrt(9 + 6*sqrt(3)) + 2848*sqrt(3) + 4933)/(4*(-277*sqrt(3)*sqrt(3 + 2*sqrt(3)) - 478*sqrt(3 + 2*sqrt(3)) + 1218 + 704*sqrt(3))), R*(-427*sqrt(9 + 6*sqrt(3)) - 736*sqrt(3 + 2*sqrt(3)) + 1081*sqrt(3) + 1877)/(4*(-39*sqrt(3)*sqrt(3 + 2*sqrt(3)) - 67*sqrt(3 + 2*sqrt(3)) + 98*sqrt(3) + 171)))

# r2
@syms d
ans_r2 = apart(res[1]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_r2 |> println

    R*(-37*sqrt(2)*3^(3/4)/552 - 21*sqrt(2)*3^(1/4)/184 + 13/46 + 4*sqrt(3)/23)

# z2
ans_z2 = apart(res[2]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_z2 |> println

    R*(-67*sqrt(2)*3^(1/4)/184 - 83*sqrt(2)*3^(3/4)/552 + 4*sqrt(3)/23 + 59/46)

eq15 = eq5(r2 => ans_r2, z2 => ans_z2) |> simplify
eq14 = eq4 |> expand
res2 = solve([eq14, eq15], (r3, z3))[2]  # 2 of 2

    ((-192*sqrt(3)*R - 312*R + 120*sqrt(2)*3^(3/4)*R + 264*sqrt(2)*3^(1/4)*R + (3*R*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3))/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)) + 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)))*(-201*sqrt(2)*3^(1/4) - 83*sqrt(2)*3^(3/4) + 156 + 96*sqrt(3)))/(63*sqrt(2)*3^(1/4) + 37*sqrt(2)*3^(3/4) + 396 + 272*sqrt(3)), 3*R*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3))/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)) + 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675)))

# r3
@syms d
ans_r3 = apart(res2[1]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_r3 |> println

    R*(-1457609*sqrt(2)*3^(3/4) - 2524500*sqrt(2)*3^(1/4) + 3^(1/4)*(-201 - 83*sqrt(3))*sqrt(-267390592*sqrt(2)*3^(3/4) - 463133856*sqrt(2)*3^(1/4) + 995984138 + 575032436*sqrt(3))/2 + 78*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 48*sqrt(-401085888*sqrt(2)*3^(3/4) - 694700784*sqrt(2)*3^(1/4) + 1493976207 + 862548654*sqrt(3)) + 6137076 + 3544118*sqrt(3))/(2618169*sqrt(2)*3^(1/4) + 1512127*sqrt(2)*3^(3/4) + 9518764 + 5497344*sqrt(3))

# z3
ans_z3 = apart(res2[2]/R, d) |> simplify |> sympy.sqrtdenest |> simplify |> x -> x*R
ans_z3 |> println

    3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675))

R = 1 のとき，r1 = 0.164191;  r2 = 0.155332;  z2 = 0.421387;  r3 = 0.0357628;  z3 = 0.607967 となる。

r3/r1 = 0.2178121821537112 である。
答えは，乙球径 = 1.32678 である。なぜこのような答えになるのか？？

function draw(R)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = R*(-sqrt(9 + 6*sqrt(3)) + sqrt(3) + 3)/2
    r2 = R*(-37*sqrt(2)*3^(3/4)/552 - 21*sqrt(2)*3^(1/4)/184 + 13/46 + 4*sqrt(3)/23)
    z2 = R*(-67*sqrt(2)*3^(1/4)/184 - 83*sqrt(2)*3^(3/4)/552 + 4*sqrt(3)/23 + 59/46)
    r3 = R*(-1457609*sqrt(2)*3^(3/4) - 2524500*sqrt(2)*3^(1/4) + 3^(1/4)*(-201 - 83*sqrt(3))*sqrt(-267390592*sqrt(2)*3^(3/4) - 463133856*sqrt(2)*3^(1/4) + 995984138 + 575032436*sqrt(3))/2 + 78*sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 48*sqrt(-401085888*sqrt(2)*3^(3/4) - 694700784*sqrt(2)*3^(1/4) + 1493976207 + 862548654*sqrt(3)) + 6137076 + 3544118*sqrt(3))/(2618169*sqrt(2)*3^(1/4) + 1512127*sqrt(2)*3^(3/4) + 9518764 + 5497344*sqrt(3))
    z3 = 3*R*(-5537*sqrt(2)*3^(3/4) - 9587*sqrt(2)*3^(1/4) + sqrt(-133695296*sqrt(2)*3^(3/4) - 231566928*sqrt(2)*3^(1/4) + 497992069 + 287516218*sqrt(3)) + 16095*sqrt(3) + 27957)/(2*(2630*sqrt(2)*3^(3/4) + 4590*sqrt(2)*3^(1/4) + 17642*sqrt(3) + 30675))
    @printf("R = %g;  r1 = %g;  r2 = %g;  z2 = %g;  r3 = %g;  z3 = %g\n", R, r1, r2, z2, r3, z3)
    p1 = plot(xlabel="x-axis", ylabel="z-axis")
    circle2(R, R, R)
    circle(2r1/√3, r1, r1, :blue)
    circle(-r1/√3, r1, r1, :blue)
    circle(0, z2, r2, :green)
    circle(2r3/√3, z3, r3, :magenta)
    circle(-r3/√3, z3, r3, :magenta)
    point(0.8, 0.9, "甲円:r2,(2r1/√3,0,r1)", :blue, :right, :bottom, delta=0.1R, mark=false)
    point(0.8, 1.1, "乙円:r2,(0,0,z2)", :green, :right, :bottom, delta=0.1R, mark=false)
    point(0.8, 1.3, "丙円:r3,(2r3/√3,0,z3)", :magenta, :right, :bottom, delta=0.08R, mark=false)
    point(R, R, "大円:R,(R,R)", :red, :right, delta=-0.1R)
    plot!(xlims=(-1, 1), ylims = (0, 1))
p2 = plot(xlabel="x-axis", ylabel="y-axis")
    circle2(R, 0, R)
    rotate(2r1/√3, 0, r1, :blue)
    circle(0, 0, r2, :green)
    rotate(2r3/√3, 0, r3, :magenta)
    plot!(xlims=(-0.5, 0.5), ylims = (-0.5, 0.5))
    plot(p1, p2)
end;

draw(1)

讃岐うどん ぼっこ屋

高松市三谷町 ぼっこ屋 三谷店
やや太麺。
特筆すべきは，「しっかり噛もうと思わなければ噛めない太麺」...　今までなかった（少ない経験ながら）

 

算額（その1670） 改訂版

算額（その1670）

直接連立方程式を解くようにした

長野県小諸市八幡町 八幡神社 寛政11年(1799)
中村信弥「改訂増補 長野県の算額」
http://www.wasan.jp/zoho/zoho.html
キーワード：等差数列，数列の一般項，部分和
#Julia, #SymPy, #算額, #和算，#数学

n 個の立方体がある。各立方体の 1 辺は等差級数をなし，公差が r とする。
n 個の立方体の 1 辺の長さの和が A，体積の和が B のとき，立方体の個数を求める術を述べよ。

注：A, B だけではなく r（または a） も与えないと解は求まらない。

using SymPy

@syms x, r, a, n, A, B

    (x, r, a, n, A, B)

各立方体の 1 辺の長さは，初項 a, 公差 r の等差数列である。
summation() は，一般項と，変数，範囲を与えると和の公式を返す。
和が A であるという方程式 eq1 を立てる。

eq1 = A - summation(a + r*(x - 1), (x, 1, n));

各立方体の体積の一般項は (a + r*(x - 1))^3 である。
体積の和が B であるという方程式 eq2 を立てる。

eq2 = B - summation((a + (x - 1)r)^3, (x, 1, n));

連立方程式を解き，n, a を求める。

res = solve([eq1, eq2], (n, a))[4]

    (-sqrt(A*(-4*A^2*r + A*r^2 + 4*B))/(2*A*r) + sqrt(A*(4*A^2*r + A*r^2 + 4*B))/(2*A*r), r/2 + sqrt(A*(-4*A^2*r + A*r^2 + 4*B))/(2*A))

# n
res[1] |> simplify

res[1](A => 33, B => 2775.96, r => 2.6).evalf() |> println

    5.00000000000000

# a
res[2] |> simplify

res[2](A => 33, B => 2775.96, r => 2.6).evalf() |> println

    1.40000000000000

末尾に示すテストデータを作るプログラムにより，初項 a = 1.4, 公差 r = 2.6 での結果，第 5 項までで A = 33, B = 2775.96 を代入すると確かに n = 5 になる。
また，初項は a = 1.4 であることもわかる。

A, B, a を与えて n, r を求める

連立方程式による方法は，何を未知数とするかを選ぶことができることである。

同じ連立方程式で，A, B, a を与えて n, r を求める解を求めることも簡単に指定できる。

res2 = solve([eq1, eq2], (n, r))[1]

    (-(-2*A^2*a + A*a^2 + B)/(2*(A*a^2 - B)) - sqrt(-4*A^4*a^2 + 8*A^3*B + 4*A^3*a^3 - 12*A^2*B*a + A^2*a^4 + 2*A*B*a^2 + B^2)/(2*(A*a^2 - B)), (A*a^2 - B)/(A*(-A + a)))

# n
res2[1] |> simplify

# r
res2[1](A => 33, B => 2775.96, a => 1.4).evalf() |> println

    5.00000000000000

res2[2] |> simplify

res2[2](A => 33, B => 2775.96, a => 1.4).evalf() |> println

    2.60000000000000

初項 a = 1.4, 公差 r = 2.6 での結果，第 5 項までで A = 33, B = 2775.96 を代入すると確かに n = 5，a = 1.4 であることがわかる。

テストデータ生成プログラム（チェック付き）

function example(a, r, n = 10)
    setprecision(BigFloat, 120)
    println("a = $a\nr = $r\n")
    A = big"0"
    B = big"0"
    for i = 1:n
        ai = a + (i - 1)r
        A += ai
        B += ai^3
        個数 = sqrt(big"2")*sqrt(1 + 4*B/(A*r^2) - sqrt(-16*A^4*r^2 + A^2*r^4 + 8*A*B*r^2 + 16*B^2)/(A*r^2))/2
        println("i = $i\nai = $ai\nA = $A\nB = $B\n個数 = $個数\n")
    end
    setprecision(BigFloat, 256)
    return  # return nothing
end;
example(big"1.4", big"2.6")

    a = 1.399999999999999999999999999999999999999999999999999999999999999999999999999997
    r = 2.599999999999999999999999999999999999999999999999999999999999999999999999999986

    i = 1
    ai = 1.3999999999999999999999999999999999997
    A = 1.3999999999999999999999999999999999997
    B = 2.7439999999999999999999999999999999996
    個数 = 0.9999999999999999999999999999999999985

    i = 2
    ai = 4.0
    A = 5.4000000000000000000000000000000000012
    B = 66.744000000000000000000000000000000015
    個数 = 2.000000000000000000000000000000000003

    i = 3
    ai = 6.5999999999999999999999999999999999988
    A = 12.0
    B = 354.2399999999999999999999999999999998
    個数 = 3.0000000000000000000000000000000000451

    i = 4
    ai = 9.1999999999999999999999999999999999976
    A = 21.20000000000000000000000000000000001
    B = 1132.9279999999999999999999999999999998
    個数 = 4.0000000000000000000000000000000002106

    i = 5
    ai = 11.79999999999999999999999999999999999
    A = 33.0
    B = 2775.9599999999999999999999999999999946
    個数 = 5.0000000000000000000000000000000002407

    i = 6
    ai = 14.399999999999999999999999999999999995
    A = 47.400000000000000000000000000000000019
    B = 5761.9439999999999999999999999999999949
    個数 = 6.0000000000000000000000000000000004815

    i = 7
    ai = 17.0
    A = 64.400000000000000000000000000000000019
    B = 10674.943999999999999999999999999999995
    個数 = 7.0000000000000000000000000000000005657

    i = 8
    ai = 19.600000000000000000000000000000000005
    A = 84.0
    B = 18204.480000000000000000000000000000016
    個数 = 7.9999999999999999999999999999999995727

    i = 9
    ai = 22.199999999999999999999999999999999986
    A = 106.19999999999999999999999999999999996
    B = 29145.527999999999999999999999999999978
    個数 = 9.0

    i = 10
    ai = 24.79999999999999999999999999999999999
    A = 131.0
    B = 44398.519999999999999999999999999999984
    個数 = 10.000000000000000000000000000000000373

 

算額（その1669）

算額（その1669）

百十一 高崎市倉賀野町 倉賀野神社 慶応3年(1867)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.

群馬の和算 111(2/2) 倉賀野神社
http://takasakiwasan.web.fc2.com/gunnsann/g111_1.html

キーワード：円弧，弦，矢，弧背
#Julia, #SymPy, #算額, #和算，#数学

矢が 1 寸，弦が 4 寸の円弧がある。円弧の周長（弧背）を求めよ。

矢，弦はそのまま変数名とする。
円弧の半径と中心座標を R, (0, 0) とおき，以下の方程式を解く。

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

using SymPy

@syms 矢, 弦, R
eq = (R - 矢)^2 + (弦/2)^2 - R^2
ans_R = solve(eq, R)[1]
ans_R |> println

    弦^2/(8*矢) + 矢/2

円弧の中心角を θ とすると θ = acosd((R - 1)/R) である。
弧背は直径 2R の円の円周の 2θ/360 である。

弧背 = 2(弦^2/(8*矢) + 矢/2)*pi*2acosd((ans_R - 1)/ans_R)/360 |> simplify
弧背 |> println

    (弦^2 + 4*矢^2)*acos((弦^2 + 4*矢*(矢 - 2))/(弦^2 + 4*矢^2))/(4*矢)

矢が 1，弦が 4 のとき，弧背は 4.63647609000806 である。

注：復元算額の答の中の「□□□」は欠損文字を表したつもりかもしれないが「000」である。欠損文字の表記「◯◯◯」と見間違えたのかもしれない。

弧背(矢 => 1, 弦 => 4).evalf() |> println

    4.63647609000806

function draw(矢, 弦, more=false)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = 弦^2/(8矢) + 矢/2
    θ = acosd((R - 1)/R)
    弧背 = 2R*pi*2θ/360
    @printf("矢が %g，弦が %g のとき，弧背は %.15g である。\n", 矢, 弦, 弧背)
    plot()
    circle(0, 0, R, beginangle=90 - θ, endangle=90 + θ)
    y = R - 1
    x = sqrt(R^2 - y^2)
    segment(-x, y, x, y)
    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)
        dimension_line(0, R - 1, 0, R, "矢", :green, :left, length=10delta)
        dimension_line(-x, y, x, y, "弦", :blue, :center, delta=-5delta, length=10delta)
        point(0, R, "R", :red, :center, :bottom, delta=2delta)
        plot!([-x, 0, x], [y, 0, y], color=:black, lw=0.5)
        circle(0, 0, 0.08R, :magenta, beginangle=90 - θ, endangle=90 + θ)
        circle(0, 0, 0.09R, :magenta, beginangle=90 - θ, endangle=90 + θ)
        point(0, 0.1R, "2θ", :magenta, :center, :bottom, delta=delta, mark=false)
        
    end
end;

draw(1, 4, true)

 

算額（その1668）

算額（その1668）

栃木県佐野市 星宮神社 天和3年(1683) 現存するもので一番古い算額
山口正義：やまぶき 第 47 号
https://yamabukiwasan.sakura.ne.jp/ymbk47.pdf
山口正義：やまぶき 第 48 号
https://yamabukiwasan.sakura.ne.jp/ymbk48.pdf
キーワード：直角三角形，鈎，股，弦，内接円，外積
#Julia #SymPy #算額 #和算 #数学

直角三角形において，
(1) A = ((弦 - 鈎)^(1//7) + (弦 - 股)^(1//5) + 円径^(1//3))^9 + 外積
(2) B = 股 - 鈎
の 2 つが与えられたときに，各長さ，面積（鈎，股，弦，円径，外積）はいかほどか。

以下は山口による説明

これを解析的に解くのは相当に難しい。

数値的に解こう。与えられる A, B から「鈎」，「股」を求める。
なお，普通の直角三角形（？）において，A は B の数万〜数十万倍である。

A には，「円径」，「外積」が含まれるが，いずれも「鈎」と「股」で計算されるものである。
「鈎」，「股」の初期値は実験により得られるデータ対の解析に基づいて，A, B の関数として与える。

まず，A, B の連立方程式を記述する。

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

using SymPy

@syms A, B, 鈎, 股, 弦, 円径, 外積
弦 = sqrt(鈎^2 + 股^2)
円径 = 鈎 + 股 - 弦
外積 = 鈎*股/2 - π*(円径/2)^2
eq1 = A - (((弦 - 鈎)^(1//7) + (弦 - 股)^(1//5) + 円径^(1//3))^9 + 外積)
eq2 = B - (股 - 鈎);

A, B から「鈎」，「股」を求める関数を定義する。

「鈎」，「股」の初期値「pred_鈎」，「pred_股」については後述する。

function driver(A, B)
    function H(u)
        (鈎, 股) = u
        return [
            A - 股*鈎/2 + pi*(股/2 + 鈎/2 - sqrt(股^2 + 鈎^2)/2)^2 - ((-股 + sqrt(股^2 + 鈎^2))^(1/5) + (-鈎 + sqrt(股^2 + 鈎^2))^(1/7) + (股 + 鈎 - sqrt(股^2 + 鈎^2))^(1/3))^9,  # eq1
            B - 股 + 鈎,  # eq2
        ]
    end;
    pred_鈎 = 1.497 + 2.337e-05*A - 4.006e-11*A^2
    pred_股 = 4.107397 + 0.934185*B - 0.003609*B^2
    iniv = BigFloat[pred_鈎, pred_股]  # 適当な初期値でも収束することが多い
    return nls(H, ini=iniv)
end;

A = 60000, B = 4 の場合の「鈎」，「股」を求める。

res = driver(60000, 4)

    ([2.6226253736776437, 6.622625373677644], true)

出力項目の最後が　true ならば，数値解は正常に収束したものであることを示している。

検算は不要であるが，念の為に検算してみる。
まず，「弦」，「円径」，「外積」を求める。

鈎 = 2.6226253736776437
股 = 6.622625373677644
弦 = sqrt(鈎^2 + 股^2)
円径 = 鈎 + 股 - 弦
外積 = 鈎*股/2 - π*(円径/2)^2
(弦, 円径, 外積)

    (7.123014157695936, 2.1222365896593516, 5.146987197425817)

次いで，A, B が最初に与えた値と同じになることを確認する。

# A
((弦 - 鈎)^(1//7) + (弦 - 股)^(1//5) + 円径^(1//3))^9 + 外積

    60000.0

# B
股 - 鈎

    4.0

「鈎」，「股」の初期値「pred_鈎」，「pred_股」について

以下のようなプログラムで，A, B に対して得られる「鈎」，「股」，「弦」，「円径」，「外積」を記録する。

普通の直角三角形（？）において，A は B の数万〜数十万倍である。

println("A B 鈎 股 弦 円径 外積")
for i = 1:20
    for j = (i+1:20) .* 10000
        res = driver(j, i)
        if res[2]
            鈎 = res[1][1]
            股 = res[1][2]
            弦 = sqrt(鈎^2 + 股^2)
            円径 = 鈎 + 股 - 弦
            外積 = 鈎*股/2 - π*(円径/2)^2
            @printf("%d %d %g %g %g %g %g\n", j, i, 鈎, 股, 弦, 円径, 外積)
        #else
        #    @printf("%d %d NA NA NA NA NA\n", j, i)
        end
    end
end

    A B 鈎 股 弦 円径 外積
    20000 1 1.73963 2.73963 3.24529 1.23398 1.18706
    30000 1 2.16956 3.16956 3.84098 1.49815 1.67551
    40000 1 2.53117 3.53117 4.34464 1.71769 2.1517
    50000 1 2.84865 3.84865 4.7882 1.90909 2.61923
    60000 1 3.13462 4.13462 5.18854 2.0807 3.07999
    70000 1 3.39665 4.39665 5.55588 2.23742 3.53519
    80000 1 3.6397 4.6397 5.89697 2.38243 3.98566
    90000 1 3.86723 4.86723 6.21654 2.51792 4.43199
    100000 1 4.08176 5.08176 6.51806 2.64546 4.87467
 　　　　：
    200000 17 4.09938 21.0994 21.4939 3.70483 32.4669
    190000 18 3.94969 21.9497 22.3022 3.59716 33.1845
    200000 18 4.06385 22.0639 22.435 3.69272 34.1223
    200000 19 4.03036 23.0304 23.3804 3.68036 35.7721

変数間の相関係数は以下のようになる。

A は「鈎」と相関が強い。
B は「股」と相関が強い。

相関係数行列

             A           B          鈎        股        弦      円径      外積
A    1.0000000  0.50000000  0.82693238 0.6276680 0.6762419 0.9899246 0.7061528
B    0.5000000  1.00000000 -0.03699511 0.9873807 0.9732902 0.4963393 0.9586976
鈎   0.8269324 -0.03699511  1.00000000 0.1217277 0.1920133 0.8254176 0.2270238
股   0.6276680  0.98738074  0.12172773 1.0000000 0.9971431 0.6237921 0.9881973
弦   0.6762419  0.97329023  0.19201330 0.9971431 1.0000000 0.6706341 0.9940134
円径 0.9899246  0.49633933  0.82541765 0.6237921 0.6706341 1.0000000 0.6897654
外積 0.7061528  0.95869763  0.22702381 0.9881973 0.9940134 0.6897654 1.0000000

「鈎」を A の二次関数で予測する

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.497e+00  2.100e-01   7.127 2.17e-11
A            2.337e-05  3.497e-06   6.682 2.63e-10
I(A^2)      -4.006e-11  1.354e-11  -2.960  0.00348

「股」を B の二次関数で予測する

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.107397   0.148959  27.574   <2e-16
B           0.934185   0.042244  22.114   <2e-16
I(B^2)      0.003609   0.002447   1.475    0.142

予測式は次のようになる。

pred_鈎 = 1.497 + 2.337e-05*A - 4.006e-11*A^2
pred_股 = 4.107397 + 0.934185*B - 0.003609*B^2

算額（その1667）

算額（その1667）

栃木県佐野市 星宮神社 天和3年(1683) 現存するもので一番古い算額
山口正義：やまぶき 第 47 号
https://yamabukiwasan.sakura.ne.jp/ymbk47.pdf
山口正義：やまぶき 第 48 号
https://yamabukiwasan.sakura.ne.jp/ymbk48.pdf
キーワード：直角三角形，鈎，股，弦，中鈎（弦を底辺としたときの高さ），方面（正方形の一辺）
#Julia #SymPy #算額 #和算 #数学

直角三角形において，
(1) 和 = (鈎 + 弦)^(1/7)
(2) 相乗 = 方面^3 * 中鈎^(1/4)
の 2 つが与えられたときに，各長さ（鈎，股，弦，中鈎，方面）はいかほどか。

以下は山口による説明

これを解析的に解くのは相当に難しい。

数値的に解こう。与えられる「和」と「相乗」から「鈎」，「股」を求める。
「和」と「相乗」には，「弦」，「方面」，「中鈎」が含まれるが，いずれも「鈎」と「股」で計算されるものである。
「鈎」，「股」の初期値は実験により得られるデータ対の解析に基づいて，「和」と「相乗」の関数として与える。

まず，「和」，「相乗」の連立方程式を記述する。

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

using SymPy

@syms 和, 相乗, 鈎, 股, 弦, 方面, 中鈎
弦 = sqrt(鈎^2 + 股^2)
中鈎 = 鈎*股/弦
方面 = 鈎*股/(鈎 + 股)
eq1 = 和 - ((鈎 + 弦)^(1//7) + (股^9)^(1//5))
eq2 = 相乗 - 方面^3 * 中鈎^(1//4);

「和」，「相乗」から「鈎」，「股」を求める関数を定義する。

「鈎」，「股」の初期値「pred_鈎」，「pred_股」については後述する。

function driver(和, 相乗)
    function H(u)
        (鈎, 股) = u
        return [
            和 - (鈎 + sqrt(股^2 + 鈎^2))^(1/7) - (股^9)^(1/5),  # eq1
            相乗 - 股^3*鈎^3*(股*鈎/sqrt(股^2 + 鈎^2))^(1/4)/(股 + 鈎)^3,  # eq2
        ]
    end;
    pred_鈎 = 1.494191 + 0.331818*相乗 - 0.009864*相乗^2
    pred_股 = 0.6073778 + 0.3156569*和 - 0.0046705*和^2
    iniv = BigFloat[pred_鈎, pred_股]  # 適当な初期値でも収束することが多い
    return nls(H, ini=iniv)
end;

和 = 7, 相乗 = 6 の場合の「鈎」，「股」を求める。

res = driver(7, 6)

    ([4.856731010389481, 2.6049699305730556], true)

出力項目の最後が　true ならば，数値解は正常に収束したものであることを示している。

検算は不要であるが，念の為に検算してみる。
まず，「弦」，「中鈎」，「方面」を求める。

鈎 = 4.856731010389481
股 = 2.6049699305730556
弦 = sqrt(鈎^2 + 股^2)
中鈎 = 鈎*股/弦
方面 = 鈎*股/(鈎 + 股)
(弦, 中鈎, 方面)

    (5.511234385005651, 2.2956088163057355, 1.6955434616110827)

次いで，「和」，「相乗」が最初に与えた値と同じになることを確認する。

# 和
((鈎 + 弦)^(1//7) + (股^9)^(1//5))

    7.0

# 相乗
方面^3 * 中鈎^(1//4)

    6.0

「鈎」，「股」の初期値「pred_鈎」，「pred_股」について

以下のようなプログラムで，「和」，「相乗」に対して得られる「鈎」，「股」，「弦」，「中鈎」，「方面」を記録する。

println("和 相乗 鈎 股 弦 中鈎 方面")
for i = 1:20
    for j = i+1:20
        res = driver(j, i)
        if res[2]
            鈎 = res[1][1]
            股 = res[1][2]
            弦 = sqrt(鈎^2 + 股^2)
            中鈎 = 鈎*股/弦
            方面 = 鈎*股/(鈎 + 股)
            @printf("%d %d %g %g %g %g %g\n", j, i, 鈎, 股, 弦, 中鈎, 方面)
        #else
        #    @printf("%d %d NA NA NA NA NA\n", j, i)
        end
    end
end

    和 相乗 鈎 股 弦 中鈎 方面
    3 1 3.78904 1.32471 4.01393 1.25049 0.981545
    4 1 2.19561 1.75111 2.8084 1.36902 0.974165
    5 1 1.8266 2.08603 2.77272 1.37422 0.973857
    6 1 1.65184 2.37886 2.89613 1.35681 0.974892
    7 1 1.54747 2.64416 3.0637 1.33556 0.976176
    8 1 1.47713 2.88926 3.24496 1.31522 0.977425
　　　　　：
    19 17 4.42182 4.91842 6.61388 3.2883 2.32846
    20 17 4.30577 5.07205 6.65323 3.28248 2.3288
    19 18 4.57333 4.91771 6.71559 3.34897 2.36964
    20 18 4.44892 5.0714 6.74626 3.34441 2.3699
    20 19 4.59103 5.07076 6.84034 3.40335 2.40949

変数間の相関係数は以下のようになる。

「和」は「股」と相関が強い。
「相乗」は「鈎」と相関が強い。

以下は相関係数行列である。

              和      相乗          鈎         股        弦      中鈎      方面
和    1.00000000 0.4881569 -0.08782835  0.9937689 0.6306035 0.4246172 0.4682552
相乗  0.48815686 1.0000000  0.71425887  0.4793127 0.8203931 0.9639511 0.9671184
鈎   -0.08782835 0.7142589  1.00000000 -0.1157736 0.6964043 0.7449299 0.7503883
股    0.99376890 0.4793127 -0.11577358  1.0000000 0.6031550 0.4265652 0.4664547
弦    0.63060351 0.8203931  0.69640434  0.6031550 1.0000000 0.7885511 0.8339510
中鈎  0.42461720 0.9639511  0.74492991  0.4265652 0.7885511 1.0000000 0.9966916
方面  0.46825522 0.9671184  0.75038826  0.4664547 0.8339510 0.9966916 1.0000000

「股」を「和」の二次関数で予測する

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.6073778  0.0161933   37.51   <2e-16
和           0.3156569  0.0026444  119.37   <2e-16
I(和^2)     -0.0046705  0.0001006  -46.45   <2e-16

縦軸が予測値 「pred_股」，横軸は「股」の実値

「鈎」を「相乗」の二次関数で予測する

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.494191   0.150165   9.950  < 2e-16
相乗         0.331818   0.042257   7.852 3.21e-13
I(相乗^2)   -0.009864   0.002438  -4.045 7.67e-05

縦軸が予測値 「pred_鈎」，横軸は「鈎」の実値

予測式は次のようになる。これらを数値解を求める関数 driver 内で初期値として用いる。

pred_鈎 = 1.494191 + 0.331818*相乗 - 0.009864*相乗^2
pred_股 = 0.6073778 + 0.3156569*和 - 0.0046705*和^2

手打 うつ海うどん

高松市下田井町 手打 うつ海うどん
シッカリした太麺

算額（その1666）

算額（その1666）

栃木県佐野市 星宮神社 天和3年(1683) 現存するもので一番古い算額
山口正義：やまぶき 第 47 号
https://yamabukiwasan.sakura.ne.jp/ymbk47.pdf
山口正義：やまぶき 第 48 号
https://yamabukiwasan.sakura.ne.jp/ymbk48.pdf
キーワード：長方形，面積，分割
#Julia #SymPy #算額 #和算 #数学

縦 120 間，横 21.5 間の土地がある。幅 3 間の道を付けて三箇所（前，左，右）の面積を等しくするとき，それぞれの縦，横はいかほどか。

注：縦横は当時は今と逆であった。すなわち，現代では横が120 間，縦が 21.5 間で，前，左，右の長方形の縦横も同じである。本問においては，縦，横は当時のとおりとする。

当時は難問として知られていたらしいが，連立方程式を解けばいとも簡単。
前，左，右の面積をそのまま変数名とする
土地の縦，横を a, b
左の縦，横を c, d，道幅を w
として以下の連理方程式を解く。

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

using SymPy

@syms a, b, w, c, d, 左, 右, 前
左 = c*d
右 = (a - c - w)*(d + w)
前 = a*(b - d - w)
eq1 = 左 - 右
eq2 = 左 - 前;

res = solve([eq1, eq2], (c, d))[1];  # 1 of 2

# c 左の縦
c = res[1] |> simplify
c |> println

    -a*b/w + a - w/2 + sqrt(4*a^2*b^2 - 4*a^2*b*w + 4*a^2*w^2 - 4*a*w^3 + w^4)/(2*w)

# d 左の横
d = res[2] |> simplify
d |> println

    (2*a*b - 4*a*w + w^2 + sqrt(4*a^2*b^2 - 4*a^2*b*w + 4*a^2*w^2 - 4*a*w^3 + w^4))/(2*(3*a - w))

c(a => 120, b => 21.5, w => 3).evalf() |> println

    64.9999999999999

d(a => 120, b => 21.5, w => 3).evalf() |> println

    12.0000000000000

c = 65, d = 12
左の縦，横は 65, 12 である。
このとき，
右の縦，横は (a - c - w) = 52, (d + w) = 15
前の縦，横は a = 120, (b - d - w) = 6.5

例題は無数に作れる。もう少しきれいな数値例は
縦 = 24，横 = 17，道幅 = 2 のとき
右 = 10 × 12，前 = 24 × 5，左 = 12 × 10 である。

function draw(a, b, w, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    t = sqrt(4a^2*b^2 - 4a^2*b*w + 4a^2*w^2 - 4a*w^3 + w^4)
    c = -a*b/w + a - w/2 + t/(2w)
    d = (2a*b - 4a*w + w^2 + t)/(2(3a - w))
    @printf("縦 = %g，横 = %g，道幅 = %g のとき\n右 = %g × %g，前 = %g × %g，左 = %g × %g である\n",
        a, b, w,
        a - c - w, d + w,
        a, b - d - w,
        c, d)
    plot(showaxis=false)
    rect(0, 0, a, b, :black)
    rect(0, 0, c, d, :aquamarine, fill=true)
    rect(0, d + w, a, b, :salmon, fill=true)
    rect(c + w, 0, a, d + w, :powderblue, fill=true)
    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, 0, "(0,0)", :black, :center, delta=-delta)
        point(a, b, "(a,b)", :black, :right, :bottom, delta=delta)
        point(c, 0, "(c,0)", :black, :right, delta=-delta)
        point(c + w, 0, "(c+w,0)", :black, :left, delta=-delta)
        point(0, d, "(0,d)", :black, :right, :vcenter)
        point(0, d + w, "(0,d+w)", :black, :right, :vcenter)
        point(c/2, d/2, "左", :black, :center, :vcenter, mark=false)
        point((a + c + w)/2, (d + w)/2, "右", :black, :center, :vcenter, mark=false)
        point(a/2, (b + d + w)/2, "前", :black, :center, :vcenter, mark=false)
        xlims!(-12delta, a + 2delta)
    end
end;

draw(120, 21.5, 3, true)

算額（その1665）

算額（その1665）

和国知恵較（わこくちえくらべ） 享保12年(1727)
和算の里みしま_パンフレット_HP.pdf
https://e-mishima.info/wp-content/uploads/2020/02/74a075f30cb806eef4a9a25c4a71608e.pdf
キーワード：魔方陣の一種，順列
#Julia #SymPy #算額 #和算 #数学

図の「い」から「ち」に 2 〜 9 の数字を当てはめ，2 つの円周に書かれている数を足しても，また，中心の 1 を除くどの直径に書かれている数を足しても，皆同じになるように数を並べなさい。

2 〜 9 の数字の並べ替え（順列）を行い，4 方向の和が等しくなるものを列挙する。

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

using Combinatorics

count = 0
a = 2:9
p = collect(1:16)
for i in 1:factorial(length(a))
    (い, ろ, は, に, ほ, へ, と, ち) = b = nthperm(a, i)
    x = い + ろ + は + に
    if x == ほ + へ + と + ち && x == い + ほ + と + は && x == ろ + へ + ち + に
        count += 1
        mod(count, 48) == 1 && println(b)
    end
end

768 通りもある。

function draw(x)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    pos = [0 2; 2 0; 0 -2; -2 0; 0 1; 1 0; 0 -1; -1 0]
    r = 0.4
    p = plot(showaxis=false)#, fontsize=24)
    delta = (fontheight = (ylims()[2]- ylims()[1]) / 500 * 10 * 2) /3  # size[2] * fontsize * 2
    circle(0, 0, 1, :gray80)
    circle(0, 0, 2, :gray80)
    segment(0, -2, 0, 2, :gray80)
    segment(-2, 0, 2, 0, :gray80)
    circle(0, 0, r)
    rotate(0, 2, r, angle=90)
    rotate(0, 1, r, angle=90)
    println(x)
    annotate!(0, 0, text("1", 12, :blue, :center, :vcenter))
    for i = 1:8
        annotate!(pos[i, 1], pos[i, 2], text(string(x[i]), 12, :blue, :center, :vcenter))
    end
    return p
end;

draw([2, 3, 8, 9, 5, 4, 7, 6])