1. 程式人生 > >[數學技巧 等比數列] 斐波那契k次冪和

[數學技巧 等比數列] 斐波那契k次冪和

題意:給定,其中,求的值。

首先k小的情況可以構造矩陣和向量[fi-1^k,fi-1^(k-1)*fi,fi-1^(k-2)*fi^2,....,fi^k,Sum]

#include<cstdio>
#include<cstdlib>
#include<algorithm>
#include<cstring>
using namespace std;
typedef long long ll;

inline char nc(){
  static char buf[100000],*p1=buf,*p2=buf;
  if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
  return *p1++;
}

inline void read(int &x){
  char c=nc(),b=1;
  for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
  for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int M=55;
const int P=1e9+9;

struct Matrix{
  int n;
  ll a[M][M];
  Matrix(){ }
  Matrix(int in,int ii=0){
    n=in; memset(a,0,sizeof(a));
    if (ii) for (int i=0;i<=n;i++) a[i][i]=1;
  }
  ll *operator[](int x){
    return a[x];
  }
  friend Matrix operator* (Matrix A,Matrix B){
    int n=A.n; Matrix ret(n);
    for (int i=0;i<=n;i++)
      for (int j=0;j<=n;j++)
        for (int k=0;k<=n;k++)
	  (ret[i][j]+=A[i][k]*B[k][j]%P)%=P;
    return ret;
  }
}A;


Matrix Pow(Matrix A,int b){
  Matrix ret(A.n,1);
  for (;b;b>>=1,A=A*A)
    if (b&1)
      ret=ret*A;
  return ret; 
}

int n,K;
ll C[M][M];

inline void Pre(){
  C[0][0]=1;
  for (int i=1;i<=K;i++){
    C[i][0]=1;
    for (int j=1;j<=i;j++)
      C[i][j]=(C[i-1][j-1]+C[i-1][j])%P;
  }
}

ll Ans;

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  read(n); read(K); Pre();
  A=Matrix(K+1);
  for (int a=0;a<=K;a++)
    for (int t=0;t<=a;t++)
      A[a][K-a+t]=C[a][t];
  A[K+1][K+1]=A[K+1][K]=1;
  A=Pow(A,n-1);
  for (int i=0;i<=K+1;i++)
    (Ans+=A[K+1][i])%=P;
  printf("%lld\n",Ans);
  return 0;
}

然後對於這道變態的題

以下來自參考博文

分析:嗯,這道題貌似有難度,如果比較小的話我們可以構造矩陣,實際上這樣做也挺麻煩的。

     以前我們做一個大Fibonacci數列模一個大素數都是用矩陣,當然這裡素數滿足條件:5是模這個素數的二

     次剩餘,那麼現在要求不要用矩陣來計算這個結果呢?那就是今天我要討論的問題,本題也是基於這種思路。

     本題我們可以直接利用Fibonacci數列的公式進行計算。因為我們知道Fibonacci數列的公式為:

      

     雖然公式中含有根號,但是我們知道是一個整數,而且5是模1000000009的二次剩餘。那麼可以通過逆元

     和二次剩餘的轉化來做。比如

就可以用中的來代替。

     為了方便表示,我們令:

     

     那麼得到

     

     把按照二項式展開得到

     

     對於每一個都這樣表示,那麼相同的合併後是等比數列,比如對於所有係數為合併後為

     , 其中

     所以到了這裡本題就明確了,列舉每一個0,依次計算即可。

     當然對於,我們可以階乘預處理然後求逆元,至於同樣預處理,這樣時間少很多。

     可以看出本題條件好在1000000009是素數,而且5是模1000000009的二次剩餘。

     經計算2模1000000009的逆元是500000005,而的一個解為383008016

     也就是說對於

     

     同理,對於

     

     嗯,貌似還有一個問題沒有處理,t = 1時咋辦? 這個很簡單啦,不用等比求和公式即可。其實吧,這裡面

     我們還可以注意到,也可以進一步優化,讀者自己體會,就說到這裡。

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int P=1000000009;
const int INV2=500000005;
const int SQRT5=383008016;
const int INVSQRT5=276601605;
const int A=691504013;
const int B=308495997;

const int N=100005;

ll n,K;
ll fac[N],inv[N];
ll pa[N],pb[N];

inline void Pre(int n){
  fac[0]=1; for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
  inv[1]=1; for (int i=2;i<=n;i++) inv[i]=(P-P/i)*inv[P%i]%P;
  inv[0]=1; for (int i=1;i<=n;i++) inv[i]=inv[i]*inv[i-1]%P;
  pa[0]=1; for (int i=1;i<=n;i++) pa[i]=pa[i-1]*A%P;
  pb[0]=1; for (int i=1;i<=n;i++) pb[i]=pb[i-1]*B%P;
}

inline ll C(int n,int m){
  return fac[n]*inv[m]%P*inv[n-m]%P;
}

inline ll Pow(ll a,ll b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%P)
    if (b&1)
      ret=ret*a%P;
  return ret;
}

inline ll Inv(ll x){
  return Pow(x,P-2);
}

inline void Solve(){
  ll Ans=0;
  for (int j=0;j<=K;j++){
    ll t=pa[K-j]*pb[j]%P,tem;
    tem=t==1?n%P:t*(Pow(t,n)-1+P)%P*Inv(t-1)%P;
    if (~j&1)
      Ans+=C(K,j)*tem%P,Ans%=P;
    else
      Ans+=P-C(K,j)*tem%P,Ans%=P;
  }
  Ans=Ans*Pow(INVSQRT5,K)%P;
  printf("%lld\n",Ans);
}

int main(){
  int T;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre(100000);
  scanf("%d",&T);
  while (T--){
    scanf("%lld%lld",&n,&K);
    Solve();
  }
}