この記事は夏のブログリレー1日目の記事です

目次

• 目次
• 初めに
• TL;DR
• 用語の定義
• 分類
• 高次元累積和
  • 一次元累積和のイメージ
  • 累積和の方向について
  • 高次元累積和の一次元への圧縮
  • 高次元累積和の一次元への圧縮: 特殊な場合
• 高次元階差
  • 一次元階差
  • 高次元階差
  • 高次元階差: 圧縮済みの形
• 前半のまとめ
• 各オプションに関する解説
  • 演算：ゼータ変換 Z
  • 演算：メビウス変換 M
  • インデックス：整数 I
  • インデックス：bitset B
  • インデックス：素因数分解 F
  • インデックス：タプル T
  • 方向：Prefix P
  • 方向：Suffix S
• 個別の実装と解説
  • PZB：Prefix Zeta transform for Bitset
  • PMB：Prefix Mobius transform for Bitset
  • SZB：Suffix Zeta transform for Bitset
  • SMB：Suffix Mobius transform for Bitset
  • PZF：Prefix Zeta transform for Factorization form
  • PMF：Prefix Mobius transform for Factorization form
  • SZF：Suffix Zeta transform for Factorization form
  • SMF：Suffix Mobius transform for Factorization form
• 添え字が連番でない場合
• 各変換同士の関係
  • PZB ↔ SZB, PMB ↔ SMB
• AND / OR convolution や GCD / LCM convolution との関係
  • AND convolution
  • OR convolution
  • GCD convolution
  • LCM convolution
• 全体を通して
• おわりに
• 参考文献

初めに

競技プログラミング文脈における「約数系包除原理」が、問題によって異なる意味で使われている、あるいは複数のアルゴリズムの総称として使われていて混乱が生じているように感じたので、私的なメモを兼ねて用語と実装を整理しました。メビウス変換、ゼータ変換、約数系包除原理といった概念を統一的に説明することを目的としています。

基本原理から全ての説明を書いたため、初学者でも頑張れば読めるものになったと思います。ゼータ変換/メビウス変換にある程度知識がある人は、後半の「各オプションに関する解説」の節から読むことをオススメします。

TL;DR

• ゼータ変換/メビウス変換は累積和/階差の一般化である
• ゼータ変換が多次元累積和のように解釈できるとき、軸ごとに累積和を取ることで高速に変換できる。これが高速ゼータ変換である。
• 累積和/階差を取る方向によって、PrefixとかSuffixという呼び名を付けた
• ゼータ変換/メビウス変換系の分類を整理した
• 以下の実装について紹介している
  • 集合の包含関係に関するゼータ/メビウス変換
    • $O(3^n)=O(N^{1.58})$ のもの
    • $O(n{2}^n)=O(N\log N)$ のもの
  • 約数倍数に関するゼータ/メビウス変換
    • $O(N\log N)$ のもの
    • $O(N\log \log N)$ のもの
• ゼータ変換/メビウス変換 をすることで LCM/GCD/AND/OR convolution ができる

用語の定義

まず、ゼータ変換とは何かを述べておきます。

ゼータ変換
数列の添え字をまとめた添え字集合 $J$ を考える。このとき、この集合 $J$ で半順序関係 $⪯$ が成立し、かつ関数 $f(x)$ が任意の $x\in J$ で定義されている時、$f$ をゼータ変換した関数 $F$ を次のように定義する。

$$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)$$

半順序関係とは？
ある集合 $S$ について, 次の $3$ つの性質を満たしているとき、関係 $⪯$ は 半順序関係であるという。
反射律 : 任意の $a\in S$ について、$a⪯a$
反対称律 : 任意の $a,b\in S$ について、$a⪯b\wedge b⪯a\leftrightarrow a=b$
推移律 : 任意の $a,b,c\in S$ について、$a⪯b\wedge b⪯c⟹a⪯c$
例 : $\mathbb{N}$(自然数) 上で, 関係 $\le$ は半順序関係である。

これだけではよく分からないので、競技プログラミングでも頻繁に登場する一次元配列 $a=\{a_0,a_1,...,a_{n-1}\}$ の累積和に対応させて考えます。
実は、この配列 $a$ を用いて定義した関数 $f(i)=a_i(i=0,1,...,n-1)$ を適切にゼータ変換すると、$a$ の累積和が得られます。
つまり、この$f$を適切にゼータ変換すると $F(i)=\sum _{k=0}^i{a}_k$ となります。
この適切なゼータ変換とは何かを言い換えていきます。
ゼータ変換の定義は、

数列の添え字をまとめた添え字集合 $J$ を考える。このとき、この集合 $J$ で半順序関係 $⪯$ が成立し、かつ関数 $f(x)$ が任意の $x\in J$ で定義されている時、$f$ をゼータ変換した関数 $F$ を次のように定義する。
$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)$

でした。
一次元配列の場合、添え字集合 $J$ は、数列が定義されている部分である $\{0,1,2,...,n-1\}$ となります。要は、数列や関数の定義域を表します。
次に、この集合 $J$ 上での半順序関係を決めます。半順序というと難しいですが、簡単に言うと、「何が自身より前なのか？」を決めるルールを用意します。$a⪯b$ は、「 $a$ は、$b$ と等しいか、$b$ より前にある」を表します。
半順序関係の定義を言い換えると次のようになります。

ある集合 $S$ について, 次の $3$ つの性質を満たしているとき、関係 $⪯$ は半順序関係であるという。
反射律 : 任意の $a\in S$ について、$a$ は、$a$ 自身と等しいか、$a$ より前にある。
反対称律 : 任意の $a,b\in S$ について、「$a$ が、$b$と等しいか $b$ より前にある」かつ「$b$が、$a$と等しいか $a$ より前にある」であるとき、$a$ と $b$ は等しい。
推移律 : 任意の $a,b,c\in S$ について、「$a$ が、$b$と等しいか $b$ より前にある」かつ「$b$が、$c$と等しいか $c$ より前にある」であるとき、「$a$ は、$c$と等しいか $c$ より前にある」

このルールとして、今回は整数の大小関係$\le$を採用します。これは、反射律も反対称律も推移律もすべて満たします。
つまり今回は、二つの添え字 $i,j$ があるとき、$i⪯j$ は、$i$ が $j$ 以下であることを表します。

これで集合 $J$、半順序関係$\le$、変換前の関数 $f$ がすべてそろったのでゼータ変換をすることができます。
変換の結果 $F$ は、次のようになります。
$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)=\sum _{i=0}^{n}f(i)$
これは、通常の累積和です。

このように、ゼータ変換は一次元配列の累積和の拡張となっていることが分かります。

また、このゼータ変換には逆変換が存在し、それを メビウス変換 と呼びます。
ゼータ変換の例として、「一次元配列の累積和の作成」がありましたが、その逆のメビウス変換の例は、「一次元配列の階差数列の作成」となります。

以上より、
ゼータ変換は累積和の一般化
メビウス変換は階差の一般化
ということが分かりました。

分類

競技プログラミングにおいては、配列に対する変換として、主に以下のようなオプションがあります。

• 和を取るか差分を取るか
  • 和: ゼータ変換
  • 差: メビウス変換
• インデックスが何か
  • 添え字集合は何か、また、それ上での半順序関係は何か
• 和や差をとる方向はどちらか
  • Prefix: 添え字の小さい方を前とする
  • Suffix: 添え字の大きい方を前とする

「和や差の方向」は後で説明します。インデックスが何か？が割と面倒な存在です。
インデックスを特徴づけるのは、添え字集合とそれ上の順序なので、それらを表にして整理します。
よくあるものには具体例を書いておきました。ほとんど扱われなかったり、そもそも定義するのが難しかったりするものは「???」にしています。

色々ありますが、ベースとなる考え方は、高次元累積和です。次の章で詳しく扱います。

高次元累積和

高次元累積和とは、次のようなゼータ変換の結果の一種です。例として三次元を採用していますが、何次元でも使えます。

三次元配列 $a_{i,j,k}(0\le i<a,0\le j<b,0\le k<c)$ がある。
これを用いて、添え字集合 $J=\{(i,j,k)|0\le i<a,0\le j<b,0\le k<c\}$ 上で定義された関数 $f((i,j,k))=a_{i,j,k}$ を用意する。
この添え字集合の順序 $(i_1,j_1,k_1)⪯(i_2,j_2,k_2)$ は、$i_1\le i_2\wedge j_1\le j_2\wedge k_1\le k_2$ と定める。
これらを用いて次のように定義した関数 $F$ が、$f$ のゼータ変換である。
$F((i,j,k))=\sum _{\begin{array}{c}(x,y,z)\in J\\ (x,y,z)⪯(i,j,k)\end{array}}f((x,y,z))=\sum _{x=0}^i\sum _{y=0}^j\sum _{z=0}^{k}f((x,y,z))$
つまり、$F((i,j,k))$ は、原点 $(0,0,0)$ と点 $(i,j,k)$ を対角とする直方体内の添え字に関して $f$ を全部足したものとなります。

一次元累積和のイメージ

一次元累積和は、前から順に伸ばしていくことで簡単に構築できます。
高次元の場合も、各次元の方向について一つずつ伸ばしていくことで構成できます。

一次元累積和のコード(元の配列をそのまま使う場合)

N = 6
A = [1, 3, 5, 7, 9, 11]
for i in range(1, N):
    A[i] += A[i-1]
print(A) #[1, 4, 9, 16, 25, 36] と出力される。

コード2(新しく配列を作る場合)

N = 6
A = [1, 3, 5, 7, 9, 11]
A_cumulative = [0]*N
for i in range(0, N):
    if i == 0:
        A_cumulative[i] = A[i]
    else:
        A_cumulative[i] = A_cumulative[i-1]+A[i]
print(A_cumulative) #[1, 4, 9, 16, 25, 36] と出力される。

競技プログラミング的には実装が楽なことが多いので、元の配列をそのまま累積和に変換することが多いです。元の配列を残したい場合も、コピーを作ってそれを変換すればよいです。

高(3)次元累積和のイメージ

初期状態

各地点は、まだその地点の値だけを表します。

一次元目の累積後

各地点は、たった今累積した方向の一次元累積和に変化します。

二次元目の累積後

各地点は、今まで累積した2つの方向が作る平面内での二次元累積和に変化します。

三次元目の累積後(完成形)

各地点は、今まで累積した3つの方向が作る空間内での三次元累積和に変化します。これで目標の三次元累積和ができました。

このように、各次元方向について一回ずつ累積する手法のことを、 高速ゼータ変換 といいます。
この高速ゼータ変換にかかる時間計算量は、次元数×値の個数となります。今回の3次元累積和なら$3\times 27=81$となる。

三次元累積和のコード(元の配列をそのまま累積和に変換)

A = [
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ]
    ]

for d in range(3): #次元方向のループ
    prev = [0,0,0]
    prev[d] = 1 # 対応する方向の「前」がどちらであるかを決める。
    
    for x in range(3): #各方向について、累積を行う
        for y in range(3):
            for z in range(3):
                if x-prev[0] >= 0 and y-prev[1] >= 0 and z-prev[2] >= 0:
                    A[x][y][z] += A[x-prev[0]][y-prev[1]][z-prev[2]]

print(A)
# 出力結果 (最初、全部の数値が1なので、原点からその地点までの直方体の体積が結果になる)
"""
[
    [
        [1, 2, 3],
        [2, 4, 6],
        [3, 6, 9]
    ],
    [
        [2, 4, 6],
        [4, 8, 12],
        [6, 12, 18]
    ],
    [
        [3, 6, 9],
        [6, 12, 18],
        [9, 18, 27]
    ]
]
"""

累積する部分はもっとごり押しでも可能です。

#x方向
for x in range(1, 3):
    for y in range(3):
        for z in range(3):
            A[x][y][z] += A[x-1][y][z]
#y方向
for x in range(3):
    for y in range(1, 3):
        for z in range(3):
            A[x][y][z] += A[x][y-1][z]
#z方向
for x in range(3):
    for y in range(3):
        for z in range(1, 3):
            A[x][y][z] += A[x][y][z-1]

累積和の方向について

先ほど扱った三次元累積和は、「原点」と「ある地点」を対角とする直方体内の和でしたが、これを「ある地点」と「終点$(a-1,b-1,c-1)$」を対角とする直方体内の和などと定義することもできます。
これは、添え字の順序を、その双対順序に変更した場合であると考えられます。双対順序とは、大小関係を逆転させることで得られる順序のことです。

つまり、$(i_1,j_1,k_1)⪯(i_2,j_2,k_2)$ という順序関係を、

• (i) ... $i_1\le i_2\wedge j_1\le j_2\wedge k_1\le k_2$ と定義した場合、ある地点と原点の直方体
• (ii) ... $i_1\ge i_2\wedge j_1\ge j_2\wedge k_1\ge k_2$ と定義した場合、ある地点と終点の直方体

の中の添え字に関して和を取ることになります。

この記事では、数値としての大小関係に基づき、
(i) の方を Prefix ゼータ変換
(ii) の方を Suffix ゼータ変換
と呼ぶことにします。

一般の順序関係に対しては、

• Prefix ゼータ変換は添え字が数値として小さい方に関して和を取るもの、
• Suffix ゼータ変換は添え字が数値として大きい方に関して和を取るもの、

と区別されます。他の具体例としては、

• 「2次元配列においてある場所と原点を対角とする長方形領域内の和を取る」場合、これは和を取る対象となる添え字のインデックスは、どちらの成分ももとの点より小さいので、Prefix である
• 「1次元配列においてあるインデックスの倍数となるようなインデックス全てに対しての和を取る」場合、倍数はもとの数より数値として大きいので、Suffix である

などのように定めます。

Suffix三次元累積和のコード(元の配列をそのまま累積和に変換)

A = [
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ]
    ]

for d in range(3): #次元方向のループ
    prev = [0,0,0]
    prev[d] = -1 # 対応する方向の「前」がどちらであるかを決める。今回は終点(2,2,2)の方が前なので、向きが反転する。
    
    for x in range(2,-1,-1): #各方向について、累積を行う。ループが逆順であることに注意
        for y in range(2,-1,-1):
            for z in range(2,-1,-1):
                if x-prev[0] < 3 and y-prev[1] < 3 and z-prev[2] < 3:
                    A[x][y][z] += A[x-prev[0]][y-prev[1]][z-prev[2]]

print(A)
# 出力結果 (最初、全部の数値が1なので、終点からその地点までの直方体の体積が結果になる)
"""
[
    [
        [27, 18, 9],
        [18, 12, 6],
        [9, 6, 3]
    ],
    [
        [18, 12, 6],
        [12, 8, 4],
        [6, 4, 2]
    ],
    [
        [9, 6, 3],
        [6, 4, 2],
        [3, 2, 1]
    ]
]
"""

高次元累積和の一次元への圧縮

先ほどのように3次元なら、工夫してループを書けばそこまで面倒ではないですが、これが $d$ 次元となるとそうはいきません。
そこで、複数の添え字を1つの整数にまとめます。

範囲が指定された $d$ 個の添え字 $i_0,i_1,...,i_{d-1}(0\le i_0<l_0,0\le i_1<l_1,...,0\le i_{d-1}<l_{d-1})$ があるとします。この添え字は、次のような形で一つの整数 $D$ に圧縮することができます。

$$D=i_0+i_1{l}_0+i_2{l}_0{l}_1+i_3{l}_0{l}_1{l}_2+...+i_{d-1}{l}_0{l}_1...l_{d-2}=\sum _{k=0}^{d-1}{i}_k\prod _{r=0}^{k-1}{l}_r$$

この圧縮法には二つの利点があります。
一つ目は、可逆であることです。

$$i_k=\lfloor \frac{D\%(l_0{l}_1...l_k)}{l_0{l}_1...l_{k-1}}\rfloor$$

とすれば簡単に復元できます。(x%yは、xをyで割った余りを0以上|y|未満で表したものです。)

二つ目は、無駄がないことです。
整数 $D$ は、最小で $0$、最大で $l_0{l}_1{l}_2...l_{d-1}-1$ になります。また、$d$ 個の添え字の組の総数は$l_0{l}_1{l}_2...l_{d-1}$ です。これは $D$ の取りうる種類数と一致します。変換が可逆であることより、元の添え字と整数 $D$ との間には1対1の対応が取れていて、まったく無駄がないです。
さらに、あり得る $D$ の値は連番の非負整数になるので、これは一次元配列に載せることができます。

これを利用して、Prefix三次元累積和を構築するコードを書き直すと次のようになります。
Prefix三次元累積和は、前から順に計算する必要がありますが、添え字の組 $a,b$ が $a⪯b$ であるとき、それぞれを整数$D_a,D_b$に変換した時には必ず $D_a\le D_b$ が成り立つので、単に $D$ の昇順に計算すればよいです。

L = [3,3,3] #各次元方向の長さ
#↑これを3次元以外に変えたり、数値を変えたりしても動きます。

N = 1
L_cum_prod = [] #Lの累積積を取っておく
for i in L:
    N *= i
    L_cum_prod.append(N)


A = [1]*N #圧縮された三次元配列を用意

#整数Dから添え字i_kを求める関数
def get_i_k(D, k):
    return D%L_cum_prod[k]//(1 if k == 0 else L_cum_prod[k-1])

#整数Dがあるとき、添え字i_kを1だけ減らした時のDを求める関数
def decrement_i_k(D, k):
    return D-(1 if k == 0 else L_cum_prod[k-1])


for d in range(len(L)): #次元の数だけループして累積する。
    for i in range(N): #全要素をDの昇順に見る
        if get_i_k(i, d) == 0: #端は飛ばす
            continue
        A[i] += A[decrement_i_k(i, d)]

print(A)
#[1, 2, 3, 2, 4, 6, 3, 6, 9, 2, 4, 6, 4, 8, 12, 6, 12, 18, 3, 6, 9, 6, 12, 18, 9, 18, 27]

高次元累積和の一次元への圧縮: 特殊な場合

基本的には先ほどのように圧縮をすることで簡潔に累積和を構築できますが、各添え字の取りうる範囲がすべて一定、すなわち、以下のような場合は、別の捉え方をすることができます。

全添え字で共通の範囲が指定された $d$ 個の添え字 $i_0,i_1,...,i_{d-1}(0\le i_0,i_1,...,i_{d-1}<l)$ がある。

このとき、先ほどのように圧縮をすると次のようになります。

$$D=i_0+i_{1}l+i_2{l}^2+i_3{l}^3+...+i_{d-1}{l}^{d-1}=\sum _{k=0}^{d-1}{i}_k{l}^k$$

これは、$l$ 進法を $10$ 進法に直していることに他なりません。
復元も、基数の変換と同じ形になります。

$$i_k=\lfloor \frac{D\%(l^{k+1})}{l^k}\rfloor$$

先ほどとは異なり、圧縮の逆変換に累乗と割り算の計算しか使わないので、累乗が高速にできれば累積積などの構築を省くことができます。
後で扱いますが、$l=2$ がかなり頻出で、累乗もビットシフトで非常に高速に計算できるので、このようにして逆変換を行うことが多いです。(除算は遅いので、実際はビット演算で書くことがほとんどです。)

$l=2$ の場合のコード
各次元方向の長さが2に固定されただけで、それ以外は先ほどの三次元累積和とそれほど実装は変わりません。

L = 3 #次元数の指定(3以外でも動きます)

N = 1<<L

A = [1]*N #圧縮された配列を用意

for d in range(L):
    for i in range(N):
        if i&(1<<d) == 0: #iの二進法表記の2^dの位が0かどうかを判定
            continue
        A[i] += A[i-(1<<d)]

print(A)
#[1, 2, 2, 4, 2, 4, 4, 8]

高次元階差

今までは累積和の作成、つまり、ゼータ変換に着目してきましたが、ほぼ同じようにしてその逆変換であるメビウス変換も実装することができます。

また、ゼータ変換にPrefix、Suffixと名前を付けたように、

• Prefixゼータ変換の逆変換のことをPrefixメビウス変換
• Suffixゼータ変換の逆変換のことをSuffixメビウス変換

と名付けておきます。

一次元階差

まず、一番簡単な例として、まずは一次元配列のPrefixゼータ変換(いわゆる普通の累積和)をPrefixメビウス変換し、元の配列を復元することを考えます。
元の配列を直接操作し、その中で完結させる方法と、新しく配列を作る方法がありますが、通常は元の配列内で完結させることが多いので、そちらの実装をします。

累積和配列があるとき、階差、すなわち、ある項とその一個前との差を取れば元の項を復元できます。しかし、Prefixゼータ変換のように添え字の昇順に行ってしまうと、階差を取りたいのに差分を取りたい項がすでに破壊されてしまっているということが起こります。
よって、添え字の降順で行います。

コード

N = 6
A = [1, 4, 9, 16, 25, 36]
for i in range(N-1, -1, 0): #降順
    A[i] -= A[i-1] #和ではなく、差を取る

print(A) #[1, 3, 5, 7, 9, 11]

Prefix累積和の復元は降順ですが、Suffix累積和の復元は逆に昇順になります。

N = 6
A = [36, 25, 16, 9, 4, 1] #Suffix累積和
for i in range(0, N-1): #(右端以外で)昇順
    A[i] -= A[i+1] #和ではなく、差を取る。Suffix累積和は、「前」は、添え字の大きい方を表すということに注意する。

print(A) #[11, 9, 7, 5, 3, 1]

高次元階差

次の例として、三次元累積和を取り上げます。
このコードで作られたPrefix三次元累積和を復元することを考えます。

A = [
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ],
        [
            [1,1,1],
            [1,1,1],
            [1,1,1]
        ]
    ]

for d in range(3): #次元方向のループ
    prev = [0,0,0]
    prev[d] = 1 # 対応する方向の「前」がどちらであるかを決める。
    
    for x in range(3): #各方向について、累積を行う
        for y in range(3):
            for z in range(3):
                if x-prev[0] >= 0 and y-prev[1] >= 0 and z-prev[2] >= 0:
                    A[x][y][z] += A[x-prev[0]][y-prev[1]][z-prev[2]]

print(A)
# 出力結果 (最初、全部の数値が1なので、原点からその地点までの直方体の体積が結果になる)
"""
[
    [
        [1, 2, 3],
        [2, 4, 6],
        [3, 6, 9]
    ],
    [
        [2, 4, 6],
        [4, 8, 12],
        [6, 12, 18]
    ],
    [
        [3, 6, 9],
        [6, 12, 18],
        [9, 18, 27]
    ]
]
"""

復元を行うコード
階差を取る場合は降順に処理する必要があることに注意してください。

A = [
    	[
    	    [1, 2, 3],
    	    [2, 4, 6],
    	    [3, 6, 9]
    	],
    	[
    	    [2, 4, 6],
    	    [4, 8, 12],
    	    [6, 12, 18]
    	],
        [
            [3, 6, 9],
            [6, 12, 18],
            [9, 18, 27]
        ]
    ]

for d in range(3):
    prev = [0,0,0]
    prev[d] = 1
    for x in range(2, -1, -1): #降順であることに注意する。
        for y in range(2, -1 ,-1):
            for z in range(2 ,-1 ,-1):
                if x-prev[0] >= 0 and y-prev[1] >= 0 and z-prev[2] >= 0:
                    A[x][y][z] -= A[x-prev[0]][y-prev[1]][z-prev[2]]

print(A)
"""
[
    [
        [1,1,1],
        [1,1,1],
        [1,1,1]
    ],
    [
        [1,1,1],
        [1,1,1],
        [1,1,1]
    ],
    [
        [1,1,1],
        [1,1,1],
        [1,1,1]
    ]
]
"""

これも、Suffix累積和の復元時は昇順に処理する必要があることに注意してください。

高次元階差: 圧縮済みの形

最後に、一次元圧縮された高次元累積和の復元について考えます。

まず、以下のコードで作成された一次元圧縮済みの高次元累積和を用意します。

L = [3,3,3] #各次元方向の長さ
#↑これを3次元以外に変えたり、数値を変えたりしても動きます。

N = 1
L_cum_prod = [] #Lの累積積を取っておく
for i in L:
    N *= i
    L_cum_prod.append(N)


A = [1]*N #圧縮された三次元配列を用意

#整数Dから添え字i_kを求める関数
def get_i_k(D, k):
    return D%L_cum_prod[k]//(1 if k == 0 else L_cum_prod[k-1])

#整数Dがあるとき、添え字i_kを1だけ減らした時のDを求める関数
def decrement_i_k(D, k):
    return D-(1 if k == 0 else L_cum_prod[k-1])


for d in range(len(L)): #次元の数だけループして累積する。
    for i in range(N): #全要素をDの昇順に見る
        if get_i_k(i, d) == 0: #端は飛ばす
            continue
        A[i] += A[decrement_i_k(i, d)]

print(A)
#[1, 2, 3, 2, 4, 6, 3, 6, 9, 2, 4, 6, 4, 8, 12, 6, 12, 18, 3, 6, 9, 6, 12, 18, 9, 18, 27]

これを復元するコードは次のようになります。

L = [3,3,3] #各次元方向の長さ
#↑これを3次元以外に変えたり、数値を変えたりしても動きます。ただし、累積和を作成した時と同じにしてください。

N = 1
L_cum_prod = [] #Lの累積積を取っておく
for i in L:
    N *= i
    L_cum_prod.append(N)

#累積和を用意
A = [1, 2, 3, 2, 4, 6, 3, 6, 9, 2, 4, 6, 4, 8, 12, 6, 12, 18, 3, 6, 9, 6, 12, 18, 9, 18, 27]

a = [0]*N
for i in range(N):
    a[i] = A[i]

#整数Dから添え字i_kを求める関数
def get_i_k(D, k):
    return D%L_cum_prod[k]//(1 if k == 0 else L_cum_prod[k-1])

#整数Dがあるとき、添え字i_kを1だけ減らした時のDを求める関数
def decrement_i_k(D, k):
    return D-(1 if k == 0 else L_cum_prod[k-1])


for d in range(len(L)): #次元の数だけループして階差を取る。
    for i in range(N-1, -1, -1): #全要素をDの降順に見る
        if get_i_k(i, d) == 0: #端は飛ばす
            continue
        a[i] -= a[decrement_i_k(i, d)]

print(a)
#[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

各次元方向の長さが一定であり、さらにそれが $2$ であるときに特殊化した実装は次のようになります。

L = 3 #次元数の指定(3以外でも動きます)

N = 1<<L

#累積和を用意
A = [1, 2, 2, 4, 2, 4, 4, 8]

a = [0]*N
for i in range(N):
    a[i] = A[i]

for d in range(L):
    for i in range(N-1, -1, -1): #降順
        if i&(1<<d) == 0: #iの二進法表記の2^dの位が0かどうかを判定
            continue
        a[i] -= a[i-(1<<d)]

print(a) #[1, 1, 1, 1, 1, 1, 1, 1]

前半のまとめ

ここまでのまとめです。

• ゼータ変換は累積和の一般化である。
• メビウス変換はゼータ変換の逆変換である。
• 適当な順序を取ったときの和を取る対象が、小さい方と定義するのがPrefix、大きい方と定義するのがSuffixである。
• 添え字集合やその順序によって様々な変換が考えられるが、高次元累積和/階差がベースとなっている。

ここからは様々な変換がどのようにして高次元累積和/階差に対応するのかを見ていきます。

各オプションに関する解説

最初に行った分類により、主に次のようなオプションがあることが分かりました。

• 「ゼータ変換(累積和の作成)」か「メビウス変換(階差数列の作成)」か？
• 添え字集合は何か？また、添え字集合での順序はどうなっているか？

これらのオプションで分類した先にも詳細を決める追加オプションがあります。

• 前はどちらか？
  • 添え字が小さい方が前ならPrefix
  • 添え字が大きい方が前ならSuffix
• 値を格納する形式は？
  • 連番の整数なら通常の一次元配列
  • 連番でない整数なら連想配列など
  • 整数の組であれば一次元圧縮して一次元配列を利用

これらのオプションのうち、よく使うものを紹介していきます。

オプションを分類するにあたって、アルファベット1文字を割り当てておきます。
イメージとしてはMBTIの 外向的(E)/内向的(I) と似た感じです。

• 前はどちらか？
  • 添え字が小さい方が前 ... $\text{P}$refix
  • 添え字が大きい方が前 ... $\text{S}$uffix
• 「ゼータ変換(累積和の作成)」か「メビウス変換(階差数列の作成)」か？
  • ゼータ変換 ... $\text{Z}$eta Transform
  • メビウス変換 ... $\text{M}$obius Transform
• 添え字集合は何か？また、添え字集合での順序はどうなっているか？
  • 整数の大小(一次元配列) ... $\text{I}$nteger
  • 二進整数のビットの包含関係 ... $\text{B}$itset
  • 約数、倍数 ... $\text{F}$actorization
  • 整数の組の大小(多次元配列かその一次元圧縮版) ... $\text{T}$uple
• 値を格納する形式
  • 上3つのオプションを決めれば自動で決定します。特に文字は割り当てません。

この文字を用いて、各パターンを (P|S)(Z|M)(I|B|F|T) というフォーマットで表すことにします。

それぞれのオプションについて、解説します

演算：ゼータ変換 Z

繰り返しになるので簡潔に述べます。
ある添え字集合 $J$ に対して定義された関数 $f(n)$ に対して、それから定まる新たな関数 $F(n)$ 、

$$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)$$

を求めることをゼータ変換といいます。つまり、「$i⪯n$ となる全ての $i\in J$ に対する $f(i)$ の総和」です。疑似コードで書くと以下のようになります。

N = 8
f = [1,1,1,1,1,1,1,1]

F = [0]*N

for a in J:
  for b in J:
    if a ⪯ b:
      F[b] += f[a]

print(F)

もちろんこのような愚直解は $O(N^2)$ ですが、実装上はいい感じの順番によりすることでこれよりも早く行うことができます。

演算：メビウス変換 M

ある添え字集合 $J$ に対して定義された関数 $F(n)$ に対して、ゼータ変換をすることでその関数になるような関数 $f(n)$ 、

$$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)$$

を求めることをメビウス変換といいます。これは、$f(n)$ に関して陽に書くと、

$$f(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}\mu (i,n)F(i)$$

ただしここで $\mu (i,n)$ は

$$\sum _{\begin{array}{c}i\in J\\ m⪯i⪯n\end{array}}\mu (i,n)=\delta (m,n)$$

を満たす関数です。この式から読み取れることは、「$f(n)$ は $i⪯n$ なる $i$ に対する $F(i)$ の線形結合で求めることができる」ということです。まあそもそもゼータ変換が、$N$ 次元ベクトル $f$ から $N$ 次元ベクトル $F$ への線形変換なので、逆向きも明らかに線形変換になります。

愚直に考えると、逆行列の計算を考えると $O(N^3)$、$J$ のトポロジカルソートを考えると以下のように $O(N^2)$ で計算できます。

疑似コード

N = 8
F = [1,2,2,4,2,4,4,8]

f = [0]*N

assert J はトポロジカルソート順になっている

for a in J:
  f[a] = F[a]  # とりあえず f(a) を項として含むものを取ってくる
  for b in J:
    if b ⪯ a and b != a:
      f[a] -= f[b]  # 余計な項を取り除く、トポソされているので f[b] は既に正しい値が入っている

print(f)

よくある設定の場合、$J$ のトポロジカルソートは解析的に容易に得られます。

インデックス：整数 I

普通の配列に、普通に $0\sim N-1$ を添え字付けて、半順序として $\le$ を取ったものです。つまり、一般的な1次元の累積和/階差の場合です。
そのため、個別の実装は省略します。

インデックス：bitset B

添え字集合 $J$ が、$n$要素からなる集合 $\{0,1,2,...,n-1\}$ の部分集合をすべて集めた集合(冪集合)の場合です。また集合の内包関係 $\subset$ を半順序とします。

添え字が集合だと分かりづらいので、添え字を$n$次元超立方体の頂点に対応させます。

添え字集合内の要素は、$0$があるかないか、$1$があるかないか、...、$N-1$があるかないかによって、$2^N$通りに分類できます。この$n$個の独立な選択を、$n$ 次元空間での座標の独立な選択に対応させて考えます。

3次元の場合を考えると分かりやすいです。
3次元空間の$8$頂点$\{(x,y,z)|x,y,z\in \{0,1\}\}$からなる立方体を用意します。
次のルールで部分集合を点に対応させます。

• $0$があるかないか ... $x$座標が$1$か$0$か
• $1$があるかないか ... $y$座標が$1$か$0$か
• $2$があるかないか ... $z$座標が$1$か$0$か

このルールで、各部分集合が立方体の各点に1対1で対応します。
これを$n$次元でも同様に行うと、各部分集合は$0,1$からなる$n$次元の座標に1対1で対応します。

次に、集合の包含関係について考えます。
$J$ 内の要素を適当に $2$ つ選んで $a,b$ とします。
目標は、$a$ が $b$ の部分集合であるときに $a⪯b$ となるような順序関係を定義することです。
$a$ が $b$ の部分集合であることは、点の座標で考えると次のように表現できます。

$a$ に対応する点が、「 $b$ に対応する点と原点 $(0,0,...,0)$ を対角とする超直方体」の中(境界含む)にある
または、次のように解釈することもできます。
$b$ に対応する点が、「 $a$ に対応する点と終点 $(1,1,...,1)$ を対角とする超直方体」の中(境界含む)にある

したがって、これは$n$次元累積和のときの順序関係の考え方と同じであることが分かります。
部分集合が対応する座標は、0以上1以下の値をとるインデックスが$n$個あるとみなせるので、これは添え字の一次元圧縮を利用することで長さ $2^n$ の配列で管理することができます。

インデックス：素因数分解 F

添え字集合を適当な自然数の集合、半順序として $a⪯b⟺(aはbの約数)$ を考えます。GCDやLCM系の問題でよく使います。

また、添え字集合に採用する自然数の集合は、

• $n\in J\wedge d⪯n⟹d\in J$

となるように取ることを前提とします。つまり、ある数が入っているならば、その約数も全て入っている必要があるということです。

このような集合の選び方として、

• $1\sim N$ の集合
• ある数 $N$ の約数全体の集合
• ある数 $A,B$ の約数全体の集合の和集合

などが挙げられます。特に、連番の場合かそうでないか、で実装が大きく変わります。実装に関しては後述します。

また、素因数分解では、素因数の個数が複数持てる ↔ 多重集合に対するゼータ変換に対応する、ということもできます。ゲーデル数的な考えであって、あくまで理論上の話ですが。

インデックス：タプル T

添字集合をタプルとする場合、つまり arr[(a,b,c)] みたいな状況を考えますが、これは事実上 arr[a][b][c] のような多次元配列を対応します。つまりタプルを添え字とする場合は、高次元累積和に対応します。
そのため、個別の実装は省略します。

方向：Prefix P

$i⪯n$ (半順序)であるときに、$i\le n$ (値としての大小関係)となるような半順序に対する変換を Prefix 的であると言います。

例えば、ゼータ変換の式、

$$F(n)=\sum _{\begin{array}{c}i\in J\\ i⪯n\end{array}}f(i)$$

に対して、$d⪯n⟺(dはnの約数)$ を考えると、これは明らかに $d\le n$ であるので Prefix ゼータ変換になります。

競プロで出てくるような半順序の場合、添え字が数値として大きかどうかを考えることでそれが Prefix 的であるか Suffix 的であるかが自然に定まります。

方向：Suffix S

Prefix の逆です。

$i⪯n$ であるときに、$n\le i$ となるような半順序に対する変換を Suffix 的であると言います。

また、半順序の向きを取り替えることで、半順序の向きと値の大小を一致させた場合、

$$F(n)=\sum _{\begin{array}{c}i\in J\\ n⪯i\end{array}}f(i)$$

といった表記にもできます。この方が Suffix っぽいですね。

個別の実装と解説

これらの全てのパターンに対して、その実装方法について見ていきます。

PZB：$\text{P}$refix $\text{Z}$eta transform for $\text{B}$itset

別名：sum over subsets, subset zeta transform
一番基礎的なゼータ変換です。ゼータ変換とググると大体上位に出てきます。
以前に @oxojo さんが このブログ にて丁寧に解説しています。

やりたいこと

$n$要素からなる集合 $\{0,1,2,...,n-1\}$ の各部分集合 $T$ について、$T$ を引数に取る関数 $f(T)$ が定義されているとする。
この $f$ から、次のような $F$ を高速に求めたい。
ただし、$F$ は $f$ と同様に、集合を引数にとる関数である。

$$F(S)=\sum _{\begin{array}{c}T\subset \{0,1,...,n-1\}\\ T\subset S\end{array}}f(T)$$

つまり、$F(S)$ は、$S$ の部分集合すべてに対して $f$ を足したものである。

実装

$N=2^n$ としています。$N$は配列長、$n$は添え字に使われるビット列の長さです。
$n$ 次元空間における $2\times 2\times 2\times ...\times 2$ の立方体の頂点に値があるとして、それの各軸に関して累積和を取ることで求めています。

コード
計算量は $O(n\cdot 2^n)=O(N\log N)$ です。

n = 3
f = [1,2,3,4,5,6,7,8]

F = [0]*(1<<n)

for i in range(1<<n):
    F[i] = f[i]

# 各次元方向について和を取る。
for i in range(n):
    for s in range(1<<n):
        if (s & (1<<i)) == 0: # 端は飛ばす
            continue
        F[s] += F[s ^ (1<<i)]
print(F)
# F = [1, 3, 4, 10, 6, 14, 16, 36]

また、計算量は悪いですが、直接部分集合を列挙をすることでも計算できます。
サイズ $|S|$ の集合の部分集合の個数は $s^{|S|}$ なので、足し算が行われる回数は、

$$\sum _{S\in P(\{0,1,...,n-1\})}{2}^{|S|}=\sum _{k=0}^n{}_n{C}_k\cdot 2^k=3^n$$

になります。よって、計算量 $O(3^n)=O(N^{1.58})$ で求めることもできます。これ以外にも、〇〇B系の変換であれば、同様の手法が使えますが、計算量が悪く使用頻度が低いため割愛します。

n = 3
f = [1,2,3,4,5,6,7,8]

F = [0]*(1<<n)
for i in range(1<<n):
  j = i
  while True:
    F[i] += f[j]
    if j == 0:
      break
    j = (j-1) & i

print(F)

例題

ARC205 E - Subset Product Problem
平方分割と組み合わせて解くことができます。
https://atcoder.jp/contests/arc205/tasks/arc205_e

PMB：$\text{P}$refix $\text{M}$obius transform for $\text{B}$itset

やりたいこと

$n$要素からなる集合 $\{0,1,2,...,n-1\}$ の各部分集合 $T$ について、$T$ を引数に取る関数 $f(T)$ が定義されているとする。
次のような関数$F$を定義する。
$F(S)=\sum _{\begin{array}{c}T\subset \{0,1,...,n-1\}\\ T\subset S\end{array}}f(T)$
つまり、$F(S)$ は、$S$ の部分集合すべてに対して $f$ を足したものである。

逆にこのように作成された$F$が与えられる。もとの$f$を復元したい。

実装

これは、先程の PZB の逆変換であるので、高次元累積和のメビウス変換となります。

コード
計算量は、$N=2^n$として、$O(n\cdot 2^n)=O(N\log N)$ です。

n = 3
F = [1,3,4,10,6,14,16,36]

f = [0]*(1<<n)

for i in range(1<<n):
    f[i] = F[i]

# 各次元方向について差を取る。
for i in range(n):
    for s in range((1<<n)-1, -1, -1): # 降順
        if (s & (1<<i)) == 0: # 端は飛ばす
            continue
        f[s] -= f[s ^ (1<<i)]
print(f)
# f = [1, 2, 3, 4, 5, 6, 7, 8]

SZB：$\text{S}$uffix $\text{Z}$eta transform for $\text{B}$itset

別名：sum over supersets

やりたいこと

$n$要素からなる集合 $\{0,1,2,...,n-1\}$ の各上位集合 $T$ について、$T$ を引数に取る関数 $f(T)$ が定義されているとする。
この $f$ から、次のような $F$ を高速に求めたい。
ただし、$F$ は $f$ と同様に、集合を引数にとる関数である。

$$F(S)=\sum _{\begin{array}{c}T\subset \{0,1,...,n-1\}\\ T\supset S\end{array}}f(T)$$

つまり、$F(S)$ は、$S$ を部分集合にもつ集合すべてに対して $f$ を足したものである。

実装

$T$ が $S$ を部分集合にもつとは、$T$ と $S$ を超立方体上の頂点に対応させたときに以下の関係が成り立つことと同値です。

$T$ に対応する点が、「 $S$ に対応する点と終点 $(1,1,...,1)$ を対角とする超直方体」の中(境界含む)にある

よって、これも同様に$n$次元累積和となる。ただし、終点からの和を取るので、Suffix $n$次元累積和になります。

コード
計算量は、$N=2^n$として、$O(n\cdot 2^n)=O(N\log N)$ です。

n = 3
f = [8,7,6,5,4,3,2,1]

F = [0]*(1<<n)

for i in range(1<<n):
    F[i] = f[i]

# 各次元方向について和を取る。
for i in range(n):
    for s in range((1<<n)-1, -1, -1): # 降順
        if (s & (1<<i)) != 0: # 端は飛ばす
            continue
        F[s] += F[s ^ (1<<i)]
print(F)
# F = [36, 16, 14, 6, 10, 4, 3, 1]

SMB：$\text{S}$uffix $\text{M}$obius transform for $\text{B}$itset

やりたいこと

SZBの逆変換です。
ただの逆変換ですが、包除原理と関連する重要なメビウス変換です。

包除原理は、$n$個の条件$A_0,A_1,...,A_{n-1}$があるとき、それらを適切に組み合わせて $A_0\cup A_1\cup ...\cup A_{n-1}$ を満たすものを数え上げるときに使う道具です。
式で書くと、次のようになります。

$$\begin{gathered} A_0\cup A_1\cup A_2\cup ...\cup A_{n-1}=\sum _{s=1}^n(-1{)}^{s-1}\sum _{\begin{array}{c}T\in \{0,1,2,...,n-1\}\\ |T|=s\end{array}}\bigcap _{i\in T}{A}_i \\ =\sum _{T\in \{0,1,\mathrm{2...},n-1\}}(-1{)}^{|T|-1}\bigcap _{i\in T}{A}_i \end{gathered}$$

上がよく紹介される形で、下は二つのシグマをくっつけた形になります。
メビウス変換と関連付けるために、次のような関数を定義します。

$$F(T)=\bigcap _{i\in T}{A}_i$$

$$f(T)=\left(\bigcap _{i\in T}{A}_i\right)\cap \left(\bigcap _{i\notin T}{\bar{A}}_i\right)$$

関数 $F(T)$ は、「集合 $T$ 内の要素に対応する条件はすべて満たすことを保証し、それ以外の条件については何でもいい」を表します。
一方で、関数 $f(T)$ は、「集合 $T$ 内の要素に対応する条件はすべて満たすことを保証し、さらに、それ以外の条件は全て満たさないことを保証する」を表します。

関数 $F(T)$ の定義に、「それ以外の条件については何でもいい」とありますが、その条件のうちどれが満たされていて、どれが満たされていないかを全て場合分けして足すことを考えます。
すると、これは $T$ を部分集合にもつ集合( $T$ の上位集合)すべてにわたって $f$ を足すことと等価であると分かります。

したがって、$f$ と $F$ には次の関係が成り立ちます。

$$F(T)=\sum _{\begin{array}{c}t\subset \{0,1,2,...,n-1\}\\ t\supset T\end{array}}f(t)$$

これは、$F$ が $f$ のSZBの結果であることを示しています。よって、逆変換のSMBを施すことで $f$ を求めることができます。

$f$さえ求まれば、$A_0\cup A_1\cup A_2\cup ...\cup A_{n-1}$ は次のように表せるため、SMBは包除原理よりも一般的な手法であると言えます。

$$A_0\cup A_1\cup A_2\cup ...\cup A_{n-1}=\sum _{\begin{array}{c}t\subset \{0,1,2,...,n-1\}\\ t\ne \varnothing \end{array}}f(t)$$

実装

計算量は、$N=2^n$ として、$O(n\cdot 2^n)=O(N\log N)$ です。

n = 3
F = [100, 50, 33, 16, 20, 10, 6, 3]

f = [0]*(1<<n)
for i in range(1<<n):
    f[i] = F[i]

for i in range(n):
    for s in range(1<<n): #昇順
        if (s & (1<<i)) != 0:
            continue
        f[s] -= f[s ^ (1<<i)]

print(f)
# f = [26, 27, 14, 13, 7, 7, 3, 3]

余談

上記のコードは、二進数の各桁を次のような条件に対応させています。
$2^0$ の位 ... $100$ 以下の整数かつ、$2$ の倍数である
$2^1$ の位 ... $100$ 以下の整数かつ、$3$ の倍数である
$2^2$ の位 ... $100$ 以下の整数かつ、$5$ の倍数である

したがって、変換後の配列の各インデックスに対応する値は次のような意味を持ちます。
$0_{(10)}=000_{(2)}$ ... 100以下の整数かつ、2の倍数でなく、3の倍数でなく、5の倍数でない ... 26個
$5_{(10)}=101_{(2)}$ ... 100以下の整数かつ、2の倍数であり、3の倍数でなく、5の倍数である ... 7個
$6_{(10)}=110_{(2)}$ ... 100以下の整数かつ、2の倍数でなく、3の倍数であり、5の倍数である ... 3個
$7_{(10)}=111_{(2)}$ ... 100以下の整数かつ、2の倍数であり、3の倍数であり、5の倍数である ... 3個

例題

ABC423 F - Loud Cicada
すでにゼータ変換された形の値なら簡単に求められます。したがって、SMBをすることで元の値を復元することができ、簡単に答えを求められます。メビウス変換を利用しなくても、包除原理の一般化を用いて解くことができます。
https://atcoder.jp/contests/abc423/tasks/abc423_f

PZF：$\text{P}$refix $\text{Z}$eta transform for $\text{F}$actorization form

別名：約数に対する累積和

やりたいこと

$1\sim N$ の全ての $f(n)$ が与えられたとき、全ての $n$ に対して $F(n)$

$$F(n)=\sum _{\begin{array}{c}d|n\end{array}}f(d)$$

を計算する。ここで $d|n$ は、$d$ は $n$ の約数であることを表す。

実装

貰うDPで実装しようとすると、添え字の約数列挙をする必要がありますが、これを配るDPで実装することでこれを回避することができます。計算量は調和級数で $\log$ のオーダーになります。
この実装は、全ての寄与を逐次足しているため、PZB における $O(3^n)=O(N^{1.58})$ 解法に相当します。

計算量：$O(N\log N)$

N = 10
f = [0,1,2,3,4,5,6,7,8,9,10] # len(f) == N + 1

F = [0]*(N+1)
for i in range(N+1):
    F[i] = f[i]

for i in range(N, 0, -1):
    for j in range(2*i, N+1, i):
        F[j] += F[i]

print(F)

また、素因数分解された結果を多次元のタプルだと思うと、高次元累積和の軸と素因数が一対一対応します。よって、エラトステネスの篩と同時に処理することで、素数を検知しつつ(軸を見つけつつ)それごとに累積和を構築することで、エラトステネスの定数倍の計算量で変換ができます。
この実装は、軸ごとに累積和を行っているので、PZB における $O(n{2}^n)=O(N\log N)$ 解法に相当します。

$O(N\log \log N)$

N = 10
f = [0,1,2,3,4,5,6,7,8,9,10] # len(f) == N + 1

F = [0]*(N+1)
for i in range(N+1):
  F[i] = f[i]

sieve = [True]*(N+1)
for p in range(2,N+1):
  if sieve[p]: # 素数1つが次元1つに対応する、多次元累積和の処理と同じ
    for k in range(1, N//p+1):
      sieve[k*p] = False
      F[k*p] += F[k]

print(F)

PMF：$\text{P}$refix $\text{M}$obius transform for $\text{F}$actorization form

別名：約数に対する階差、約数系包除

やりたいこと

PZFの逆変換、つまり、
$1\sim N$ の全ての $F(n)$ が与えられたとき、全ての $n$ に対して $f(n)$

$$F(n)=\sum _{\begin{array}{c}d|n\end{array}}f(d)$$

を計算する。ここで $d|n$ は、$d$ は $n$ の約数であることを表す。

実装

PZF の処理を逆順にし、累積和を構築する部分を階差に置き換えてやれば実装できます。

配るDPの方法に場合は、
F[1] = f[1] から f[1] を確定
F[2] = f[1] + f[2] から f[2] を確定
F[3] = f[1] + f[3] から f[3] を確定
F[4] = f[1] + f[2] + f[4] から f[4] を確定
...
とすればよいです。

実装をより詳細に言うと、
f[1] は全てに足されているので自身以外の全てから引く
f[2] は2ずつ飛んで足されているので自身以外の偶数番目から引く
f[3] は3ずつ飛んで足されているので自身以外の3の倍数番目から引く
...
となります。

計算量：$O(N\log N)$

N = 10
F = [0,1,3,4,7,6,12,8,15,13,18] # len(f) == N + 1

f = [0]*(N+1)
for i in range(N+1):
  f[i] = F[i]

for i in range(1,N+1):
    for j in range(2*i, N+1, i):
        f[j] -= f[i]

print(f)

エラトステネスと同時に行う方に関しては、素数が前からしか検知できなから完全に逆の処理をしている訳では無いではないですが、別にどの軸から先に累積和を取る/階差を取るとしても最終結果は変わらないので、結局は軸ごとに階差を取るだけでよいです。

計算量：$O(N\log \log N)$

N = 10
F = [0,1,3,4,7,6,12,8,15,13,18] # len(f) == N + 1

f = [0]*(N+1)
for i in range(N+1):
  f[i] = F[i]

sieve = [True]*(N+1)
for p in range(2,N+1):
  if sieve[p]:
    for k in range(N//p,0,-1): # 降順
      sieve[k*p] = False
      F[k*p] -= F[k] # 階差

print(F)

SZF：$\text{S}$uffix $\text{Z}$eta transform for $\text{F}$actorization form

別名：倍数に対する累積和

やりたいこと

$1\sim N$ の全ての $f(n)$ が与えられたとき、全ての $n$ に対して $F(n)$

$$F(n)=\sum _{\begin{array}{c}n|d\\ d\le N\end{array}}f(d)$$

を計算する。ここで $n|d$ は、$d$ は $n$ の倍数であることを表す。

実装

シンプルに貰うDPで実装すればよいです。エラトステネスを使うこともできます。

計算量：$O(N\log N)$

N = 10
f = [0,1,2,3,4,5,6,7,8,9,10] # len(f) == N + 1

F = [0]*(N+1)
for i in range(N+1):
    F[i] = f[i]

for i in range(1,N+1):
    for j in range(2*i,N+1,i):
        F[i] += F[j]

print(F)

PZF の処理とほぼ同じです。配るDP/貰うDPの方式よりも対応が明確だと思います。
計算量：$O(N\log \log N)$

N = 10
f = [0,1,2,3,4,5,6,7,8,9,10] # len(f) == N + 1

F = [0]*(N+1)
for i in range(N+1):
    F[i] = f[i]

sieve = [True]*(N+1)
for p in range(2,N+1):
  if sieve[p]:
    for k in range(N//p, 0, -1):  # PZF は昇順に処理
      sieve[k*p] = False
      F[k] += F[k*p]              # PZF は F[k*p] += F[k]

print(F)

SMF：$\text{S}$uffix $\text{M}$obius transform for $\text{F}$actorization form

別名：倍数に対する階差、約数系包除

やりたいこと

SZFの逆変換、つまり、
$1\sim N$ の全ての $F(n)$ が与えられたとき、全ての $n$ に対して $f(n)$

$$F(n)=\sum _{\begin{array}{c}n|d\\ d\le N\end{array}}f(d)$$

を計算する。ここで $n|d$ は、$n$ は $d$ の約数であることを表す。

実装

貰うDPの方針
上のほうから確定していきます。
F[N] = f[N] から f[N] が確定する
F[N-1] = f[N-1] から f[N-1] が確定する
...
F[N//2] = f[N//2] + f[2*(N//2)] から f[N//2] が確定する
...
F[3] = f[3] + f[6] + f[9] + ... + f[3*(N//3)] からf[3] が確定する
F[2] = f[2] + f[4] + f[6] + ... + f[2*(N//2)] からf[2] が確定する
F[1] = f[1] + f[2] + f[3] + ... + f[1*(N//1)] からf[1] が確定する

計算量: $O(N\log N)$

N = 10
F = [0,55,30,18,12,15,6,7,8,9,10] # len(F) == N + 1

f = [0]*(N+1)
for i in range(N+1):
    f[i] = F[i]

for i in range(N,0,-1):
    for j in range(2*i, N+1,i):
        f[i] -= f[j]
        
print(f)

エラトステネスも同様です。篩の二重ループの内側のループで、昇順/降順なのか、累積和/階差なのか、しか変わりません。

計算量: $O(N\log \log N)$

N = 10
F = [0,55,30,18,12,15,6,7,8,9,10] # len(F) == N + 1

f = [0]*(N+1)
for i in range(N+1):
    f[i] = F[i]

sieve = [True]*(N+1)
for p in range(2, N+1):
    if sieve[p]:
        for k in range(1, N//p+1): #昇順
          sieve[k*p] = False
          f[k] -= f[k*p] #階差

print(f)

例題

$N$ 以下の整数の組 $(x,y)$ であって、$x$ と $y$ が互いに素であるようなものはいくつあるか求めよ。
制約 $1\le N\le 10^6$

添え字が連番でない場合

72の約数全体といった、連番でない整数が添え字集合になる場合があります。

特に、$735134400$ などのように、大きい数の約数全体を考える場合、普通の配列では空間が足りなくなってしまいます。(それを確保する時間も膨大) そのため、このような場合は連想配列を用います。値の順序は重要でないので、ハッシュテーブルによる連想配列で十分です。

配列が連想配列になったこと以外は実装はほぼ同じです。
Prefix Zeta transform for Factorization form の場合の実装例だけ載せておきます。その他の場合もほとんど同じです。

ある数の約数全体の集合を考える場合は、その数が持つ素因数が次元に対応するので、予め素因数分解を行い、高次元累積和として考えると楽です。

from collections import defaultdict

def makediv_from_primes(p):
    d = [1]
    if p[0][0] == 1:
        return d
    
    for v in p:
        t = len(d)
        temp = 1
        for w in range(v[1]):
            temp *= v[0]
            for i in range(t):
                d.append(d[i] * temp)
                
    d.sort()
    return d

def factorize(n):
    if n == 1:
        return [[1, 0]]
    ret = []
    i = 2
    while i * i <= n:
        if n % i == 0:
            a = 0
            while n % i == 0:
                a += 1
                n //= i
            ret.append([i, a])
        i += 1
    if n > 1:
        ret.append([n, 1])
    return ret

N = 72
pfact = factorize(72)
divs = makediv_from_primes(pfact)

f = defaultdict(int)


for d in divs:
    f[d] = 1

"""
f = {
    1 : 1,
    2 : 1,
    3 : 1,
    4 : 1,
    6 : 1,
    8 : 1,
    9 : 1,
    12: 1,
    18: 1,
    24: 1,
    36: 1,
    72: 1
}
"""
# ここまで : fの定義
# ここから : PZF

F = f.copy()

for p in pfact:
    for d in divs:
        if d%p[0] != 0: #端を飛ばす
            continue
        F[d] += F[d//p[0]]

print(F)
"""
{
    1 : 1,
    2 : 2,
    3 : 2,
    4 : 3,
    6 : 4,
    8 : 4, 
    9 : 3,
    12: 6,
    18: 6,
    24: 8,
    36: 9,
    72: 12
}
"""

各変換同士の関係

ゼータ変換とメビウス変換が互いに逆変換の関係にあることは最初にも触れました。
同様に、いくつかの変換の間には深い関連性があります。

PZB $\leftrightarrow$ SZB, PMB $\leftrightarrow$ SMB

これらは、原点と終点を入れ替える、すなわち、添え字となっている01列を反転(0は1に、1は0に置換)することで、互いに入れ替わることが分かります。
また、添え字の01列の反転は、配列を逆順に並び替えることと等価です。

AND / OR convolution や GCD / LCM convolution との関係

これらのconvolutionは、適切に式変形をすることでゼータ変換と各点積とメビウス変換に分解することができます。

説明のために、以下のように演算子の意味を新しく定義しておきます。
非負整数 $i,j$ に対して、$i\subset j$ であるとは、次の条件と同値である。

$$i\&j=i$$

別の言い方をすると、「$i,j$ をそれぞれ二進法表記したとき、$i$ の方で $1$ が立っている場所の集合が、$j$ の方で $1$ が立っている場所の集合の部分集合となっていること」になります。

同様に、$i\supset j$ も $i\&j=j$ と定義します。

AND convolution

AND convolutionとは、次のような問題です。

長さが $2^n$ の数列 $A=\{A_0,A_1,...,A_{2^n-1}\}$ と $B=\{B_0,B_1,...,B_{2^n-1}\}$ が与えられる。
以下のように定義される長さ $2^n$ の数列 $C=\{C_0,C_1,...,C_{2^n-1}\}$ を求めよ。
ただし、a&b はビット単位AND演算を表す。
$C_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i\&j=k\end{array}}{A}_i{B}_j$

まずは、これを別の数列 $D$ を求める問題とメビウス変換に分解します。
次のような長さ $2^n$ の数列 $D=\{D_0,D_1,...,D_{2^n-1}\}$ を求める問題を考えます。
$D_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i\&j\supset k\end{array}}{A}_i{B}_j$

$D$ の定義から、$D$ と $C$ の関係性は次のように表せます。

$D_k=\sum _{\begin{array}{c}0\le l<2^n\\ l\supset k\end{array}}{C}_l$
したがって、$D$ が求まればそれをSMBすることで $C$ を得ることができます。

次に、$D$ を求める問題をゼータ変換と各点積に分解します。

$i\&j\supset k$ という条件は、$i\supset k\wedge j\supset k$ に分解できます。
したがって、$D$ の定義を次のように変更できます。
$D_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i\supset k\wedge j\supset k\end{array}}{A}_i{B}_j=\left(\sum _{\begin{array}{c}0\le i<2^n\\ i\supset k\end{array}}{A}_i\right)\left(\sum _{\begin{array}{c}0\le j<2^n\\ j\supset k\end{array}}{B}_j\right)$
これは、数列 $A,B$ をそれぞれSZBし、そのあとに各点積を取ることで $D$ が得られるということを意味します。

したがって、AND convolutionは、次の手順で計算することができます。

• 数列 $A,B$ をSZBする。
• 変換後の数列の各点積をとって、数列 $D$ を作成する。
• $D$ をSMBして目的の $C$ を得る。

計算量は、$O(n\cdot 2^n)$ です。

OR convolution

OR convolutionとは、次のような問題です。

長さが $2^n$ の数列 $A=\{A_0,A_1,...,A_{2^n-1}\}$ と $B=\{B_0,B_1,...,B_{2^n-1}\}$ が与えられる。
以下のように定義される長さ $2^n$ の数列 $C=\{C_0,C_1,...,C_{2^n-1}\}$ を求めよ。
$a|b$ はビット単位OR演算を表す。
$C_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i|j=k\end{array}}{A}_i{B}_j$

まずは、これを別の数列 $D$ を求める問題とメビウス変換に分解します。
次のような長さ $2^n$ の数列 $D=\{D_0,D_1,...,D_{2^n-1}\}$ を求める問題を考えます。
$D_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i|j\subset k\end{array}}{A}_i{B}_j$

$D$ の定義から、$D$ と $C$ の関係性は次のように表せます。

$D_k=\sum _{\begin{array}{c}0\le l<2^n\\ l\subset k\end{array}}{C}_l$
したがって、$D$ が求まればそれをPMBすることで $C$ を得ることができます。

次に、$D$ を求める問題をゼータ変換と各点積に分解します。

$i|j\subset k$ という条件は、$i\subset k\wedge j\subset k$ に分解できます。
したがって、$D$ の定義を次のように変更できます。
$D_k=\sum _{\begin{array}{c}0\le i,j<2^n\\ \\ i\subset k\wedge j\subset k\end{array}}{A}_i{B}_j=\left(\sum _{\begin{array}{c}0\le i<2^n\\ i\subset k\end{array}}{A}_i\right)\left(\sum _{\begin{array}{c}0\le j<2^n\\ j\subset k\end{array}}{B}_j\right)$
これは、数列 $A,B$ をそれぞれPZBし、そのあとに各点積を取ることで $D$ が得られるということを意味します。

したがって、OR convolutionは、次の手順で計算することができます。

• 数列 $A,B$ をPZBする。
• 変換後の数列の各点積をとって、数列 $D$ を作成する。
• $D$ をPMBして目的の $C$ を得る。

計算量は、$O(n\cdot 2^n)$ です。

GCD convolution

GCD convolutionとは、次のような問題です。

長さが $N$ の数列 $A=\{A_1,A_2,...,A_N\}$ と B=\{B_\red{1}, B_2,...,B_{N}\} が与えられる。
以下のように定義される長さ $N$ の数列 C=\{C_\red1,C_2,...,C_{N}\} を求めよ。
$C_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \gcd (i,j)=k\end{array}}{A}_i{B}_j$

これも、AND/OR convolutionのように、分解します。
次のような数列 D = \{D_\red1, D_2,...,D_N\} を考えます。
ただし、$a|b$ は $a$ が $b$ の約数であることを表します。(ビット単位ORではありません。)
$D_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \\ k|\gcd (i,j)\end{array}}{A}_i{B}_j$
$D$ と $C$ の関係は次のようになります。
$D_k=\sum _{\begin{array}{c}1\le l\le N\\ \\ k|l\end{array}}{C}_l$
したがって、$D$ をSMFすることで $C$ を求められます。

次に、$D$ を求める問題をゼータ変換と各点積に分解します。
$k|\gcd (i,j)$ は $k|i\wedge k|j$ に分解できます。よって、次のように変形できます。
$D_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \\ k|i\wedge k|j\end{array}}{A}_i{B}_j=\left(\sum _{\begin{array}{c}1\le i\le N\\ k|i\end{array}}{A}_i\right)\left(\sum _{\begin{array}{c}1\le j\le N\\ k|j\end{array}}{B}_j\right)$
これは、数列 $A,B$ をそれぞれSZFし、そのあとに各点積を取ることで $D$ が得られるということを意味します。

したがって、GCD convolutionは、次の手順で計算することができます。

• 数列 $A,B$ をSZFする。
• 変換後の数列の各点積をとって、数列 $D$ を作成する。
• $D$ をSMFして目的の $C$ を得る。

計算量は、$O(N\log N)$ または $O(N\log \log N)$ です。SZF, ZMFの計算量によります。

LCM convolution

LCM convolutionとは、次のような問題です。

長さが $N$ の数列 $A=\{A_1,A_2,...,A_N\}$ と B=\{B_\red{1}, B_2,...,B_{N}\} が与えられる。
以下のように定義される長さ $N$ の数列 C=\{C_\red1,C_2,...,C_{N}\} を求めよ。
$C_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \text{lcm}(i,j)=k\end{array}}{A}_i{B}_j$

これも、AND/OR convolutionのように、分解します。
次のような数列 D = \{D_\red1, D_2,...,D_N\} を考えます。
ただし、$a|b$ は $a$ が $b$ の約数であることを表します。(ビット単位ORではありません。)
$D_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \\ \text{lcm}(i,j)|k\end{array}}{A}_i{B}_j$
$D$ と $C$ の関係は次のようになります。
$D_k=\sum _{\begin{array}{c}1\le l\le N\\ \\ l|k\end{array}}{C}_l$
したがって、$D$ をPMFすることで $C$ を求められます。

次に、$D$ を求める問題をゼータ変換と各点積に分解します。
$\text{lcm}(i,j)|k$ は $i|k\wedge j|k$ に分解できます。よって、次のように変形できます。
$D_k=\sum _{\begin{array}{c}1\le i,j\le N\\ \\ i|k\wedge j|k\end{array}}{A}_i{B}_j=\left(\sum _{\begin{array}{c}1\le i\le N\\ i|k\end{array}}{A}_i\right)\left(\sum _{\begin{array}{c}1\le j\le N\\ j|k\end{array}}{B}_j\right)$
これは、数列 $A,B$ をそれぞれPZFし、そのあとに各点積を取ることで $D$ が得られるということを意味します。

したがって、LCM convolutionは、次の手順で計算することができます。

• 数列 $A,B$ をPZFする。
• 変換後の数列の各点積をとって、数列 $D$ を作成する。
• $D$ をPMFして目的の $C$ を得る。

計算量は、$O(N\log N)$ または $O(N\log \log N)$ です。PZF, PMFの計算量によります。

全体を通して

これらの convolution は、それぞれ

• AND convolution：SZB → 各点積 → SMB
• OR convolution：PZB → 各点積 → PMB
• GCD convolution：SZF → 各点積 → SMF
• LCM convolution：PZF → 各点積 → PMF

という構造をしています。これは

• 通常(添え字の和に関する) convolution：離散フーリエ変換 → 各点積 → 離散逆フーリエ変換
• XOR convolution：アダマール変換 → 各点積 → 逆アダマール変換

にも対応します。

おわりに

高速ゼータ変換/高速メビウス変換について、約数系包除原理との並べて解説しました。混同されがちな概念ですが、結局は多次元累積和をいい感じに取るテクニックなんだな～ということさえ分かればOKな気がします。ここまでの内容を読んで理解できたら、ゼータ変換/メビウス変換がいざ必要になったときでもスラスラ書けると思います。

というか、書いてる途中に気づいたんですが、AtCoder Lecture でゼータ変換/メビウス変換が紹介されてますね。だいぶ説明に被りが多いな...

まあ、このブログを読んで、「ゼータ変換なんかちょっと複雑な多次元累積和だな～」と思えるようになって頂けたら、AtCoder Lecture と独立な理解を得たということで喜ばしく思います。

よく分からなくなったら、三次元の場合を考えてみましょう。四次元以上は想像が難しいので、まずは簡単な具体例で理解するのがいいと思います。
実装も、「三次元直方体の内部の総和を求めている」というイメージを持つことでバグを防ぐことができると思います。

参考文献

https://naoyat.hatenablog.jp/entry/zeta-moebius
https://codeforces.com/blog/entry/45223
https://noshi91.hatenablog.com/entry/2018/12/27/121649
https://blog.ngtkana.com/entry/2021/04/25/135636

明日の投稿者は@sumire0517です