算額（その1417）

算額（その1417）

百三十八 利根郡月夜野町上津 八幡神社 明治22年(1889)
群馬県和算研究会:群馬の算額，上武印刷株式会社，高崎市，1987年3月31日.
キーワード：数列，和
#Julia, #SymPy, #算額, #和算

蕎麦の実のように（正四面体）積み上げた玉が全部で 1769億7205万8701 個ある。一番上が 1 個のとき，底面（一番下）には何個の玉があるか。

参照：算額（その1402），算額（その1023）

ちょっとした落とし穴がある。単純に4つの球が互いに接し合うようにする場合は，配置法は 1 通りしかないが，2層以上に積み上げる場合には，何通りかあるようだ。

根本的には「六方最密構造」と「面心立方格子」の違いだそうである。
「化学の質問から自分の思い込みに気がついた話(結晶構造)」

算額（その1402）や算額（その1023）のような配置法は，「六方最密構造」である。上から満たずで黄色で示したものが最上位の球である。

1 層だけの場合：1 層目にある玉の数は 1，総数は 1
2 層ある場合：2 層目にある玉の数は 2+1 = 3，総数は 4
3 層ある場合：3 層目にある玉の数は 4+3+2+1 = 10，総数は 14
4 層ある場合：4 層目にある玉の数は 5+4+3+2+1 = 15，総数は 29
5 層ある場合：5 層目にある玉の数は 6+5+4+3+2+1 = 21，総数は 50
のようになる。

もう一つの配置法は，「面心立方格子」である。

1 層だけの場合：1 層目にある玉の数は 1，総数は 1
2 層ある場合：2 層目にある玉の数は 2+1 = 3，総数は 4
3 層ある場合：3 層目にある玉の数は 3+2+1 = 6，総数は 10
4 層ある場合：4 層目にある玉の数は 4+3+2+1 = 10，総数は 20
5 層ある場合：5 層目にある玉の数は 5+4+3+2+1 = 15，総数は 35
である。後者の場合は各層にある玉の数はいわゆる三角数である。六方最密構造と比べて隙間が少なく，密度が高いことがわかる。

算額は三角数が前提のようで，面心立方格子を指しているようだ。
ということで以下では面心立方格子を採用する。

1.　ブルートフォース

ブルートフォース（brute force）とは，問題を解決するために可能な全ての選択肢や組み合わせを網羅的に試す方法のことである。この手法は非常に単純で，特定の問題に対して確実に解を見つけることができるが，効率的ではないことが多い。

以下のようなプログラムにより，「上から 10201 番目の層の玉の数は 52035301 個で，それまでの層にある玉の総数が 176972058701 個である」ことが見つかる。

いい加減にプログラムしたり，長精度数をサポートしていない言語を使用すると，整数オーバーフローで解が求まらないことがあるかもしれないので注意。

# Julia によるプログラム
using Printf
function s(n = 99999)
    t = 0
    for i in 1:n
        a = Int(i*(i + 1)/2)
        t += a
        if t == 176972058701
            println("第 $i 層の玉の数は $a 個，そこまでの玉の総数は $t 個")
            return
        end
    end
end;
s()

    第 10201 層の玉の数は 52035301 個，そこまでの玉の総数は 176972058701 個

Python の場合には以下のようにする。

    >>> def s():
    ...     t = 0
    ...     for i in range(1, 999999):
    ...         a = i*(i + 1)/2
    ...         t += a
    ...         if t == 176972058701:
    ...             print(i, a, t)
    ...             return
    ...    
    ... 
    >>> s()
    10201 52035301.0 176972058701.0

この問題ではブルートフォースでも難なく解を得られるが，多くの場合はブルートフォースだと計算時間がかかりすぎて解が求めにくいということもある。

そんなときには，何らかの法則性を利用することも一法である。

2. 賢い方法

今の場合，上から n 番目の層にある玉の数は a_n = n(n + 1)/2 個，n 番目の層までにある玉の総数は t＿n = n(n + 1)(n + 2)/6 個である。
 n   a_n   t_n
 1     1     1
 2     3     4
 3     6    10
 4    10    20
 5    15    35
 6    21    56
 7    28    84
 8    36   120
 9    45   165
10    55   220
 :     :     :

t_n = 176972058701 のときの a_n を求めればよい。
以下の連立方程式を解く。

using SymPy
@syms n, a
eq1 = n*(n + 1)*(n + 2)/6 - 176972058701
eq2 = n*(n + 1)/2 - a
solve([eq1, eq2], (a, n))[1]

    (52035301, 10201)

10201 番目の層には 52035301 個の玉があり，1番目から 10201 番目の層にある玉の総数は 176972058701 個である。

また，SymPy には数列の部分和を求める関数 summation があるので，一般項の公式 a_n があれば部分和の公式 t_n は不要で，以下のようにすることもできる。

@syms i
eq1 = summation(i*(i + 1)/2, (i, 1, n)) - 176972058701
eq2 = n*(n + 1)/2 - a
solve([eq1, eq2], (a, n))[1]

    (52035301, 10201)

summation(i*(i + 1)/2, (i, 1, n)) |> factor |> println

    n*(n + 1)*(n + 2)/6

Python の sympy では以下のようにする。

    >>> from sympy import var, solve, summation
    >>> var('a, n, i')
    (a, n, i)
    >>> eq1 = n*(n + 1)*(n + 2)/6 - 176972058701
    >>> eq2 = n*(n + 1)/2 - a
    >>> solve([eq1, eq2], (a, n))[0]
    (52035301, 10201)

    >>> eq1 = summation(i*(i + 1)/2, (i, 1, n)) - 176972058701
    >>> eq2 = n*(n + 1)/2 - a
    >>> solve([eq1, eq2], (a, n))[0]
    (52035301, 10201)