「BZOJ 5104」Fib数列-BSGS+二次剩余

求在 $\bmod 10 ^ 9 + 9$ 的意义下,数字 $C$ 在 $\text{Fib}$ 数列的哪个位置,无解输出 $-1$。

链接

BZOJ 5104

题解

由于模数是 $10 ^ 9 + 9$,其使得我们能找到 $5$ 的二次剩余,所以我们可以先考虑通项公式。
$$\text{Fib}_n = \frac {\sqrt{5}} {5}((\frac {1 + \sqrt{5}} {2}) ^ n - (\frac {1 - \sqrt{5}} {2}) ^ n)$$

我们知道 $\sqrt{5} \equiv 383008016 \pmod {10 ^ 9 + 9}$,所以我们是可以很容易地根据通项公式 + 快速幂求出 $\text{Fib}$ 数列的第 $n$ 项。

下面令 $x = (\frac {1 + \sqrt{5}} {2}) ^ n$,则 $\text{Fib}_n = \frac {\sqrt{5}} {5}(x + \frac {1}{x})$。

那么根据题意
$$\frac {\sqrt{5}} {5}(x + \frac {1}{x}) \equiv C \pmod {10 ^ 9 + 9}$$

两边同乘 $x$ 后变为一个一元二次方程,方程有两根(无解输出 $-1$),我们解出方程的根设其为 $T$,问题变为

$$x \equiv T \pmod {10 ^ 9 + 9}$$

而 $x$ 为 $(\frac {1 + \sqrt{5}} {2}) ^ n$ 或 $(\frac {1 - \sqrt{5}} {2}) ^ n$,直接用 BSGS 求出最小的 $n$ 即可。

这里还有一个问题,我们如何求出这个一元二次方程在 $\bmod 10 ^ 9 + 9$ 意义下的根,由于求根公式为 $\frac {-b \pm \sqrt{\Delta}} {2a}$,只要能求出 $\sqrt{\Delta}$ 其他都能方便的计算。

这里用 cipolla 算法求解:

我们要求 $x ^ 2 \equiv n \pmod {P}$ 的解 $x$(设 $n$ 在模 $P$ 意义下二次剩余),cipolla 算法一开始先不断的随机,找到一个数 $a$ 使得 $\left (\frac{a ^ 2 - n} {p} \right) = -1$ (这里是勒让德符号),此时即这个数无法开根,于是我们扩域,给出一个域 $\mathbb{F}_{p ^ 2}$,通过令 $\omega ^ 2 = a ^ 2 - n$,使得其在 $\mathbb{F}_{p ^ 2}$ 下可以开根,那么
$$x \equiv (a + \omega) ^ \frac {P + 1}{2} \pmod P$$

时间复杂度 $O(\sqrt {P} \log P)$。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/**
* Copyright (c) 2017, xehoth
* All rights reserved.
* 「BZOJ 5104」Fib数列 07-12-2017
* BSGS + 二次剩余
* @author xehoth
*/
#include <bits/stdc++.h>
namespace {
const int MOD = 1e9 + 9;
const int SQRT_5 = 383008016;
const int INV_5 = 200000002;
const int FRAC_SQRT_5_5 = 276601605;
const int INV_FRAC_SQRT_5_5 = 383008016;
const int INV_2 = 500000005;
const int X = 691504013;
const int FRAC_SQRT_5_5_2_MUL_4 = 800000008;
typedef unsigned long long ulong;
typedef unsigned int uint;
inline uint nextUint() {
static uint x = 495;
return x ^= x << 13, x ^= x >> 17, x ^= x << 5, x;
}
inline int modPow(register int a, register int b) {
register int ret = 1;
for (; b; b >>= 1, a = (ulong)a * a % MOD)
(b & 1) ? ret = (ulong)ret * a % MOD : 0;
return ret;
}
struct Complex {
int r, i;
static int w;
Complex(int r, int i) : r(r), i(i) {}
inline Complex operator*(const Complex &x) const {
return Complex(((ulong)r * x.r + (ulong)i * x.i % MOD * w) % MOD,
((ulong)r * x.i + (ulong)i * x.r) % MOD);
}
};
inline Complex modPow(register Complex a, register int b) {
register Complex ret(1, 0);
for (; b; b >>= 1, a = a * a)
if (b & 1) ret = ret * a;
return ret;
}
int Complex::w;
inline int legendre(const register int n) {
return modPow(n, MOD - 1 >> 1) == 1 ? 1 : (-(n != 0));
}
inline int cipolla(register int n) {
if (n == 0) return n;
if (legendre(n) == -1) return -1;
register int a;
for (;;) {
a = nextUint() % MOD;
Complex::w = ((ulong)a * a + MOD - n) % MOD;
if (legendre(Complex::w) == -1) break;
}
return modPow(Complex(a, 1), MOD + 1 >> 1).r;
}
const int H = 42899;
struct HashMap {
std::vector<std::pair<int, int> > d[H];
inline std::pair<int, int> *find(int i) {
register int idx = i % H;
for (register int j = 0; j < d[idx].size(); j++)
if (d[idx][j].first == i) return &d[idx][j];
return NULL;
}
inline void put(const int i, const int v) {
register int idx = i % H;
for (register int j = 0; j < d[idx].size(); j++)
if (d[idx][j].first == i) return (void)(d[idx][j].second = v);
d[idx].push_back(std::pair<int, int>(i, v));
}
} map[2];
int S;
inline int bsgs(int pos, int b) {
register int min = INT_MAX, least = modPow(modPow(X, S), MOD - 2);
register std::pair<int, int> *it;
for (register int i = 0, p, tmp = b; i <= MOD / S;
i++, tmp = (ulong)tmp * least % MOD) {
p = i * S, it = map[(pos - (p & 1) + 2) % 2].find(tmp);
if (!it) continue;
min = std::min(min, it->second + p);
}
return min;
}
inline void solve() {
register long long n;
scanf("%lld", &n);
if (n >= MOD) {
puts("-1");
return;
}
n = (ulong)SQRT_5 * n % MOD;
S = sqrt(MOD);
for (register int i = 0, tmp = 1; i < S; i++, tmp = (ulong)tmp * X % MOD)
map[i & 1].put(tmp, i);
register int delta = cipolla(((ulong)n * n + 4ull) % MOD);
register int ans = INT_MAX;
if (delta != -1) {
ans = std::min(ans, bsgs(0, ((ulong)n + delta) * INV_2 % MOD));
ans = std::min(ans, bsgs(0, ((ulong)n + MOD - delta) * INV_2 % MOD));
}
delta = cipolla(((ulong)n * n + MOD - 4ull) % MOD);
if (delta != -1) {
ans = std::min(ans, bsgs(1, ((ulong)n + delta) * INV_2 % MOD));
ans = std::min(ans, bsgs(1, ((ulong)n + MOD - delta) * INV_2 % MOD));
}
if (ans == INT_MAX) {
puts("-1");
return;
}
printf("%d", ans);
}
} // namespace
int main() {
solve();
return 0;
}
分享到