题意
已知
x
0
,
x
1
,
a
,
b
(
1
≤
1
0
9
)
x_0,x_1,a,b(1\leq10^9)
x0,x1,a,b(1≤109),求
x
n
=
a
x
n
−
1
+
b
x
n
−
2
(
m
o
d
m
)
x_n=ax_{n-1}+bx_{n-2}(\ mod \ m)
xn=axn−1+bxn−2( mod m)
1
≤
n
≤
1
0
1
0
6
,
1
0
9
<
m
≤
2
×
1
0
9
1\leq n\leq 10^{10^6},10^9<m\leq2\times10^9
1≤n≤10106,109<m≤2×109
思路
由于题目中给的
n
n
n特别大,所以直接做矩阵快速幂不行,当然可以用10进制快速幂做。这里介绍一下找循环节的办法。
显然的在模意义下广义斐波那契必有循环节,那么循环节和什么有关呢,对于广义斐波那契数列
{
f
[
1
]
=
c
f
[
2
]
=
d
f
[
n
]
=
a
f
[
n
−
1
]
+
b
f
[
n
−
2
]
,
n
>
2
\left\{\begin{matrix} f[1]=c\\ f[2]=d\\ f[n]=af[n-1]+bf[n-2],n>2 \end{matrix}\right.
⎩⎨⎧f[1]=cf[2]=df[n]=af[n−1]+bf[n−2],n>2
我们可以尝试固定
a
,
b
a,b
a,b改变
c
,
d
c,d
c,d,会发现数列的循环节其实没有变化的,要说证明的话我们可以构造一下矩阵
[
f
n
f
n
−
1
]
=
[
a
b
1
0
]
[
f
n
−
1
f
n
−
2
]
(
m
o
d
m
)
\begin{bmatrix} f_n\\ f_{n-1} \end{bmatrix}=\begin{bmatrix} a&b \\ 1&0 \end{bmatrix}\begin{bmatrix} f_{n-1}\\ f_{n-2} \end{bmatrix}(\ mod \ m)
[fnfn−1]=[a1b0][fn−1fn−2]( mod m)
而循环节就是在模意义下
[
a
b
1
0
]
k
=
[
1
0
0
1
]
(
m
o
d
m
)
\begin{bmatrix} a&b \\ 1&0 \end{bmatrix}^k=\begin{bmatrix} 1&0 \\ 0&1 \end{bmatrix}(\ mod \ m )
[a1b0]k=[1001]( mod m),这里的
k
k
k就是循环节,所以实际上只会和
a
,
b
a,b
a,b有关和
c
,
d
c,d
c,d是无关的
接下来有这么一个定理若模数为
m
=
p
1
a
1
p
2
a
2
⋯
p
k
a
k
m=p_{1}^{a_1}p_{2}^{a_2}\cdots p_{k}^{a_k}
m=p1a1p2a2⋯pkak,
g
(
m
)
g(m)
g(m)表示模
m
m
m下循环节,那么就有
g
(
m
)
=
l
c
m
(
g
(
p
1
a
1
)
,
g
(
p
2
a
2
)
,
⋯
 
,
g
(
p
k
a
k
)
)
g(m)=lcm(g(p_{1}^{a_1}),g(p_{2}^{a_2}),\cdots, g(p_{k}^{a_k}))
g(m)=lcm(g(p1a1),g(p2a2),⋯,g(pkak))
l
c
m
lcm
lcm为最小公倍数
并且对于
g
(
p
k
)
g(p^k)
g(pk),
p
p
p为质数,有
g
(
p
k
)
=
g
(
p
)
∗
p
k
−
1
g(p^k)=g(p)*p^{k-1}
g(pk)=g(p)∗pk−1
那么我们现在就可以将给的模数
m
m
m唯一分解了然后对于其每一个质因数求解循环节,现在考虑对于质因数
p
p
p怎么求解循环节
我们可以对其的通项下手进行分析,对于这种2阶线性常系数齐次递推关系,我们可以用特征方程的方式来求解其的通项即
x
2
−
a
x
−
b
=
0
x^2-ax-b=0
x2−ax−b=0
考虑用求根公式来求解
Δ
=
a
2
+
4
b
\Delta=a^2+4b
Δ=a2+4b,考虑
Δ
\Delta
Δ在
(
m
o
d
p
)
(\ mod\ p)
( mod p)的情况
p
>
2
p>2
p>2
-
Δ
≡
0
(
m
o
d
p
)
\Delta\equiv0(\ mod \ p)
Δ≡0( mod p),即有重根的情况,就有
f
[
n
]
=
A
n
(
B
n
+
C
)
f[n]=A^n(Bn+C)
f[n]=An(Bn+C)
根据求根公式 A = a + Δ 2 = a 2 A=\frac{a+\sqrt{\Delta}}{2}=\frac{a}{2} A=2a+Δ=2a
有 A = a ∗ i n v 2 ( m o d p ) A=a*inv2(\ mod\ p) A=a∗inv2( mod p)那我们就可以用费马小定理来得到循环节了
A p − 1 ≡ A 0 ≡ 1 ( m o d p ) A^{p-1}\equiv A^0\equiv 1(\ mod \ p) Ap−1≡A0≡1( mod p)
又有 p ≡ 0 ( m o d p ) p\equiv 0(\ mod \ p) p≡0( mod p)
所以循环节必是 l c m ( p − 1 , p ) = p ( p − 1 ) lcm(p-1,p)=p(p-1) lcm(p−1,p)=p(p−1)的倍数 -
Δ
̸
≡
0
(
m
o
d
p
)
\Delta\not\equiv 0(\ mod \ p)
Δ̸≡0( mod p),即有两个不同的根的情况有
f
[
n
]
=
A
B
n
+
C
D
n
f[n]=AB^n+CD^n
f[n]=ABn+CDn
根据求根公式我们求的是 Δ \sqrt{\Delta} Δ的所以我们得考虑 Δ \Delta Δ是否是模 p p p意义下的二次剩余
若 Δ \Delta Δ为二次剩余即 x 2 ≡ Δ ( m o d p ) x^2\equiv\Delta(\ mod \ p) x2≡Δ( mod p), x x x有解,那么表示在模 p p p意义下 Δ \sqrt{\Delta} Δ可以表示,那么
B = a + Δ 2 , D = a + Δ 2 B=\frac{a+\sqrt{\Delta}}{2},D=\frac{a+\sqrt{\Delta}}{2} B=2a+Δ,D=2a+Δ在模意义下都是整数,根据费马小定理或者是扩展欧拉降幂, B n , C n B^n,C^n Bn,Cn在模 p p p的情况下循环节就是 p − 1 p-1 p−1的倍数
若 Δ \Delta Δ不为二次剩余即 x 2 ≡ Δ ( m o d p ) x^2\equiv\Delta(\ mod \ p) x2≡Δ( mod p), x x x无解,那么我们可以考虑将 B p + 1 B^{p+1} Bp+1展开
B p + 1 = ( a + Δ 2 ) p + 1 = 1 2 p + 1 ∑ i = 0 p + 1 C p + 1 i a p + 1 − i Δ i B^{p+1}=(\frac{a+\sqrt{\Delta}}{2})^{p+1}=\frac{1}{2^{p+1}}\sum_{i=0}^{p+1}C_{p+1}^ia^{p+1-i}\sqrt{\Delta}^i Bp+1=(2a+Δ)p+1=2p+11i=0∑p+1Cp+1iap+1−iΔi
≡ 1 4 ( Δ p + 1 + a Δ p + 0 + ⋯ + 0 + a Δ + a 2 ) ( m o d p ) \equiv\frac{1}{4}(\sqrt{\Delta}^{p+1}+a\sqrt{\Delta}^p+0+\cdots+0+a\sqrt{\Delta}+a^2)(\ mod \ p) ≡41(Δp+1+aΔp+0+⋯+0+aΔ+a2)( mod p)
上式根据组合数 C p + 1 i = ( p + 1 ) ( p ) ( p − 1 ) ⋯ 1 ∗ 2 ∗ 3 ∗ ⋯ C_{p+1}^{i}=\frac{(p+1)(p)(p-1)\cdots}{1*2*3*\cdots} Cp+1i=1∗2∗3∗⋯(p+1)(p)(p−1)⋯,所以中间部分都为 0 0 0, 1 2 p + 1 \frac{1}{2^{p+1}} 2p+11根据逆元可以变为 1 4 \frac{1}{4} 41
再根据欧拉准则
Δ p − 1 2 ≡ − 1 ( m o d p ) \Delta^{\frac{p-1}{2}}\equiv -1(\ mod\ p) Δ2p−1≡−1( mod p)
带入上式有
B p + 1 ≡ 1 4 ( − Δ − a Δ + a Δ + a 2 ) ( m o d p ) B^{p+1}\equiv \frac{1}{4}(-\Delta-a\sqrt{\Delta}+a\sqrt{\Delta}+a^2)(\ mod \ p) Bp+1≡41(−Δ−aΔ+aΔ+a2)( mod p)
有 Δ = a 2 + 4 b \Delta=a^2+4b Δ=a2+4b,所以有
B p + 1 ≡ − b ( m o d p ) B^{p+1}\equiv -b(\ mod \ p) Bp+1≡−b( mod p)
再根据费马小定理
( B p + 1 ) p − 1 ≡ ( − b ) p − 1 ≡ 1 ( m o d p ) (B^{p+1})^{p-1}\equiv (-b)^{p-1}\equiv 1(\ mod \ p) (Bp+1)p−1≡(−b)p−1≡1( mod p)
所以
B ( p + 1 ) ( p − 1 ) ≡ 1 ( m o d p ) B^{(p+1)(p-1)}\equiv 1(\ mod \ p) B(p+1)(p−1)≡1( mod p)
同理
D ( p + 1 ) ( p − 1 ) ≡ 1 ( m o d p ) D^{(p+1)(p-1)}\equiv 1(\ mod \ p) D(p+1)(p−1)≡1( mod p)
也就是说循环节必是 ( p + 1 ) ( p − 1 ) (p+1)(p-1) (p+1)(p−1)的倍数
解决完
p
p
p为奇素数的情况,再考虑
p
=
2
p=2
p=2的情况,我们可以通过打表打法发现对于
p
=
2
p=2
p=2时任意的
a
,
b
a,b
a,b可以循环节可以为
3
3
3的倍数,但是
p
=
2
k
p=2^k
p=2k时,循环节并不遵从
g
(
2
k
)
=
g
(
2
)
∗
2
k
−
1
g(2^k)=g(2)*2^{k-1}
g(2k)=g(2)∗2k−1,其实本质在于求根公式中有一个除以
2
2
2的操作,所以做不了同余式
对于
p
=
2
k
p=2^k
p=2k时,例如
p
=
4
,
a
=
2
,
b
=
1
p=4,a=2,b=1
p=4,a=2,b=1时循环节为
4
4
4,所以循环节并不符合
g
(
2
k
)
=
g
(
2
)
∗
2
k
−
1
=
6
g(2^k)=g(2)*2^{k-1}=6
g(2k)=g(2)∗2k−1=6,而是4,那么我们其实可以取这两者的最小公倍数,那么其实可以发现对于
p
=
2
k
p=2^k
p=2k的情况循环节满足是
g
(
2
k
)
=
3
∗
2
k
g(2^k)=3*2^k
g(2k)=3∗2k.这样就可以找到循环节了。
顺便提一下如何判定一个数是否是模
p
p
p的二次剩余
a
p
−
1
2
≡
1
(
m
o
d
p
)
a^{\frac{p-1}{2}}\equiv1(\ mod \ p)
a2p−1≡1( mod p)时
a
a
a就是模
p
p
p下的二次剩余
以上题解可以用一个例子加深理解
f
[
n
]
=
f
[
n
−
1
]
+
f
[
n
−
2
]
%
1000000009
f[n]=f[n-1]+f[n-2]\% 1000000009
f[n]=f[n−1]+f[n−2]%1000000009
f
n
=
5
5
[
(
1
+
5
2
)
n
−
(
1
−
5
2
)
n
]
f_n=\frac{\sqrt{5}}{5}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n]
fn=55[(21+5)n−(21−5)n]
38300801
6
2
≡
5
(
m
o
d
1000000009
)
383008016^2\equiv 5(\ mod \ 1000000009)
3830080162≡5( mod 1000000009)
383008016
≡
5
(
m
o
d
1000000009
)
383008016\equiv \sqrt{5}(\ mod \ 1000000009)
383008016≡5( mod 1000000009)
5
5
≡
276601605
(
m
o
d
1000000009
)
\frac{\sqrt{5}}{5}\equiv 276601605 (\ mod \ 1000000009)
55≡276601605( mod 1000000009)
1
+
5
2
≡
691504013
(
m
o
d
1000000009
)
\frac{1+\sqrt{5}}{2}\equiv 691504013 (\ mod \ 1000000009)
21+5≡691504013( mod 1000000009)
1
−
5
2
≡
308495997
(
m
o
d
1000000009
)
\frac{1-\sqrt{5}}{2}\equiv 308495997 (\ mod \ 1000000009)
21−5≡308495997( mod 1000000009)
f
n
≡
276601605
(
69150401
3
n
−
30849599
7
n
)
(
m
o
d
1000000009
)
f_n\equiv 276601605(691504013^n-308495997^n)(\ mod \ 1000000009)
fn≡276601605(691504013n−308495997n)( mod 1000000009)
本题由于求得的循环节可能很大所以在对输入的数进行乘10取模的时候可能会爆LL,所以用快速乘代替
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2;
LL mod = 1000000007;
int cnt,ct;
LL pri[100005];
LL num[100005];
struct Matrix
{
LL m[N][N];
} ;
Matrix A;
Matrix I = {1, 0, 0, 1};
char str[1000005];
map<int,int>y;
int flag=0;
Matrix multi(Matrix a,Matrix b)
{
Matrix c;
for(int i=0; i<N; i++)
{
for(int j=0; j<N; j++)
{
c.m[i][j] =0;
for(int k=0; k<N; k++)
{
c.m[i][j] =(c.m[i][j]+a.m[i][k] * b.m[k][j]%mod)%mod;
}
}
}
return c;
}
Matrix power(Matrix A,LL n)
{
Matrix ans = I, p = A;
while(n>0)
{
if(n & 1)
{
ans = multi(ans,p);
n--;
}
n >>= 1;
p = multi(p,p);
}
return ans;
}
LL quick_mod(LL a,LL b,LL MOD)
{
LL ans = 1;
a %= MOD;
while(b>0)
{
if(b & 1)
{
ans = ans * a % MOD;
b--;
}
b >>= 1;
a = a * a % MOD;
}
return ans;
}
LL quick_mul(LL a,LL b,LL MOD)
{
LL ans=0;
while(b)
{
if(b&1)
{
ans=ans+a;
if(ans>MOD)
{
ans%=MOD;
flag=1;
}
}
a=a+a;
if(a>MOD)
{
a%=MOD;
flag=1;
}
b=b/2;
}
return ans;
}
LL Legendre(LL a,LL p,LL MOD)
{
LL t = quick_mod(a,(p-1)>>1,MOD);
if(t == 1) return 1;
return -1;
}
void yinzi(LL x)
{
for(int i=2; i*i<=x; i++)
{
if(x%i==0)
{
while(x%i==0)
{
y[i]++;
x=x/i;
}
}
}
if(x>1)
{
y[x]++;
}
}
LL lcm(LL a,LL b)
{
return a/__gcd(a,b)*b;
}
LL get_len(LL a,LL b,LL mod)
{
yinzi(mod);
for(map<int,int>::iterator i=y.begin(); i!=y.end(); i++)
{
pri[cnt]=i->first;
num[cnt]=i->second;
cnt++;
}
LL ans=1;
LL c=a*a+4*b;
for(int i=0; i<cnt; i++)
{
LL p=1;
if(pri[i]==2)
{
p*=3;
for(int j=0; j<=num[i]; j++)
p*=pri[i];
}
else if(c%pri[i]==0)
{
p*=pri[i]*(pri[i]-1);
for(int j=1; j<num[i]; j++)
p*=pri[i];
}
else if(Legendre(c,pri[i],pri[i])==1)
{
p*=(pri[i]-1);
for(int j=1; j<num[i]; j++)
p*=pri[i];
}
else
{
p*=(pri[i]-1)*(pri[i]+1);
for(int j=1; j<num[i]; j++)
p*=pri[i];
}
ans=lcm(ans,p);
}
return ans;
}
int main()
{
LL a,b,c,d;
cin>>c>>d>>a>>b;
cin>>str>>mod;
LL p=get_len(a,b,mod);
LL n=0;
int len=strlen(str);
for(int i=0; i<len; i++)
{
n=quick_mul(n,10,p)+str[i]-'0';
}
if(flag)
n+=p;
A.m[0][0]=a;
A.m[0][1]=b;
A.m[1][0]=1;
A.m[1][1]=0;
Matrix T=power(A,n-1);
long long ans=(T.m[0][0]*d%mod+T.m[0][1]*c%mod)%mod;
printf("%lld\n",ans);
return 0;
}