题解：P10663 BZOJ4833 最小公倍佩尔数

· · 题解

Solution

容易发现 f(n) = \dfrac{(1+\sqrt 2)^n - (1 - \sqrt 2)^n}{2 \sqrt 2}。

经常学数竞的同学容易发现，上面两个幂的底数都是代数整数，我们很容易整理出一个关于 f 的二阶齐次线性递推：f(n) = 2 f(n-1) + f(n-2)，其中 n \ge 3。

有一个常用的 \rm GCD - LCM 反演（实际上就是 \rm Min-Max 容斥）：

\text{lcm}(a_1,a_2,\dots,a_n) = \prod_{S \subseteq \{1,2,\dots,n\}} (\gcd_{u \in S} a_u)^{(-1)^{|S|-1}}

对每个质数分开考虑就是 \rm Min-Max 容斥。这还有一道例题。

考虑如何刻画 \gcd_{u \in S} f(u)？

斐波那契数列的递推式和 f 很像，而且有 fib(\gcd(a,b)) = \gcd(fib(a),fib(b)) 这一经典重要性质。我们推测 f 也有这个性质。

考虑递推的转移矩阵有 \left[\begin{matrix}2 &1 \\1 &0 \end{matrix} \right]^n = \left[\begin{matrix}f(n+1) & f(n) \\f(n) & f(n-1) \end{matrix} \right]，那么应当有

\left[\begin{matrix} f(a) & f(a-1) \end{matrix}\right]\left[\begin{matrix}2 &1 \\1 &0 \end{matrix} \right]^{b-a} = \left[ \begin{matrix} f(b) & f(b-1)\end{matrix}\right]

即

f(b) = f(a)f(b-a-1) + f(a-1) f(b-a)

则

\gcd(f(a),f(b)) = \gcd(f(a),f(a-1)f(b-a))

根据归纳法易得 \gcd(f(a),f(b))=1，有 \gcd(f(a),f(b)) = \gcd(f(a),f(b-a))。

根据辗转相除法，即得。

一般的，形如 f(n) = k f(n-1) + f(1)，k \in \mathbb Z 的递推，若 f(1) = 1，f(2) = k，则 \gcd(f(a),f(b)) = f(\gcd(a,b)) 就成立。

下面开始劲爆莫比乌斯反演。

套用上述公式和莫比乌斯反演的公式易得

g(n)=\prod_{d=1}^n f(d)^{\sum_{t=1}^{\lfloor\frac{n}{d}\rfloor} \mu(t)}

考虑 \dfrac{g(n)}{g(n-1)}，显然就等于

\prod_{d \mid n} f(d)^{\mu(\frac{n}{d})}

枚举所有使得 \mu(\dfrac{n}{d}) 有值的 d 即可（其实可以枚举 n 的所有因子，但是这样做会使得常数更大）

复杂度 O(n \ln n)，常数略大。

#include<bits/stdc++.h>
#define int long long
#define ffor(i,a,b) for(int i=(a);i<=(b);i++)
#define roff(i,a,b) for(int i=(a);i>=(b);i--)
using namespace std;
const int MAXN=1e6+10;
int tc,n,p,pw[MAXN],f[MAXN],inv[MAXN],mu[MAXN],flg[MAXN],cnt[MAXN],g[MAXN];
vector<int> vc;
int qpow(int base,int p,int MOD) {
    int ans=1;
    while(p) {
        if(p&1) ans=ans*base%MOD;
        base=base*base%MOD,p>>=1;
    }
    return ans;
}
void init(int mx) {
    mu[1]=1;
    ffor(i,2,mx) {
        if(!flg[i]) vc.push_back(i),mu[i]=-1;
        for(auto v:vc) {
            if(i*v>mx) break ;
            flg[i*v]=1;
            if(i%v==0) break ;
            mu[i*v]=-mu[i];
        }
    }
    return ;
}
signed main() {
    ios::sync_with_stdio(false),cin.tie(0),cout.tie(0);
    init(1000000);
    vc.clear();
    ffor(i,1,1000000) if(mu[i]) vc.push_back(i);
    cin>>tc;
    while(tc--) {
        cin>>n>>p;
        f[1]=1,f[2]=2;
        ffor(i,3,n) f[i]=(2*f[i-1]+f[i-2])%p;
        ffor(i,1,n) inv[i]=qpow(f[i],p-2,p);
        ffor(i,0,n) g[i]=1;
        ffor(d,1,n) {
            int lim=n/d;
            for(auto id:vc) {
                if(id>lim) break ;
                int Id=id*d;
                if(mu[id]==1) g[Id]=g[Id]*f[d]%p;
                else g[Id]=g[Id]*inv[d]%p;  
            }
        }
        ffor(i,1,n) g[i]=g[i-1]*g[i]%p;
        int ans=0;
        ffor(i,1,n) ans=(ans+i*g[i])%p;
        cout<<ans<<'\n';
    }
    return 0;
}

本题完整的输入，可以用来调试：Link