HDU 4959 / BC 5D Poor Akagi

题意比较简单，就是求数列Li的k次方的前n项和 

看了之后感觉无从下手，看了这个题解后才知道http://blog.csdn.net/ahm001/article/details/38724607 这个数列是有规律的 准确的说是叫卢卡斯数，通项公式为

L(i)=((1+√5) ^ n + (1-√5) ^ n)/ 2^n  然后再把k次方打开，然后枚举第i项，定义二次域 然后 用快速幂和分治求前缀和 来计算每项的值。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<vector>
#include<list>
#include<set>
#include<map>
#include<stack>
#include<queue>
#include<deque>
#define mem(x,y) memset(x,y,sizeof(x))
#define pb push_back
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,double> pii;
#define bug puts("===========");
#define zjc puts("");
const double pi=(acos(-1.0));
const double eps=1e-8;
const ll INF=1e18+10;
const ll inf=1e9+10;
const int mod=1e9+7;
const int maxn=100000+10;
/*=======================================*/
const ll c=5;
struct eq{///数值为x+y*sqrt(c)
    ll x,y;
    eq(const ll &x = 0, const ll &y = 0): x(x), y(y) {}
    friend eq operator +  (const eq  &a, const eq  &b){ return eq((a.x + b.x)%mod, (a.y + b.y)%mod); }
    friend eq operator *  (const eq  &a, const ll &b) { return eq(a.x * b, a.y * b); }
    friend eq operator *  (const eq &a, const eq  &b) { return eq((a.x*b.x+c*a.y*b.y)%mod,(a.x*b.y+a.y*b.x)%mod); }
}r1[maxn],r2[maxn];
inline eq eqpow(eq a,ll b){
    eq ans(1,0);
    for(;b;b>>=1,a=a*a) if(b&1) ans=ans*a;
    return ans;
}
inline eq powsum(eq A,ll k) {
    eq ans(1,0), I(1,0), tmp(1,0);
    while(k > 0) {
        if (k & 1) ans= ans*A+tmp;
        tmp= tmp * (I + A);
        A = A * A;
        k >>= 1;
    }
    return ans;
}
inline ll powmod(ll a, ll b){
    ll ans=1;
    for(;b;b>>=1,a=a*a%mod) if(b&1) ans=ans*a%mod;
    return ans;
}
ll inv[maxn];
int main()
{
    r2[0]=r1[0]=eq(1,0);
    r1[1]=eq(1,1);
    r2[1]=eq(1,-1);
    for(int i=1;i<=100000;i++){
        r1[i]=r1[i-1]*r1[1];
        r2[i]=r2[i-1]*r2[1];
        inv[i]=powmod(i,mod-2);
    }
    int T_T;
    scanf("%d",&T_T);
    while(T_T--){
        ll n,k;
        scanf("%I64d%I64d",&n,&k);
        eq ans(0,0),fm(powmod(inv[2],k),0);
        ll Cki=1;
        for(int i=0;i<=k;i++){
            eq tmp=r1[i]*r2[k-i]*fm;
            ans=ans+powsum(tmp,n)*eq(Cki,0);
            Cki=Cki*(k-i)%mod*inv[i+1]%mod;
        }
        printf("%I64d\n",(ans.x+mod)%mod);
    }
    return 0;
}