传送门
思路
唉,我太弱了,T1 都做不来,唉,我太弱啦!
显然这个题可以用矩阵快速幂在
O(logn)
O
(
log
n
)
的时间复杂度内解决。这里的
n
n
为 106,
logn
log
n
为
3.3×106
3.3
×
10
6
,但是不要忘了我们只能用高精度来存储这个数,这样我们的时间复杂度就变成了
O(n18logn)
O
(
n
18
log
n
)
,显然不能这么做。
观察发现,最后的结果仅要求给出在模
109+7
10
9
+
7
的意义下的值,又观察递推式,发现
f(x)
f
(
x
)
仅由
f(x−1)
f
(
x
−
1
)
决定,这意味着
f(x)
f
(
x
)
在模意义下存在循环,且循环节长度至多为
109+7
10
9
+
7
。另外,
f(1)
f
(
1
)
一定包含在循环内。
这意味着,如果我们知道了循环节长度,那么这道题就做出来了。当
a=1
a
=
1
时(
c
c
和 d 同理,这里我们用
a
a
和 b 来讨论),根据学习 Pólya 定理时的经验,我们知道循环节长度为
109+7gcd(a,109+7)
10
9
+
7
gcd
(
a
,
10
9
+
7
)
,由于模数是个质数,因此循环节长度就是
109+7
10
9
+
7
。
当
a≠1
a
≠
1
时,又该怎么做呢?我太弱了,什么都不会,于是随便选了几个数打了一个表(用 Folyd 判圈法),发现循环节长度要么为
109+6
10
9
+
6
,要么为它的一半。由于计算一次需要花费数秒时间,因此这不大可能是全部情况。另外,由于
f(1)
f
(
1
)
一定在循环内,因此这里的循环节指的是除去一开始的
1
1
,当前数又一次变成 1 所需的步数。
然后因为我太弱了,所以其它的都不会做了。唉,我太弱啦!
正解
按照上面的思路,显然我们可以猜想一个做法:对于
a=1
a
=
1
的情况,我们将
n
n
对 109+7 取模;对于
a≠1
a
≠
1
的情况,我们将
n
n
对 109+6 取模(因为观察到循环节是
109+6
10
9
+
6
的因数)。这样做过后,我们就能通过此题了。唉,发现了因数的关系,却没有想到以
109+6
10
9
+
6
为模数,我太弱啦!
粗略看了一下,网上居然没有针对指数取模的矩阵快速幂的题解,要不就是神仙们说矩阵也满足费马小定理,他们也不用证明,太强啦!唉,我太弱啦!
关键是,我又看错题了,所以上面的白写了,唉,我太弱啦!
先不管看错题的问题。现在的问题是,如何证明递推式
f(n)=af(n−1)+b
f
(
n
)
=
a
f
(
n
−
1
)
+
b
当系数
a
a
不为 1 时循环节的长度是
109+6
10
9
+
6
的因数。(当
a=1
a
=
1
时根据经验,其实我们已经解决了,所以我们先只考虑
a≠1
a
≠
1
的情况)
我们来看递推公式
f(n)=af(n−1)+b
f
(
n
)
=
a
f
(
n
−
1
)
+
b
。我们还是来推导一下通项公式:
f(n)+1a−1b=af(n−1)+aa−1b
f
(
n
)
+
1
a
−
1
b
=
a
f
(
n
−
1
)
+
a
a
−
1
b
f(n)+1a−1bf(n−1)+1a−1b=a
f
(
n
)
+
1
a
−
1
b
f
(
n
−
1
)
+
1
a
−
1
b
=
a
f(n)+1a−1bf(1)+1a−1b=an−1
f
(
n
)
+
1
a
−
1
b
f
(
1
)
+
1
a
−
1
b
=
a
n
−
1
f(n)=an−1+(an−1−1)ba−1
f
(
n
)
=
a
n
−
1
+
(
a
n
−
1
−
1
)
b
a
−
1
观察与
n
n
有关的东西,只存在于 a 的指数。由费马小定理,
a109+6≡1(mod109+7)
a
10
9
+
6
≡
1
(
mod
10
9
+
7
)
,因此每
109+6
10
9
+
6
次操作一定会回到起点。故我们可以将矩阵乘法的指数模上
109+6
10
9
+
6
。
同样我们用数列的知识来看看为什么当
a=1
a
=
1
时指数要模
109+7
10
9
+
7
。这时原递推公式对应的数列就从一个类似于等比数列的东西变成了一个等差数列,其通项公式为:
f(n)=(n−1)×b+f(1)
f
(
n
)
=
(
n
−
1
)
×
b
+
f
(
1
)
与
n
n
有关的东西只存在于 b 的系数。由模意义下的乘法,
ab≡(amod(109+7))×b(mod109+7)
a
b
≡
(
a
mod
(
10
9
+
7
)
)
×
b
(
mod
10
9
+
7
)
。故在
a=1
a
=
1
时循环节变成了
109+7
10
9
+
7
。
我之前看成了下面的位置是由正上方的位置推导出来的,结果是第一个是由上一行的最后一个推导出来的。不过这并不影响循环节的性质:我们现在可以构造一个可计算的矩阵使得从某一行的第一个数转移到最后一个数,把它再乘上
cd
c
d
的转移矩阵,就得到了从某一行的第一个数推导出下一行的第一个数的矩阵。显然的是(你手算一下,发现第二行确实如此),这个矩阵的形状如下:
[x0y1]
[
x
y
0
1
]
因此它仍然代表着一个
f(i)=xf(i−1)+y
f
(
i
)
=
x
f
(
i
−
1
)
+
y
,循环节仍然存在。不过需要注意的是,第二步就不能通过
c
c
特判 1,而要通过上面的
x
x
特判 1。
参考代码
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <cassert>
#include <cctype>
#include <climits>
#include <ctime>
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <list>
#include <functional>
typedef long long LL;
typedef unsigned long long ULL;
using std::cin;
using std::cout;
using std::endl;
typedef LL INT_PUT;
INT_PUT readIn()
{
INT_PUT a = 0; bool positive = true;
char ch = getchar();
while (!(ch == '-' || std::isdigit(ch))) ch = getchar();
if (ch == '-') { positive = false; ch = getchar(); }
while (std::isdigit(ch)) { a = a * 10 - (ch - '0'); ch = getchar(); }
return positive ? -a : a;
}
void printOut(INT_PUT x)
{
char buffer[20]; int length = 0;
if (x < 0) putchar('-'); else x = -x;
do buffer[length++] = -(x % 10) + '0'; while (x /= 10);
do putchar(buffer[--length]); while (length);
}
const int mod = int(1e9) + 7;
const int maxn = int(1e6) + 5;
char strn[maxn];
char strm[maxn];
LL n, m;
int a, b, c, d;
struct Matrix
{
int c[2][2];
Matrix() : c() {}
Matrix(bool unit) : c()
{
c[0][0] = c[1][1] = 1;
}
int* operator[](int x) { return c[x]; }
const int* operator[](int x) const { return c[x]; }
Matrix operator*(const Matrix& b) const
{
Matrix ret;
for (int i = 0; i < 2; i++)
for (int k = 0; k < 2; k++) if (c[i][k])
for (int j = 0; j < 2; j++)
ret[i][j] = (ret[i][j] + (LL)c[i][k] * b[k][j]) % mod;
return ret;
}
Matrix operator^(int y) const
{
Matrix ret(true);
Matrix x(*this);
while (y)
{
if (y & 1) ret = ret * x;
x = x * x;
y >>= 1;
}
return ret;
}
} matrix, nextline;
void run()
{
scanf("%s%s", strn, strm);
a = readIn();
b = readIn();
c = readIn();
d = readIn();
int tempMod;
tempMod = mod - (a != 1);
for (char* ch = strm; *ch; ch++)
m = (((m << 3) + (m << 1)) + (*ch - '0')) % tempMod;
m = (m - 1 + tempMod) % tempMod;
matrix[0][0] = a;
matrix[0][1] = b;
matrix[1][0] = 0;
matrix[1][1] = 1;
matrix = matrix ^ m;
nextline[0][0] = c;
nextline[0][1] = d;
nextline[1][0] = 0;
nextline[1][1] = 1;
nextline = matrix * nextline;
tempMod = mod - (nextline[0][0] != 1);
for (char* ch = strn; *ch; ch++)
n = (((n << 3) + (n << 1)) + (*ch - '0')) % tempMod;
n = (n - 1 + tempMod) % tempMod;
matrix = (nextline ^ n) * matrix;
printOut((matrix[0][0] + matrix[0][1]) % mod);
}
int main()
{
run();
return 0;
}
Remarks
暴力不要在二进制下进行,要在十进制下进行,这样的话就不用高精度了,时间复杂度为
O(80n)
O
(
80
n
)
,被卡常得到
0
0
<script type="math/tex" id="MathJax-Element-65">0</script> 分。
还可以用通项公式(高一数列内容),然后用逆元来做。然而我太弱了,什么都不会。神仙们都是直接十进制快速幂过了,或者用一个没有证明的矩阵费马小定理(说白了就是错的),太强啦!
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)