裏 RjpWiki

Julia ときどき R, Python によるコンピュータプログラム,コンピュータ・サイエンス,統計学

算額(その1580)

2025年02月03日 | Julia

算額(その1580)

福島県田村郡三春町御木沢 厳島神社
三春まちなか寺子屋 30.12.22
https://miharukoma.com/wp-content/uploads/2018/12/寺子屋12.22 問題解説.pdf

キーワード:円4個,半円2個,正方形
#Julia, #SymPy, #算額, #和算, #数学

正方形の中に半円 2 個,乙円 2 個,丙円 2 個を容れる。半円の直径が 96 寸のとき,乙円,丙円の直径はいかほどか。

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

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

using SymPy
@syms r1::positive, r2::positive, r3::positive
eq1 = (r1 - r3)^2 + r1^2 - (r1 + r3)^2
eq2 = (r1 - 2r3 - r2)^2 + r1^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r2, r3))[1]

    (r1/12, r1/4)

乙円,丙円の直径は,甲円の直径の 1/12, 1/4 倍である。

甲円の直径が 96 寸のとき,乙円,丙円の直径は 8 寸,24 寸である。

using Colors
function draw(r1, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r2, r3) = (r1/12, r1/4)
    @printf("半円の直径が %g のとき,乙円,丙円の直径は %g, %g である。\n", 2r1, 2r2, 2r3)
    plot()
    circlef(0, r1, r1, parse(Colorant, "#FFF565"), beginangle=180, endangle=360)
    circlef(0, -r1, r1, parse(Colorant, "#FFF565"), beginangle=0, endangle=180)
    circle2f(r1 - r3, 0, r3, parse(Colorant, "#F02F22"))
    circle2f(r1 - 2r3 - r2, 0, r2, parse(Colorant, "#757CF7"))    
    circle(0, r1, r1, :black, beginangle=180, endangle=360)
    circle(0, -r1, r1, :black, beginangle=0, endangle=180)
    circle2(r1 - r3, 0, r3, :black)
    circle2(r1 - 2r3 - r2, 0, r2, :black)    
    rect(-r1, -r1, r1, r1, :black)
    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,r1)", :black, :center, :bottom, delta=delta/2)
        point(r1 - r3, 0, "丙円:r2\n(r1-r3,0)", :black, :center, delta=-delta/2)
        point(r1 - 2r3 - r2, 0, "乙円:r2,(r1-2r3-r2,0)", :black, :right, delta=-3delta)
    end
end;

draw(96/2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1579)

2025年02月02日 | Julia

算額(その1579)

福島県田村市 日渡神社 明治21年(1888)
https://www.city.tamura.lg.jp/soshiki/30/bunkazai.html
キーワード:円9個,外円
#Julia, #SymPy, #算額, #和算, #数学

外円の中に,大円 2 個と,小円 6 個を容れる。大円の直径が与えられたときに小円の直径を求める術を述べよ。

外円の半径と中心座標を R, (0, 0)
大円の半径と中心座標を r1, (r1, 0)
小円の半径と中心座標を r2, (0, y21), (r2, y22)
とおき,以下の連立方程式を解く。

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

using SymPy
@syms R::positive, r1::positive,
      r2::positive, y21::positive, y22::positive
r1 = R/2
eq1 = r1^2 + y21^2 - (r1 + r2)^2
eq2 = r2^2 + (y22 - y21)^2 - 4r2^2
eq3 = r2^2 + y22^2 - (R - r2)^2
res = solve([eq1, eq2, eq3], (r2, y21, y22))[1]

    (R*(-5 + sqrt(33))/4, R*(-5*sqrt(2) + sqrt(66))/(2*sqrt(7 - sqrt(33))), R*sqrt(7/2 - sqrt(33)/2))

小円の半径 r2 は,大円の半径 R の (√33 - 5)/4 倍である。
大円の直径が 3954 のとき,736.0001761028412 である。

3954*(√33 - 5)/4

    736.0001761028412

using Colors
function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = R/2
    (r2, y21, y22) = (R*(-5 + sqrt(33))/4, R*(-5*sqrt(2) + sqrt(66))/(2*sqrt(7 - sqrt(33))), R*sqrt(7/2 - sqrt(33)/2))
    @printf("外円の直径が %g のとき,小円の直径は %g である。\n", 2R, 2r2)
    plot()
    circlef(0, 0, R, RGB(194/255, 214/255, 236/255))
    circle2f(r1, 0, r1, RGB(234/255, 51/255, 35/255))
    circle22f(0, y21, r2, RGB(229/255, 74/255, 233/255))
    circle4f(r2, y22, r2,  RGB(229/255, 74/255, 233/255))
    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", :green, :center, :bottom, delta=delta/2)
        point(r1, 0, "大円:r1,(r1,0)", :white, :center, delta=-delta/2)
        point(0, y21, "小円:r2,(0,y21)", :black, :left, :bottom, delta=delta)
        point(r2, y22, " 小円:r2,(r2,y22)", :black, :left, :bottom, delta=delta)
    end
end;

draw(1977, true)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1578)

2025年02月01日 | Julia

算額(その1578)

『算法續淺問答 巻之七』
terasima3.pdf
http://www.wasan.jp/terasima/terasima.html
http://www.wasan.jp/terasima/terasima3.pdf

キーワード:円環,外円
#Julia, #SymPy, #算額, #和算

図のように 8 個の等円が互いに外接し,外円に内接している。外円の直径が 10 寸のとき,等円の直径を求めよ。

外円の半径と中心座標を R, (0, 0)
等円の半径と中心座標を r, (R - r, 0)
等円の個数を n
とおき,以下の方程式を解く。

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

using SymPy
@syms R::positive, r::positive, n::positive
eq = (R - r)*sind(180/n) - r
res = solve(eq, r)[1]
res |> println

    R*sin(pi/n)/(sin(pi/n) + 1)

n = 8, R = 10/2 を代入し,得られた半径を 2 倍すれば直径は 2.76768653914155 である。

2*res(n => 8, R => Sym(10)/2) |> println
2*res(n => 8, R => Sym(10)/2).evalf() |> println

    10*sqrt(1/2 - sqrt(2)/4)/(sqrt(1/2 - sqrt(2)/4) + 1)
    2.76768653914155

function draw(R, n, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = R*sin(pi/n)/(sin(pi/n) + 1)
    p = plot()
    circle(0, 0, R, :green)
    circle(0, 0, R - r, :gray80)
    rotate(0, R - r, r, angle=360/n)
    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", :green, :center, :bottom, delta=delta/2)
        point(0, R - r, "R-r", :gray, :center, :bottom, delta=delta/2)
    end
    return p
end;

draw(10, 8, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1577)

2025年02月01日 | Julia

算額(その1577)

愛知県名古屋市中央区 大須観音(北野山真福寺寶生院) 天保8年(1837)
藤安順:所掲大須観音堂解義
山田潤:動的幾何ソフトウエアを利用した平面図形の作図についての一考察 - 和算所にある平面図形問題の利用 -
2023018-Yamada.pdf

https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2273-19.pdf

キーワード:円10個,外円
#Julia, #SymPy, #算額, #和算

外円の中に大円 4 個,甲円 3 個,乙円を 2 個容れる。甲円の直径が 1 寸のとき,乙円および外円の直径はいかほどか。

算額(その1575)は,深川が意図的に大円が交差する部分にある円を省略し,甲円と乙円の関係を述べさせるものであったようだ。

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

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

using SymPy
@syms R::positive, r1::positive, x1::positive,
      r2::positive, r3::positive
x1 = r1 - r2
eq1 = R - (2r1 + r2)
eq2 = x1^2 + r1^2 - (r1 + r2)^2  # 大円と甲円が接する
eq3 = x1^2 + (R - r3 - r1)^2 - (r1 + r3)^2;  # 大円と乙円が接する

res = solve([eq1, eq2, eq3], (R, r1, r3))[1]

    (9*r2, 4*r2, r2)

外円の直径は甲円の直径の 9 倍,大円の直径は甲円の直径の 4 倍,乙円の直径は甲円の直径に等しい。

考察:大円の交差するところに甲円が入っているのが,算額(その1575)での y1 の制約条件になっている。

function draw(r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (R, r1, r3) = (9*r2, 4*r2, r2)
    plot()
    circle(0, 0, R, :green)
    circle4(r1 - r2, r1, r1, :blue)
    circle(0, 0, r2)
    circle22(0, r1, r3)
    circle22(0, R - r3, r3, :magenta)
    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", :green, :center, :bottom, delta=delta/2)
        point(r1 - r2, r1, "大円:r1,(r1-r2,r1)", :blue, :center, delta=-delta, deltax=5delta)
        point(0, R - r3, "乙円:r3,(0,R-r3)", :magenta, :left, :bottom, delta=delta/2, deltax=4delta)
        point(0, r1, "甲円:r2,(0,r1)", :red, :right, :vcenter, deltax=-4delta)
    end
end;

draw(1/2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1576)

2025年02月01日 | Julia

算額(その1576)

佐久間纉:算法起源集,1877.
https://kokusho.nijl.ac.jp/biblio/100234582/69?ln=ja

山田潤:動的幾何ソフトウエアを利用した平面図形の作図についての一考察 - 和算所にある平面図形問題の利用 -
2023018-Yamada.pdf

https://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/2273-19.pdf

キーワード:円7個,外円
#Julia, #SymPy, #算額, #和算

外円の中に甲円,乙円,丙円を 2 個ずつ容れる。甲円,乙円,丙円の直径が与えられたとき,外円の直径はいかほどか。

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

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

using SymPy
@syms R::positive, r1::positive,
      r2::positive, x2::positive, y2::negative,
      r3::positive, x3::positive, y3::positive
eq1 = x2^2 + y2^2 - (R - r2)^2 |> expand
eq2 = x3^2 + y3^2 - (R - r3)^2 |> expand
eq3 = x2^2 + (y2 - r1 + R)^2 - (r1 + r2)^2 |> expand
eq4 = x3^2 + (R - r1 - y3)^2 - (r1 + r3)^2 |> expand
eq5 = (x3 - x2)^2 + (y3 - y2)^2 - (r2 + r3)^2 |> expand;

まず,eq1, eq3 の連立方程式から x2, y2 を求める。

res2 = solve([eq1, eq3], (x2, y2))[1]

    (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2)/(-R + r1), -(-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1))

次に,eq2, eq4 の連立方程式から x3, y3 を求める。

res3 = solve([eq2, eq4], (x3, y3))[1]

    (-2*sqrt(R)*sqrt(r1)*sqrt(r3)*sqrt(R - r1 - r3)/(-R + r1), (-R^2 + R*r1 + R*r3 + r1*r3)/(-R + r1))

eq5 の x2, y2, x3, y3 に求められた解を代入し,整理する。

eq15 = eq5(x2 => res2[1], y2 => res2[2], x3 => res3[1], y3 => res3[2])/R |> simplify |> numerator;

しかし,このままでは SymPy では解けない。

eq15 |> println

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

8*r1*sqrt(r2)*sqrt(r3)*sqrt(R - r1 - r2)*sqrt(R - r1 - r3) の項とそれ以外の項をそれぞれ二乗し,方程式を再構成する。

eq25 = (4*R^3 - 8*R^2*r1 - 4*R^2*r2 - 4*R^2*r3 + 4*R*r1^2 + 4*R*r1*r2 + 4*R*r1*r3 + 8*r1*r2*r3)^2 - ( - 8*r1*sqrt(r2)*sqrt(r3)*sqrt(R - r1 - r2)*sqrt(R - r1 - r3))^2;

ここまで来ると SymPy でも解ける。

res = solve(eq25)[2];  # 2 of 5
res[R] |> println

    r1 + r2 + r3

外円の半径は,甲円,乙円,丙円の半径の合計である。

なお,x2 = x3 であることがわかる。

# x2
res2[1](R => r1 + r2 + r3) |> println
# x3
res3[1](R => r1 + r2 + r3) |> println

    -2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(-r2 - r3)
    -2*sqrt(r1)*sqrt(r2)*sqrt(r3)*sqrt(r1 + r2 + r3)/(-r2 - r3)

function draw(r1, r2, r3, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    R = r1 + r2 + r3
    (x2, y2) = (-2*sqrt(R)*sqrt(r1)*sqrt(r2)*sqrt(R - r1 - r2)/(-R + r1), -(-R^2 + R*r1 + R*r2 + r1*r2)/(-R + r1))
    (x3, y3) = (-2*sqrt(R)*sqrt(r1)*sqrt(r3)*sqrt(R - r1 - r3)/(-R + r1), (-R^2 + R*r1 + R*r3 + r1*r3)/(-R + r1))
    @printf("r1 = %g;  r2 = %g;  r3 = %g;  R = %g;  x2 = %g;  y2 = %g;  x3 = %g;  y3 = %g\n",
        r1, r2, r3, R, x2, y2, x3, y3)
    plot()
    circle(0, 0, R, :green)
    rotate(0, R - r1, r1, angle=180)
    rotate(x2, y2, r2, :blue, angle=180)
    rotate(x3, y3, r3, :magenta, angle=180)
    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", :green, :center, :bottom, delta=delta/2)
        point(0, R - r1, "甲円:r1,(0,R-r1)", :red, :center, delta=-delta/2)
        point(x2, y2, "乙円:r2,(x2,y2)", :blue, :center, delta=-delta/2)
        point(x3, y3, "丙円:r3\n(x3,y3)", :magenta, :center, delta=-delta/2)
    end
end;

draw(0.5, 0.3, 0.2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1575)

2025年02月01日 | Julia

算額(その1575)

愛知県名古屋市中央区 大須観音(北野山真福寺寶生院) 天保8年(1837)
深川英俊,トニー・ロスマン:聖なる数学:算額,森北出版株式会社,2010年4月22日.

キーワード:円8個,外円
#Julia, #SymPy, #算額, #和算

外円の中に,大円 4 個,小円 3 個を容れる。外円と小円の直径が与えられたとき,大円の直径はいかほどか。

深川は,「中央の小円1と左右の小円2の直径の関係を示せ」という問題としている。
図を見ただけで直径は等しいとわかるが,ちゃんと式で示せということではある。
大円の中心の y 座標はどんな非負の値も取れるが,3個の小円の直径は同じであるということで,算額問題としては,「大円の直径を外円と小円の直径であらわせ」という問題にした。

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

using SymPy
@syms R::positive, r1::positive, y1::positive, r2::positive, r3::positive
eq1 = r1^2 + y1^2 - (R - r1)^2
eq2 = r1^2 + y1^2 - (r1 + r2)^2
res = solve([eq1, eq2], (r1, y1))[1]

    (-(-R + r2)/2, sqrt(R)*sqrt(r2))

大円の直径は,外円の直径から小円の直径を差し引いたものの半分である。

外円の直径が同じでも,小円の直径により印象の異なる図が得られる。

function draw(R, r2, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    (r1, y1) = ((R - r2)/2, sqrt(R)*sqrt(r2))
    r3 = r2
    p = plot()
    circle(0, 0, R, :green)
    circle4(r1, y1, r1, :magenta)
    circle2(R - r2, 0, r2, :blue)
    circle(0, 0, r3)
    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(R, 0, "R", :green, :left, :bottom, delta=delta/2, deltax=delta/2)
        point(R - r2, 0, "R-r2", :blue, :center, delta=-delta)
        point(r3, 0, "r3", :red, :right, :vcenter, deltax=-delta/2)
        point(r1, y1, "大円:r1,(r1,y1)", :magenta, :center, :bottom, delta=delta/2)
        pos = max(R, r2)
        str = @sprintf(" R=%.3g, r2=%.3g", R, r2)
        point(0, pos, str, :black, :center, :bottom, delta=delta, mark=false)
    end
    return p
end;

draw(1, 0.1, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1574)

2025年01月31日 | Julia

算額(その1574)

三重県鳥羽市 観世音 寛政12年(1800)
深川英俊,トニー・ロスマン:聖なる数学:算額,森北出版株式会社,2010年4月22日.

キーワード:円1個,二等辺三角形2個
#Julia, #SymPy, #算額, #和算

外円の中に大小の正三角形を容れる。外円の直径が与えられたとき,小正三角形の一辺の長さはいかほどか。

外円の半径と中心座標を R, (0, 0)
小正三角形の一辺の長さを q
とおき,以下の方程式を解く。

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

using SymPy
@syms R, q
eq = (q/2)^2 + (-R/2 - √Sym(3)*q/2)^2 - R^2
res = solve(eq, q)[1]
res |> println

    R*(-sqrt(3) + sqrt(15))/4

小正三角形の一辺の長さ q は,外円の半径 R の (√15 - √3)/4 倍である。
外円の直径が 100 のとき,小正三角形の一辺の長さは 26.76165673298175 である。

100/2*(√15 - √3)/4

    26.76165673298175

function draw(R, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    q = R*(√15 - √3)/4
    @printf("外円の直径が = %g のとき,小正三角形の一辺の長さは %g である。\n", 2R, q)
    plot()
    circle(0, 0, R)
    polygon(0, 0, R, 3, color=:blue)
    polygon(0, -R/2 - q/√3, q/√3, 3, color=: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, R, "R", :red, :center, :bottom, delta=delta/2)
        point(0, -R/2, "-R/2", :blue, :center, :bottom, delta=delta/2)
        point(√3R/2, -R/2, "(√3R/2,-R/2)  ", :blue, :right, :bottom, delta=delta/2)
        point(q/2, -R/2 - √3q/2, "(q/2,-R/2-√3q/2)", :blue, :left, delta=-delta/2)
    end  
end;

draw(100/2, true)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

Julia の歪度と尖度

2025年01月31日 | Julia

1. Julia の歪度と尖度は古典的定義に基づく

1.1. skewness, kurtosis は StatsBase に含まれる。

# Skewness
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    skewness(v, [wv::AbstractWeights], m=mean(v))

Compute the standardized skewness of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

# (excessive) Kurtosis
# This is Type 1 definition according to Joanes and Gill (1998)
"""
    kurtosis(v, [wv::AbstractWeights], m=mean(v))

Compute the excess kurtosis of a real-valued array `v`, optionally
specifying a weighting vector `wv` and a center `m`.
"""

1.2. R は 3 通りの定義に基づく関数値を返す

Joanes and Gill (1998) [^1]discuss three methods for estimating kurtosis and skewness.

[1]: D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189.

1.2.1. Kurtosis

Type 1: This is the typical definition used in many older textbooks.

g2 = m4/m2^2 - 3

Type 2: Used in SAS and SPSS.

G2 = ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))

Type 3: Used in MINITAB and BMDP.

b2 = m4/s^4 − 3 = (g2 + 3)*(1 − 1/n)^2 − 3

Only G2 (corresponding to type = 2) is unbiased under normality.

デフォルトは type=3 である。

1.2.2. Skewness

Type 1: This is the typical definition used in many older textbooks.

g1 = m3/m2^(3/2)

Type 2: Used in SAS and SPSS.

G1 = g1*sqrt(n*(n - 1))/(n - 2)

Type 3: Used in MINITAB and BMDP.

b1 = m3/s^3 = g1*((n - 1)/n)^(3/2)

All three skewness measures are unbiased under normality.

デフォルトは type=3 である。

1.3. Julia でも type を選べる関数を定義する

1.3.1. 尖度

using StatsBase

kurt1(x) = kurt(x, type=1)
kurt2(x) = kurt(x, type=2)
kurt3(x) = kurt(x)
function kurt(x; type=3)
    doc = """
    kurt(x; type=n)
    type: an integer between 1 and 3 selecting one of the algorithms for computing kurtosis detailed below.

    Type 1: g2. This is the typical definition used in many older textbooks.
    Type 2: G2. Used in SAS and SPSS.
    Type 3: b2. default. Used in MINITAB and BMDP.
    """
    type in 1:3 || error(doc)
    n = length(x)
    g2 = kurtosis(x)
    if type == 2
        return  ((n + 1)*g2 + 6)*(n − 1)/((n − 2)*(n − 3))
    elseif type== 3
        return (g2 + 3)*(1 − 1/n)^2 − 3
    else
        return g2
    end
end;           

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1];

kurt(x), kurt3(x)

    (-1.6335639958376689, -1.6335639958376689)

kurt(x, type=1), kurt1(x)

    (-1.3130419701699616, -1.3130419701699616)

kurt(x, type=2), kurt2(x)

    (-1.356984911550468, -1.356984911550468)

kurt(x, type=3), kurt3(x)

    (-1.6335639958376689, -1.6335639958376689)

1.3.2. 歪度

using StatsBase

skew1(x) = skew(x, type=1)
skew2(x) = skew(x, type=2)
skew3(x) = skew(x)
function skew(x; type=3)
    doc = """
    skew(x; type=n)
    type: an integer between 1 and 3 selecting one of the algorithms for computing skewness detailed below.

    Type 1: g1. This is the typical definition used in many older textbooks.
    Type 2: G1. Used in SAS and SPSS.
    Type 3: b1. default. Used in MINITAB and BMDP.
    """
    type in 1:3 || error(doc)
    n = length(x)
    g1 = skewness(x)
    if type == 2
        return g1*sqrt(n*(n - 1))/(n - 2)
    elseif type== 3
        return g1*((n - 1)/n)^(3/2)
    else
        return g1
    end
end;           

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 1]
skew(x), skew3(x)

    (0.10905343709973031, 0.10905343709973031)

skew(x, type=1), skew1(x)

    (0.12772490663150174, 0.12772490663150174)

skew(x, type=2), skew2(x)

    (0.15146310708295876, 0.15146310708295876)

skew(x, type=3), skew3(x)

    (0.10905343709973031, 0.10905343709973031)

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1572)

2025年01月31日 | Julia

算額(その1572)

四十 岩手県一関市牧沢 牧沢八幡神社 明治5年(1872)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

一関市博物館 第23回 令和6年度 和算に挑戦!
https://www.city.ichinoseki.iwate.jp/index.cfm/7,177295,c,html/177295/20241125-090413.pdf

キーワード:円,長方形,折り紙
#Julia, #SymPy, #算額, #和算

長辺,短辺が 4,3 の長方形の紙を折ってできる直角三角形に内接する円の直径はいかほどか。

山村の図は相変わらず不明瞭・不正確で(そもそも何を描いているのかさえわからない),算額の穴埋めも失敗している。「一関博物館」は流石である。

長辺,短辺の長さを a, b
折り線と長辺の交点座標を (c, 0)
元の頂点が対角線上に重なる座標を (x0, y0)
円の半径と中心座標を r, (x, y)
とおき,以下の連立方程式を解く。

円の半径だけを求めるならば eq1, eq2 だけの連立方程式を解けばよい。

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

using SymPy
@syms a, b, c, r, diag, x0, y0, x, y
eq1 = ((sqrt(a^2 + b^2) - b)^2 + c^2) - (a - c)^2
eq2 = c + sqrt(a^2 + b^2) - b - (a - c) - 2r
eq3 = y0/(a - x0) - b/a
eq4 = (x0 - c)^2 + y0^2 - c^2;

res = solve([eq1, eq2, eq3, eq4], (c, r, x0, y0))[1]

    (b*(-b + sqrt(a^2 + b^2))/a, (a*(-a - b + sqrt(a^2 + b^2))/2 - b^2 + b*sqrt(a^2 + b^2))/a, a*b/sqrt(a^2 + b^2), -b^2/sqrt(a^2 + b^2) + b)

対角線の長さを d = sqrt(a^2 + b^2) とすれば,
折り線と長方形の頂点の交点の x 座標は c = b*(d - b)/a
円の半径は (d - a - b)/2 + c である。

以下は図を描くために円の中心座標を求める連立方程式であるが,a, b を記号のままで解くと長い式になるので,eq5, eq6 の a, b に特定の値を与えて解くほうが現実的である。

eq5 = dist2(a, 0, 0, b, x, y, r)(r => res[2], x0 => res[3], y0 => res[4])(a => 4, b => 3)
eq6 = dist2(c, 0, x0, y0, x, y, r)(c => res[1], r => res[2], x0 => res[3], y0 => res[4])(a => 4, b => 3);

solve([eq5, eq6], (x, y))[3]  # 3 of 4

    (5/2, 1/2)

function draw(a, b, more)
    pyplot(size=(500, 500), grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    d = sqrt(a^2 + b^2)
    c = b*(d - b)/a
    r = (d - a - b)/2 + c
    (x0, y0) = (a*b/d, -b^2/d + b)
    (x, y) = (5/2, 1/2)
    @printf("a = %g;  b = %g;  d = %g;  c = %g;  r = %g\n", a, b, d, c, r)
    plot([0, a, a, 0, 0], [0, 0, b, b, 0], color=:green, lw=0.5, seriestype=:shape, fillcolor=:lemonchiffon)
    plot!([c, a, x0, c], [0, 0, y0, 0], color=:gray, seriestype=:shape, fillcolor=:orange)
    plot!(a .- [c, a, x0, c], b .- [0, 0, y0, 0], color=:gray, seriestype=:shape, fillcolor=:orange)
    plot!([0, c, x0, 0], [b, 0, y0, b], color=:red, seriestype=:shape, fillcolor=:white)
    plot!(a .- [0, c, x0, 0], b .- [b, 0, y0, b], color=:red, seriestype=:shape, fillcolor=:white)
    circlef(x, y, r)
    circlef(a - x, b - y, 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(0, b, " b", :blue, :left, :bottom, delta=delta/2) 
        point(a, b, "(a,b)", :blue, :right, :bottom, delta=delta/2) 
        point(a - c, b, "(a-c,b)", :blue, :center, :bottom, delta=delta/2) 
        point(c, 0, "c ", :blue, :right, :bottom, delta=delta/2) 
        point(x0, y0, "(x0,y0)", :blue, :left, :bottom, delta=delta/2)
        point(x, y, "(x,y)", :blue, :left, :bottom, delta=delta/2)
    end  
end;

draw(4, 3, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1571)

2025年01月30日 | Julia

算額(その1571)

福島県田村市 安倍文殊菩薩堂 明治10年(1877)
五輪教一:和算で遊ぼう!! 「三春まちなか寺子屋」2017 レポート

https://miharukoma.com/wp-content/uploads/2018/01/三春まちなか寺子屋2017レポート.pdf
キーワード:3次元図形,体積
#Julia, #SymPy, #算額, #和算

辺の長さが等しい正方形 6 面,正三角形 8 面からなる立体「截篭(キリコ)」がある。辺の長さが 1 寸のとき,この立体の体積はいかほどか。

截篭は立方体の頂点を含む 8 個の三角錐(O・ABC など)を切り取ってできる立体である。A,B,C は各辺の中点である。

AB = a とおく。
OA = b = a/√2
立方体の一辺の長さは 2*b = √2*a
立方体の体積は V1 = 2√2*a^3

三角錐の体積は,底辺 ABC の面積×高さ/3 であるが,そのように考えると「高さ」の計算が面倒。
底面OAB,高さOC と考えて,三角錐の体積の公式を使う。
底面積OAB = b^2/2,高さ = b
よって,三角錐の体積 V2 = (b^3/2)/3 = √2a^3/24

截篭の体積 = V1 - 8V2 = 5√2a^3/3

一辺の長さが 1 寸のとき,截篭の体積 =  5√2/3 = 2.35702260395516 である。

using SymPy
@syms a, b
b = a/√Sym(2)
V1 = (2b)^3
V1 |> println

    2*sqrt(2)*a^3

V2 = (b^3/2)/3
V2 |> println

    sqrt(2)*a^3/24

V = V1 - 8V2
V |> println

    5*sqrt(2)*a^3/3

V(a => 1).evalf() |> println

    2.35702260395516

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1570)

2025年01月29日 | Julia

算額(その1570)

東京都八王子市片倉町(片倉城城跡公園内) 住吉神社 嘉永4年(1851) (昭和62年復元)
山口正義:あきる野市有形民族文化財 二宮神社の算額

https://yamabukiwasan.sakura.ne.jp/ninomiyasasshi2.pdf
https://yamabukiwasan.sakura.ne.jp/page2.html
キーワード:円2個,面積
#Julia, #SymPy, #算額, #和算

互いに外接する 2 個の等円と共通接線によって囲まれる面積(空積)を与えたとき,等円の直径を求めよ

空積は,等円の半径を一辺とする正方形の面積から等円の面積の 1/4 を差し引いたものを 2 倍したものである。

等円の半径と中心座標を r, (r, 0), (-r, 0)
とおき,以下の方程式を解く。

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

using SymPy
@syms r::positive, S::positive, 円周率::positive, 円積率::positive

eq = 2(r^2 - 円周率*r^2/4) - S;

ans_r = solve(eq, r)[2]  # 2 of 2
ans_r |> println

    sqrt(2)*sqrt(S)*sqrt(-1/(円周率 - 4))

等円の半径は sqrt(2S/(4 - π)) である。

術では円積率 0.79 を使っているので,それに合わせるように円周率として 3.16 を使い,空積が 1 のときは以下のようになる。

2ans_r(S => 1, 円周率 => 3.16).evalf() |> println

    3.08606699924184

術は以下のようになっており,同じ答えを与える。

空積 = 1
円積率 = 0.79
天 = 1 - 円積率
sqrt(天*2*空積)/天

    3.0860669992418384

function draw(S, more=false)  # b:楕円の短半径と,r:円の半径
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r = sqrt(2S/(4 - π))
    plot()
    circle2(r, r, r)
    rect(-r, 0, r, r, :gray80)
    θ = 270:0.1:360
    x = r*cosd.(θ) .- r
    y = r*sind.(θ) .+ r
    θ = 180:0.1:270
    append!(x, r*cosd.(θ) .+ r)
    append!(y, r*sind.(θ) .+ r)
    append!(x, -r)
    append!(y, 0)
    plot!(x, y, color=:black, lw=0.5, seriestype=:shape, fillcolor= :gray90)
    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(r, r, " 等円:r,(r,r)", :red, :left, :vcenter)
        point(0, r/4, "空積:S", :black, :center, :vcenter, mark=false)
    end
end;

draw(1, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1569)

2025年01月29日 | Julia

算額(その1569)

百 大船渡市猪川町 雨宝堂 現雨宝山竜宝院 文政10年(1827)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:円14個
#Julia, #SymPy, #算額, #和算

内部にそれぞれ 3 個の中円が入った 2 個の大円が,互いに接し合う小円 6 と接している。小円の直径が 1 寸のとき,中円の直径はいかほどか。

大円の半径と中心座標を r1, (x1, 0)
中円の半径を r2
小円の半径と中心座標を r3, (0, r3), (r3, y3)
とおき,以下の連立方程式を解く。

まず,大円と小円について以下の連立方程式を解く。大円の中の中円については後で。

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

using SymPy
@syms r1::positive, x1::positive,
      r3::positive, y3::positive

eq1 = x1^2 + r3^2 - (r1 + r3)^2
eq2 = (x1 - r3)^2 + y3^2 - (r1 + r3)^2
eq3 = r3^2 + (y3 - r3)^2 - 4r3^2
res = solve([eq1, eq2, eq3], (r1, x1, y3))[1];

# r1
res[1] |> sympy.sqrtdenest |> simplify |> println

    r3*(-1 + sqrt(2) + sqrt(6))

大円の半径 r1 は,小円の半径 r3 の (√2 + √6 - 1) 倍である。

# x1
res[2] |> println

    r3*(sqrt(3) + 2)

# y3
res[3] |> println

    r3*(1 + sqrt(3))

大円とその中にある 3 個の中円の関係については,算額(その197) に示した。
中円の半径 r2 は,大円の半径 r1 の 3/(3 + 2√3) 倍である。
大円の半径は上述のように r3*(√2 + √6 - 1) なので,r2 = r3*(3 + 3√2 - 2√3 - √6) である。

r2 = 3res[1]/(3 + 2√Sym(3)) |> simplify |> sympy.sqrtdenest |> simplify
r2 |> println

    r3*(-2*sqrt(3) - sqrt(6) + 3 + 3*sqrt(2))

小円の半径が 1/2 のとき,中円の半径は 0.664524664599176,直径は 1.329049329198352 である。

r2(r3 => 1/2).evalf() |> println

    0.664524664599176

2*0.664524664599176

    1.329049329198352

術は以下のようになっている。

小円径 = 1
天 = √12 - 3
A = sqrt((天 + 7)*2) - 1
中円径 = A*天*小円径

    1.3290493291983518

小円径に掛ける倍数を 1 つの式にまとめて簡約化すれば,上で得たものと同じであることがわかる。

倍数 = (sqrt((√Sym(12) - 3 + 7)*2) - 1)*(√Sym(12) - 3)
倍数 |> expand |> simplify |> sympy.sqrtdenest |> println

    -2*sqrt(3) - sqrt(6) + 3 + 3*sqrt(2)

function draw(r3, more=false)  # b:楕円の短半径と,r:円の半径
    pyplot(size=(500, 500), showaxis=true, grid=false, aspectratio=1, label="", fontfamily="IPAMincho")
    r1 = r3*(√2 + √6 - 1)
    x1 = r3*(√3 + 2)
    y3 = r3*(1 + √3)
    r2 = r3*(3 + 3√2 - 2√3 - √6)
    plot()
    circle2(x1, 0, r1)
    circle22(0, r3, r3, :blue)
    circle4(r3, y3, r3, :blue)
    circle2(x1 - r1 + r2, 0, r2, :green)
    circle4(x1 + (r1 - r2)/2, (r1 - r2)√3/2, r2, :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(x1, 0, "大円:r1,(x1,0)", :red, :center, delta=-delta)
        point(x1 - r1 + r2, 0, "中円:r2\n(x1-r1+r2,0)", :green, :center, :bottom, delta=delta)
        point(0, r3, "小円:r3,(0,r3)", :blue, :right, delta=-3delta, deltax=5delta)
        point(r3, y3, "小円:r3,(r3,y3) ", :blue, :right, :vcenter)
    end
end;

draw(1/2, true)

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1568)

2025年01月29日 | Julia

算額(その1568)

九十七 岩手県大船渡市猪川町 雨宝堂(現雨宝山竜宝院) 文政7年(1824)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html
キーワード:球4個,3次元
#Julia, #SymPy, #算額, #和算

直方体の中に,甲球 1 個,乙球 1 個,丙球 2 個を容れる。甲球の直径は丙球の直径の 2 倍である。乙球の直径が 4 寸のとき,甲球の直径はいかほどか。

直方体の,原点と対角の頂点座標を (a, 2r1, 2r1)
甲球の半径と中心座標を r1, (r1, r1, r1)
乙球の半径と中心座標を r2, (x2, r1, 2r1 - r2)
丙球の半径と中心座標を r3, (a - r3, r3, r3), (a - r3, 3r3, r3)
とおき,以下の連立方程式を解く。

using SymPy
@syms a::positive, r1::positive, r3::positive, r2::positive, x3::positive, z2::positive
r3 = r1/2
eq1 = (a - r3 - r1)^2 + (r3 - r1)^2 + (r3 - r1)^2 - (r1 + r3)^2
eq2 = (x3 - r1)^2 + (r1 - r2)^2 - (r1 + r2)^2
eq3 = (a - r3 - x3)^2 + (r1 - r3)^2 + (2r1 - r2 - r3)^2 - (r3 + r2)^2;

res = solve([eq1, eq2, eq3], (r1, a, x3))[1];

# r1
res[1] |> println

    7*r2/4

甲球の半径 r1 は,乙球の半径 r2 の 7/4 倍である。
乙球の直径が 4 寸のとき,甲球の直径は 7 寸である。

# a
res[2] |> simplify |> println

    7*r2*(sqrt(7) + 3)/8

# x3
res[3] |> println

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

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1567)

2025年01月28日 | Julia

算額(その1567)

九十五 大船渡市猪川町長谷堂 気仙長谷寺 文政5年(1822)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

キーワード:球5個,3次元
#Julia, #SymPy, #算額, #和算

盤上に大球 1 個,小球 4 個を載せる。小球は互いに外接しあい,大球は小球 3 個と外接している。小球の直径が 1 寸のとき,大球の直径はいかほどか。

大球の半径と中心座標を r1, (0, 0, r1)
小球の半径と中心座標を r2, (x22, r2, r2), (x21, 0, r2), (x22 + r2/√3, 0, z2)
とおき,以下の連立方程式を解く。

using SymPy
@syms r1::positive, r2::positive, x21::positive, x22::positive, z2::positive
s3 = √Sym(3)
eq1 = x22^2 + r2^2 + (r1 - r2)^2 - (r1 + r2)^2  # A:B
eq2 = (x22 + r2/s3)^2 + (z2 - r1)^2 - (r1 + r2)^2  # A:D
eq3 = r2^2/3 + r2^2 + (z2 - r2)^2 - 4r2^2  # B:D
eq4 = (x22 - x21)^2 + r2^2 - 4r2^2  # BC
res = solve([eq1, eq2, eq3, eq4], (r1, x21, x22, z2))[2];  # 2 of 3

# r1
res[1] |> sympy.sqrtdenest |> simplify |> println

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

大球の半径 r1 は,小球の半径 r2 の (√6 + 3)/2 倍である。
小球の直径が 1 寸のとき,大球の直径は (√6 + 3)/2 = 2.724744871391589 寸である。

# x21
res[2] |> sympy.sqrtdenest |> simplify |> println

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

# x22
res[3] |> sympy.sqrtdenest |> simplify |> println

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

# z2
res[4] |> sympy.sqrtdenest |> simplify |> println

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

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

算額(その1566)

2025年01月28日 | Julia

算額(その1566)

八十六 岩手県室根村室根山 室根神社 明治32年(1899)
山村善夫:現存 岩手の算額,昭和52年1月30日,熊谷印刷,盛岡市.

http://www.wasan.jp/yamamura/yamamura.html

キーワード:整数方程式
#Julia, #SymPy, #算額, #和算

第 1 問

今如図三宝ニ大球ヲ載ルアリ各数ヲ不知只云大球(一十一段ト小球五段ト)和シテ一百八拾三個 又云大球ト小球数ヲ相併テ七除シ奇令ナシ大球数幾何但シ大球数小球数共ニ一位ニ止無奇令

注:「段」は段に積むという意味ではなく,「倍」ということ。

意訳:大球と小球がある。大球の個数の 11 倍と小球の個数の 5 倍を足すと 183,大球と小球をあわせた個数は 7 で割り切れる。それぞれ何個ずつか。

以下のブルートフォースで,大球 13 個,小球 8 個である。
検算:
13*11 + 8*5 = 183
(13 + 8)/7 = 3

for 大 = 1:16, 小 = 1:36
    11大 + 5小 == 183 && rem(大 + 小, 7) == 0 && (println("大 = $大, 小 = $小"); break)
end

    大 = 13, 小 = 8

第 3 問

今如図甲乙ノ数球ヲ高円台ニ重載スルアリ只云甲球一個量目五拾三匁又云乙球一個量目四拾八匁別云乙総球ノ目弱ナルコト一百九拾四匁甲乙球数幾何 但シ甲乙ノ球数壱位止無奇令

「乙総球ノ目弱ナルコト一百九拾四匁」を山村は「総量 - 乙量 = 194」と読んでいるが,答えからいうと「乙量 - 甲量 = 194」であろう。

意訳:甲球,乙球の 1 個あたりの重さは 53, 48 である。それぞれいくつあるかはわからないが,乙球の総重量から甲球の総重量を引くと 194 である。甲球,乙球はそれぞれいくつか。

ブルートフォースで以下のようにして,甲球数 = 38, 乙球数 = 46 である。

検算:
甲球の総重量 = 53*38 = 2014
乙球の総重量 = 48*46 = 2208
乙球の総重量 - 甲球の総重量 = 2208 - 2014 = 194

甲球重量 = 53
乙球重量 = 48
for 甲球数 = 1:100, 乙球数 = 1:100
    (乙球数*乙球重量) - (甲球数*甲球重量) == 194 && (println("甲球数 = $甲球数, 乙球数 = $乙球数"); break)
end

    甲球数 = 38, 乙球数 = 46

第 5 問

今如図円台ニ甲球若干数ト一大球ニ若干乙球ヲ容テ載ルアリ甲球ト乙球ハ等数ナリ只云甲球ノ総径ヲ乙球一ケ径ヲ以テ除之二拾四個ヲ得又云乙球ノ総径ヲ甲球一ケノ径ヲ以除之六ケヲ得等数幾何

意訳:甲球,乙球が同数個ある。甲球の直径の和を乙球の直径で割ると 24 になり,乙球の直径の和を甲球の直径で割ると 6 になる。何個ずつあるのか。

以下の方程式を解いて,甲球,乙球は 12 個ずつある。甲球の直径は乙球の直径の 2 倍である。

using SymPy

@syms 数::positive, 甲球径, 乙球径
eq1 = 甲球径*数//乙球径 - 24
eq2 = 乙球径*数//甲球径 - 6
solve([eq1, eq2], (数, 甲球径))[1]

    (12, 2*乙球径)

ブルートフォースで解くと以下のような結果になる。

for 数 = 1:100, 甲球径 = 1:30, 乙球径 = 1:甲球径 - 1
    甲球径*数//乙球径 == 24 && 乙球径*数//甲球径 == 6 && println("数 = $数, 甲球径 = $甲球径, 乙球径 = $乙球径")
end

    数 = 12, 甲球径 = 2, 乙球径 = 1
    数 = 12, 甲球径 = 4, 乙球径 = 2
    数 = 12, 甲球径 = 6, 乙球径 = 3
    数 = 12, 甲球径 = 8, 乙球径 = 4
    数 = 12, 甲球径 = 10, 乙球径 = 5
    数 = 12, 甲球径 = 12, 乙球径 = 6
    数 = 12, 甲球径 = 14, 乙球径 = 7
    数 = 12, 甲球径 = 16, 乙球径 = 8
    数 = 12, 甲球径 = 18, 乙球径 = 9
    数 = 12, 甲球径 = 20, 乙球径 = 10
    数 = 12, 甲球径 = 22, 乙球径 = 11
    数 = 12, 甲球径 = 24, 乙球径 = 12
    数 = 12, 甲球径 = 26, 乙球径 = 13
    数 = 12, 甲球径 = 28, 乙球径 = 14
    数 = 12, 甲球径 = 30, 乙球径 = 15

 

コメント
  • X
  • Facebookでシェアする
  • はてなブックマークに追加する
  • LINEでシェアする

PVアクセスランキング にほんブログ村

PVアクセスランキング にほんブログ村