【2020年杭电暑假第一场】6755 Fibonacci Sum

传送门:http://acm.hdu.edu.cn/showproblem.php?pid=6755

前置技能

二项式定理、二次剩余、等比数列、欧拉降幂、阶乘逆元

参考:牌王无敌

题意

计算一个式子

思路

由斐波那契通项公式得 $$F_n=\frac{1}{\sqrt 5}((\frac{1+ \sqrt 5}{2})^n-(\frac{1-\sqrt5}{2})^n)$$ 可以通过二次剩余的方法求出$\sqrt5=383008016或者616991993$,因为二次剩余的解有两个。下面顺便在给出二次剩余的板子,我也是花了一天的时间才弄懂计算过程。 我们让$A=\frac{1+ \sqrt 5}{2}、B=\frac{1- \sqrt 5}{2}、D=\frac{1}{\sqrt5}$ 进而利用逆元计算出 $A=691504013$ $B=308495997$ $D=276601605$

所以式子就变成了: $$F_n=D(A^n-B^n)(mod 1e9 + 9)$$ 如果这样暴力算的话一定会T,因为N和C都达到了1e18。 我们先尝试列出前几项: $当n=1时,F_c^k=D^k(A^c-B^c)^k$ 用二项式定理展开后变为 $$F_c^k=D^k\sum_{i=0}^{k}(-1)^iC_k^i(A^c)^{k-i}(B^c)^{i}$$

$当n=2时F_{2c}^k=D^k(A^{2c}-B^{2c})^k$ 用二项式定理展开后变为 $$F_{2c}^k=D^k\sum_{i=0}^{k}(-1)^iC_k^i(A^{2c})^{k-i}(B^{2c})^{i}$$ $$F_{2c}^k=D^k\sum_{i=0}^{k}(-1)^iC_k^i(A^{c})^{k-i}(B^{c})^{i}C_k^i(A^{c})^{k-i}(B^{c})^{i}$$ 所以当i相同时,这就是一个等比数列求和。

首项:$(-1)^i(A^{k-i}B^{i})^c$ 公比:$(A^{k-i}B^{i})^c$

最后整理一下整个式子得 $$\sum_{i=0}^{n}F_{ic}^{K}=D^k\sum_{i=0}^{K}C_k^i\frac{(-1)^i(A^{k-i}B^{i})^c(((A^{k-i}B^{i})^c)^n-1)}{(A^{k-i}B^{i})^c-1}$$

当公比为1时,$(A^{k-i}B^{i})^c=1$,所以只要计算$C_k^i*(-1)^i*(A^{k-i}B^{i})^c*n$

为了快速处理$C_m^n$,需要预处理1e5以内的阶乘和阶乘逆元。 $$C_m^n=\frac{m!}{n!*(m-n)!}=m!*inv[n!]*inv[(m-n)!]$$

公比的优化,因为随着i的增大,第$i$项和第$i+1$项相差$(\frac{B}{A})^c$,所以我们率先算出$A$的逆元$invA$,就可以$O(1)$递推求出该公比,即: $$q_{i+1} = qi *(\frac{B}{A})^c$$

对于使用快速幂求$c$次方,可以使用欧拉降幂来优化时间,因为题目给的$mod$是个质数,即: $$a^b = a^{b\%(1e9+9-1)}(mod 1e9+9)$$

Code

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
##include <bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef long double ld;
typedef pair<int, int> pdd;

##define INF 0x7f7f7f
##define mem(a,b) memset(a , b , sizeof(a))
##define FOR(i, x, n) for(int i = x;i <= n; i++)

const ll mod = 1e9 + 9;
const int maxn = 1e5;
// const double eps = 1e-6;

const ll B = 308495997;
const ll invA = 691504012;
const ll A = 691504013;
const ll D = 276601605;

ll F[maxn + 10]; // 阶乘
ll invF[maxn + 10]; // 阶乘逆元
ll invn[maxn + 10];

inline ll quick_pow(ll a, ll b) // 快速幂
{
ll ans = 1, base = a;
while (b)
{
if (b & 1)
{
ans = ans * base % mod;
}
base = base * base % mod;
b >>= 1;
}
return ans % mod;
}

inline void Init()
{
F[0] = F[1] = invF[0] = invF[1] = invn[0] = invn[1] = 1;
for (int i = 2; i <= maxn; i++)
{
F[i] = F[i - 1] * i % mod;
invn[i] = (mod - mod / i) * invn[mod % i] % mod;
invF[i] = invF[i - 1] * invn[i] % mod;
}
}

inline ll getC(int m, int n)
{
if (m < 0 n < 0 n > m)
return 0;
ll ans = F[m];
ans = ans * invF[n] % mod;
ans = ans * invF[m - n] % mod;
return ans;
}

inline ll MOD(ll a, ll b)
{
a += b;
if (a >= mod)
a -= mod;
return a;
}

void solve()
{
ll n, c, k;
cin >> n >> c >> k;
ll ans = 0;
ll a1 = quick_pow(quick_pow(A, k), c % (mod - 1)) % mod; // 首项(公比)
ll q = quick_pow(invA * B % mod, c % (mod - 1)) % mod; // 公比递推
ll n1 = n % mod;
ll n2 = n % (mod - 1); // 欧拉降幂使用
ll a1power = quick_pow(a1, n2) % mod;
ll qpower = quick_pow(q, n2) % mod;

for (int i = 0; i <= k; i++)
{
ll sum = getC(k, i) % mod;
if (i & 1)
sum = mod - sum; // 负数
if (a1 == 1) // 公比为1
{
ans = (ans + n1 * sum % mod) % mod;
}
else
{
sum = sum * (a1 * (a1power - 1 + mod) % mod) % mod;
sum = sum * quick_pow(a1 - 1, mod - 2) % mod;
ans = MOD(ans, sum);
}
a1 = (a1 * q) % mod; // 更换公比
a1power = a1power * qpower % mod;
}
cout << ans * quick_pow(D, k) % mod << endl;
}

signed main()
{
ios_base::sync_with_stdio(false);
//cin.tie(nullptr);
//cout.tie(nullptr);
##ifdef FZT_ACM_LOCAL
// freopen("in.txt", "r", stdin);
// freopen("out.txt", "w", stdout);
##else
ios::sync_with_stdio(false);
int T = 1;
cin >> T;
Init();
while (T--)
solve();
##endif
return 0;
}

二次剩余板子

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
##include <iostream>
##include <ctime>

using namespace std;

typedef long long ll;

typedef struct{
ll x, y; // 把求出来的w作为虚部,则为a + bw
}num;

num num_mul(num a, num b, ll w, ll p) // 复数乘法
{
num ans = {0, 0};
ans.x = (a.x * b.x % p + a.y * b.y % p * w % p + p) % p;
ans.y = (a.x * b.y % p + a.y * b.x % p + p) % p;
return ans;
}

ll num_pow(num a, ll b, ll w, ll p) // 复数快速幂
{
num ans = {1, 0};
while(b)
{
if(b & 1)
{
ans = num_mul(ans, a, w, p);
}
a = num_mul(a, a, w, p);
b >>= 1;
}
return ans.x % p;
}

ll quick_pow(ll a, ll b, ll p) // 快速幂求解
{
ll ans = 1, base = a;
while(b)
{
if(b & 1)
{
ans = ans * base % p;
}
base = base * base % p;
b >>= 1;
}
return ans % p;
}

ll legendre(ll a, ll p) // 勒让德符号 = {1, -1, 0}
{
return quick_pow(a, (p - 1) >> 1, p);

}

ll Cipolla(ll n, ll p) // 输入a和p,是否存在x使得x^2 = a (mod p),存在二次剩余返回x,存在二次非剩余返回-1 注意:p是奇质数
{
n %= p;
if(n == 0)
return 0;
if(p == 2)
return 1;
if(legendre(n, p) + 1 == p) // 二次非剩余
return -1;

ll a, w;

while(true) // 找出a,求出w,随机成功的概率是50%,所以数学期望是2
{
a = rand() % p;
w = ((a * a - n) % p + p) % p;
if(legendre(w, p) + 1 == p) // 找到w,非二次剩余条件
break;
}
num x = {a, 1};
return num_pow(x, (p + 1) >> 1, w, p) % p; // 计算x,一个解是x,另一个解是p-x,这里的w其实要开方,但是由拉格朗日定理可知虚部为0,所以最终答案就是对x的实部用快速幂求解
}

int main()
{
ll n, p;
cin >> n >> p;
srand((unsigned)time(NULL));
cout << Cipolla(n, p) << endl;
return 0;
}

后记

感谢牌王的亲情指导和博客

本文作者:jujimeizuo
本文地址https://blog.jujimeizuo.cn/2020/07/23/hdu-6755/
本博客所有文章除特别声明外,均采用 CC BY-SA 3.0 协议。转载请注明出处!