Two Snuke

エイシングプログラミングコンテスト 2020
maspyさんによる解説を理解したい

問題

問題条件でNが与えられる
N以下のnについて、それを5つの数に分ける
この一つに注目する、数をmとする
このmを2つの数i, jに分ける
- $i + j = m$
- ここで $0 \leq i < j$ という制約がつけられている
とりうるすべてのi,jについての差の和を求めたい
- $f (m) = \sum_{i, j} j - i$
- ここまでは自力で求めたのだが、ここから先がわからなくなった python

def f(x):
    return (x + x % 2) * ((x + 2) // 2) // 2

形式的べき級数で解く

この関数fは自然数を引数に取る関数なので、数列と解釈できる
数列にはその母関数である形式的べき級数が対応する、これをFとする
- $F = \sum_{m = 0}^{\infty} f (m)$
- 僕のfの実装では剰余や切り捨てが発生していたが、これを有理式で表すことができれば形式的べき級数の道具が使える
- 数列を観察 python

In [85]: np.array([f(x) for x in range(0, 100)])
Out[85]: 
array([   0,    1,    2,    4,    6,    9,   12,   16,   20,   25,   30,
         36,   42,   49,   56,   64,   72,   81,   90,  100,  110,  121,
        132,  144,  156,  169,  182,  196,  210,  225,  240,  256,  272,
        289,  306,  324,  342,  361,  380,  400,  420,  441,  462,  484,
        506,  529,  552,  576,  600,  625,  650,  676,  702,  729,  756,
        784,  812,  841,  870,  900,  930,  961,  992, 1024, 1056, 1089,
       1122, 1156, 1190, 1225, 1260, 1296, 1332, 1369, 1406, 1444, 1482,
       1521, 1560, 1600, 1640, 1681, 1722, 1764, 1806, 1849, 1892, 1936,
       1980, 2025, 2070, 2116, 2162, 2209, 2256, 2304, 2352, 2401, 2450,
       2500])

- 僕は一般項をコードで書いてるので偶奇が影響するのわかってるから分けてみる

python

In [94]: fs[::2]
Out[94]: 
array([   0,    2,    6,   12,   20,   30,   42,   56,   72,   90,  110,
        132,  156,  182,  210,  240,  272,  306,  342,  380,  420,  462,
        506,  552,  600,  650,  702,  756,  812,  870,  930,  992, 1056,
       1122, 1190, 1260, 1332, 1406, 1482, 1560, 1640, 1722, 1806, 1892,
       1980, 2070, 2162, 2256, 2352, 2450])

In [95]: fs[1::2]
Out[95]: 
array([   1,    4,    9,   16,   25,   36,   49,   64,   81,  100,  121,
        144,  169,  196,  225,  256,  289,  324,  361,  400,  441,  484,
        529,  576,  625,  676,  729,  784,  841,  900,  961, 1024, 1089,
       1156, 1225, 1296, 1369, 1444, 1521, 1600, 1681, 1764, 1849, 1936,
       2025, 2116, 2209, 2304, 2401, 2500])

- この観察の時点でmaspyさんは$F = \frac{P}{(1-x^2)^3}$が確定したと書いているが僕にはそれがわからなかったので噛み砕く

一つ飛ばしの数列

1, 0, 1, 0, … という数列を形式的べき級数でどう表現するのか
素朴に式を書いてみると
$G = 1 + x^{2} + x^{4} + \dots = \sum_{n = 0}^{\infty} x^{2 n}$
形式的べき級数の逆元を使った無限和圧縮を使う
$G = 1 + x^{2} + x^{4} + \dots$
$x^{2} G = x^{2} + x^{4} + \dots$
$(1 - x^{2}) G = 1$
$G = 1/ (1 - x^{2})$

一つ飛ばしの平方数

0, 1, 0, 4, 0, 9, … という数列を形式的べき級数でどう表現するのか
これは $H = \sum_{n = 0}^{\infty} ((n + 1) /2)^{2} x^{n}$ の奇数番目だけ取ったやつ
- なお数列が1-originか0-originかで混乱したので0に揃えました
引いてみた $F - H$ python

In [101]: fs = np.array([f(x) - ((x+1)/2)**2 for x in range(100)])
In [102]: fs
Out[102]: 
array([-0.25,  0.  , -0.25,  0.  , -0.25,  0.  , -0.25,  0.  , -0.25,
        0.  , -0.25,  0.  , -0.25,  0.  , -0.25,  0.  , -0.25,  0.  ,...

つまり $F = G /4 + H$
- 計算が面倒なので全部4倍して4H=Jとしておく
- $4 F = G + J$
- $J = \sum_{n = 0}^{\infty} (n + 1)^{2} x^{n}$
- 初回は計算ミスをしててこの先で正しくない式を出してしまっていた python

In [35]: reduce(np.convolve, [[1, -1]] * 0,  np.array([(x + 1) ** 2 for x in range(100)]))[:20]
Out[35]: 
array([  1,   4,   9,  16,  25,  36,  49,  64,  81, 100, 121, 144, 169,
       196, 225, 256, 289, 324, 361, 400])

In [36]: reduce(np.convolve, [[1, -1]] * 1,  np.array([(x + 1) ** 2 for x in range(100)]))[:20]
Out[36]: 
array([ 1,  3,  5,  7,  9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33,
       35, 37, 39])

In [37]: reduce(np.convolve, [[1, -1]] * 2,  np.array([(x + 1) ** 2 for x in range(100)]))[:20]
Out[37]: array([1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

In [38]: reduce(np.convolve, [[1, -1]] * 3,  np.array([(x + 1) ** 2 for x in range(100)]))[:20]
Out[38]: array([1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

- $J = \frac{1 + x}{(1-x)^3}$
- (1-x)^3は3回差を取ることに対応していたのか
- $4F = G + J = \frac{1}{1-x^2} + \frac{1 + x}{(1-x)^3} = \frac{(1 - x)^2 + (1 + x)^2}{(1 + x)(1 - x)^3} = \frac{4x}{(1 + x)(1 - x)^3}$

もっと簡単な式を求めて

maspyさんのやってる「分母に(1-x)が来そうだから掛けてみる」を試す数列を有理式にする
- 形式的べき級数の積は、数列の畳み込み
- np.convolveを使えば1行で計算できる python

In [120]: reduce(np.convolve, [[1, -1], [1, -1], [1, -1], [1, 1]], fs)
Out[120]: 
array([    0,     1,     0,     0,     0,     0,     0,     0,     0, ...

    - (1 - x)(1 - x)(1 - x)(1 + x) を掛けたらxが残った
    - $\frac{x}{(1-x)^3(1+x)}$

僕が途中で計算間違いをしてたのだろうか？後で検証する。
ちなみにmaspyさんの当初の話どおり $(1 - x^{2})^{3}$ を掛けてみたらこうなった python

In [121]: reduce(np.convolve, [[1, 0, -1], [1, 0, -1], [1, 0, -1]], fs)
Out[121]: 
array([    0,     1,     2,     1,     0,     0,     0,     0,     0, ...

- $\frac{x(1 + 2 x + x^2)}{(1-x^2)^3} = \frac{x(1 + x)^2}{(1-x)^3(1+x)^3} = \frac{x}{(1-x)^3(1 + x)}$
- 整理すれば結果は同じ

Fが求まったとする

$F = \frac{x}{( 1 - x ) ^{3} ( 1 + x )}$
Fは5つに分けたブロックの1つについてのべき級数だった
- これはmを引数とする関数と解釈できて、合計がmである時の答えだった
  - $[x^{m}] F = f (m) = \sum_{i, j} j - i$
$[x^{n}] F^{5}$ が、5個のブロックの総和がnの時の答え
問題条件は「総和がN以下の時」なので、得たい最終的な答えは
- $[x^{n}] \sum_{n = 0}^{\infty} F^{5}$
この和を取るところと(1-x)で割ることはどう関係するのだっけ？？？
- 形式的べき級数の逆元を使った無限和圧縮かと思ったけど違うな…
- (1-x)で割るということは $\sum_{n = 0}^{\infty} x^{n}$ を掛けるということ
- 形式的べき級数の積は、数列の畳み込み。すべて1の数列と畳み込んでいるので、畳み込み結果の第N項は、元の数列のNまでの和になる
- $\sum_{n = 0}^{\infty} F^{5} = F^{5} / (1 - x) = \frac{x ^{5}}{( 1 - x ) ^{16} ( 1 + x ) ^{5}}$ 求めるべき形式的べき級数が決まったので、次はそれを解く
多項式剰余を理解して実装すれば良さそう
まず形式的べき級数の偶数次元/奇数次元だけ通す関数を定義する
- $E (G (x)) = \frac{G ( x ) + G ( - x )}{2}$
- $O (G (x)) = \frac{G ( x ) - G ( - x )}{2}$
- 要するにx軸で鏡映して足し合わせると奇関数は打ち消しあって消滅するということ
- ここは説明のために(x)を明記したけど以下では省く, $G^{*} = G (- x)$ とする
$[x^{n}] F = P / Q$ を求めたい(Q(0)=1)
- n=0のとき
- $[x^{n}] F = P (0)$
  - これがループの終了条件
- nが偶数の時
- $[x^{n}] F = [x^{n}] E (F)$
- $= [x^{n}] (\frac{P}{Q} + \frac{P ^{*}}{Q ^{*}}) /2$
- $= [x^{n}] (\frac{P Q ^{*} + P ^{*} Q}{Q Q ^{*}}) /2$
- $= [x^{n}] (\frac{E ( P Q ^{*} )}{Q Q ^{*}}) /2$
- nが奇数の時
- $[x^{n}] F = [x^{n}] O (F)$
- $= [x^{n}] (\frac{P}{Q} - \frac{P ^{*}}{Q ^{*}}) /2$
- $= [x^{n}] (\frac{P Q ^{*} - P ^{*} Q}{Q Q ^{*}}) /2$
- $= [x^{n}] (\frac{O ( P Q ^{*} )}{Q Q ^{*}}) /2$
- この時 $Q Q^{*}$ は偶関数
  - $Q = q_{0} + q_{1} x + q_{2} x^{2} + ...$ とする
  - $Q^{*} = q_{0} - q_{1} x + q_{2} x^{2} + ...$ を掛けると、奇数次の項が消えるから。
  - つまり数列としては、奇数番目が0の飛び飛びの列
  - これを圧縮したQ’が作れる
    - such as $Q^{'} (x^{2}) = Q Q^{*}$
    - 形式的べき級数の積なのだから、数列の畳み込みで実装できて、飛び飛びで読めば良いだけ
  - 同じように分子も飛び飛び数列なので圧縮できる
  - $P_{e} (x^{2}) = E (P Q^{*})$
  - $P_{o} (x^{2}) = O (P Q^{*}) / x$

実装 AC

maspyさんはもっとコンパクトな分母を、おそらく手元での探索によって見つけて使っていた(勘違い)
- が、オーパーツなので僕はここまでの理屈どおり $\frac{x ^{5}}{( 1 - x ) ^{16} ( 1 + x ) ^{5}}$ を使った
- 下記のコードの[1,-2,0,2,-1]を見て勘違いしてたが、式は同じだった python

den = np.array([1,-1], np.int64)
for _ in range(5):
    den = np.convolve(den, np.array([1,-2,0,2,-1]))

- これは$(1-x)^3(1+x)$を展開したものだった

python

In [287]: reduce(np.convolve, [[1, -1]] * 3 + [[1, 1]], [1])
Out[287]: array([ 1, -2,  0,  2, -1])

maspyさんはFFTを使った高速な畳み込みを使っていたけど、疲れたのでそちらを読むのはやめた
- np.convoluteで済まそうとしたのだけど、int64をオーバーフローするので仕方なく自分で作った
Qm: $Q^{*}$
- Qの奇数次の係数を符号反転するだけ
- Qm = Q.copy() Qm[1::2] *= -1 python

def solve(N):
    Q = reduce(
        np.convolve,
        [[1, -1]] * 16 + [[1, 1]] * 5,
        np.array([1], dtype=np.int64))
    P = np.zeros(6, dtype=np.int64)
    P[5] = 1

    def conv(X, Y):
        n = X.shape[0]
        m = Y.shape[0]
        ret = np.zeros(n + m - 1, dtype=np.int64)
        for i in range(n):
            ret[i:i + m] += X[i] * Y
            ret[i:i + m] %= MOD
        return ret

    while N:
        Qm = Q.copy()
        Qm[1::2] *= -1
        QQm = conv(Q, Qm)
        PQm = conv(P, Qm)
        if N % 2:
            P = PQm[1::2]
        else:
            P = PQm[::2]
        Q = QQm[::2]
        N //= 2
    return P[0]

🪴 Quartz 4.0

Two Snuke

Graph View

Backlinks