累積和とゼータ変換

AtCoder過去問チャレンジとか単元別の解説記事を呼んだりしていた。 けど最近は解くのをそぼりがち。というのも競プロなブログを読んでいて累積和やゼータ変換って言葉の使い方で混乱したため。 これを解消するために数列の累積和や畳み込み積について整理する。まず大枠をぱっと書く。

きっかけはこの記事この記事です。

結論

累積和と呼ばれるものは色々あるが一部の亜種はゼータ変換ではない。でも多くの累積和はゼータ変換となっている。 これについて説明して関連する話題を補足としてまとめる。

  • 資料
  • 結論に向かって
  • 補足的な話題

資料

結論に向かって

数列をA Aその要素をA[n] A[n]と書く。 下のような形式の数列の積を畳み込み積と総称する。

C[k]=Pred(i,j,k)A[i]B[j] \begin{aligned} C[k] = \sum_{\it{Pred}(i,j,k)}A[i]B[j] \end{aligned}

ただし有意義な畳み込み積は数列を写像や級数などと対応づけられるようになっていることが多い。 よくある有意義な畳み込み積は添字集合を半群とし値集合を可換環として以下のように書かれている。

C[k]=ij=kA[i]B[j] \begin{aligned} C[k] = \sum_{i*j=k}A[i]B[j] \end{aligned}

添字集合が半群であることは総和の対象を指定するところに使っている。可換環の性質は総和の計算で使っている。

話題にしてるゼータ変換は隣接代数(Incidence algebra)で定義される。 そこではゼータ変換はゼータ関数との畳み込み積を指す。

この隣接代数(Incidence algebra)は添字集合が局所有限半順序(poset)を前提とした世界。 これは前述した畳み込み積よりも強い構造を持っている。

C[k]=ij=kA[i]B[j] \begin{aligned} C[k] = \sum_{i*j=k}A[i]*B[j] \end{aligned}

半群に単位元と可換性を加えた可換モノイドが導くのが前順序なので、局所有限半順序はなかなかに強い性質だ。 だけど普通の配列は添字に自然数として使っているので、自然数の順序を採用して可換モノイドとなる演算で考えている限りはたぶん大丈夫。 (本当に大丈夫か自信がない。) 数列の添字を剰余群に対応させてたりすると難しそう。

前式(半順序つき半群上の畳み込み)のB[j] B[j]をゼータ関数に置き換えるとゼータ変換の式になる。

B[k]=ikA[i]1 \begin{aligned} B[k] &= \sum_{i \le k}A[i]*1 \end{aligned}

最初に累積和の一部の亜種はゼータ変換ではないと書いた。 ゼータ変換は隣接代数に基づいた畳み込み積を前提としていた。 この畳み込み積は配列の値集合が可換環であることを前提としている。 だが一部の累積和の亜種は最小値(min(a,b) \it{min}(a,b))など可逆的でない計算を行う。 そのためゼータ変換とは呼べないと考えた。また逆変換となるメビウス変換を定義できない。 そもそも、これらの演算は加法には当たらず積なので最初から累積積と呼んだ方が良い。

補足的な話題

半群Gの形式的R係数線型結合と畳み込み

はじめに記号を整理する。 数列をA Aその要素をA[n] A[n]と書く。 配列は写像や形式的級数に対応させて考えると累積和を自然に理解できる。 添字を半群G G値を環R Rとすると配列の集合は下のようになる。

A=iGRi A = \bigotimes_{i \in G}{R}_i

数列の畳み込み積を定義するときに添字が半群であることや値が環が使われる。

畳み込み(convolution): C[k]=ij=kA[i]B[j]C=AB \begin{aligned} \text{畳み込み(convolution): } C[k] &= \sum_{i*j = k}A[i]B[j] \\ C &= A \ast B \\ \end{aligned}

半群として自然数N N に積(× \times)や和(+ +)、他に最大値(max(x,y) max(x,y))を導入したものなど色々と選べる。(bit演算ももちろん使える) G G が連続であるなら積分などを用いて同様に畳み込みを定義できる。(今回は数列の話なので連続値は扱わない) このように定義すると数列の積を「半群Gを形式的R係数線型結合(R[G])という代数に拡張したときの積」として解釈できる。

この畳み込み積は以下を満たす。

  • 積の和に対する分配法則
  • R[G] の積の結合法則
  • R[G] の和の結合法則,交換法則
  • G が交換法則を満たす場合、 R[G] も積の交換法則
  • G が単位元を持つときには、 R[G] も単位元を持つ

R[G]と対応した畳み込み積の具体例

畳み込みの具体例としてこんなものがある。他にも i,j をビットベクトルとして捉えAND,OR,XORなどを演算として採用するようなものもある。

C[k]=i+j=kA[i]B[j]C[k]=ij=kA[i]B[j]C[k]=max(i,j)=kA[i]B[j]C[k]=min(i,j)=kA[i]B[j]C[k]=gcd(i,j)=kA[i]B[j]C[k]=lcm(i,j)=kA[i]B[j] \begin{aligned} C[k] &= \sum_{i+j=k}A[i]B[j] \\ C[k] &= \sum_{ij=k}A[i]B[j] \\ C[k] &= \sum_{\max(i,j)=k}A[i]B[j] \\ C[k] &= \sum_{\min(i,j)=k}A[i]B[j] \\ C[k] &= \sum_{\gcd(i,j)=k}A[i]B[j] \\ C[k] &= \sum_{\it{lcm}(i,j)=k}A[i]B[j] \end{aligned}

隣接代数に対応した畳み込み積

前述の畳み込み積に使われている半群 G は実は全て可換モノイドである。 可換モノイドは以下の順序を導く。

ab:=c,ac=b a \le b := \exist c, a*c = b

またこれらは反対称性も満たすため半順序になっており隣接代数と対応づけられる。 具体例としてmin(i,j) \min(i,j)を使って隣接代数と紐づく数列を定義する。

min(i,j) \min(i,j)から導く隣接代数・畳み込み積・ゼータ変換

項目 定義
順序 iminj i \le_{\min} j k,min(i,k)=jij \exist k, \min(i,k) = j \Leftrightarrow i \ge j
単位元 + +\infty
隣接代数の積 (fg)(a,b)=aminxminbf(a,x)g(x,b) (f \ast g)(a,b) = \sum_{a \le_{\min} x \le_{\min} b}f(a, x)g(x, b)

これで隣接代数の積が定義できた。 注意点としてはmin \le_{\min}は自然数の順序の逆になっているところ。 また隣接代数は閉区間[a,b] [a,b] に対する関数を扱うところ。

ここから隣接代数の積と関連するゼータ変換を定義する。 まず隣接代数は閉区間に対する関数なので数列との対応づけで工夫が必要になる。 とりあえず端点の一方を固定する。 f(a,b) f(a, b)のどちらかを固定する。

それぞれでやってみる。

まずa=n a = nのとき。 隣接代数と対応する畳み込みは下の通り。 

(fg)(n,b)=nminxminbf(n,x)g(x,b) (f \ast g)(n,b) = \sum_{n \le_{\min} x \le_{\min} b}f(n, x)g(x, b) H[b]=bxnF[x]G[x][b] H[b] = \sum_{b \le x \le n}F[x]G[x][b]

つづいてb=1 b = 1のとき。 隣接代数と対応する畳み込みは下の通り。 

(fg)(a,1)=aminxmin1f(a,x)g(x,1) (f \ast g)(a,1) = \sum_{a \le_{\min} x \le_{\min} 1}f(a, x)g(x, 1) H[a]=1xaF[a][x]G[x] H'[a] = \sum_{1 \le x \le a}F[a][x]G[x]

これらの2変数関数として残していた畳み込みの相手を隣接代数のゼータ関数(ζ(a,b) \zeta(a,b))を代入する。 それによりゼータ変換が手に入る。

ζ(a,b)={1(if: aminb)0(otherwise) \zeta(a,b) = \begin{cases} 1 & (\text{if: } a \le_{\min} b) \\ 0 & (\text{otherwise}) \end{cases}

まずa=n a = nのとき。 H[b]=bxnF[x]1=bxnF[x] \begin{aligned} H[b] &= \sum_{b \le x \le n}F[x]*1 \\ &= \sum_{b \le x \le n}F[x] \end{aligned}

つづいてb=1 b = 1のとき。 H[a]=1xa1G[x]=1xaG[x] \begin{aligned} H'[a] &= \sum_{1 \le x \le a}1*G[x] \\ &= \sum_{1 \le x \le a}G[x] \end{aligned}

隣接代数のメビウス関数

組み合わせ論でのメビウス反転公式を考える。 前述の(b=1 b = 1のときを例に)ゼータ変換の逆変換を以下のように記述する。

G[a]=1xaμ(x,a)H[x] \begin{aligned} G[a] &= \sum_{1 \le x \le a}\mu(x,a)H'[x] \end{aligned}

このμ(a,b) \mu(a,b)をメビウス関数と呼び以下を満たす

x:aminxminbμ(a,x)=δ(x,b) \sum_{x: a \le_{\min} x \le_{\min} b}\mu(a, x) = \delta(x, b)

式変形して陽形式にできるっぽいけど古典的な定義じゃなく隣接代数の一般形式はかなり複雑になる。 ただプログラムを書く分には素朴にゼータ変換を戻すように書けばいいので同じ添字について逆算すればいい。 インプレースに実装したい場合はループを逆向きにする必要はある。

典型的な畳み込み

ここまで一般化された畳み込みを説明し半順序を前提に一般化されたゼータ変換と逆変換のメビウス変換を確認した。 ここからは畳み込みの具体例をいくつか確認する。

Dirichlet convolution

ディリクレ畳み込みは約数に関してたたみ込んでいる。意味づけは後述する。

(fg)(n1,n2,)=d1n1,d2n2,f(d1,d2,)g(n1/d1,n2/d2,) (f*g)(n_1, n_2,\cdots) = \sum_{d_1|n_1, d_2|n_2, \cdots}f(d_1, d_2, \cdots)g(n_1/d_1,n_2/d_2,\cdots)

1変数ではした (fg)(n)=dnf(d)g(n/d) (f*g)(n) = \sum_{d|n}f(d)g(n/d)

Unitary convolution

ユニタリー畳み込みはユニタリーな約数に関してたたみ込んでいる。

(fg)(n1,n2,)=d1n1,d2n2,f(d1,d2,)g(n1/d1,n2/d2,) (f*g)(n_1, n_2,\cdots) = \sum_{d_1||n_1, d_2||n_2, \cdots}f(d_1, d_2, \cdots)g(n_1/d_1,n_2/d_2,\cdots)

1変数ではした (fg)(n)=dnf(d)g(n/d) (f*g)(n) = \sum_{d||n}f(d)g(n/d)

ユニタリー除数 d とは dn and gcd(d,n/d)=1 d|n \text{ and } \gcd(d, n/d) = 1を満たす。

他にも Gcd,Lcm,Binomial convolution などがある。

Dirichlet級数との対応づけ

古典的な Dirichlet 畳み込み積とDirichlet級数との対応づけを追っていく。

Dirichlet級数: DA(s)=A11s+A22s+A33s+ \text{Dirichlet級数: } D_A(s) = \frac{A_1}{1^s} + \frac{A_2}{2^s} + \frac{A_3}{3^s} + \cdots

数列を A の値A[i] A[i]がDirichlet級数の係数Ai A_iに対応すると定義する。

こうすると対応する畳み込み積は下のようになる。

C[k]=ij=kA[i]B[j] C[k] = \sum_{i*j=k} A[i]B[j]

ゼータ関数

古典的なリーマンのゼータ関数は下の通り。

Riemann Zeta関数: ζ(s)=11s+12s+13s+ \text{Riemann Zeta関数: } \zeta(s) = \frac{1}{1^s} + \frac{1}{2^s} + \frac{1}{3^s} + \cdots

このゼータ関数と数列のDirichlet級数の積に対応した数列を求める変換を数列のゼータ変換と呼ぶ。 ゼータ関数の系数列をζ \zetaとするとC=ζA C = \zeta*Aとなる。

添字が半順序を持つ半群であるため隣接代数に基づいたゼータ変換だが、もともとはこっちが先にあったっぽい(要出典)。

C[k]=ikA[i] C[k] = \sum_{i|k} A[i]

ここで dn d|n は「 d dn n の約数である」と読む。

ゼータ関数はEuler積と呼ばれる標式でも書き表せます。

ζ(s)=p:primes(1+1ps+1p2s+)=p:primes(1ps)1 \begin{aligned} \zeta(s) &= \prod_{p:\text{primes}}(1+\frac{1}{p^s}+\frac{1}{p^{2s}}+ \cdots) \\ &= \prod_{p:\text{primes}}(1-p^{-s})^{-1} \\ \end{aligned}

気持ちとしては、全ての自然数は全ての素数の無限等比級数を掛け合わせれば取り出せるってこと。 このオイラー積をメビウス関数を求めるのに使う。

単位元

単位元: I=11s+02s+03s+=1 \text{単位元: } I = \frac{1}{1^s} + \frac{0}{2^s} + \frac{0}{3^s} + \cdots = 1

メビウス関数

ζ(s) \zeta(s)の逆元を考える。

1ζ(s)=μ(1)1s+μ(2)2s+μ(3)3s+ \frac{1}{\zeta (s)} = \frac{\mu(1)}{1^s} + \frac{\mu(2)}{2^s} + \frac{\mu(3)}{3^s} + \cdots

この係数μ(n) \mu(n)をメビウス関数と呼びます。 定義からζμ=I \zeta*\mu=Iとなる。

ゼータ関数の逆数の系数列がメビウス関数となる。ゼータ関数をオイラー積で表すと扱いやすくなる。

1ζ(s)=p:primes(1ps)=+μ(n)ns+ \begin{aligned} \frac{1}{\zeta(s)} &= \prod_{p:\text{primes}}(1-p^{-s}) \\ &= \cdots + \frac{\mu(n)}{n^s} + \cdots \end{aligned}

よってn=p:primespe n=\prod_{p:\text{primes}}p^eの時メビウス関数μ(n) \mu(n)

μ(n)={1(e=0)1(e=1)0(e2) \mu(n) = \begin{cases} 1 & (e = 0) \\ -1 & (e = 1) \\ 0 & (e \ge 2) \end{cases}

となる。

オイラー関数

べき乗関数

べき乗関数をNk=[1k,2k,3k,] N^k = [1^k, 2^k, 3^k, \cdots]と定義する。 べき乗関数に対応するDirichlet級数ははゼータ関数を用いてζ(sk) \zeta(s-k)と書き表わせる。

DNk=11sk+12sk+13sk+=ζ(sk) \begin{aligned} D_{N^k} &= \frac{1}{1^{s-k}} + \frac{1}{2^{s-k}} + \frac{1}{3^{s-k}} + \cdots \\ &= \zeta(s-k) \end{aligned}

約数和関数

約数和関数をσk(n)=dndk \sigma^k(n) = \sum_{d|n}d^kで定義する。 特にσ0 \sigma^0は約数の個数となる。 定義からσk=ζNk \sigma^k = \zeta*N^kを満たす。 またDσk=ζ(s)ζ(sk) D_{\sigma^k} = \zeta(s)\zeta(s-k)の級数は

Euler関数

Euler関数(φ(n) \varphi(n))は1以上n以下でnと互いに素な整数の個数。 n との最大公約数が d なる n/d と互いに素なので n 以下の数の個数は φ(n/d) \varphi(n/d) となる。

n=dnφ(n/d)=dnφ(d)N1=φζDφ=ζ(s1)ζ(s) \begin{aligned} n &= \sum_{d|n}\varphi(n/d) \\ &= \sum_{d|n}\varphi(d) \\ N^1 &= \varphi \ast \zeta \\ D_{\varphi} &= \frac{\zeta(s-1)}{\zeta(s)} \end{aligned}

となる。

乗法的関数

n, m が互いに素であるなら。乗法的関数は以下を満たす。 全ての n,m について満たす場合は完全乗法的関数と呼ぶ。

f(nm)=f(n)f(m)f(n1m1,n2m2,)=f(n1,n2,)f(m1,m2,) \begin{aligned} f(nm) &= f(n)*f(m) \\ f(n_1 m_1, n_2 m_2, \cdots) &= f(n_1, n_2, \cdots)*f(m_1, m_2, \cdots) \end{aligned}

畳み込みの定理が成り立つ畳み込み

畳み込み積(Ψ \Psi)が関数の乗法性を壊さない場合には畳み込み定理が成り立つ。 必要条件は確認してない。 具体例にあげたDirichlet畳み込み積、Unitary畳み込み積は乗法性を壊さない。

Ψ(fg)(n)=d1dr=n(fg)(d1,,dr)=d1dr=na1b1=d1,,arbr=drf(a1,,ar)g(b1,,br)=xy=na1ar=xf(a1,,ar)b1br=xg(b1,,br)=xy=nΨ(f)(x)Ψ(g)(y)=(Ψ(f)Ψ(g))(n) \begin{aligned} \Psi(f*g)(n) &= \sum_{d_1 \cdots d_r = n} (f*g)(d_1, \cdots, d_r) \\ &= \sum_{d_1 \cdots d_r = n} \sum_{a_1 b_1 = d_1, \cdots, a_r b_r = d_r} f(a_1, \cdots, a_r) g(b_1, \cdots, b_r) \\ &= \sum_{xy=n} \sum_{a_1 \cdots a_r = x} f(a_1, \cdots, a_r) \sum_{b_1 \cdots b_r = x} g(b_1, \cdots, b_r) \\ &= \sum_{xy=n}\Psi(f)(x)\Psi(g)(y) = (\Psi(f)*\Psi(g))(n) \end{aligned}

隣接代数と関連した畳み込み積の場合、添字が束だと畳み込みの定理を満たす積を定義できるらしい。 そのような畳み込みは、結び( \vee)や交わり( \wedge)を使って以下のように定義する。

C[k]=ij=kA[i]B[j]C[k]=ij=kA[i]B[j] \begin{aligned} C[k] &= \sum_{i \wedge j=k}A[i]B[j] \\ C[k] &= \sum_{i \vee j=k}A[i]B[j] \end{aligned}

高速ゼータ変換

ゼータ変換、メビウス変換を高速に実装したアルゴリズムを指してるっぽい。 日本語なネット解説記事だとDirichlet級数に紐づけた畳み込み積bit表現された集合関数に紐づけた畳み込み積についてが見つかった

畳み込み積に関わる部分が疎なので適切な探索方法を使えば高速になるよねって工夫になっている。 Dirichlet畳み込み積の法はEuler積の表現を使っていて面白かった。

メビウス変換が存在しないとき

半群Gに逆元が存在しない or 総和で加算の代わりに非可逆な演算を採用するといった状況で似た感じの問題を解く場面があると思う。 そんな時は、

  • 逆算を諦める(範囲による問い合わせの回答能力が落ちる)
  • もとの情報を失わないように管理する
    • 配列範囲全体を復元できるデータ構造を使う: セグ木, フェニック木
    • 一部を復元できるだけで十分: しゃくとり法(優先度付きキュー)

という感じで調べるモチベーションが自然にセグ木やフェニック木やしゃくとり法に繋がった。 やった!勉強を続けられる!

こんな寄り道をして鍛錬をサボっている横で息子は積み木を積み上げられるようになるなど圧倒的な成長を見せつけてくれて心強い。

comments powered by Disqus