BZOJ #4833. [Lydsy2017年4月月赛]最小公倍佩尔数

题面传送门

神仙题，做了一个下午。

Section 1

首先我们要对这个东西观察出一点性质来。

按照题目中讲的，有递推式\((e_{n-1}+f_{n-1}\sqrt 2)(\sqrt 2+1)=e_n+f_n\sqrt 2\)，于是展开可以得到\(e_n=e_{n-1}+2f_{n-1},f_n=e_{n-1}+f_{n-1}\)。

递推式的形式也不好处理，观察发现\(e_n=f_n+f_{n-1}\)，故\(f_n=2f_{n-1}+f_{n-2}\)，到此我们完成了初步转化。

Section 2

对于\(f_n=af_{n-1}+f_{n-2}\)形式，我们有以下引理：

引理一：\(f_{n+m}=f_{n}f_{m+1}+f_{m}f_{n+1}\)

引理二：\(\gcd(f_n,f_{n+1})=1\)

证明：考虑辗转相减，\(\gcd(f_{n+1},f_n)=\gcd(f_{n-1}+af_{n},f_n)=\gcd(f_{n-1},f_n)\)即可归纳。

引理三：\(\gcd(f_n,f_m)=f_{\gcd(n,m)}\)

证明：我们设\(m=n+s\)，则\(\gcd(f_n,f_m)=\gcd(f_n,f_nf_{s+1}+f_{n+1}f_s)=\gcd(f_n,f_s)=\gcd(f_n,f_{m-n})=\gcd(f_n,f_{m \bmod n})\)即可归纳。

Section 3

我们要求的式子长这样：\(g_n=\operatorname{lcm}_{i=1}^{n}{f_i}\)。

仿照min-max容斥，上gcd-lcm容斥：\(g_n=\prod\limits_{T\subset S}{\gcd(T)^{(-1)^{|T|+1}}}\)。

由引理三，我们有\(g_n=\prod\limits_{T\subset S}{f_{\gcd(T)}^{(-1)^{|T|+1}}}\)

枚举\(gcd(T)\)，有\(g_n=\prod\limits_{d=1}^{n}{f_{d}^{\sum\limits_{T\subset S}{[\gcd(T)==d](-1)^{|T|+1}}}}\)。

看到gcd容易想到上莫反：\(g_n=\prod \limits_{d=1}^{n}{f_d^{\sum\limits_{k=1}^{\frac{n}{d}}{\mu(k)\sum\limits_{T\subset S}{(-1)^{|T|+1}[kd\mid T]}}}}\)

根据二项式定理，后面所构成的集合大小至少为\(1\)，因此后面一坨为\(1\)，故原式为\(g_n=\prod \limits_{d=1}^{n}{f_d^{\sum\limits_{k=1}^{\frac{n}{d}}{\mu(k)}}}\)

这个东西差分后暴力跳调和级数即可，时间复杂度\(O(n\log n)\)

code:

#include<bits/stdc++.h>
#define I inline
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define abs(x) ((x)>0?(x):-(x))
#define ll long long
#define db double
#define lb long db
#define N (1000000+5)
#define M (6000000+5)
#define K (1500+5)
#define mod 998244353
#define Mod (mod-1)
#define eps (1e-5)
#define ull unsigned ll
#define it iterator
#define Gc() getchar() 
#define Me(x,y) memset(x,y,sizeof(x))
#define Mc(x,y) memcpy(x,y,sizeof(x))
#define d(x,y) ((k+1)*(x)+(y))
#define R(n) (1ll*rand()*rand()%(n)+1)
#define Pc(x) putchar(x)
#define LB lower_bound
#define UB upper_bound
#define PB push_back
using namespace std;const int INF=N-5;
int n,p,T,Fl[N],pr[N],ph,mu[N],F[N];ll Ans,W[N],frc[N],Inv[N];
I ll mpow(ll x,int p,int y){ll Ans=1;while(y) y&1&&(Ans=Ans*x%p),y>>=1,x=x*x%p;return Ans;}
I void Solve(){
	int i,j;scanf("%d%d",&n,&p);for(i=1;i<=n;i++) W[i]=1;F[0]=0;F[1]=1;for(i=2;i<=n;i++) F[i]=(2ll*F[i-1]+F[i-2])%p;Ans=0;for(frc[0]=i=1;i<=n;i++) frc[i]=frc[i-1]*F[i]%p;
	Inv[n]=mpow(frc[n],p,p-2);for(i=n-1;~i;i--) Inv[i]=Inv[i+1]*F[i+1]%p;for(i=1;i<=n;i++) Inv[i]=frc[i-1]*Inv[i]%p;for(i=1;i<=n;i++) for(j=1;j*i<=n;j++) mu[j]&&(W[i*j]=(~mu[j]?W[i*j]*F[i]%p:W[i*j]*Inv[i]%p));
	Ans=0;W[0]=1;for(i=1;i<=n;i++) W[i]=W[i]*W[i-1]%p,Ans+=W[i]*i%p;printf("%lld\n",Ans%p);
}
int main(){
	freopen("1.in","r",stdin);
	int i,j;mu[1]=1;for(i=2;i<=INF;i++){!Fl[i]&&(pr[++ph]=i,mu[i]=-1);for(j=1;j<=ph&&i*pr[j]<=INF;j++) {Fl[i*pr[j]]=1;if(i%pr[j]==0) break;mu[i*pr[j]]=-mu[i];}}
	scanf("%d",&T);while(T--) Solve();
}