Two Snuke
maspyさんによる解説を理解したい
問題
問題条件でNが与えられる
N以下のnについて、それを5つの数に分ける
この一つに注目する、数をmとする
このmを2つの数i, jに分ける
i + j = m
ここで 0 \le i < jという制約がつけられている
とりうるすべてのi,jについての差の和を求めたい
f(m) = \sum_{i,j} j - i
ここまでは自力で求めたのだが、ここから先がわからなくなった
pythondef f(x):
return (x + x % 2) * ((x + 2) // 2) // 2形式的べき級数で解く
この関数fは自然数を引数に取る関数なので、数列と解釈できる
F = \sum_{m=0}^\infty f(m)
僕のfの実装では剰余や切り捨てが発生していたが、これを有理式で表すことができれば形式的べき級数の道具が使える
数列を観察
pythonIn [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])僕は一般項をコードで書いてるので偶奇が影響するのわかってるから分けてみる
pythonIn [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 + \ldots = \sum_{n=0}^\infty x^{2n}
G = 1 + x^2 + x^4 + \ldots
x^2G = x^2 + x^4 + \ldots
(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
pythonIn [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としておく
4F = G + J
J = \sum_{n=0}^\infty (n+1)^2 x^n
初回は計算ミス をしててこの先で正しくない式を出してしまっていた
pythonIn [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行で計算できる
pythonIn [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を掛けてみたらこうなった
pythonIn [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}
求めるべき形式的べき級数が決まったので、次はそれを解く
多項式剰余を理解して実装すれば良さそう
まず形式的べき級数の偶数次元/奇数次元だけ通す関数を定義する
\mathcal{E}(G(x)) = \frac{G(x) + G(-x)}{2}
\mathcal{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]\mathcal{E}(F)
= [x^n](\frac{P}{Q} + \frac{P^*}{Q^*})/2
= [x^n](\frac{PQ^* + P^*Q}{QQ^*})/2
= [x^n](\frac{\mathcal{E}(PQ^*)}{QQ^*})/2
nが奇数の時
[x^n]F = [x^n]\mathcal{O}(F)
= [x^n](\frac{P}{Q} - \frac{P^*}{Q^*})/2
= [x^n](\frac{PQ^* - P^*Q}{QQ^*})/2
= [x^n](\frac{\mathcal{O}(PQ^*)}{QQ^*})/2
この時 QQ^*は偶関数
Q = q_0 + q_1x + q_2 x^2 + ...とする
Q^* = q_0 - q_1x + q_2 x^2 + ...を掛けると、奇数次の項が消えるから。
つまり数列としては、奇数番目が0の飛び飛びの列
これを圧縮したQ'が作れる
such as Q'(x^2) = QQ^*
形式的べき級数の積なのだから、数列の畳み込みで実装できて、飛び飛びで読めば良いだけ
同じように分子も飛び飛び数列なので圧縮できる
P_e(x^2) = \mathcal{E}(PQ^*)
P_o(x^2) = \mathcal{O}(PQ^*) / x
実装 AC
maspyさんはもっとコンパクトな分母を、おそらく手元での探索によって見つけて使っていた(勘違い)
が、オーパーツなので僕はここまでの理屈どおり\frac{x^5}{(1-x)^{16}(1+x)^5} を使った
下記のコードの
[1,-2,0,2,-1] を見て勘違いしてたが、式は同じだったpythonden = 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)を展開したものだった
pythonIn [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 pythondef 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]Related Pages
- →オイラー路×グリッド上の幅優先探索×ハンガリアン法×二部グラフの最大マッチング×二重辺連結成分分解×二重頂点連結成分分解×全点対間最短路×単一始点最短路×強連結成分分解×彩色数×最大クリーク×最大流×最大独立集合×最小全域有向木×最小全域木×最小流量制限付き最大流×最小費用流×橋/関節点×bit×binary-trie×convex-hull-trick-add-monotone×li-chao-tree×link-cut木_部分木クエリ×link-cut木×ウェーブレット行列×スパーステーブル×スライド区間の昇順k個の和×セグメント木×トライ木×マージ可能ヒープ×列の平方分割×平衡二分探索木×永続配列×素集合データ構造×unionfind×ローリングハッシュ×接尾辞配列×最長共通接頭辞×最長回文×回文×複数文字列検索×hl分解×全方位木dp×最小共通祖先×木の直径×木の重心分解×根付き木に変換×mod-pow×オイラーのφ関数×オイラーのφ関数テーブル×ベル数×ラグランジュ補間×二項係数×二項係数テーブル×任意mod畳み込み×分割数×分割数テーブル×商列挙×形式的冪級数×形式的べき級数×拡張ユークリッドの互除法×拡張ユークリッド互除法×第2種スターリング数×約数列挙×素因数分解×素数テーブル×素数判定×組合せ×行列演算×進数変換×階乗×離散対数問題×高速フーリエ変換×divide-and-conquer-optimization×monotone-minima×スライド最小値×一次元累積和×二次元累積和×個数制限付きナップサック×最大長方形×最適二分探索木×最長増加部分列×ダブリング×包除原理×燃やす埋める問題×燃やす埋める×牛ゲー×mo’s_algorithm×offline-dynamic-connectivity×座標圧縮×アルゴリズム→
- →ホッケースティック恒等式×二項係数の公式×vandermondeの恒等式×巨大なnについての二項係数×形式的べき級数×ARC110D_FPS×oeis×重複組合せの畳み込み×arc110d_nonfps×arc110→
- →地図読み込み時に番兵をつける×数列を有理式にする×ラグランジュ補間×対称性で次元削減×余事象を引く×未ac×足し算の順序を変える×xごとのf(x,y)の和はyごとのf(x,y)の和×最大値の期待値×順位統計量×maxの不等号は不等号のand×0,1→
"Engineer's way of creating knowledge" the English version of my book is now available on [Engineer's way of creating knowledge]
(C)NISHIO Hirokazu / Converted from [Scrapbox] at [Edit]