题目地址ZOJ3707 Caluculate Prime S
题目中
S(n)
S
(
n
)
表示在
1,2,3,...,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)
S
(
n
)
=
F
i
b
(
n
+
2
)
又有一个重要的性质:
GCD(Fib(n),Fib(m))=Fib(GCD(n,m))
G
C
D
(
F
i
b
(
n
)
,
F
i
b
(
m
)
)
=
F
i
b
(
G
C
D
(
n
,
m
)
)
其中,
GCD(x,y)
G
C
D
(
x
,
y
)
表示
x,y
x
,
y
的最大公约数
由此可以知道如果一个
S(n)
S
(
n
)
与之前的所有
S(i)
S
(
i
)
都互质的话,那么
GCD(n,i)=1或GCD(n,i)=2
G
C
D
(
n
,
i
)
=
1
或
G
C
D
(
n
,
i
)
=
2
,因为只有这样
Fib(GCD(i,m))=1
F
i
b
(
G
C
D
(
i
,
m
)
)
=
1
所以就可以构造一个特殊的质数数组,分别是{3,4,5,7…},总个数为1000000,所以用埃式筛法数组的大小大约为16000000.其中4是特殊的数,因为
4−2=2
4
−
2
=
2
为一个质数.
然后根据以上的定理和推论,我们可以利用矩阵快速幂的方法求得大于 Fib(P(K)) F i b ( 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;
}