2021 计蒜之道 精英组 决赛 Day2 类斐波那契数据分析
设\(f_0=0,f_1=1,f_i=f_{i-1}c+f_{i-2}\)。
求\(\sum_{i=1}^n\sum_{j=1}^i [i\perp j]\gcd(f_i,f_{j+d})\)。
\(n\le 10^6,d,c\le 2^{30}\)，模数不固定也不保证是质数。
首先有结论：\(\gcd(f_x,f_y)=f_{\gcd(x,y)}\)。
证明：
首先可证\(\gcd(f_i,f_{i+1})=1\)，归纳易证。
由\(f_{x}=f_{x-y}f_{y+1}+f_{x-y+1}f_y\)得\(原式=\gcd(f_{x-y}f_{y+1}+f_{x-y+1}f_y,f_y)=\gcd(f_{x-y}f_{y+1},f_y)=\gcd(f_{x-y},f_y)\)。
令\(\sum_{d|n}h(d)=f(n)\)。推式子：
发现只有\(k\perp d\)有贡献。
\(\phi(i)=\sum_{j=1}^i [i\perp j]=\sum_{d=0}^{k-1}[k\perp d]\sum_{j=1}^i[i\perp j][k|j+d]\)
接下来证明：\(k|i\)，对于每个满足\(k\perp d\)的\(d\)，\(\sum_{j=1}^i[i\perp j][k|j+d]\)相等。
把\(i\)写成\(kx\)。上式化为：\(\sum_{j=1}^x [kx\perp kj-d]=\sum_{j=1}^x [x\perp kj-d]\)。
以下可以证明：对于不同的满足\(k\perp d\)的\(d\)，\(j=1..x,\gcd(x,kj-d)\)形成的可重集相同。
记下每个\(\gcd\)的出现次数，反演一下变成记下每个约数的出现次数。
于是要证，不同的\(d\)之间对于每个约数\(t|x\)，都满足\(\sum_{j=1}^x[t|kj-d]\)相同。
因为\(t|x\)不妨比较\(\sum_{j=0}^{t-1}[kj\equiv d\pmod t]\)。
看\(kj\equiv d\pmod t\)是否有解。拆开得\(kj-ty=d\)。
如果有解一定有\(\gcd(k,t)|d\)，又因\(k\perp d,\gcd(k,t)|k\)，则\(\gcd(k,t)=1\)。
于是如果有解，\(k\)在模\(t\)意义下有逆元。\(j\equiv\frac{d}{k}\pmod t\)，于是对于不同的\(d\)都会有个确定的解。
以上结论证毕。
既然\(\gcd\)的可重集相同了，那么\(\sum_{j=1}^x [x\perp kj-d]\)也自然相同。
因此，对于某个\(d\)，贡献是\(\frac{\phi(i)}{\phi(d)}\)。
后面的计算就没有什么难度了。
using namespace std;
#include <bits/stdc++.h>
#define N 1000005
#define ll long long
int gcd(int a,int b){
int k;
while (b)
k=a%b,a=b,b=k;
return a;
}
int n,d,c,mo;
int p[N],np;
bool inp[N];
int phi[N];
void initp(int n=1000000){
phi[1]=1;
for (int i=2;i<=n;++i){
if (!inp[i])
p[++np]=i,phi[i]=i-1;
for (int j=1;j<=np && i*p[j]<=n;++j){
inp[i*p[j]]=1;
if (i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
}
}
}
ll f[N];
int main(){
// freopen("in.txt","r",stdin);
freopen("analyze.in","r",stdin);
freopen("analyze.out","w",stdout);
initp();
int T;
scanf("%d",&T);
while (T--){
scanf("%d%d%d%d",&n,&d,&c,&mo);
f[0]=0,f[1]=1;
for (int i=2;i<=n;++i)
f[i]=(f[i-1]*c+f[i-2])%mo;
for (int i=1;i<=n;++i){
f[i]%=mo;
for (int j=i+i;j<=n;j+=i)
f[j]-=f[i];
}
ll ans=0;
for (int k=1;k<=n;++k)
if (gcd(k,d)==1){
ll s=0;
for (int i=k;i<=n;i+=k)
s+=phi[i];
(ans+=s/phi[k]%mo*f[k])%=mo;
}
ans=(ans+mo)%mo;
printf("%lld\n",ans);
}
return 0;
}
