【BZOJ4833】最小公倍佩尔数（min-max容斥）
题面
题解
首先考虑怎么求(f(n))，考虑递推这个东西
((1+sqrt 2)(e(n-1)+f(n-1)sqrt 2)=e(n)+f(n)sqrt 2)
拆开之后可以得到：(e(n)=e(n-1)+2f(n-1),f(n)=f(n-1)+e(n-1))。
把每一层的(e)都给展开，得到：(displaystyle f(n)=1+f(n-1)+2sum_{i=1}^{n-2}f(i))
然后差分搞搞，(displaystyle f(n)-f(n-1)=f(n-1)-f(n-2)+2*f(n-2))。
得到(f(n)=2f(n-1)+f(n-2))，特殊的(f(0)=0,f(1)=1)。
然后我们发现要求(lcm)，那么就先考虑(f(a))和(f(b))的(gcd)是什么。
这个东西显然可以类似斐波那契数列那样子利用辗转相减得到(gcd(f(a),f(b))=f(gcd(a,b)))。
接下来就可以考虑怎么求答案了。
然后(lcm)的式子是对于每个质因子，考虑其(max)。
考虑(min-max)容斥，把(max)变成(min)，那么就可以从(lcm)变成(gcd)。
然后把(min-max)容斥的式子给写出来：
[max(S)=sum_{Tsubset S}(-1)^{|T|+1}min(T)
]
套到(lcm)上就是：
[lcm(S)=prod_{Tsubset S}gcd(T)^{(-1)^{|T|+1}}
]
那么就有
[g(n)=prod_{Tsubset S}f_{gcd(T)}^{(-1)^{|T|+1}}=prod_{i=1}^n f_i^{sum_{Tsubset S}[gcd(T)=i](-1)^{|T|+1}}
]
上面那个指数看着就可以莫比乌斯反演一下之类的，然后令上面那一堆东西是(a[i])，然后令(b[i]=sum_{i|d}a[d])这个系数稍微推一下，得到：
[b[i]=sum_{i|d}a[d]=sum_{Tsubset S}[i|gcd(T)](-1)^{|T|+1}
]
这个值显然之和是否存在(i)倍数的数相关，存在就是(1)，没有就是(0)。
而莫比乌斯反演可以得到
[a[i]=sum_{i|d}mu(frac{d}{i})b[d]
]
再把这个东西带回去
[egin{aligned}
g[n]&=prod_{i=1}^n f_i^{a[i]}\
&=prod_{i=1}^n f_i^{sum_{i|d}mu(frac{d}{i})b[d]}\
&=prod_{i=1}^nprod_{i|d}f_i^{mu(frac{d}{i})b[d]}
end{aligned}]
因为(d)的范围在(n)以内，所以必定存在(d)的倍数，所以(b[d]=1)，那么只需要提前一个(log)预处理后面一半就行了。
#include<iostream>
#include<cstdio>
using namespace std;
#define MAX 1000100
inline int read()
{
int x=0;bool t=false;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')t=true,ch=getchar();
while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
return t?-x:x;
}
int n,MOD;
bool zs[MAX];
int pri[MAX],mu[MAX],tot;
int f[MAX],g[MAX],s[MAX],inv[MAX];
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
void Sieve(int n)
{
mu[1]=1;
for(int i=2;i<=n;++i)
{
if(!zs[i])pri[++tot]=i,mu[i]=-1;
for(int j=1;j<=tot&&i*pri[j]<=n;++j)
{
zs[i*pri[j]]=true;
if(i%pri[j])mu[i*pri[j]]=-mu[i];
else{mu[i*pri[j]]=0;break;}
}
}
}
int main()
{
Sieve(MAX-1);
int T=read();
while(T--)
{
n=read();MOD=read();
f[1]=1;for(int i=2;i<=n;++i)f[i]=(2ll*f[i-1]+f[i-2])%MOD;
for(int i=1;i<=n;++i)s[i]=1,inv[i]=fpow(f[i],MOD-2);
for(int i=1;i<=n;++i)
for(int j=i;j<=n;j+=i)
if(mu[j/i]==1)s[j]=1ll*s[j]*f[i]%MOD;
else if(mu[j/i]==-1)s[j]=1ll*s[j]*inv[i]%MOD;
g[0]=1;for(int i=1;i<=n;++i)g[i]=1ll*g[i-1]*s[i]%MOD;
int ans=0;for(int i=1;i<=n;++i)ans=(ans+1ll*g[i]*i)%MOD;
printf("%d
",ans);
}
}