BZOJ4833: [Lydsy1704月赛]最小公倍佩尔数（min-max容斥&莫比乌斯反演）（线性多项式多个数求LCM）

4833: [Lydsy1704月赛]最小公倍佩尔数

Time Limit: 8 Sec  Memory Limit: 128 MB
Submit: 240  Solved: 118
[Submit][Status][Discuss]

Description

 令(1+sqrt(2))^n=e(n)+f(n)*sqrt(2),其中e(n),f(n)都是整数，显然有(1-sqrt(2))^n=e(n)-f(n)*sqrt(2)。令g(

n)表示f(1),f(2)…f(n)的最小公倍数，给定两个正整数n和p,其中p是质数，并且保证f(1),f(2)…f(n)在模p意义

下均不为0,请计算sigma(i*g(i)),1<=i<=n.其在模p的值。

Input

第一行包含一个正整数 T ，表示有 T 组数据，满足 T≤210 。接下来是测试数据。每组测试数据只占一行，包含

两个正整数 n 和 p ，满足 1≤n≤10^6,2≤p≤10^9+7 。保证所有测试数据的 n 之和不超过 3×10^6  。

Output

对于每组测试数据，输出一行一个非负整数，表示这组数据的答案。

Sample Input

5
1 233
2 233
3 233
4 233
5 233

Sample Output

1
5
35
42
121

思路：C表示LCM，可以得到暴力：

        scanf("%d%d",&N,&P);
        A=B=C=ans=1;
        rep(i,2,N){
            int tA=A,tB=B;
            A=(tA+2*tB%P)%P; B=(tA+tB)%P; C=(ll)C/__gcd(B,C)*B;
            ans=(ans+(ll)i*C%P)%P;
        }
        printf("%d
",ans);
   

但是最小公倍数C会越来越大，而且LCM不能去%P，所以会出错。

由于A和B是线性递推的，应该会有通项公式，我们最后得到F[n]=2*F[n-1]+F[n-2]；

后面的就是参考的，证明可以看其他人的，这里只说代码需要什么，简单的说，就是：

  1，我们构造数论g[]，满足F[n]=∏g[d]（d是n的因子，即所有因子对应的g之积，注意不是之和）。

  2，Ci=g1*g2*...*gi。

所以我们只需要求g就可以了。 F[N]=g[N]*∏g[d]（d是小于N的因子），则g[N]=F[N]/∏g[d](d是小于N的因子)；所以我们可以用筛法，O(NlgN)求出g，顺便求出C。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define ll long long
using namespace std;
const int maxn=1000010;
int T,N,P,g[maxn],f[maxn];
int qpow(int a,int x)
{
    int res=1; while(x){
        if(x&1) res=1LL*res*a%P;
        a=1LL*a*a%P; x>>=1;
    } return res;
}
int main()
{
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&N,&P);
        f[0]=0;  f[1]=g[1]=1;
        rep(i,2,N) g[i]=f[i]=(2LL*f[i-1]+f[i-2])%P;
        rep(i,2,N){
            int nw=qpow(g[i],P-2);
            for(int j=i+i;j<=N;j+=i)  g[j]=1LL*g[j]*nw%P;
        }
        int ans=0,C=1;
        rep(i,1,N)
            C=(ll)C*g[i]%P,ans=(1LL*ans+1LL*i*C)%P;
        printf("%d
",ans);
    }
    return 0;
}

至于此题用到的结论。 gcd(F[x],F[y])=F[gcd(x,y)]；以及其他过程，可以参考知乎：https://www.zhihu.com/question/61218881