[BZOJ3328]PYXFIB

description

BZOJ
给定(n,k,mod),求$$sum_{i=0}^{lfloorfrac{n}{k} floor}inom{n}{ki}f_{ki}$$
其中(f_0=1,f_1=1,f_n=f_{n-1}+f_{n-2}(n>1))
多组数据。(nle 10^{18},kle 20000,Tle 20),(modle 10^9)且(k|(mod-1))

solution

单位根反演的一个常见应用。
记(F(x)=sum_{i=0}^{n}inom{n}{i}f_ix^i)
原式即$$sum_{i=0}^{n}[i mod k=0]inom{n}{i}f_i$$
由单位根反演的公式我们知道(frac{1}{k}sum_{j=0}^{k-1}w_k^{ij}=[i mod k=0])
原式可转化为

[egin{aligned} &sum_{i=0}^{n}frac{1}{k}sum_{j=0}^{k-1}w_k^{ij}inom{n}{i}f_i\ =&frac{1}{k}sum_{j=0}^{k-1}sum_{i=0}^{n}inom{n}{i}f_i(w_k^j)^i\ =&frac{1}{k}sum_{j=0}^{k-1}F(w_k^j)\ end{aligned}]

于是我们通过单位根反演把多项式间隔k位的系数求和问题变成了多点求值问题。
要是(n)不大,常见的多点求值即可。但是这里(nle 10^{18})。
在下面的表述中,简写(w_k^j)为(w)。

令(g_i=w^if_i),那么有(g_i=wg_{i-1}+w^2g_{i-2}),转移矩阵为(A=egin{equation}left(egin{array}& w &w^2 \ 0 &w \end{array} ight)end{equation}),
设初始矩阵为(S=(1 w)),点值表示即为

[egin{aligned} F(w)&=sum_{i=0}^{n}inom{n}{i}g_i\ &=sum_{i=0}^{n}inom{n}{i}(SA^i)[0][0]\ &=(Ssum_{i=0}^{n}inom{n}{i}A^i)[0][0]\ &=(S(I+A)^n)[0][0]\ end{aligned} ]

由二项式定理可得。
可以发现这样的矩阵构造适宜一系列式子(sum_{i=0}^ninom{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
",solve(n,k));
  }
  return 0;
}