BZOJ
給定\(n,k,mod\),求\[\sum_{i=0}^{\lfloor\frac{n}{k}\rfloor}\binom{n}{ki}f_{ki}\]
其中\(f_0=1,f_1=1,f_n=f_{n-1}+f_{n-2}(n>1)\)
多組數據。\(n\le 10^{18},k\le 20000,T\le 20\),\(mod\le 10^9\)且\(k|(mod-1)\)php
單位根反演的一個常見應用。
記\(F(x)=\sum_{i=0}^{n}\binom{n}{i}f_ix^i\)
原式即\[\sum_{i=0}^{n}[i\ mod\ k=0]\binom{n}{i}f_i\]
由單位根反演的公式咱們知道\(\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}=[i\ mod\ k=0]\)
原式可轉化爲
\[\begin{aligned} &\sum_{i=0}^{n}\frac{1}{k}\sum_{j=0}^{k-1}w_k^{ij}\binom{n}{i}f_i\\ =&\frac{1}{k}\sum_{j=0}^{k-1}\sum_{i=0}^{n}\binom{n}{i}f_i(w_k^j)^i\\ =&\frac{1}{k}\sum_{j=0}^{k-1}F(w_k^j)\\ \end{aligned}\]
因而咱們經過單位根反演把多項式間隔k位的係數求和問題變成了多點求值問題。
要是\(n\)不大,常見的多點求值便可。可是這裏\(n\le 10^{18}\)。
在下面的表述中,簡寫\(w_k^j\)爲\(w\)。c++
令\(g_i=w^if_i\),那麼有\(g_i=wg_{i-1}+w^2g_{i-2}\),轉移矩陣爲\(A=\begin{equation}\left(\begin{array}& w &w^2 \\ 0 &w \\\end{array}\right)\end{equation}\),
設初始矩陣爲\(S=(1\ w)\),點值表示即爲
\[ \begin{aligned} F(w)&=\sum_{i=0}^{n}\binom{n}{i}g_i\\ &=\sum_{i=0}^{n}\binom{n}{i}(SA^i)[0][0]\\ &=(S\sum_{i=0}^{n}\binom{n}{i}A^i)[0][0]\\ &=(S(I+A)^n)[0][0]\\ \end{aligned} \]
由二項式定理可得。
能夠發現這樣的矩陣構造適宜一系列式子\(\sum_{i=0}^n\binom{n}{i}a_i\)的求解,
僅需\(a_i\)存在線性遞推關係。spa
最後矩陣快速冪便可。code
#include<bits/stdc++.h> #define mp make_pair #define pb push_back #define fi first #define se second #define FL "a" using namespace std; typedef long long ll; const int N=10001; inline ll read(){ ll data=0,w=1;char ch=getchar(); while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar(); if(ch=='-')w=-1,ch=getchar(); while(ch<='9'&&ch>='0')data=data*10+ch-48,ch=getchar(); return data*w; } inline void file(){ freopen(FL".in","r",stdin); freopen(FL".out","w",stdout); } ll n;int k,mod,wn; inline void upd(int &a,int b){a+=b;if(a>=mod)a-=mod;} inline void dec(int &a,int b){a-=b;if(a<0)a+=mod;} inline int poww(int a,int b){ int res=1; for(;b;b>>=1,a=1ll*a*a%mod) if(b&1)res=1ll*res*a%mod; return res; } inline int getroot(int x){ static int cnt,pri[N];x--;cnt=0; for(int i=2;1ll*i*i<=x;i++) if(x%i==0){pri[++cnt]=i;if(i*i!=x)pri[++cnt]=x/i;} for(int i=2,b;i<=x;i++){ b=1; for(int j=1;j<=cnt;j++) if(poww(i,x/pri[j])==1){b=0;break;} if(b)return i; } return 0; } struct matrix{int a[2][2];inline int* operator[](int x){return a[x];}}S,T; inline matrix operator *(matrix a,matrix b){ return (matrix){ (int)((1ll*a[0][0]*b[0][0]+1ll*a[0][1]*b[1][0])%mod), (int)((1ll*a[0][0]*b[0][1]+1ll*a[0][1]*b[1][1])%mod), (int)((1ll*a[1][0]*b[0][0]+1ll*a[1][1]*b[1][0])%mod), (int)((1ll*a[1][0]*b[0][1]+1ll*a[1][1]*b[1][1])%mod)}; } inline int solve(ll n,int k){ int ans=0;wn=getroot(mod);wn=poww(wn,(mod-1)/k); for(int i=0,w=1;i<k;i++,w=1ll*w*wn%mod){ S=(matrix){1,w,0,0}; T=(matrix){1,(int)(1ll*w*w%mod),1,(w+1)%mod}; ll b=n;while(b){if(b&1)S=S*T;T=T*T;b>>=1;}upd(ans,S[0][0]); } return 1ll*ans*poww(k,mod-2)%mod; } int main() { int T=read(); while(T--){ n=read();k=read();mod=read(); printf("%d\n",solve(n,k)); } return 0; }