[BZOJ3328]PYXFIB

description

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)\)

solution

单位根反演的一个常见应用。
记\(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\)。

令\(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\)存在线性递推关系。

最后矩阵快速幂即可。

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;
}