# 原创OI题目 GCD卷积 Problem and Solution

## Problem

GCD卷积 (gcdconv.cpp/.in/.out) (1s,512MB)

Description

$h_k = \sum_{\gcd(i,j) = k} f_i g_j$

Input

Output

Sample 1

Sample 1 Input

3
5 1 4
2 3 3


Sample 1 Output

78


Sample 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 1


Sample 2 Output

158


Sample 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

Hints

Source

sun123zxy

• 样例2比较暴力。

## 记号说明

$h(n) = (f \circ g)(n) = \sum_{\gcd(i,j) = n} f(i) g(j)$

## Solution

$\hat f(n) \hat g(n) = \widehat {(f \circ g)}(n)$

$\hat f$ 即对函数 $f$ 进行该变换后得到的函数。

$\hat f(n) = \sum_{n|d} f(d)$

$f(n) = \sum_{n|d} \mu(\frac n d) \hat f(d)$

（这实际上是标准莫比乌斯反演的另一种形式，详见后文 [Further Thoughts](#Further Thoughts) ）

$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}

## 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)$

$h_k = \sum_{\mathrm{lcm}(i,j) = k} f_i g_j$

## Further Further Thoughts (update 2020/10/28)

• 离散傅里叶变换用的是单位根反演——其偏序集是定义在自然数集上的小于等于关系；
• 子集和/超集和变换用的是容斥原理——其偏序集是定义在集合上的包含关系；
• 而约数和/倍数和变换则运用了莫比乌斯反演——其偏序集是定义在自然数集上的整除关系。

## 后记 & 致谢

Idea是在语文课上摸出来的，Solution是在随后的数学课上想出来的

And you...

——sun123zxy

Sep. 2020 初稿完成

Nov. 2020 最后更新

posted @ 2020-12-06 15:12  sun123zxy  阅读(364)  评论(0编辑  收藏  举报