Fibonacci数列的幂和 zoj 3774
题目大意:
求斐波那契数列前n项的k次幂和 Mod 1000000009。 n<=1e18, k<=1e5
这题的k比较大,所以不能用矩阵乘法来递推。学到了新姿势... http://blog.csdn.net/acdreamers/article/details/23039571
基本思想就是求出通项公式,把里面的$\sqrt{5}$ 用 $x$ 代替, 其中 $x^2\equiv 5\pmod{1000000009}$
然后二项式展开求和就好了。
一个合法的$x$是383008016。 正好前几天做了个高次剩余方程的题,直接拉过去跑出来。。
为什么这样替代是合法的呢? 网络上好像没找到严谨的证明。
下面给出我个人的不严谨证明:
将原式子合并之后最终会得到$\frac{P(t)}{Q(t)}$ 这样的形式。 其中$t=\sqrt{5}$ $P(t)$和$Q(t)$是关于$t$的多项式.
把$P(t)$ $Q(t)$中次数>=2的可以利用$t^2=5$降幂成次数不超过1的, 最后变成$\frac{P'(t)}{Q'(t)}$ 这样的形式.
由于最终答案一定是整数,所以 $Q'(t)$ 能整除 $P'(t)$ 。 设 $\frac{P'(t)}{Q'(t)}=K$
那么把$t$ 换成 上面求出的 $x$ 得到的 $\frac{P'(x)}{Q'(x)}$ 也是 $K$. 不影响答案。
AC代码:
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <iostream>
2 #include <cstdio>
3 #include <cmath>
4 #include <cstring>
5 #include <algorithm>
6 #include <vector>
7 #include <map>
8 #include <cstdlib>
9 #include <set>
10 #include <queue>
11 using namespace std;
12
13 #define X first
14 #define Y second
15 #define Mod 1000000009
16 #define N 100010
17 #define M 101
18
19 typedef long long ll;
20
21 const int INF=1<<30;
22 const int sqrt5=383008016;
23 int P,Q;
24 int fac[N],fac_inv[N];
25
26 int Power(int x,ll p)
27 {
28 int res=1;
29 for (;p;p>>=1)
30 {
31 if (p&1) res=1ll*res*x%Mod;
32 x=1ll*x*x%Mod;
33 }
34 return res;
35 }
36
37 int C(int n,int m)
38 {
39 if (m==0) return 1;
40 int res=1ll*fac[n]*fac_inv[m]%Mod;
41 return 1ll*res*fac_inv[n-m]%Mod;
42 }
43
44 int main()
45 {
46 //freopen("in.in","r",stdin);
47 //freopen("out.out","w",stdout);
48
49 P=1ll*(1+sqrt5)*Power(2,Mod-2)%Mod;
50 Q=1ll*(1-sqrt5)*Power(2,Mod-2)%Mod;
51 if (Q<0) Q+=Mod;
52 fac[0]=fac_inv[0]=1;
53 for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod,fac_inv[i]=Power(fac[i],Mod-2);
54 int T,K,ans; ll n;
55 scanf("%d",&T);
56 while (T--)
57 {
58 scanf("%lld%d",&n,&K);
59 ans=0;
60 for (int i=K,op=1;i>=0;i--,op=-op)
61 {
62 int tmp1=op*C(K,i),tmp2,x=1ll*Power(P,i)*Power(Q,K-i)%Mod;
63 if (x==1) tmp2=n%Mod;
64 else tmp2=1ll*x*(Power(x,n)-1)%Mod,tmp2=1ll*tmp2*Power(x-1,Mod-2)%Mod;
65 ans+=1ll*tmp1*tmp2%Mod;
66 ans%=Mod;
67 }
68 if (ans<0) ans+=Mod;
69 ans=1ll*ans*Power(Power(sqrt5,Mod-2),K)%Mod;
70 printf("%d\n",ans);
71 }
72 return 0;
73 }
我的求高次剩余方程的代码(目前只会做模数是质数的情况...):
http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1038
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
1 #include <cstdio>
2 #include <iostream>
3 #include <queue>
4 #include <cmath>
5 #include <map>
6 #include <algorithm>
7 #include <cstring>
8 #include <set>
9 using namespace std;
10
11
12 #define X first
13 #define Y second
14 #define N 100010
15 #define M 5783
16 typedef long long ll;
17 typedef unsigned long long ull;
18
19 const int Mod=1000000007;
20
21 int prime[N];
22 bool flag[N];
23 vector<int> ans;
24
25 int Power(int x,int p,int Mod)
26 {
27 int res=1;
28 for (;p;p>>=1)
29 {
30 if (p&1) res=1ll*res*x%Mod;
31 x=1ll*x*x%Mod;
32 }
33 return res;
34 }
35
36 bool Check(int g,int P,int a[])
37 {
38 for (int i=1;i<=a[0];i++) if (Power(g,(P-1)/a[i],P)==1) return false;
39 return true;
40 }
41
42 int Find_Root(int P)
43 {
44 int x=P-1,a[100]; a[0]=0;
45 for (int i=1;i<=prime[0] && prime[i]*prime[i]<=x;i++)
46 {
47 if (x%prime[i]==0)
48 {
49 a[++a[0]]=prime[i];
50 while (x%prime[i]==0) x/=prime[i];
51 }
52 }
53 if (x!=1) a[++a[0]]=x;
54 for (int i=2;;i++) if (Check(i,P,a)) return i;
55 }
56
57 vector<pair<int,int> > hash[M];
58
59 int Calc_Exp(int a,int b,int P) // 求解a^x=b (mod P)
60 {
61 int m=ceil(sqrt(P+0.5));
62 for (int i=0;i<M;i++) hash[i].clear();
63
64 int v=Power(a,m,P),tmp=1; v=Power(v,P-2,P);
65 hash[1].push_back(make_pair(1,0));
66 for (int i=1;i<m;i++)
67 {
68 tmp=1ll*tmp*a%P;
69 hash[tmp%M].push_back(make_pair(tmp,i));
70 }
71
72 for (int i=0;i<m;i++)
73 {
74 int t=b%M;
75 for (int j=0;j<hash[t].size();j++)
76 {
77 if (hash[t][j].X==b) return i*m+hash[t][j].Y;
78 }
79 b=1ll*b*v%P;
80 }
81 }
82
83 void ex_gcd(ll a,ll b,ll &x,ll &y,ll &d)
84 {
85 if (!b)
86 {
87 d=a;
88 x=1;
89 y=0;
90 }
91 else
92 {
93 ex_gcd(b,a%b,y,x,d);
94 y-=a/b*x;
95 }
96 }
97
98 //X^A = B (mod P)
99 void Solve_Equation(int A,int B,int P)
100 {
101 ans.clear();
102 int g=Find_Root(P); //求出P的原根
103 int b=Calc_Exp(g,B,P); // B=g^b
104 ll x,y,d;
105 ex_gcd(A,P-1,x,y,d);
106 if (b%d==0)
107 {
108 int del=(P-1)/d,tmp,t;
109 x*=b/d;
110 x%=del;
111 if (x<0) x+=del;
112 tmp=Power(g,(int)x,P);
113 t=Power(g,del,P);
114 while (x<=P-2) ans.push_back(tmp),x+=del,tmp=1ll*tmp*t%P;
115 sort(ans.begin(),ans.end());
116 for (int i=0;i<ans.size();i++) printf("%d%c",ans[i],i==ans.size()-1? '\n':' ');
117 }
118 if (ans.size()==0) printf("No Solution\n");
119 }
120
121 int main()
122 {
123 //freopen("in.in","r",stdin);
124 //freopen("out.out","w",stdout);
125
126 for (int i=2;i<N;i++)
127 {
128 if (!flag[i]) prime[++prime[0]]=i;
129 for (int j=1;j<=prime[0] && i*prime[j]<N;j++)
130 {
131 flag[i*prime[j]]=true;
132 if (i%prime[j]==0) break;
133 }
134 }
135
136 int T,A,B,P; scanf("%d",&T);
137 while (T--)
138 {
139 scanf("%d%d%d",&P,&A,&B);
140 Solve_Equation(A,B,P);
141 }
142
143 return 0;
144 }
可以扩展到模数不是质数的情况,等我学会了在写篇文章 记录具体做法吧。
Every day is meaningful, keeping learning!