[数论]【ZOJ3707】Calculate Prime S

题目地址ZOJ3707 Caluculate Prime S
题目中$S(n)$表示在${1,2,3,...,n}$这个集合中不存在连续项的子集的个数，其中包括空集.
所以可以列举几个，比如
S(1)=2（{空集}，{1}）
S(2)=3（{空集}，{1}，{2}）
S(3)=5（{空集}，{1}，{2}，{3}，{1,3}）
……
依次列举之后可以发现这是斐波那契数列，且$S(n)=Fib(n+2)$
又有一个重要的性质：
$GCD(Fib(n),Fib(m))=Fib(GCD(n,m))$
其中，$GCD(x,y)$表示$x,y$的最大公约数

由此可以知道如果一个$S(n)$与之前的所有$S(i)$都互质的话，那么$GCD(n,i)=1或GCD(n,i)=2$，因为只有这样$Fib(GCD(i,m))=1$
所以就可以构造一个特殊的质数数组，分别是{3,4,5,7…}，总个数为1000000，所以用埃式筛法数组的大小大约为16000000.其中4是特殊的数，因为$4-2=2$为一个质数.

然后根据以上的定理和推论，我们可以利用矩阵快速幂的方法求得大于$Fib(P(K))$中可以整除X的那一项，由于X的大小为【3，100】，所以只需要处理至多100次即可。其中 P P <script type="math/tex" id="MathJax-Element-26">P</script>为之前的质数数组，K为题意中的第K大。

找到以后还要利用同余公式转化最后的求解式
由同余公式(a/b)%c=（a%（b*c)）/b可以把最后的结果转化为S(i)/X%M=S(i)%(X*M)/X

比赛的时候能够做出来的也绝对是大佬了

#include<bits/stdc++.h>
#define CLR(x) memset(x,0,sizeof(x))
#define ll long long
using namespace std;

const int maxn = 1.6e7+100;

int p[2000010];
int pnum;
bool vis[maxn];
int K,X,M;

struct matrix{
    ll mat[3][3];
    matrix(){CLR(mat);}
    ll *operator [](int x){return mat[x];}
};

matrix mul(matrix &a,matrix &b,int mod){
    matrix c;
    for(int i=1;i<=2;i++){
        for(int j=1;j<=2;j++){
            for(int k=1;k<=2;k++){
                c[i][j] += a[i][k]*b[k][j];
                c[i][j] %= mod;
            }
        }
    }
    return c;
}

void get_prime(){
    pnum=1;
    memset(vis,false,sizeof(vis));
    for(int i=2;i<16000001;i++){
        if(!vis[i]){
            p[pnum++]=i;
            for(int j=i*2;j<16000001;j+=i){
                vis[j]=true;
            }
        }
    }
    p[1]=3;
    p[2]=4;
}

matrix qmp(matrix a,int n,int mod){
    matrix c;
    for(int i=1;i<=2;i++)c[i][i]=1;
    while(n){
        if(n&1) c = mul(c,a,mod);
        a=mul(a,a,mod);
        n/=2;
    }
    return c;
}

int main()
{
    get_prime();
    int T;
    scanf("%d",&T);
    while(T--){
        int ansi;
        int temp;
        matrix a,c;
        a[1][1]=a[1][2]=a[2][1]=1;
        scanf("%d%d%d",&K,&X,&M);
        for(int i=p[K];;i++){
            c=qmp(a,i-2,X);
            temp=(c[1][1]+c[1][2])%X;
            if(temp==0){
                ansi=i;break;
            }
        }
        c=qmp(a,ansi-2,M*X);
        temp=(c[1][1]+c[1][2])%(M*X);
        printf("%d\n",temp/X);
    }
    return 0;
}