[BZOJ4833] 最小公倍佩尔数 题解

在这篇题解中，我会将各个部分的证明分成不同的推导过程，以达到逐一击破的效果。

引理 1：\(f(n)=2f(n-1)+f(n-2)\)

我的证明挺繁琐的，过程如下：

\[(1+\sqrt 2)^{n-2}=e(n-2)+f(n-2)\sqrt 2 \]

\[(1+\sqrt 2)^{n-1}=e(n-1)+f(n-1)\sqrt 2 \]

\[(1+\sqrt 2)^{n-1}=(1+\sqrt 2)^{n-2}(1+\sqrt 2) \]

\[=(e(n-2)+2f(n-2))+(e(n-2)+f(n-2))\sqrt 2 \]

\[e(n-1)=e(n-2)+2f(n-2),f(n-1)=e(n-2)+f(n-2) \]

\[f(n)=e(n-1)+f(n-1)=2(e(n-2)+f(n-2))+f(n-2) \]

\[f(n)=2f(n-1)+f(n-2) \]

那么递推式有了，如何和 \(\operatorname{lcm}\) 扯上关系呢？

引理 2：\(f(n+m)=f(n+1)f(m)+bf(n)f(m-1)\)

这个引理的完整描述过程如下：

如果 \(f(0)=0,f(1)=1,f(n)=af(n-1)+bf(n-2)\)，且 \(\gcd(a,b)=1\)，则有 \(\gcd(f(x),f(y))=f(\gcd(x,y))\)。

考虑这个东西长得和个兔子数列似的，于是给他弄个转移矩阵。容易发现有：

\[\begin{bmatrix}f(n+1)\ \ f(n)\end{bmatrix}=\begin{bmatrix}a\ \ 1\\b\ \ 0\end{bmatrix}\begin{bmatrix}f(n)\ \ f(n-1)\end{bmatrix} \]

这个转移矩阵是不平凡的：

\[\begin{bmatrix}f(n+1)\ \ f(n)\\bf(n)\ \ bf(n-1)\end{bmatrix}\begin{bmatrix}a\ \ 1\\b\ \ 0\end{bmatrix}=\begin{bmatrix}af(n+1)+bf(n)\ \ f(n+1)\\bf(n+1)\ \ bf(n)\end{bmatrix}=\begin{bmatrix}f(n+2)\ \ f(n+1)\\bf(n+1)\ \ bf(n)\end{bmatrix} \]

\[\begin{bmatrix}a\ \ 1\\b\ \ 0\end{bmatrix}^1=\begin{bmatrix}f(2)\ \ f(1)\\bf(1)\ \ bf(0)\end{bmatrix},\begin{bmatrix}a\ \ 1\\b\ \ 0\end{bmatrix}^n=\begin{bmatrix}f(n+1)\ \ f(n)\\bf(n)\ \ bf(n-1)\end{bmatrix} \]

那么就有：

\[\begin{bmatrix}f(n+m+1)\ \ f(n+m)\end{bmatrix}=\begin{bmatrix}a\ \ 1\\b\ \ 0\end{bmatrix}^m\begin{bmatrix}f(n+1)\ \ f(n)\end{bmatrix} \]

\[=\begin{bmatrix}f(m+1)\ \ f(m)\\bf(m)\ \ bf(m-1)\end{bmatrix}\begin{bmatrix}f(n+1)\ \ f(n)\end{bmatrix} \]

\[=\begin{bmatrix}f(n+1)f(m+1)+bf(m)f(n)\ \ f(n+1)f(m)+bf(n)f(m-1)\end{bmatrix} \]

引理 1 得证。

引理 3：\(\gcd(f(x),f(y))=f(\gcd(x,y))\)

显然有：

\[\gcd(f(n+m),f(n))=\gcd(f(n+1)f(m)+bf(n)f(m-1),f(n)) \]

\[\gcd(f(n+m),f(n))=\gcd(f(n+1)f(m),f(n)) \]

若 \(\gcd(f(n-1),f(n-2))=1\)，则有：

\[\gcd(f(n),f(n-1))=\gcd(af(i-1)+bf(n-2),f(n-1))=1 \]

由于 \(\gcd(f(1),f(0))=1\)，所以 \(\gcd(f(n),f(n-1))=1\)。

所以 \(\gcd(f(n+m),f(n))=\gcd(f(m),f(n))\)。

我们从辗转相除的角度化简，可以得到：

\[\gcd(f(n+m),f(n))=\gcd(f(\gcd(n+m,n)),f(0)) \]

\[=f(\gcd(n+m,n))=f(\gcd(m,n))=\gcd(f(m),f(n)) \]

引理 4：\(\operatorname{lcm}(S)=\prod\limits_{T\subset S}\gcd(T)^{(-1)^{|T|-1}}\)

特殊说明一下，\(S,T\) 为非空集合。

注意力惊人的注意到：

\[\operatorname{lcm}(S)=\prod_{p\in prime}p^{\max\limits_{i\in S}cnt_{i,p}} \]

\[\gcd(S)=\prod_{p\in prime}p^{\min\limits_{i\in S}cnt_{i,p}} \]

其中 \(cnt_{i,p}\) 表示 \(i\) 中 \(p\) 这个质因子的个数。

想到最值反演。那么有：

\[\operatorname{lcm}(S)=\prod_{p\in prime}p^{\max\limits_{i\in S}cnt_{i,p}} \]

\[\operatorname{lcm}(S)=\prod_{p\in prime}p^{\sum\limits_{T\subseteq S}\min\limits_{i\in T}cnt_{i,p}(-1)^{|T|-1}} \]

\[\operatorname{lcm}(S)=\prod_{T\subseteq S}(\prod_{p\in prime}p^{\min\limits_{i\in T}cnt_{i,p}})^{(-1)^{|T|-1}} \]

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

\(g(n)\) 的第一推导

设 \(U=\{1,2,\dots,n\}\)，则有：

\[g(n)=\prod_{T\subseteq U}(\gcd\limits_{i\in T}f(i))^{(-1)^{|T|-1}} \]

\[g(n)=\prod_{T\subseteq U}f(\gcd(T))^{(-1)^{|T|-1}} \]

看到 \(\gcd\)，想到莫比乌斯反演：

\[g(n)=\prod_{i=1}^nf(i)^{\sum\limits_{T\subseteq U}[\gcd(T)=i](-1)^{|T|-1}} \]

记 \(h(i)=\sum\limits_{T\subseteq U}[\gcd(T)=i](-1)^{|T|-1},H(i)=\sum\limits_{i|d}h(d)\)，我们就要开启下一部分了。

引理 5：\(H(i)=1\)

设 \(S_i=\{1,2,\dots,\lfloor\dfrac ni\rfloor\}\)（显然 \(|S|>0\)），\(T\) 仍然是非空集合，则有：

\[H(i)=\sum_{i|d}\sum_{T\subseteq U}[\gcd(T)=d](-1)^{|T|-1} \]

\[=\sum_{T\subseteq U}[i|\gcd(T)](-1)^{|T|-1} \]

\[=\sum_{T\subseteq S_i}(-1)^{|T|-1}=-\sum_{T\subseteq S_i}(-1)^{|T|} \]

\[=-(\sum_{i=0}^{|S_i|}\binom{|S|}i(-1)^i-1) \]

\[=-(1-1)^{|S_i|}+1=1 \]

\(g(n)\) 的最终推导

根据莫比乌斯反演，有：

\[h(i)=\sum_{i|d}\mu(\frac di)H(d)=\sum_{i|d}\mu(\frac di) \]

带回 \(g(n)\) 中，得：

\[g(n)=\prod_{i=1}^nf(i)^{\sum\limits_{i|d}\mu(\frac di)} \]

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

很好算。时间复杂度 \(O(\sum n\log n)\)。

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e6+5;
int t,n,p,mu[N],f[N],s[N];
int qpow(int x,int y){
    int re=1;
    while(y){
        if(y&1) re=re*x%p;
        x=x*x%p,y>>=1;
    }return re;
}signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>t,mu[1]=1;
    for(int i=1;i<N;i++)
        for(int j=i*2;j<N;j+=i) mu[j]-=mu[i];
    while(t--){
        cin>>n>>p,f[1]=s[1]=1;
        for(int i=2;i<=n;i++)
            f[i]=(f[i-1]*2+f[i-2])%p,s[i]=1;
        for(int i=1;i<=n;i++){
            int x=f[i],y=qpow(x,p-2);
            for(int j=i,k=1;j<=n;j+=i,k++)
                s[j]=s[j]*(!mu[k]?1:mu[k]>0?x:y)%p;
        }int re=0,g=s[1];
        for(int i=1;i<=n;)
            re=(re+g*i)%p,g=g*s[++i]%p;
        cout<<re<<"\n";
    }return 0;
}