4833: [Lydsy2017年4月月赛]最小公倍佩尔数数论变换-优快云博客

先吐槽一下：唐老师下标全崩了，最重要的后面的一坨推导全都挂了qwq，我推到倒数第二步然后和题解一对发现我好想全推错了233333？最后还是找的A了的zyz才弄明白他写的是什么

这题和一道斐波那契公倍数比较像，先说这个吧

先考虑给你的那个式子，你上面减下面再移个项就变成通项公式了

然后通过特征根/观察/打表你get了这个递推式是 $f(i) = 2f(i-1) + f(i-2)$

观察这个递推式的性质，首先我们能够得到 $\gcd(f(i),f(i-1))=1$ ,原因是 $f(i)\% f(i-1)=f(i-2)$ 然后无限递归一下就知道他们等于 $\gcd(f(1),f(2))=1$

然后考虑 $f(n+m) = f(n-1)f(m)+f(n)f(m+1)$ ，我证这个的时候非常傻，用通项弄得，实际上考虑矩阵乘法的递推式就可以了，注意到结合律之后拆出中间一个就好了

然后我们get了一个非常关键的结论是 $\gcd(f(i),f(j))=f(\gcd(i,j)))$

证明的话考虑一种让 $i>j$ 那么考虑令 $j=m$ , $i=n+m$ ，

此时 $\gcd(f(n+m),f(m)) = \gcd(f(n-1)f(m)+f(n)f(m+1), f(m))$

前面那项显然可以直接去掉了，因为是 $f(m)$ 的倍数

那么剩下的是 $\gcd(f(n)f(m+1), f(m))$

注意到 $\gcd(f(m+1),f(m))=1$ 所以乘不乘他都一样了

所以他等于 $\gcd(f(n),f(m))$ 递归即可

现在可以考虑 $g(i)=lcm(f(1),f(2),...,f(i))$

考虑 $g(S) = lcm(f(S))$ ，其中S是一个集合从1到n

然后我们子集反演一下就得到了 $g(S)=\prod_{T\subseteq S}\gcd(f(T))^{(-1)^{|T|+1}}$

注意到 $\gcd(f(T)) = f(\gcd(T))$

一个经典思路考虑构造数列 ${h}$ , $f(n)=\prod_{d|n}h(d)$

把这个带进去就有

$g(S)=\prod_{T\subseteq S}\prod_{d|\gcd(T)}h(d)^{(-1)^{|T|+1}}$

考虑直接枚举一个 $h(d)$ 的贡献

$g(S)=\prod_{d=1}^n h(d)^{\prod_{T\subseteq S}(-1) \prod_{i\in T}(-1)([d|i])}$

$g(S)=\prod_{d=1}^n h(d)^{1-\prod_{i\in S}(1-[d|i])}$

这个是怎么来的呢？考虑这样一个事情，对于一个元素来说他被选择而且是 $d$ 的倍数的话他的贡献是 $-1$ 否则的话贡献就是1

所以一个元素总的贡献就是 $[1-[d|i]]$ ，我们枚举了所有子集，由于乘法原理每个数字之间一点关系都没有所以把所有的贡献乘在一起就好啦

但是注意到前面是符号所以是减去这个贡献还有空集，空集的贡献就是1，所以就有了上面的那个式子了

然后考虑后面那个东西在 $n=1$ 的时候等于0，此时 $g(1)=1$ ,当 $n!=1$ 的时候总存在一个 $d$ 不能作为所有 $i$ 的约数那么就等于0，此时 $g(S)=\prod_{d=1}^n h(d)$

那么 $h(d)$ 等于啥呢？

通过广义莫比乌斯反演我们可以知道他等于同等运算在mu意义下的逆运算

$f(n)=\prod_{d|n}h(d), h(n)=\prod_{d|n}f(d)^{\mu{\frac{n}{d}}}$

问题完美解决，时间复杂度 $n\log(n)$

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
using namespace std;
const int N = 1e6+5;
int prime[N], cnt, mu[N];
bool F[N];
typedef long long LL;

inline int pow(int a,int b,int mod) {
    LL res = 1;
    for(;b;b>>=1,a=(LL)a*a%mod)
        if(b&1)
            res = res * a % mod;
    return res;
}

inline void init() {
    const int Max = 1e6;
    mu[1] = 1;
    for(int i=2;i<=Max;++i) {
        if(!F[i]) mu[i] = -1, prime[++cnt] = i;
        for(int j=1;prime[j]*i<=Max;++j) {
            F[i*prime[j]] = 1;
            if(i%prime[j]==0) break;
            mu[i*prime[j]] = -mu[i];
        } 
    }
}

int f[N], h[N], g[N], inv[N];

int main() {
    int T;
    cin >> T;
    init();
    while(T--)  {
        int Max, mod;
        cin >> Max >> mod;
        f[1] = 1, f[2] = 2;
        register int i; int x;
        for(i=3;i<=Max;++i) f[i] = (2ll * f[i-1] + f[i-2]) % mod;
        for(i=1;i<=Max;++i) h[i] = 1;
        for(i=1;i<=Max;++i) inv[i] = pow(f[i], mod-2,mod);
        for(int d=1;d<=Max;++d) {
            for(i=1;i*d<=Max;++i) {
                x = i * d;
                if(mu[i] == 1) 
                    h[x] = (LL) h[x] * f[d] % mod;
                else if(mu[i] == -1) 
                    h[x] = (LL) h[x] * inv[d] % mod;
            }
        }
        for(i=2;i<=Max;++i) h[i] = (LL) h[i] * h[i-1] % mod;
        for(i=1;i<=Max;++i) g[i] = h[i];
        LL ans = 0;
        for(i=1;i<=Max;++i) ans = (ans + (LL) g[i] * i % mod) % mod;
        cout << ans << endl;
    }
}