解题思路
太难了搞了我整整一天,数论咋这么难
参考题解
1.累加法化简 g ( n ) g(n) g(n):
在看看g(n)的定义 g(n) = f(1)^2 + f(2)^2 + … + f(n)^2;
因为
所以:
多写几个就可以看出累加法:
所以:
2.矩阵快速幂求 f ( n ) f(n) f(n):
由
[
f
(
1
)
,
f
(
0
)
]
×
A
n
=
[
f
(
n
+
1
)
,
f
(
n
)
]
[f(1),f(0)] \times A^{n} = [f(n+1),f(n)]
[f(1),f(0)]×An=[f(n+1),f(n)]
其中
f
(
1
)
=
1
,
f
(
0
)
=
0
,
f(1)=1, f(0)=0,
f(1)=1,f(0)=0,
A
=
[
2
1
1
0
]
A= \left[ \begin{matrix} 2 & 1\\ 1 & 0 \end{matrix} \right]
A=[2110]
设结果矩阵为
r
e
s
res
res,则
f
(
n
+
1
)
=
r
e
s
[
0
]
[
0
]
,
f
(
n
)
=
r
e
s
[
0
]
[
1
]
f(n+1)=res[0][0],f(n)=res[0][1]
f(n+1)=res[0][0],f(n)=res[0][1]。
3.模运算的处理。
现在的式子化为: x g ( n ) % ( s + 1 ) = x f ( n ) × f ( n + 1 ) 2 % ( s + 1 ) x^{g(n)}\%(s+1)=x^{\frac {f(n) \times f(n+1)}{2}}\%(s+1) xg(n)%(s+1)=x2f(n)×f(n+1)%(s+1)
涉及幂运算的取模,运用欧拉降幂公式的第三种情况
a
b
%
p
=
{
a
b
%
ϕ
(
p
)
%
p
g
c
d
(
a
,
p
)
=
1
a
b
%
p
g
c
d
(
a
,
p
)
≠
1
,
b
<
ϕ
(
p
)
a
b
%
ϕ
(
p
)
+
ϕ
(
p
)
%
p
g
c
d
(
a
,
p
)
≠
1
,
b
≥
ϕ
(
p
)
a^b\%p=\left\{ \begin{aligned} &a^{b\%\phi(p)}\%p& \quad gcd(a,p)=1\\ &a^b\%p& \quad gcd(a,p) \neq1,b<\phi(p)\\ &a^{b\%\phi(p)+\phi(p)}\%p& \quad gcd(a,p) \neq1,b\ge\phi(p)\\ \end{aligned} \right .
ab%p=⎩⎪⎪⎨⎪⎪⎧ab%ϕ(p)%pab%pab%ϕ(p)+ϕ(p)%pgcd(a,p)=1gcd(a,p)=1,b<ϕ(p)gcd(a,p)=1,b≥ϕ(p)
化为:
x
f
(
n
)
×
f
(
n
+
1
)
2
%
ϕ
(
s
+
1
)
+
ϕ
(
s
+
1
)
%
(
s
+
1
)
x^{\frac {f(n) \times f(n+1)}{2}\%\phi(s+1)+\phi(s+1)}\%(s+1)
x2f(n)×f(n+1)%ϕ(s+1)+ϕ(s+1)%(s+1)
单独求一个欧拉函数的板子:
LL phi(LL n) {
LL res = n;
for(int i = 2; i * i <= n; i++) {
if (n % i == 0) {
res = res - res / i;
while (n % i == 0) n /= i;
}
}
if (n > 1) res = res / n * (n - 1);
return res;
}
不知道啥是欧拉函数戳这里,很详细~
单独拎出 x x x 的指数部分的前半部分看,要求 f ( n ) × f ( n + 1 ) 2 % ϕ ( s + 1 ) \frac {f(n) \times f(n+1)}{2}\%\phi(s+1) 2f(n)×f(n+1)%ϕ(s+1)。
因为 2 2 2 和 ϕ ( s + 1 ) \phi(s+1) ϕ(s+1) 不一定互质,因此不能用乘法逆元求解。
正确做法是应用这个公式: a b % m o d = a % m o d × b b \frac ab\%mod=\frac {a\%mod\times b}{b} ba%mod=ba%mod×b
- 证明:
a b m o d k = d a b = k x + d a = k b x + b d a m o d k b = b d a m o d k b b = d \frac ab\mod k = d \\ \frac ab=kx+d\\ a=kbx+bd\\ a\mod kb=bd\\ \frac {a \mod kb}{b}=d bamodk=dba=kx+da=kbx+bdamodkb=bdbamodkb=d
因此 f ( n ) × f ( n + 1 ) 2 % ϕ ( s + 1 ) = f ( n ) × f ( n + 1 ) % 2 × ϕ ( s + 1 ) 2 \frac {f(n) \times f(n+1)}{2}\%\phi(s+1)= \frac {f(n) \times f(n+1) \%2 \times\phi(s+1)}{2} 2f(n)×f(n+1)%ϕ(s+1)=2f(n)×f(n+1)%2×ϕ(s+1),可劲儿求吧。
注意事项
求 f ( n ) f(n) f(n) 的时候运用矩阵快速幂取模,注意此时的模不是 s + 1 s+1 s+1,而是 2 × ϕ ( s + 1 ) 2 \times\phi(s+1) 2×ϕ(s+1)。
参考代码
#include<stdio.h>
#include<iostream>
#include<vector>
#include<cstring>
#include<cstdio>
#include<climits>
#include<cmath>
#include<algorithm>
#include<queue>
#include<deque>
#include<map>
#include<set>
#include<stack>
//#define LOCAL //提交时一定注释
#define VI vector<int>
#define io ios::sync_with_stdio(false); cin.tie(0); cout.tie(0)
using namespace std;
typedef long long LL;
typedef double db;
const int inf = 0x3f3f3f3f;
const LL INF = 1e18;
typedef vector<VI> mat;
inline int readint() {int x; scanf("%d", &x); return x;}
LL mod;
mat multi(mat &A, mat &B) {
mat C(A.size(), VI(B[0].size()));
for(int i = 0; i < A.size(); i++) {
for(int j = 0; j < B.size(); j++) {
for(int k = 0; k < B.size(); k++) {
C[i][j] = (C[i][j] + 1LL * A[i][k] * B[k][j] % mod) % mod;
}
}
}
return C;
}
mat fm(mat A, LL n) {
mat res(A.size(), VI(A.size()));
for(int i = 0; i < A.size(); i++)
res[i][i] = 1;
while (n) {
if (n & 1) res = multi(res, A);
A = multi(A, A);
n >>= 1;
}
return res;
}
LL phi(LL n) {
LL res = n;
for(int i = 2; i * i <= n; i++) {
if (n % i == 0) {
res = res - res / i;
while (n % i == 0) n /= i;
}
}
if (n > 1) res = res / n * (n - 1);
return res;
}
LL qm(LL a, LL b, LL c) {
LL res = 1;
while (b) {
if (b & 1) res = res * a % c;
a = a * a % c;
b >>= 1;
}
return res;
}
int main() {
#ifdef LOCAL
freopen("input.txt", "r", stdin);
// freopen("output.txt", "w", stdout);
#endif
int t = readint();
mat A(2, VI(2)), B(1, VI(2));
A[0][0] = 2;
A[0][1] = A[1][0] = 1;
B[0][0] = 1, B[0][1] = 0;
LL n, y, x, s;
mat res, A_n;
while (t--) {
scanf("%lld%lld%lld%lld", &n, &y, &x, &s);
mod = 2 * phi(s + 1);
A_n = fm(A, n * y); //求A^n
res = multi(B, A_n);
LL mi = (res[0][0] % mod * res[0][1] % mod) % mod;
mi = mi / 2 + mod / 2;
printf("%lld\n", qm(x, mi, s + 1));
}
return 0;
}