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代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <vector>
#include <map>
#include <cstdlib>
#include <set>
#include <queue>
using namespace std;

#define X first
#define Y second
#define Mod 1000000009
#define N 100010
#define M 101

typedef long long ll;

const int INF=1<<30;
const int sqrt5=383008016;
int P,Q;
int fac[N],fac_inv[N];

int Power(int x,ll p)
{
    int res=1;
    for (;p;p>>=1)
    {
        if (p&1) res=1ll*res*x%Mod;
        x=1ll*x*x%Mod;
    }
    return res;
}

int C(int n,int m)
{
    if (m==0) return 1;
    int res=1ll*fac[n]*fac_inv[m]%Mod;
    return 1ll*res*fac_inv[n-m]%Mod;
}

int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);
    
    P=1ll*(1+sqrt5)*Power(2,Mod-2)%Mod;
    Q=1ll*(1-sqrt5)*Power(2,Mod-2)%Mod;
    if (Q<0) Q+=Mod;
    fac[0]=fac_inv[0]=1;
    for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod,fac_inv[i]=Power(fac[i],Mod-2);
    int T,K,ans; ll n;  
    scanf("%d",&T);
    while (T--)
    {
        scanf("%lld%d",&n,&K);
        ans=0;
        for (int i=K,op=1;i>=0;i--,op=-op)
        {
            int tmp1=op*C(K,i),tmp2,x=1ll*Power(P,i)*Power(Q,K-i)%Mod;
            if (x==1) tmp2=n%Mod;
            else tmp2=1ll*x*(Power(x,n)-1)%Mod,tmp2=1ll*tmp2*Power(x-1,Mod-2)%Mod;
            ans+=1ll*tmp1*tmp2%Mod;
            ans%=Mod;
        }
        if (ans<0) ans+=Mod;
        ans=1ll*ans*Power(Power(sqrt5,Mod-2),K)%Mod;
        printf("%d\n",ans);
    }
    return 0;
}

 

我的求高次剩余方程的代码(目前只会做模数是质数的情况...)： 

http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1038

#include <cstdio>
#include <iostream>
#include <queue>
#include <cmath>
#include <map>
#include <algorithm>
#include <cstring>
#include <set>
using namespace std;


#define X first
#define Y second
#define N 100010
#define M 5783
typedef long long ll;
typedef unsigned long long ull;

const int Mod=1000000007;

int prime[N];
bool flag[N];
vector<int> ans;

int Power(int x,int p,int Mod)
{
    int res=1;
    for (;p;p>>=1)
    {
        if (p&1) res=1ll*res*x%Mod;
        x=1ll*x*x%Mod;
    }
    return res;
}

bool Check(int g,int P,int a[])
{
    for (int i=1;i<=a[0];i++) if (Power(g,(P-1)/a[i],P)==1) return false;
    return true;
}

int Find_Root(int P)
{
    int x=P-1,a[100]; a[0]=0;
    for (int i=1;i<=prime[0] && prime[i]*prime[i]<=x;i++)
    {
        if (x%prime[i]==0)
        {
            a[++a[0]]=prime[i];
            while (x%prime[i]==0) x/=prime[i];
        }
    }
    if (x!=1) a[++a[0]]=x;
    for (int i=2;;i++) if (Check(i,P,a)) return i;
}

vector<pair<int,int> > hash[M];

int Calc_Exp(int a,int b,int P)  // 求解a^x=b (mod P)
{
    int m=ceil(sqrt(P+0.5)); 
    for (int i=0;i<M;i++) hash[i].clear();
    
    int v=Power(a,m,P),tmp=1; v=Power(v,P-2,P);
    hash[1].push_back(make_pair(1,0));
    for (int i=1;i<m;i++) 
    {
        tmp=1ll*tmp*a%P;
        hash[tmp%M].push_back(make_pair(tmp,i));
    }
    
    for (int i=0;i<m;i++)
    {
        int t=b%M;
        for (int j=0;j<hash[t].size();j++)
        {
            if (hash[t][j].X==b) return i*m+hash[t][j].Y;
        }
        b=1ll*b*v%P;
    }
}

void ex_gcd(ll a,ll b,ll &x,ll &y,ll &d)
{
    if (!b)
    {
        d=a;
        x=1;
        y=0;
    }
    else
    {
        ex_gcd(b,a%b,y,x,d);
        y-=a/b*x;
    }
}

//X^A = B (mod P)
void Solve_Equation(int A,int B,int P)
{
    ans.clear();
    int g=Find_Root(P); //求出P的原根 
    int b=Calc_Exp(g,B,P);  // B=g^b
    ll x,y,d;
    ex_gcd(A,P-1,x,y,d);
    if (b%d==0)
    {
        int del=(P-1)/d,tmp,t;
        x*=b/d;
        x%=del;
        if (x<0) x+=del;
        tmp=Power(g,(int)x,P);
        t=Power(g,del,P);
        while (x<=P-2) ans.push_back(tmp),x+=del,tmp=1ll*tmp*t%P;
        sort(ans.begin(),ans.end());
        for (int i=0;i<ans.size();i++) printf("%d%c",ans[i],i==ans.size()-1? '\n':' ');
    }
    if (ans.size()==0) printf("No Solution\n");
}

int main()
{
    //freopen("in.in","r",stdin);
    //freopen("out.out","w",stdout);
    
    for (int i=2;i<N;i++)
    {
        if (!flag[i]) prime[++prime[0]]=i;
        for (int j=1;j<=prime[0] && i*prime[j]<N;j++)
        {
            flag[i*prime[j]]=true;
            if (i%prime[j]==0) break;
        }
    }
    
    int T,A,B,P; scanf("%d",&T);
    while (T--)
    {
        scanf("%d%d%d",&P,&A,&B);
        Solve_Equation(A,B,P);
    }
    
    return 0;
}

可以扩展到模数不是质数的情况,等我学会了在写篇文章 记录具体做法吧。