原創題目 GCD卷積 Problem and Solution
比賽用題面、題解、標程和資料生成器均已掛在 [email protected]:sun123zxy/gcdconv.git 上。
Problem
GCD卷積 (gcdconv.cpp/.in/.out) (1s,512MB)
Description
定義一種新的卷積 —— GCD卷積，其接受兩個長度為 \(n\) 的序列 \(f,g\) ，依據下式生成長度為 \(n\) 的序列 \(h\) ：
\[h_k = \sum_{\gcd(i,j) = k} f_i g_j \]現給定序列 \(f,g\) ，求各 \(h_i\) 對 \(998244353\) 取模後的值。
Input
第一行輸入一個正整數 \(n\)
，表示 \(f,g\) 的長度。第二行輸入 \(n\) 個整數 \(f_i\) 。
第三行輸入 \(n\) 個整數 \(g_i\) 。
Output
為減少輸出量，只需輸出1個整數，表示各 \(h_i\) 對 \(998244353\) 取模後的異或和。
Sample 1
Sample 1 Input
3 5 1 4 2 3 3Sample 1 Output
78Sample 1 Explanation
\[\begin{aligned} h_1 &= f_1 ( g_1 + g_2 + g_3 ) + g_1 (f_2 + f_3) + f_2 g_3 + f_3 g_2 = 65 \\ h_2 &= f_2 g_2 = 3 \\ h_3 &= f_3 g_3 = 12 \end{aligned} \]\[65 \oplus 3 \oplus 12 = 78 \]Sample 2
Sample 2 Input
4 7 1 8 0 6 2 9 1Sample 2 Output
158Sample 2 Explanation
\[\begin{aligned} h_1 &= 213 \\ h_2 &= 3 \\ h_3 &= 72 \\ h_4 &= 0 \end{aligned} \]\[213 \oplus 3 \oplus 72 \oplus 0 = 158 \]Sample 3
見
sample目錄下gcdconv3.in/.ans。Constraints
對20%的資料，\(1 \le n \le 2000\) 。
對100%的資料， \(1 \le n \le 4 \times 10^5\) ， \(0 \le f_i, g_i \le 998244352\) 。
Hints
時限在std的1.5倍左右。std沒有卡常，資料有一定梯度，請放心食用。
Source
sun123zxy
???
- 樣例2比較暴力。
記號說明
預設諸如 \(n,d,i,j,k\) 的下標變數的最大值為題目中的給出的序列長度，並把序列更換為數論函式來表示。
注意本題解中的 \(n\) 通常是一個變數，和題目中定義的序列長度 \(n\) 不同。
另外，用 \(\circ\) 代表 \(\gcd\) 卷積，即
\[h(n) = (f \circ g)(n) = \sum_{\gcd(i,j) = n} f(i) g(j) \]Solution
我們按照快速傅立葉變換（FFT）、快速莫比烏斯變換（FMT）解決卷積的思路來解決該問題——構造一種變換來滿足卷積定理：
\[\hat f(n) \hat g(n) = \widehat {(f \circ g)}(n) \]\(\hat f\) 即對函式 \(f\) 進行該變換後得到的函式。
通過一些敏銳的直覺，我們能感受到 \(\gcd\) 和列舉約數或者倍數有一些關係。
容易想到構造出一種變換，稱之為“倍數和變換”：
\[\hat f(n) = \sum_{n|d} f(d) \]並根據莫比烏斯反演得到它的逆變換
\[f(n) = \sum_{n|d} \mu(\frac n d) \hat f(d) \]（這實際上是標準莫比烏斯反演的另一種形式，詳見後文 [FurtherThoughts](#Further Thoughts) ）
這個變換對 \(\gcd\) 卷積滿足卷積定理，證明如下。
首先，寫出 \(\gcd\) 卷積
\[h(k) = (f \circ g)(k) = \sum_{i,j} [\gcd(i,j)=k] f(i) g(j) \]左右兩邊做倍數和變換：
\[\begin{aligned} \hat h(n) &= \sum_{n|k} \sum_{i,j} [\gcd(i,j)=k] f(i) g(j) \\ &= \sum_{i,j} \left( \sum_{n|k} [\gcd(i,j)=k] \right) f(i) g(j) \quad \\ &= \sum_{i,j} [n|i][n|j] f(i) g(j) \quad \\ &= \sum_{n|i} f(i) \sum_{n|j} g(j) \quad \\ &= \hat f(n) \hat g(n) \end{aligned} \]得證。上述證明的核心在於 \(\sum_{n|k} [\gcd(i,j)=k] = [n|i] [n|j]\) 。
於是，先對 \(f,g\) 做倍數和變換，然後直接 \(O(n)\) 點值相乘，再逆變換回來，就能得到 \(f \circ g\) 。
那麼剩下的問題在於如何快速做倍數和變換及其逆變換。這是非常simple的，直接暴力就好了。複雜度為 \(O(n H(n))\) ，其中調和級數 \(H(n)= \sum_{k=1}^{n} \frac 1 k\) ，有 \(\lim_{n \to \infty} H(n) = \ln(n) + c\) 。尤拉常數 \(c \approx 0.57721566490153286060651209\) 。
可以稱之為“快速倍數和變換”。
總時間複雜度約為 \(O(n \ln n)\) 。
Code
/*
gcd卷積 (gcdconv) std
by sun123zxy
*/
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<ctime>
#include<cstdlib>
#include<queue>
#include<vector>
#include<map>
#include<set>
using namespace std;
typedef long long ll;
ll Rd(){
ll ans=0;bool fh=0;char c=getchar();
while(c<'0'||c>'9'){if(c=='-') fh=1; c=getchar();}
while(c>='0'&&c<='9') ans=ans*10+c-'0',c=getchar();
if(fh) ans=-ans;
return ans;
}
const ll MOD=998244353;
#define _ %MOD
ll PMod(ll x){
if(x<0) return x+MOD;
else if(x>=MOD) return x-MOD;
else return x;
}
const ll MXN=5E5+5;
ll P[MXN],mu[MXN];ll pN;
bool notP[MXN];
void LinearSieve(ll n){
notP[1]=1;for(ll i=2;i<=n;i++) notP[i]=0;
P[1]=0;mu[1]=1;
pN=0;
for(ll i=2;i<=n;i++){
if(!notP[i]){
P[++pN]=i;
mu[i]=-1;
}
for(ll j=1;i*P[j]<=n;j++){
notP[i*P[j]]=1;
if(i%P[j]==0){
mu[i*P[j]]=0;
break;
}
mu[i*P[j]]=mu[i]*mu[P[j]];
}
}
}
class Poly{public:
ll& operator [] (ll idx){return cof[idx];}
ll n;vector<ll> cof;
Poly(){}
Poly(ll n){Resize(n);}
void Resize(ll n){
this->n = n;
cof.resize(n+1,0);
}
};
namespace PC{//PolyCalc
void FMT(Poly& A,ll typ){
ll n=A.n;
Poly B(n);
for(ll i=1;i<=n;i++){
for(ll j=1;i*j<=n;j++){
ll t=A[i*j];if(typ==-1) t=t*mu[j]_;
B[i]=PMod(B[i]+t);
}
}
A=B;
}
Poly GcdConv(Poly A,Poly B){
ll n=min(A.n,B.n);
Poly C(n);
FMT(A,1);FMT(B,1);
for(ll i=1;i<=n;i++) C[i]=A[i]*B[i]_;
FMT(C,-1);
return C;
}
}
ll N;
int main(){
freopen("gcdconv.in","r",stdin);
freopen("gcdconv_std.out","w",stdout);
N=Rd();
LinearSieve(N);
Poly A(N),B(N);
for(ll i=1;i<=N;i++) A[i]=Rd();
for(ll i=1;i<=N;i++) B[i]=Rd();
Poly C=PC::GcdConv(A,B);
ll ans=0;
for(ll i=1;i<=N;i++) ans^=C[i];
printf("%lld",ans);
return 0;
}
Further Thoughts
本題用到的是莫比烏斯反演的一種變形
\[g(n) = \sum_{n|d} f(n) \iff f(d) = \sum_{n|d} \mu(\frac n d) g(d) \]原版莫比烏斯反演長這樣
\[g(n) = \sum_{d|n} f(n) \iff f(d) = \sum_{d|n} \mu(\frac n d) g(d) \]若把題目改成 \(\mathrm{lcm}\) 卷積，即
\[h_k = \sum_{\mathrm{lcm}(i,j) = k} f_i g_j \]用到的就是原版莫比烏斯反演。相應的構造一個“約數和變換”即可解決。而約數和變換可以用“快速約數和變換”（實際上就是個埃篩）實現，有興趣的同學可以試試。
Further Further Thoughts (update 2020/10/28)
我們來探討一些更加接近本質的東西。
普通卷積（多項式乘法）、集合並/交卷積、還有 \(\mathrm{lcm}\) / \(\gcd\) 卷積，這三個問題我們是如何解決的？
關鍵在於——我們分別構造出了離散傅立葉變換、子集和/超集和變換（莫比烏斯變換）和約數和/倍數和變換來滿足卷積定理。
而這三個變換的共性在哪裡呢？
他們都是偏序集上的具有優美性質的廣義反演。
- 離散傅立葉變換用的是單位根反演——其偏序集是定義在自然數集上的小於等於關係；
- 子集和/超集和變換用的是容斥原理——其偏序集是定義在集合上的包含關係；
- 而約數和/倍數和變換則運用了莫比烏斯反演——其偏序集是定義在自然數集上的整除關係。
反演的意義就在於——它為將原來的“係數”轉化為“點值”， \(O(n)\) 乘起來後又轉回“係數”提供了可能。而這三個特殊的反演之所以優美，就是因為我們發現了他們的“快速 (Fast) ”變換方式，讓我們能通過這個Trick，免去了 \(O(n^2)\) 的暴力運算。
還有更重要的一點——子集和/超集和變換與約數和/倍數和變換有著更加緊密的聯絡——這也是筆者由集合的FMT想到約數和變換的關鍵所在——以全體素數為基底可以張成一個包含全體正整數的“素數空間”（其良定義依賴於整數惟一分解定理）。
約數和變換就是“素數空間”中的高維字首和，對應著“集合空間”中的子集和變換。換句話說，子集和變換、約數和變換分別是莫比烏斯變換在“素數空間”、“集合空間”的具體體現，它們的本質都是高維字首和。
如果理解到上面這一點，就不難發現集合並卷積和 \(\mathrm{lcm}\) 卷積的共性了——他們都把數放到了對兩個輸入元素各維度座標取 \(\max\) 後的位置；而集合交卷積和 \(\gcd\) 卷積則是取 \(\min\) 。
綜上所述，有了上面的認識，可以說 \(\gcd\) 卷積的發現是非常自然的。
後記 & 致謝
Idea是在語文課上摸出來的，Solution是在隨後的數學課上想出來的
很早以前就覺得FMT和莫比烏斯反演有些說不清道不明的關係，出這道題也讓筆者對其理解更加深入了。
其實早有大佬對 \(\mathrm{lcm}\) / \(\gcd\) 卷積展開過研究。Google一波可以發現國外數學社群有這方面的討論，最晚2013年arXiv上就有討論其性質的論文了（筆者駑鈍，確實無力理解，感興趣的大佬可以直接Google "GCD Convolution"瞭解）。只可惜如此自然而美妙的 \(\gcd\) 卷積，竟然沒有隨著那篇集合冪級數的論文引入國內，讓筆者感到有些遺憾。
總之，這道題深入的考察對FFT、FMT原理的理解以及對 \(\gcd\) 、莫比烏斯反演的敏銳直覺，是不可多得的數論好題（
下面是戰術感謝環節。
感謝keke學長為我們教授集合冪級數。
感謝TbYangZ、diong神、changruinian2020幾位神仙的點撥。
And you...
——sun123zxy
Sep. 2020 初稿完成
Nov. 2020 最後更新