算額（その1534）
奈良県大和郡山市 庚申堂 明治13年(1880)
http://www.wasan.jp/nara/kosindo.html
（五）奈良県 (5)庚申堂
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：等差数列，部分和
#Julia, #SymPy, #算額, #和算
一辺の長さが a[i], i=1,2,...,n, a[k] > a[k+1] の等差数列に従う n 個の立方体がある。
それらの立体の体積の和，表面積の和および公差 K = a[k] - a[k+i] が与えられたとき，末尾の（最小の）立方体の一辺の長さを求めよ。
まず最初に，体積の和と表面積の和を一辺の長さの初項 a，公差 K, と項数 n で表す式を決定する。
1. n 個の立方体の体積の和
初項は a^3, 第 2 項までは a^3 + (a - K)^3, 第 3 項までは a^3 + (a - K)^3 + (a - 2K)^3 ,... なので，以下のようにプログラムすれば n ごとの式が得られる。
using SymPy
@syms s, a, K
s = 0
for n = 0:10
s += (a - n*K)^3
s |> expand |> println
end
a^3
-K^3 + 3*K^2*a - 3*K*a^2 + 2*a^3
-9*K^3 + 15*K^2*a - 9*K*a^2 + 3*a^3
-36*K^3 + 42*K^2*a - 18*K*a^2 + 4*a^3
-100*K^3 + 90*K^2*a - 30*K*a^2 + 5*a^3
-225*K^3 + 165*K^2*a - 45*K*a^2 + 6*a^3
-441*K^3 + 273*K^2*a - 63*K*a^2 + 7*a^3
-784*K^3 + 420*K^2*a - 84*K*a^2 + 8*a^3
-1296*K^3 + 612*K^2*a - 108*K*a^2 + 9*a^3
-2025*K^3 + 855*K^2*a - 135*K*a^2 + 10*a^3
-3025*K^3 + 1155*K^2*a - 165*K*a^2 + 11*a^3
ついで，これらの係数の一般項を求める。
ai^3 の項の係数は (n + 1) * K^0
ai^2 の項の係数は -3*n*(n + 1)/2 * K^1
ai^1 の項の係数は n*(n+1)*(2*n+1)/2 * K^2
ai^0 の項の係数は -(n*(n+1)/2)^2 * K^3
である。
まとめると，体積の和は以下のようになる。
@syms n::positive, K
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
V |> println
-K^3*n^2*(n + 1)^2/4 + K^2*a*n*(n + 1)*(2*n + 1)/2 - 3*K*a^2*n*(n + 1)/2 + a^3*(n + 1)
n = 10 のときの式は以下のようになり，前述の式のとおりになっていることが確認できる。
V(n => 10) |> println
-3025*K^3 + 1155*K^2*a - 165*K*a^2 + 11*a^3
2. n 個の立方体の表面積の和
体積の和の一般項を求めるのも前節と同様にして以下のようになる。
using SymPy
@syms s, a, K
s = 0
for i = 0:10
s += 6(a - i*K)^2
s |> expand |> println
end
6*a^2
6*K^2 - 12*K*a + 12*a^2
30*K^2 - 36*K*a + 18*a^2
84*K^2 - 72*K*a + 24*a^2
180*K^2 - 120*K*a + 30*a^2
330*K^2 - 180*K*a + 36*a^2
546*K^2 - 252*K*a + 42*a^2
840*K^2 - 336*K*a + 48*a^2
1224*K^2 - 432*K*a + 54*a^2
1710*K^2 - 540*K*a + 60*a^2
2310*K^2 - 660*K*a + 66*a^2
ai^2 の項の係数は 6(n + 1)
ai^1 の項の係数は -6n*(n+1) * K
ai^0 の項の係数は n*(n+1)*(2n+1) * K^2
@syms n::positive, K
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2
S |> println
K^2*n*(n + 1)*(2*n + 1) - 6*K*a*n*(n + 1) + a^2*(6*n + 6)
S(n => 10) |> println
2310*K^2 - 660*K*a + 66*a^2
3. 算額問題を解く
1, 2 節で得られた一般項の式と具体的な体積と表面積の値から立方体の個数を求める。
テストデータとして，初項 a = 20，公差 K = 3（実際は -3）のときの n = 4 個の立方体の一辺の長さが [20, 17, 14, 11] のとき，体積の和は 16988，表面積の和は 6036 である。
a = [20, 17, 14, 11]
sum(a.^3), sum(6a.^2)
(16988, 6036)
@syms n::positive, a::positive, K::positive, 末尾::positive
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2;
eq1 = V(K => 3) - 16988
eq2 = S(K => 3) - 6036
eq3 = (a - n*K)(K => 3) - 末尾
solve([eq1, eq2, eq3], (末尾, n, a))[1]
(11, 3, 20)
一番小さい立方体の一辺の長さは 11 である。
n = 3 であるが，n は 0 から数えるので，立方体の個数は n + 1 であることに注意。
一番大きい立方体の一辺の長さ（初項 a） は 20 である。
もう一つテストデータとして，初項 a = 100，公差 K = 5（実際は -5）のときの n = 16 個の立方体の一辺の長さは，
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
である。
体積の和は 5500000，表面積の和は 426000 である。
a = collect(range(100, step=-5, length=16));
a |> println
(sum(a.^3), sum(6a.^2)) |> println
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
(5500000, 426000)
@syms n::positive, a::positive, K::positive, 末尾::positive
V = -(n*(n+1)/2)^2*K^3 + n*(n+1)*(2*n+1)/2*K^2*a - 3*n*(n + 1)/2*K*a^2 + (n + 1)*a^3
S = n*(n+1)*(2n+1)*K^2 - 6n*(n+1)*K*a + 6(n + 1)*a^2;
eq1 = V(K => 5) - 5500000
eq2 = S(K => 5) - 426000
eq3 = (a - n*K)(K => 5) - 末尾
solve([eq1, eq2, eq3], (末尾, n, a))[1]
(25, 15, 100)
最小の立方体の一辺の長さは 25，最大の立方体の一辺の長さは 100, 立方体の個数は n + 1 = 16 である。
