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\n",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\n",ans);
    }
    return 0;
}

至於此題用到的結論。 gcd(F[x],F[y])=F[gcd(x,y)]；可以參考知乎：https://www.zhihu.com/question/61218881