原题题面
The Fibonacci numbers are defined as below:
F
0
=
0
,
F
1
=
1
F_0=0,F_1=1
F0=0,F1=1
F
n
=
F
n
−
2
+
F
n
−
1
(
n
>
1
)
F_n=F_{n-2}+F_{n-1}(n>1)
Fn=Fn−2+Fn−1(n>1)
Given three integers
N
N
N,
C
C
C and
K
(
1
≤
N
,
C
≤
1
e
18
,
1
≤
K
≤
1
e
5
)
K(1\leq N,C \leq 1e18,1 \leq K \leq 1e5)
K(1≤N,C≤1e18,1≤K≤1e5), calculate the following summation:
(
F
0
)
K
+
(
F
C
)
K
+
(
F
2
C
)
K
+
(
F
3
C
)
K
.
.
.
+
(
F
N
C
)
K
(F_0)^K+(F_{C})^K+(F_{2C})^K+(F_{3C})^K...+(F_{NC})^K
(F0)K+(FC)K+(F2C)K+(F3C)K...+(FNC)K
Since the answer can be huge, output it modulo
1000000009
1000000009
1000000009 (
1
0
9
+
9
10^9+9
109+9).
输入样例
2
5 1 1
2 1 2
输出样例
12
2
题面分析
众所周知, 斐波那契数列的通项公式是
F
n
=
1
5
(
(
1
+
5
2
)
n
−
(
1
−
5
2
)
n
)
F_n=\frac{1}{\sqrt 5}((\frac{1+\sqrt 5}{2})^n-(\frac{1-\sqrt 5}{2})^n)
Fn=51((21+5)n−(21−5)n)
因为通式带根号,所以要想办法用二次剩余处理,
具体过程见我去年牛客的这篇博客->2019南昌H题
我们可以轻松 求出
5
=
383008016
\sqrt 5= 383008016
5=383008016或者
616991993
616991993
616991993。
以
616991993
616991993
616991993为例,我们可以得到
1
+
5
2
=
691504013
\frac{1+\sqrt 5}{2}=691504013
21+5=691504013
1
−
5
2
=
308495997
\frac{1-\sqrt 5}{2}=308495997
21−5=308495997
1
5
=
276601605
\frac{1}{\sqrt 5}=276601605
51=276601605
我们记
A
=
1
+
5
2
,
B
=
1
−
5
2
,
D
=
1
5
A=\frac{1+\sqrt 5}{2},B=\frac{1-\sqrt 5}{2},D=\frac{1}{\sqrt 5}
A=21+5,B=21−5,D=51,则
F
n
=
D
(
A
n
−
B
n
)
(
m
o
d
1
0
9
+
9
)
F_n=D(A^n-B^n)(mod\ 10^9+9)
Fn=D(An−Bn)(mod 109+9)
虽然可以在
O
(
l
o
g
n
)
O(logn)
O(logn)的时间内求出单个
F
n
F_n
Fn,但
N
,
C
N,C
N,C都是
1
e
18
1e18
1e18,因此考虑别的办法。
注意到
n
=
1
n=1
n=1时,
F
c
k
=
[
D
(
A
c
−
B
c
)
]
k
F_c^k=[D(A^c-B^c)]^k
Fck=[D(Ac−Bc)]k,将其展开,由二项式定理可以得到
F
c
k
=
D
k
∑
i
=
0
k
C
k
i
(
A
c
)
k
−
i
(
−
B
c
)
i
=
D
k
∑
i
=
0
k
(
−
1
)
i
C
k
i
(
A
c
)
k
−
i
(
B
c
)
i
F_c^k=D^k\sum_{i=0}^{k}C_{k}^{i}(A^c)^{k-i}(-B^c)^i=D^k\sum_{i=0}^{k}(-1)^iC_{k}^{i}(A^c)^{k-i}(B^c)^i
Fck=Dk∑i=0kCki(Ac)k−i(−Bc)i=Dk∑i=0k(−1)iCki(Ac)k−i(Bc)i
当
n
=
2
n=2
n=2时,把
c
c
c改为
2
c
2c
2c,可以得到
F
2
c
k
=
D
k
∑
i
=
0
k
(
−
1
)
i
C
k
i
(
A
2
c
)
k
−
i
(
B
2
c
)
i
F_{2c}^k=D^k\sum_{i=0}^{k}(-1)^iC_{k}^{i}(A^{2c})^{k-i}(B^{2c})^i
F2ck=Dk∑i=0k(−1)iCki(A2c)k−i(B2c)i
=
D
k
∑
i
=
0
k
(
−
1
)
i
C
k
i
(
A
c
)
k
−
i
(
B
c
)
i
∗
(
A
c
)
k
−
i
(
B
c
)
i
=D^k\sum_{i=0}^{k}(-1)^iC_{k}^{i}(A^{c})^{k-i}(B^{c})^i*(A^{c})^{k-i}(B^{c})^i
=Dk∑i=0k(−1)iCki(Ac)k−i(Bc)i∗(Ac)k−i(Bc)i
注意到,当列举
i
i
i相同的项时,
∑
\sum
∑内的项是个等比。
首项是
(
−
1
)
i
(
A
k
−
i
B
i
)
c
(-1)^i(A^{k-i}B^{i})^c
(−1)i(Ak−iBi)c,公比是
(
A
k
−
i
B
i
)
c
(A^{k-i}B^{i})^c
(Ak−iBi)c。
因此我们可以考虑枚举
k
k
k,每一次去计算一个等比数列的和即可。
所以我们可以整理式子,得到
∑
i
=
0
n
F
i
c
K
\sum_{i=0}^{n}F_{ic}^K
∑i=0nFicK
=
D
k
∑
i
=
0
K
C
k
i
(
−
1
)
i
(
A
k
−
i
B
i
)
c
∗
(
(
(
A
k
−
i
B
i
)
c
)
n
−
1
)
(
A
k
−
i
B
i
)
c
−
1
=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}
=Dk∑i=0KCki(Ak−iBi)c−1(−1)i(Ak−iBi)c∗(((Ak−iBi)c)n−1)
当然要处理一下公比是1的情况,即
(
A
k
−
i
B
i
)
c
=
1
(A^{k-i}B^{i})^c=1
(Ak−iBi)c=1
此时
∑
\sum
∑内的项为
C
k
i
∗
n
∗
(
−
1
)
i
(
A
k
−
i
B
i
)
c
C_k^i*n*(-1)^i(A^{k-i}B^{i})^c
Cki∗n∗(−1)i(Ak−iBi)c
接下来预处理
1
e
5
1e5
1e5的阶乘,
1
e
5
1e5
1e5的阶乘逆元(为了快速处理
C
n
m
C_{n}^{m}
Cnm)。
当然算到这步可能还会超时。所以依旧需要优化。
优化1:关于计算首项
(
−
1
)
i
(
A
k
−
i
B
i
)
c
(-1)^i(A^{k-i}B^{i})^c
(−1)i(Ak−iBi)c和公比的
n
n
n次方,从
i
i
i项变为
i
+
1
i+1
i+1项时其实相当于乘上了
B
A
c
\frac{B}{A}^c
ABc,这样每次就可以递推
O
(
1
)
O(1)
O(1)求出首项(以及公比),同理公比的
n
n
n次也可以用类似方法处理。
优化2:欧拉降幂。本题的
m
o
d
mod
mod是个质数,所以可以直接用欧拉降幂
a
b
≡
a
b
%
(
1
e
9
+
9
−
1
)
(
m
o
d
1
e
9
+
9
)
a^b\equiv a^{b\%(1e9+9-1)}(mod\ 1e9+9)
ab≡ab%(1e9+9−1)(mod 1e9+9)
AC代码(873ms)
#include<bits/stdc++.h>
using namespace std;
const long long mod=1e9+9;
const long long MAXN=1e5;
long long factor[MAXN+10];
long long invFactor[MAXN+10];
long long invn[MAXN+10];
const long long A=691504013;
const long long invA=691504012;
const long long B=308495997;
const long long D=276601605;
//F(n)=D(A^n-B^n)
inline long long quick_pow(long long a, long long b)
{
long long ans=1, base=a;
while(b!=0)
{
if (b&1)
ans=(long long) ans*base%mod;
base=(long long) base*base%mod;
b>>=1;
}
return ans;
}
inline long long MOD(long long a, long long b)
{
a+=b;
if (a>=mod)
a-=mod;
return a;
}
inline void init()
{
factor[0]=invFactor[0]=invn[0]=factor[1]=invFactor[1]=invn[1]=1;
for(int i=2; i<=MAXN; i++)
{
factor[i]=factor[i-1]*i%mod;
invn[i]=(long long) (mod-mod/i)*invn[mod%i]%mod;
invFactor[i]=invFactor[i-1]*invn[i]%mod;
}
}
inline long long getC(long long m, long long n)
{
if (n<0 || m<0 || m>n)
return 0;
long long ans=factor[n];
ans=(long long) ans*invFactor[m]%mod;
ans=(long long) ans*invFactor[n-m]%mod;
return ans;
}
int main()
{
init();
int t;
scanf("%d", &t);
while(t--)
{
long long n, c, k;
scanf("%lld%lld%lld", &n, &c, &k);
long long ans=0;
long long a1=quick_pow(quick_pow(A, k), c%(mod-1));
long long q=quick_pow((long long) invA*B%mod, c%(mod-1));
long long n1=n%mod;
long long n2=n%(mod-1);
long long a1power=quick_pow(a1, n2);
long long qpower=quick_pow(q, n2);
for(int i=0; i<=k; i++)
{
long long sum=getC(i, k);
if (i&1)
{
sum=mod-sum;
}
if (a1==1)
{
ans=(ans+(long long) sum*n1%mod)%mod;
}
else
{
sum=(long long) sum*((long long) a1*(a1power-1+mod)%mod)%mod;
sum=(long long) sum*quick_pow(a1-1, mod-2)%mod;
ans=MOD(ans, sum);
}
a1=(long long) a1*q%mod;
a1power=(long long) a1power*qpower%mod;
}
printf("%lld\n", (long long) ans*quick_pow(D, k)%mod);
}
}
后记
再也不用long double的快速乘了,我是傻逼。
DrGilbert 2020.7.22