算額（その1548）
神壁算法 皇都 安井聖天堂 文化2年()
藤田貞資門人 皇都粟田 町野左一郎好謙
藤田貞資(1807)：続神壁算法
http://www.wasan.jp/jinpeki/zokujinpekisanpo.pdf
キーワード：数列
#Julia, #SymPy, #算額, #和算，#数学
最上段に 1 個，その下に 3 個，更にその下に 6 個，次いで，10，15...のように玉が積まれている。玉の重さは一番上が最も重く，以下順に 2 ずつ軽くなっている。玉の重さは全部で 6860 である。玉の個数は全部で 56 個である。一番上の玉の重さはいかほどか。
一般項を求めるのが肝要である。
一番上の玉の重さを w, n 段あるとする。
段ごとの玉の数は，1, 3, 6, 10, ... で，n 段目にある玉の数は n*(n + 1)/2 である。
最上段から n 段目までにある玉の重さの合計は n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n + 1)*(n + 2)/4 である。
最上段から各段までの玉の数の合計は，1, 4, 10, 20, ... で，n*(n + 1)*(n + 2)/6 である。
最上段から n 段目までにある玉の重さの合計を K1
最上段から各段までの玉の数の合計を K2
とおき，以下の連立方程式を解く。
K1, K2 を実地数とすれば容易に解ける。
問の場合は，最上段の重さは 130 である。
include("julia-source.txt");
# # julia-source.txt ソース https://blog.goo.ne.jp/r-de-r/e/ad3a427b84bb416c4f5b73089ae813cf
@syms w::positive, n::positive, K1::positive, K2::positive
eq1 =n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n+1)*(n+2)/4 - 6860
eq2 = n*(n + 1)*(n + 2)/6 - 56;
res = solve([eq1, eq2], (n, w))[1]
(6, 130)
K1, K2 を記号のまま解くと虚数解になるが，虚部はほとんど 0 なので，実部だけをとればよい。
eq1 =n*(n + 1)*(n + 2)/6 * w - (n - 1)*n*(n+1)*(n+2)/4 - K1
eq2 = n*(n + 1)*(n + 2)/6 - K2;
res = solve([eq1, eq2], (n, w))[1] # 1 of 3
((-(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/3 - (-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) - 3^(1/6)*I*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(2/3) + 3^(5/6)*I*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 2)/((3^(1/3) - 3^(5/6)*I)*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3)), (I*K1*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3^(5/6)*K1*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) - I*K2*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/2 - 3*I*K2*(-81*K2 + 3*sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3*3^(1/6)*K2*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(2/3)/2 - 3*3^(5/6)*K2*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3) + 3*I*K2)/(K2*(3^(5/6) + 3^(1/3)*I)*(-27*K2 + sqrt(3)*sqrt(243*K2^2 - 1))^(1/3)))
res[1](K1 => 6860, K2 => 56).evalf() |> println
res[2](K1 => 6860, K2 => 56).evalf() |> println
6.0 + 5.47542772767019e-26*I
130.0 + 3.28284985518336e-24*I
15 段目までの球の総数は 680 個，重さの合計は 74120 である。
res[1](K1 => 74120, K2 => 680).evalf() |> println
res[2](K1 => 74120, K2 => 680).evalf() |> println
15.0 - 2.03553346724214e-19*I
130.0 - 1.17148914145633e-18*I
