题解 最小公倍佩尔数

传送门

又一道神题……

先扔一个坑在这里

• 类菲波那切数列(满足 f(n+m) = f(n+1) * f(m) + f(n) * f(m-1)) 都是满足 gcd(f(x) , f(y)) = f( gcd(x,y) )的

dalao题解
首先打表可以找到 \(f\) 的递推式（当然像上面题解一样直接推也是可以的只是我没推出来）
于是有

\[pell_i=2\times pell_{i-1}+pell_{i-2} \]

这里要求lcm，\(n\) 又大的离谱，猜测 \(\gcd(pell_i, pell_j)\) 有什么性质
继续猜测+打表验证可得 \(\gcd(pell_i, pell_j)=pell_{\gcd(i, j)}\)
于是发现现在gcd很好求，考虑min-max容斥出lcm
那就是

\[lcm(S)=\prod_{T\subset S}gcd(T)^{(-1)^{|T|+1}} \]

我在这里脑残了两个地方：一个是 \(\prod\) 写着写着写成 \(\sum\) 了，一个是试图枚举集合大小将指数展开
但其实集合gcd的话也可以不展开，直接带着集合反演

然后后面的式子见上面题解
这里解释一下题解里那句 「这个值显然之和是否存在i倍数的数相关，存在就是1，没有就是0」
证明在这里
大意是若只有 \(i\) 的1倍，贡献是1
此后每个可能的 \(k*i\) 选或不选造成的两种贡献互相抵消了，所以总贡献仍为1

然后考虑怎么方便地统计每个 \(g(i)\)
枚举每个 \(pell\) 数的倍数，开桶存一下就好
总复杂度 \(O(Tn\ln n)\)

点击查看代码

#include <bits/stdc++.h>
using namespace std;
#define INF 0x3f3f3f3f
#define N 1000010
#define ll long long
//#define int long long

char buf[1<<21], *p1=buf, *p2=buf;
#define getchar() (p1==p2&&(p2=(p1=buf)+fread(buf, 1, 1<<21, stdin)), p1==p2?EOF:*p1++)
inline int read() {
	int ans=0, f=1; char c=getchar();
	while (!isdigit(c)) {if (c=='-') f=-f; c=getchar();}
	while (isdigit(c)) {ans=(ans<<3)+(ans<<1)+(c^48); c=getchar();}
	return ans*f;
}

int n;
bool npri[N];
int pri[N], pcnt, mu[N];
ll mod, pell[N], buc[N], g[N];
inline ll qpow(ll a, int b) {ll ans=1; for (; b; a=a*a%mod,b>>=1) if (b&1) ans=ans*a%mod; return ans;}

signed main()
{
	mu[1]=1;
	for (int i=2; i<N; ++i) {
		if (!npri[i]) pri[++pcnt]=i, mu[i]=-1;
		for (int j=1,x; j<=pcnt&&i*pri[j]<N; ++j) {
			npri[x=i*pri[j]]=1;
			if (!(i%pri[j])) break;
			else mu[x]=-mu[i];
		}
	}

	int T=read();
	while (T--) {
		n=read(); mod=read();
		pell[1]=g[0]=1;
		for (int i=2; i<=n; ++i) pell[i]=(2*pell[i-1]+pell[i-2])%mod;
		for (int i=1; i<=n; ++i) buc[i]=1;
		for (int i=1; i<=n; ++i) {
			ll inv=qpow(pell[i], mod-2);
			for (int j=i; j<=n; j+=i)
				if (mu[j/i]==1) buc[j]=buc[j]*pell[i]%mod;
				else if (mu[j/i]==-1) buc[j]=buc[j]*inv%mod;
		}
		for (int i=1; i<=n; ++i) g[i]=g[i-1]*buc[i]%mod;
		for (int i=2; i<=n; ++i) g[i]=(g[i-1]+i*g[i])%mod;
		printf("%lld\n", g[n]);
	}
	
	return 0;
}