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)=∏T⊆Sgcd(f(T))(−1)|T|+1

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

一个经典思路考虑构造数列h,f(n)=∏d|nh(d)

把这个带进去就有

g(S)=∏T⊆S∏d|gcd(T)h(d)(−1)|T|+1

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

g(S)=∏nd=1h(d)∏T⊆S(−1)∏i∈T(−1)([d|i])

g(S)=∏nd=1h(d)1−∏i∈S(1−[d|i])

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

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

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

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

那么h(d)等于啥呢？

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

f(n)=∏d|nh(d),h(n)=∏d|nf(d)μnd

问题完美解决，时间复杂度nlog(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;    }}