算額（その1549）
百三 群馬県高崎市八幡町 八幡宮 安政7年(1860)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：数列，部分和
#Julia, #SymPy, #算額, #和算，#数学
逐増する「平方垜」がある，第 1 項（底子 1）では 1 個，第 2 項までの合計は 6 個（5 個増えたことになる。この増え方が平方垜になる），第 3 項までの合計は 20個（14 個増えたことになる）であるとき，第 n 項までの合計は何個か。
注：垜（だ）とは，球，土嚢などを一番下から順次規則的に積み上げた構造物である。
規則はいろいろあるが
1. 三角数の積み上げ
三角数は、上から順に**1、2、3…**と増加する層で構成される。Tn = n(n + 1)/2
2. 四角数の積み上げ
四角数は正方形の形を形成する。Qn = n^2
3. 五角数の積み上げ
五角数は星型や角形を意識した積み方Pn = n(3n - 1)/2
4. 六角数（蜂の巣型） Hn = n(2n - 1)
など
平方垜は四角数の積み上げと思いきや，合計が 1, 6, 20, ... となり，差分が（問題文にもあるが）5, 14 となる。例が少ないが前を補うと 0, 1, 5, 14, ... で，更に差分を取ると 1, 4, 9, ...となるこれは平方数の数列であろうと推測できる。
当時の人は，平方垜といえばすぐにこの数列を書き下せるほど，常識だったのであろう。
ひとまずその仮定に基づくと，平方垜の第 n 項までの合計の数列は，
0, 1, 6, 20, 50, 105, ... 合計
0, 1, 5, 14, 30, 55, ... 各段の個数
1, 4, 9, 16, 25, ... 差分
「オンライン整数列大辞典」を参照すれば以下のものが見つかる。
A002415 4-dimensional pyramidal numbers: a(n) = n^2*(n^2-1)/12.
0, 0, 1, 6, 20, 50, 105, 196, 336, ...
ただし，この数式は n = 0 から始まっているので，その調整をする必要がある。
using SymPy
@syms n
an = n^2*(n^2-1)/12 # A002415
an |> println
n^2*(n^2 - 1)/12
an(n => n + 1) |> factor|> println # 項の調整
n*(n + 1)^2*(n + 2)/12
各項を計算して，推定どおりであることを確認する。
for n = 1:10
println("n = $n, an = $(n*(n + 1)^2*(n + 2)/12)")
end
n = 1, an = 1.0
n = 2, an = 6.0
n = 3, an = 20.0
n = 4, an = 50.0
n = 5, an = 105.0
n = 6, an = 196.0
n = 7, an = 336.0
n = 8, an = 540.0
n = 9, an = 825.0
n = 10, an = 1210.0
「術」は n*(n^3 +4n^2 + 5n + 2)/12 であるが，因数分解すると同じであることがわかる。
n*(n^3 +4n^2 + 5n + 2)/12 |> factor |> println
n*(n + 1)^2*(n + 2)/12
整数列大辞典に見つからないような場合は，以下のように多項式の係数を推定するというアプローチを取る。
using SymPy
@syms a1, a2, a3, a4
f(n) = a1*n + a2*n^2 + a3*n^3 + a4*n^4;
eq1 = f(1) - 1
eq2 = f(2) - 6
eq3 = f(3) - 20
eq4 = f(4) - 50
res = solve([eq1, eq2, eq3, eq4], (a1, a2, a3, a4))
res |> println
Dict{Sym{PyCall.PyObject}, Sym{PyCall.PyObject}}(a3 => 1/3, a1 => 1/6, a2 => 5/12, a4 => 1/12)
一般項を表す，同じ多項式が得られた。
@syms n
res[a1]*n + res[a2]*n^2 + res[a3]*n^3 + res[a4]*n^4 |> factor |> println
n*(n + 1)^2*(n + 2)/12
一般項がわかっている場合の和の数列は summatin 関数で求めることができる。
using SymPy
@syms n, m
# 一般項
# an = n*(n + 1)*(2n + 1)/6
# 和の一般項
summation(n*(n + 1)*(2n + 1)/6, (n, 1, m)) |> factor |> println
m*(m + 1)^2*(m + 2)/12
