算額（その1535）
滋賀県新旭町太田 大田神社 (1868)
（四）滋賀県 (8)大田神社
近畿数学史学会：近畿の算額「数学の絵馬を訪ねて」，平成4年5月16日 初版第一刷，大阪教育図書株式会社，大阪市.
キーワード：等差数列，部分和
#Julia, #SymPy, #算額, #和算
一辺の長さが a[i], i=1,2,...,n, a[k] > a[k+1] の等差数列に従う n 個の正方形がある。
それらの正方形の面積の和，一辺の長さの和および公差 K = a[k] - a[k+i] が与えられたとき，正方形の総数を求めよ。
算額（その1534）の類題である。こちらのほうが易しい。
まず最初に，面積の和と一辺の長さの和を一辺の長さの初項 a，公差 K, と項数 n で表す式を決定する。
1. n 個の正方形の面積の和
初項は a^2, 第 2 項までは a^2 + (a - K)^2, 第 3 項までは a^2 + (a - K)^2 + (a - 2K)^2 ,... なので，以下のようにプログラムすれば n ごとの式が得られる。
using SymPy
@syms s, a, K
s = 0
for n = 0:10
s += (a - n*K)^2
s |> expand |> println
end
a^2
K^2 - 2*K*a + 2*a^2
5*K^2 - 6*K*a + 3*a^2
14*K^2 - 12*K*a + 4*a^2
30*K^2 - 20*K*a + 5*a^2
55*K^2 - 30*K*a + 6*a^2
91*K^2 - 42*K*a + 7*a^2
140*K^2 - 56*K*a + 8*a^2
204*K^2 - 72*K*a + 9*a^2
285*K^2 - 90*K*a + 10*a^2
385*K^2 - 110*K*a + 11*a^2
ついで，これらの係数の一般項を求める。
ai^2 の項の係数は (n + 1)
ai^1 の項の係数は -n*(n+1) * K
ai^0 の項の係数は n*(n+1)*(2*n+1)/6 * K^2
である。
まとめると，体積の和は以下のようになる。
@syms n::positive, K
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
S |> println
K^2*n*(n + 1)*(2*n + 1)/6 - K*a*n*(n + 1) + a^2*(n + 1)
n = 10 のときの式は以下のようになり，前述の式のとおりになっていることが確認できる。
S(n => 10) |> println
385*K^2 - 110*K*a + 11*a^2
2. n 個の正方形の一辺の長さの和
一辺の長さの和の一般項を求めるのも前節と同様にして以下のようになる。
using SymPy
@syms s, a, K
s = 0
for i = 0:10
s += (a - i*K)
s |> expand |> println
end
a
-K + 2*a
-3*K + 3*a
-6*K + 4*a
-10*K + 5*a
-15*K + 6*a
-21*K + 7*a
-28*K + 8*a
-36*K + 9*a
-45*K + 10*a
-55*K + 11*a
ai^1 の項の係数は (n+1)
ai^0 の項の係数は n*(n+1)/2 * K
@syms n::positive, K
L = -n*(n+1)/2*K + (n + 1)*a
L |> println
-K*n*(n + 1)/2 + a*(n + 1)
L(n => 10) |> println
-55*K + 11*a
3. 算額問題を解く
1, 2 節で得られた一般項の式と具体的な面積の和と一辺の長さの和から正方形の個数を求める。
テストデータとして，初項 a = 20，公差 K = 3（実際は -3）のときの n = 4 個の正方形の一辺の長さが [20, 17, 14, 11] のとき，面積の和は 1006，一辺の長さの和は 62 である。
a = [20, 17, 14, 11]
sum(a.^2), sum(a)
(1006, 62)
@syms n::positive, a::positive, K::positive, 末尾::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => 3) - 1006
eq2 = L(K => 3) - 62
solve([eq1, eq2], (n, a))[1]
(3, 20)
n = 3 であるが，n は 0 から数えるので，正方形の個数は n + 1 = 4 個である。
もう一つのテストデータとして，初項 a = 100，公差 K = 5（実際は -5）のときの n = 16 個の正方形の一辺の長さは，
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
である。
面積の和は 71000，一辺の長さの和は 1000 である。
a = collect(range(100, step=-5, length=16));
a |> println
(sum(a.^2), sum(a)) |> println
[100, 95, 90, 85, 80, 75, 70, 65, 60, 55, 50, 45, 40, 35, 30, 25]
(71000, 1000)
@syms n::positive, a::positive, K::positive, 末尾::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => 5) - 71000
eq2 = L(K => 5) - 1000
solve([eq1, eq2], (n, a))[1]
(15, 100)
正方形の個数は n + 1 = 16 個である。
4. 一般解を求める
面積の和 K1，一辺の長さの和 K2，公差 K3 を記号定数として一般解を求める。
一度には解けないので，段階を追って解く。
using SymPy
@syms n::positive, a::positive, K::positive, K1::positive, K2::positive, K3::positive
S = n*(n+1)*(2*n+1)/6*K^2 - n*(n+1)*K*a + (n + 1)*a^2
L = -n*(n+1)/2*K + (n + 1)*a
eq1 = S(K => K3) - K1
eq2 = L(K => K3) - K2;
ans_a = solve(eq2, a)[1];
eq11 = eq1(a => ans_a) |> simplify
ans_n = solve(eq11, n)[3]; # 3 of 4
ans_n |> println
(-6*K3^(5/2) + sqrt(3)*sqrt(K3^(13/3)*(144*K2^2 + K3^2)/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + K3^(11/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + 2*K3^5) - sqrt(3)*sqrt(72*sqrt(3)*K1*K3^(22/3)/sqrt((2*K3^(26/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + K3^(22/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(2/3) + K3^8*(144*K2^2 + K3^2))/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3)) - K3^(13/3)*(144*K2^2 + K3^2)/(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) - K3^(11/3)*(1944*K1^2 + 432*K2^2*K3^2 - K3^4 + sqrt(-K3^2*(144*K2^2 + K3^2)^3 + (1944*K1^2 + 432*K2^2*K3^2 - K3^4)^2))^(1/3) + 4*K3^5))/(6*K3^(5/2))
大変長い式であるが，ちゃんと数式解は求まる。
ans_n(K1 => 71000, K2 => 1000, K3 => 5).evalf() |> println
15.0000000000000
