写在前面
需要学会的前置技能:
(快速幂+二次剩余+逆元+扩展欧几里得) 或者 (矩阵快速幂+二次剩余)
笔者是一名十八线蒟蒻
A
C
M
e
r
ACMer
ACMer,如果文章有误,请在评论区下留言,我会尽快处理,非常感谢!
原题题面
For a series
F
F
F:
F
(
0
)
=
0
,
F
(
1
)
=
1
F(0) = 0,F(1) = 1
F(0)=0,F(1)=1
F
(
n
)
=
3
∗
F
(
n
−
1
)
+
2
∗
F
(
n
−
2
)
,
(
n
≥
2
)
F(n) = 3*F(n-1)+2*F(n-2),(n \geq 2)
F(n)=3∗F(n−1)+2∗F(n−2),(n≥2)
We have some queries. For each query
N
N
N, the answer
A
A
A is the value
F
(
N
)
F(N)
F(N) modulo
998244353
998244353
998244353.
Moreover, the input data is given in the form of encryption, only the number of queries
Q
Q
Q and the first query
N
1
N_1
N1 are given. For the others, the query
N
i
(
2
≤
i
≤
Q
)
N_i(2\leq i\leq Q)
Ni(2≤i≤Q) is defined as the xor of the previous
N
i
−
1
N_{i-1}
Ni−1 and the square of the previous answer
A
i
−
1
A_{i-1}
Ai−1.For example, if the first query
N
1
N_1
N1 is
2
2
2, the answer
A
1
A_1
A1 is
3
3
3, then the second query
N
2
N_2
N2 is
2
x
o
r
(
3
∗
3
)
=
11
2\ xor\ (3*3)=11
2 xor (3∗3)=11.
Finally, you don’t need to output all the answers for every query, you just need to output the xor of each query’s answer
A
1
x
o
r
A
2
.
.
.
x
o
r
A
Q
A_1\ xor\ A_2 ... xor\ A_Q
A1 xor A2...xor AQ
Input:
The input contains two integers,
Q
,
N
,
(
1
≤
Q
≤
1
e
7
,
0
≤
N
≤
1
e
18
)
,
Q
Q, N,(1\leq Q\leq1e7,0\leq N\leq 1e18),Q
Q,N,(1≤Q≤1e7,0≤N≤1e18),Q representing the number of queries and
N
N
N representing the first query.
Output:
An integer representing the final answer.
输入
17 473844410
输出
903193081
题面分析
已知
F
(
n
)
F(n)
F(n),
N
1
N_1
N1求
A
1
x
o
r
A
2
.
.
.
x
o
r
A
Q
A_1\ xor\ A_2 ... xor\ A_Q
A1 xor A2...xor AQ,其中
A
i
=
{
N
i
−
1
x
o
r
(
A
i
−
1
2
)
(
i
≥
2
)
F
(
N
)
(
i
=
1
)
A_i=\begin{cases} N_{i-1}\ xor \ (A_{i-1}^2)\ \ \ (i\geq2)\\ F(N)\ \ \ (i=1)\\ \end{cases}
Ai={Ni−1 xor (Ai−12) (i≥2)F(N) (i=1)
因为异或没有快速或者简便算法,所以我们把重心放在求
F
(
n
)
F(n)
F(n)上。
本文将介绍两种方法供读者选择:
1.二阶线性数列递推+快速幂优化
2.矩阵快速幂+广义斐波那契循环节
二阶线性数列递推
首先我们直接给不喜欢看证明过程的读者写结论。
对于已知
F
(
1
)
,
F
(
2
)
F(1),F(2)
F(1),F(2)的二阶线性数列
F
(
n
)
=
a
∗
F
(
n
−
1
)
+
b
∗
F
(
n
−
2
)
,
(
n
≥
3
)
F(n) = a*F(n-1)+b*F(n-2),(n \geq 3)
F(n)=a∗F(n−1)+b∗F(n−2),(n≥3)
我们先计算方程
x
2
−
a
x
−
b
=
0
x^2-ax-b=0
x2−ax−b=0的解,
如果该方程有二实根
p
,
q
p,q
p,q,则
F
(
n
)
=
A
p
n
+
B
q
n
,
A
=
F
(
2
)
−
q
F
(
1
)
p
(
p
−
q
)
,
B
=
p
F
(
1
)
−
F
(
2
)
q
(
p
−
q
)
F(n)=Ap^n+Bq^n,A=\frac{F(2)-qF(1)}{p(p-q)},B=\frac{pF(1)-F(2)}{q(p-q)}
F(n)=Apn+Bqn,A=p(p−q)F(2)−qF(1),B=q(p−q)pF(1)−F(2)------------①
如果该方程只有一实根
p
p
p,则
F
(
n
)
=
(
A
+
b
n
)
p
n
,
A
=
p
F
(
1
)
−
F
(
2
)
p
2
,
B
=
F
(
2
)
−
p
F
(
1
)
p
2
F(n)=(A+bn)p^n,A=\frac{pF(1)-F(2)}{p^2},B=\frac{F(2)-pF(1)}{p^2}
F(n)=(A+bn)pn,A=p2pF(1)−F(2),B=p2F(2)−pF(1)------------②
所以我们设
a
=
3
,
b
=
2
a=3,b=2
a=3,b=2,代入得到
F
(
n
)
=
1
17
(
(
3
+
17
2
)
n
−
(
3
−
17
2
)
n
)
F(n)=\frac{1}{\sqrt{17}}((\frac{3+\sqrt{17}}{2})^n-(\frac{3-\sqrt{17}}{2})^n)
F(n)=171((23+17)n−(23−17)n)
下面给出这种方法的证明,可以不需要掌握,建议有一定数学基础的读者观看。
证明
出于个人偏好,笔者把
F
(
n
)
F(n)
F(n)记为
a
n
=
A
a
n
−
1
+
B
a
n
−
2
a_{n}=Aa_{n-1}+Ba_{n-2}
an=Aan−1+Ban−2。
为了化简该方程,我们努力构造一种
b
n
+
1
=
a
n
+
1
−
x
a
n
(
x
为
待
定
系
数
)
b_{n+1}=a_{n+1}-xa_{n}(x为待定系数)
bn+1=an+1−xan(x为待定系数),使得原式可以化简为
b
n
+
1
=
k
b
n
(
k
为
待
定
系
数
)
b_{n+1}=kb_{n}(k为待定系数)
bn+1=kbn(k为待定系数),从而将二阶递推转化为一阶线性递推。
那我们展开
b
n
b_n
bn,得到
a
n
+
1
=
(
k
+
x
)
a
n
−
k
x
a
n
−
1
a_{n+1}=(k+x)a_{n}-kxa_{n-1}
an+1=(k+x)an−kxan−1
那就得到
{
k
+
x
=
A
−
k
x
=
B
\begin{cases} k+x=A\\ -kx=B\\ \end{cases}
{k+x=A−kx=B
那我们很容易就会想到韦达定理,所以我们得到一个方程:
x
2
−
A
x
−
B
=
0
x^2-Ax-B=0
x2−Ax−B=0。
这个叫做这个数列的 特征方程 ,我们可以用求根公式解出答案。
1.当方程只有一实根时(记为
p
p
p):
首先我们可以简单地求出
b
1
b_1
b1,那显然
b
n
=
b
1
p
n
b_n=b_1p^n
bn=b1pn,
所以我们得到
a
n
+
1
−
p
a
n
=
b
1
p
n
a_{n+1}-pa_{n}=b_1p^n
an+1−pan=b1pn
再利用和刚才的方法一样,求它的特征方程解出来即可,②式得证。
2.当方程有两实根时(记为
p
,
q
p,q
p,q)
我们有
a
n
+
1
−
q
a
n
=
b
1
p
n
a_{n+1}-qa_{n}=b_1p^n
an+1−qan=b1pn,同上的方法求解即可,①式得证。
二次剩余
你可能会问:
1
17
(
(
3
+
17
2
)
n
−
(
3
−
17
2
)
n
)
\frac{1}{\sqrt{17}}((\frac{3+\sqrt{17}}{2})^n-(\frac{3-\sqrt{17}}{2})^n)
171((23+17)n−(23−17)n)这个式子有无理数和分数存在,会丢失精度,为什么我们要求出它?
所以,我们需要引入二次剩余。
通俗来说,就是已知
n
,
p
n,p
n,p,求
x
2
≡
n
(
m
o
d
p
)
x^2\equiv n(mod\ p)
x2≡n(mod p)的解,具体内容本文不赘述。
我们可以令
n
=
17
,
p
=
998244353
n=17,p=998244353
n=17,p=998244353,代入方程得到其中一个解
x
=
524399943
x=524399943
x=524399943。
这个数特殊在哪里呢?
注意
x
2
≡
17
(
m
o
d
998244353
)
x^2\equiv 17(mod\ 998244353)
x2≡17(mod 998244353)
也就是说
x
≡
17
(
m
o
d
998244353
)
x\equiv \sqrt{17}(mod\ 998244353)
x≡17(mod 998244353),我们记为
i
n
v
s
q
r
t
(
17
)
invsqrt(17)
invsqrt(17)
所以
524399943
524399943
524399943就是在模998244353的情况下,等价于
17
\sqrt{17}
17的数字!!!
这是我们今天解出这题的第一个关键点!
有了这个数字,我们可以得到
3
+
17
2
=
(
3
+
524399943
)
∗
i
n
v
2
m
o
d
998244353
=
262199973
\frac{3+\sqrt{17}}{2}=(3+524399943)*inv2\ mod\ 998244353=262199973
23+17=(3+524399943)∗inv2 mod 998244353=262199973
其中,
i
n
v
2
inv2
inv2是2的逆元。利用费马小定理求解。
同理求出
17
\sqrt{17}
17的逆元---->
52439994
3
998244353
−
2
m
o
d
998244353
524399943^{998244353-2}\ mod\ 998244353
524399943998244353−2 mod 998244353,记为
i
n
v
s
q
r
t
(
17
)
invsqrt(17)
invsqrt(17)
所以也可以得到
3
−
17
2
=
736044383
\frac{3-\sqrt{17}}{2}=736044383
23−17=736044383
因此
F
(
n
)
=
(
26219997
3
n
−
73604438
3
n
)
∗
i
n
v
s
q
r
t
(
17
)
m
o
d
998244353
F(n)=(262199973^n-736044383^n)*invsqrt(17)\ mod\ 998244353
F(n)=(262199973n−736044383n)∗invsqrt(17) mod 998244353。
n次幂的部分我们可以用快速幂求解,复杂度为
O
(
l
o
g
2
n
)
O(log_2n)
O(log2n)。
优化快速幂
按理来说,这题到这里也就结束了,但因为本题
Q
Q
Q最大可以是
1
e
7
1e7
1e7,
O
(
l
o
g
2
n
)
O(log_2n)
O(log2n)的复杂度会超时!
那应该怎么改呢?
我们来考虑考虑优化快速幂。
我们原本使用的快速幂其实是基于2进制的快速幂。
举个例子,求解
2
23
2^{23}
223时,我们有
23
=
16
+
4
+
2
+
1
23=16+4+2+1
23=16+4+2+1(拆成2进制上的权值)
所以
2
23
=
2
16
2
4
2
2
2
1
2^{23}=2^{16}2^42^22^1
223=216242221。
这样计算不仅快捷,在模
p
p
p时还会避免溢出。
尽管它已经是如此优秀,但依旧达不到我们的要求。
因此我们提出对快速幂的优化-------->k进制快速幂。
k进制快速幂的原理与普通快速幂相同,就是先把一个数拆成k进制上的权值,再相乘求和。
对于一些题而言,
O
(
l
o
g
10
n
)
O(log_{10}n)
O(log10n)已经可以应付大部分,甚至大数的情况了,但我们此处对其要求更高,于是可以考虑
k
≥
1
e
4
k\geq1e4
k≥1e4。
笔者在写题时考虑到
n
n
n可能为
1
e
18
1e18
1e18,就设
k
=
2
16
=
65536
k=2^{16}=65536
k=216=65536,这样权值位只有四种:
1
,
65536
,
6553
6
2
,
6553
6
3
1,65536,65536^2,65536^3
1,65536,655362,655363。
那以
F
(
n
)
F(n)
F(n)中的
26219997
3
n
262199973^n
262199973n的优化为例,我们有
26219997
3
n
=
x
1
n
/
6553
6
3
+
x
2
(
n
m
o
d
6553
6
3
)
/
6553
6
2
+
x
3
(
n
m
o
d
6553
6
2
)
/
65536
+
x
4
n
m
o
d
65536
262199973^n=x_1^{n/65536^3}+x_2^{(n\ mod\ 65536^3)/65536^2}+x_3^{(n\ mod\ 65536^2)/65536}+x_4^{n\ mod\ 65536}
262199973n=x1n/655363+x2(n mod 655363)/655362+x3(n mod 655362)/65536+x4n mod 65536。
其中
x
1
=
26219997
3
281474976710656
m
o
d
998244353
x_1=262199973^{281474976710656}\ mod\ 998244353
x1=262199973281474976710656 mod 998244353
x
2
=
26219997
3
4294967296
m
o
d
998244353
x_2=262199973^{4294967296}\ mod\ 998244353
x2=2621999734294967296 mod 998244353
x
3
=
26219997
3
65536
m
o
d
998244353
x_3=262199973^{65536}\ mod\ 998244353
x3=26219997365536 mod 998244353
x
4
=
262199973
x_4=262199973
x4=262199973
这四个数可以用普通快速幂提前处理出来。
同理也是可以优化
73604438
3
n
736044383^n
736044383n。
这样复杂度就从原来的
O
(
l
o
g
2
n
)
O(log_2n)
O(log2n)降到
O
(
l
o
g
65536
n
)
O(log_{65536}n)
O(log65536n)左右了,即使
n
=
1
e
18
n=1e18
n=1e18,复杂度大约也在
O
(
96
)
O(96)
O(96)不到。当然
k
k
k可以取更大的或者稍小一些的数字,有兴趣的读者可以试一试。
需要注意的是,在求
A
i
A_i
Ai的时候可能会产生循环节,因此把每一次的结果都要记忆化。 好像群里大佬说只有最小循环节期望只有几万…
事实上,这个方法可以再优化:利用欧拉定理降幂,可以把n降到
(
p
−
1
)
(p-1)
(p−1)以下后,再用上述方法取
k
=
65536
k=65536
k=65536求解,就只有两个权值:
1
1
1和
65536
65536
65536了。复杂度可以降到原来的
1
2
\frac{1}{2}
21左右。这部分留给读者自己去探索。
AC代码(305ms)
//代码很丑,大家见谅TAT
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=998244353;
map<long long,long long>P;
long long qsm(long long a,long long b,long long c)//普通快速幂
{
long long ans=1,base=a;
while(b!=0)
{
if(b&1)
ans=ans*base%c;
base=base*base%c;
b>>=1;
}
return ans%c;
}
long long inv17=524399943;//根号17
long long ans1=262199973,ans11=133071688,ans12=296255346,ans13=253151293;
long long ans2=736044383,ans21=19600266,ans22=683624219,ans23=530734902;
//都是预处理用的数字
long long f1(long long ans,long long n)//优化快速幂
{
long long x0=n/281474976710656LL;
long long x1=x1=(n%281474976710656LL)/4294967296LL;
long long x2=(n%4294967296LL)/65536LL;
long long x3=n%65536LL;
long long sum=(((((qsm(ans12,x1,mod)%mod)*qsm(ans11,x2,mod)%mod)%mod)*qsm(ans1,x3,mod)%mod)*qsm(ans13,x0,mod)%mod)%mod;
return sum;
}
long long f2(long long ans,long long n)
{
long long x0=n/281474976710656LL;
long long x1=x1=(n%281474976710656LL)/4294967296LL;
long long x2=(n%4294967296LL)/65536LL;
long long x3=n%65536LL;
long long sum=((((qsm(ans22,x1%mod,mod)%mod)*(qsm(ans21,x2,mod)%mod)%mod)*qsm(ans2,x3,mod)%mod)*qsm(ans23,x0,mod)%mod)%mod;
return sum;
}
long long f(long long n)//求F(n)
{
if(P.find(n)!=P.end())
return P[n];
long long ans=(f1(ans1,n)-f2(ans2,n)+mod)%mod;
ans=ans*qsm(inv17,mod-2,mod)%mod;
P[n]=ans%mod;
return ans%mod;
}
int main()
{
long long q,n;
scanf("%lld%lld",&q,&n);
long long ans=f(n),sum=f(n);
for(int i=1;i<=q-1;i++)
{
n=n^(ans*ans);
ans=f(n);
sum=sum^ans;
}
printf("%lld\n",sum);
}
真是道一言难尽的题…
但在赛后听了群里大佬们的做法之后,我仔细想了想,发现 这个做法缺乏一般性。
了解二次剩余的读者可以发现,对于
x
2
≡
n
(
m
o
d
p
)
x^2\equiv n(mod\ p)
x2≡n(mod p),
x
x
x不一定有解!,这也就意味着今天如果
F
(
N
)
F(N)
F(N)下的系数如果被设计到没有二次剩余(比如改成
F
(
n
)
=
3
∗
F
(
n
−
1
)
+
F
(
n
−
2
)
F(n) = 3*F(n-1)+F(n-2)
F(n)=3∗F(n−1)+F(n−2),根据上述方法做的话,你会发现根号
x
2
≡
13
(
m
o
d
998244353
)
x^2\equiv 13(mod\ 998244353)
x2≡13(mod 998244353)无解!),该方法便无法使用!
那咋整呢?那咋整呢?那咋整呢?
这里引出第二种方法(通用):广义斐波那契循环节
广义斐波那契循环节
对于数列
F
(
n
)
=
a
∗
F
(
n
−
1
)
+
b
∗
F
(
n
−
2
)
F(n) = a*F(n-1)+b*F(n-2)
F(n)=a∗F(n−1)+b∗F(n−2),我们称其为广义斐波那契数列。
在数学上可以证明:
F
(
n
)
m
o
d
p
F(n)\ mod\ p
F(n) mod p是存在循环节的。这里不提供证明,可以说下结论。
记
c
=
a
2
+
4
b
c=a^2+4b
c=a2+4b,若
c
c
c是模
p
p
p的二次剩余时,枚举
p
−
1
p-1
p−1的所有因子,找到最小的因子
a
n
s
ans
ans,使得
[
a
b
1
0
]
a
n
s
(
m
o
d
p
)
≡
[
1
0
0
1
]
\begin{bmatrix} a & b \\ 1 & 0 \end{bmatrix}^{ans}(mod\ p)\equiv \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
[a1b0]ans(mod p)≡[1001]
a
n
s
ans
ans就是所求的循环节。
若
c
c
c不是模
p
p
p的二次剩余,则枚举
p
2
−
1
p^2-1
p2−1的所有因子,方法同上。
找出循环节
T
T
T之后(笔者之前算出循环节为
499122476
499122476
499122476(即
p
−
1
2
)
\frac{p-1}{2})
2p−1)),
F
(
N
)
=
F
(
N
m
o
d
T
)
F(N)=F(N\ mod\ T)
F(N)=F(N mod T)。因为这个循环节还是很大,所以要再用矩阵快速幂去得到解循环节以内的内容。
当然也可以用类似优化快速幂的方法去优化矩阵快速幂。
这里给出我队友的代码。
AC代码(683ms)
#pragma warning(disable:4996)
#include<bits/stdc++.h>
using namespace std;
#define per(a,b,c) for(int a=b;a<c;++a)
#define rep(a,b,c) for(int a=b;a>c;--a)
#define ll long long
const int maxn(2e2 + 10);
const int mod(998244353);
const int G(3);
const int xunhuanjie(499122176);
#define MAXN maxn
ll a[3][3],b[3][3],c[3][3],dabiao[100000][3][3];
inline void qmatrix(ll m) {
ll d[3][3];
per(i, 1, 3) per(j, 1, 3) b[i][j] = a[i][j];
c[1][1] = 1;
c[1][2] = 0;
c[2][1] = 0;
c[2][2] = 1;
while(m){
if (m & 1) {
d[1][1] = c[1][1] * b[1][1] + c[1][2] * b[2][1];
d[1][2] = c[1][1] * b[1][2] + c[1][2] * b[2][2];
d[2][1] = c[2][1] * b[1][1] + c[2][2] * b[2][1];
d[2][2] = c[2][1] * b[1][2] + c[2][2] * b[2][2];
per(i, 1, 3)per(j, 1, 3)c[i][j] = d[i][j] % mod;
per(i,1, 3)per(j,1, 3)c[i][j] = d[i][j] % mod;
}
d[1][1] = b[1][1] * b[1][1] + b[1][2] * b[2][1];
d[1][2] = b[1][1] * b[1][2] + b[1][2] * b[2][2];
d[2][1] = b[2][1] * b[1][1] + b[2][2] * b[2][1];
d[2][2] = b[2][1] * b[1][2] + b[2][2] * b[2][2];
per(i,1, 3)per(j,1, 3)b[i][j] = d[i][j] % mod;
m >>= 1;
}
}
inline void tenmatrix(ll m) {
ll d[3][3];
per(i,1, 3) per(j,1, 3) b[i][j] = a[i][j];
c[1][1] = 1;
c[1][2] = 0;
c[2][1] = 0;
c[2][2] = 1;
int i = 0;
while (m) {
ll tmp = m % 100;
m /= 100;
if (tmp == 0) {
i += 100; continue;
}
d[1][1] = c[1][1] * dabiao[i+tmp][1][1] + c[1][2] * dabiao[i+tmp][2][1];
d[1][2] = c[1][1] * dabiao[i+tmp][1][2] + c[1][2] * dabiao[i+tmp][2][2];
d[2][1] = c[2][1] * dabiao[i+tmp][1][1] + c[2][2] * dabiao[i+tmp][2][1];
d[2][2] = c[2][1] * dabiao[i+tmp][1][2] + c[2][2] * dabiao[i+tmp][2][2];
per(k, 1, 3)per(j,1, 3)c[k][j] = d[k][j] % mod;
i += 100;
}
}
ll Q, N,ten[20];
int main() {
#ifndef ONLINE_JUDGE
freopen("./Acm.in", "r", stdin);
freopen("./Acm.out", "w", stdout);
#endif // !ONLINE_JUDGE
ios_base::sync_with_stdio(false);
cin.tie(0);
cout.tie(0);
a[1][1] = 3;
a[1][2] = 2;
a[2][1] = 1;
a[2][2] = 0;
ten[0] = 1;
per(i, 1, 20)ten[i] = ten[i - 1] * 100 % xunhuanjie;
per(i, 1, 1001) {
if (i % 100 == 0) {
continue;
}
else {
qmatrix((ten[i/100] * (i % 100))%xunhuanjie);
per(k, 1, 3)per(j,1, 3)dabiao[i][k][j] = c[k][j];
}
}
cin >> Q >> N;
ll ans, finans;
tenmatrix((N-1)%xunhuanjie);
ans = c[1][1] % mod;
finans = ans;
//cout << ans << endl;
per(i, 2, Q+1) {
N = N ^ (ans * ans);
tenmatrix((N-1)%xunhuanjie);
//N = M; qmu(N - 1);
ans = c[1][1]%mod;
finans ^= ans;
//cout << ans << endl;
}
cout << finans << "\n";
}
后记
感谢一波今天的银川网络赛,让我享受赛前起不了床,赛时交不了题,赛后莫名AC的待遇。
祝主办方 提乾日日涉经,令堂天天下葬。
今天和A题斗智斗勇两个小时,赛后一看防AK的TAT…
DrGilbert 2019.9.8