BZOJ 4833: [Lydsy1704月赛]最小公倍佩尔数

好诡异的一道题的说，不过我连那个\(\operatorname{lcm}\)和\(\gcd\)的互化都不知道，真是太屑了

首先题目中给出的\(f(n)\)是通项公式，我们可以用待定系数法解特征根方程然后求出递推式，求出来是

\[\begin{cases}0&n=0\\1&n=1\\f(n)=2f(n-1)+f(n-2)&n\ge 2\\\end{cases} \]

然后我们知道形如\(f(n)=af(n-1)+bf(n-2)\)的递推方程都有性质\(\gcd(f(n),f(m))=f(\gcd(n,m))\)，这个和斐波那契数列的证明类似，直接大力归纳即可

然后接下来由于我们要求\(g(n)=\operatorname{lcm}_{i=1}^n f(i)\)，我们考虑经典的\(\operatorname{lcm}\)和\(\gcd\)的互化（以下\(S=\{1,2,\cdots,n\}\)）：

\[\operatorname{lcm}(S)=\prod_{T\subseteq S,T\not=\empty} \gcd(T)^{(-1)^{|T|+1}} \]

什么，你说这很新鲜，怎么推出来的？其实就是把min-max容斥放在了质数的指数上，然后指数的加法相当于乘法，于是就有了这个式子

然后我们就有：

\[g(n)=\prod_{T\subseteq S,T\not=\empty}f(\gcd(T))^{(-1)^{|T|+1}} \]

然后是不是很摸不着头脑，然后这道题最精髓的一步来了（绝妙！），由于题中保证了\(f(n)\not =0\)，并且是在模意义下求答案，因此我们一定可以构造一个函数\(h(n)\)满足\(f(n)=\prod_{d|n} h(d)\)

然后代进去就有：

\[g(n)=\prod_{T\subseteq S,T\not=\empty}f(\gcd(T))^{(-1)^{|T|+1}}\\=\prod_{T\subseteq S,T\not=\empty}[\prod_{d|\gcd(T)} h(d)]^{(-1)^{|T|+1}}\\=\prod_{d=1}^n h(d)^{\sum_{T\subseteq S',T\not=\empty} (-1)^{|T|+1}} \]

其中\(S'=\{1,2,\cdots,\lfloor\frac{n}{d}\rfloor\}\)，最后一步转化是提前\(h(d)\)考虑它的贡献，注意关于\(\prod\)的式子转化时本来的乘法到指数上就要变成加法

看起来并没有让问题更简单是么，我们来看一下\(h(d)\)的指数：

\[\sum_{T\subseteq S',T\not=\empty} (-1)^{|T|+1}=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} (-1)^{i+1} C_{\lfloor\frac{n}{d}\rfloor}^i\\=1+\sum_{i=0}^{\lfloor\frac{n}{d}\rfloor} (-1)^iC_{\lfloor\frac{n}{d}\rfloor}^i=1 \]

woc全体起立，思路清奇！这样就有\(g(n)=\prod_{i=1}^n h(i)\)

然后考虑已知\(f(n)\)的情况下怎么推出\(h(n)\)，很简单，我们都知道莫比乌斯反演：

\[f(n)=\sum_{d|n} g(d)\Leftrightarrow g(n)=\sum_{d|n} \mu(d)\times f(\frac{n}{d}) \]

然后把它的求和变成乘积即可（考虑反演公式的本质是关于质因数的容斥，而质数指数的加减就是原式子的乘除）：

\[f(n)=\prod_{d|n} g(d)\Leftrightarrow g(n)=\prod_{d|n} f(\frac{n}{d})^{\mu(d)} \]

然后就做完了，说实话那一步转化真是太妙了，佩服出题人的脑洞

#include<cstdio>
#define RI register int
#define CI const int&
using namespace std;
const int N=3e6+5;
int t,n,mod,prime[N],cnt,mu[N],f[N],ivf[N],h[N],ans; bool vis[N];
#define P(x) prime[j]
inline void init(CI n)
{
	RI i,j; for (vis[1]=mu[1]=1,i=2;i<=n;++i)
	{
		if (!vis[i]) prime[++cnt]=i,mu[i]=-1;
		int tp; for (j=1;j<=cnt&&(tp=i*P(j))<=n;++j)
		{
			vis[tp]=1; if (i%P(j)) mu[tp]=-mu[i]; else break;
		}
	}
}
#undef P
inline int sum(CI x,CI y)
{
	int t=x+y; return t>=mod?t-mod:t;
}
inline int quick_pow(int x,int p=mod-2,int mul=1)
{
	for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
int main()
{
	for (scanf("%d",&t),init(3000000);t;--t)
	{
		RI i,j; for (scanf("%d%d",&n,&mod),f[1]=ivf[1]=1,i=2;i<=n;++i)
		f[i]=sum(2*f[i-1]%mod,f[i-2]),ivf[i]=quick_pow(f[i]);
		for (i=0;i<=n;++i) h[i]=1; for (i=1;i<=n;++i) for (j=i;j<=n;j+=i)
		if (mu[j/i]==1) h[j]=1LL*h[j]*f[i]%mod; else
		if (mu[j/i]==-1) h[j]=1LL*h[j]*ivf[i]%mod;
		for (ans=0,i=1;i<=n;++i) h[i]=1LL*h[i]*h[i-1]%mod,ans=sum(ans,1LL*i*h[i]%mod);
		printf("%d\n",ans);
	}
	return 0;
}