bzoj 3328 PYXFIB——单位根反演

题目：https://www.lydsy.com/JudgeOnline/problem.php?id=3328

单位根反演主要就是有

( [k|n] = frac{1}{k}sumlimits_{i=0}^{k-1}w_{k}^{i*n} )

如果 k | n ，转 n 下就会是 1 ；不然用等比数列求和公式可知是 0 。

一般是构造一个 ( f(x) = ( 1+x )^n ) 之类的，来求含有组合数的式子。比如

　　　　　　　　　　　　　( sumlimits_{i=0}^{n}C_{n}^{i*k} = sumlimits_{i=0}^{n}C_{n}^{i}[ i | k ] )

　　　　　　　　　　　　　　　( = frac{1}{k}sumlimits_{i=0}^{n}C_{n}^{i}sumlimits_{j=0}^{k-1}w_{k}^{j*i} )

　　　　　　　　　　　　　　　( = frac{1}{k}sumlimits_{j=0}^{k-1}sumlimits_{i=0}^{n}C_{n}^{i}w_{k}^{j*i} )

　　　　　　　　　　　　　　　( = frac{1}{k}sumlimits_{j=0}^{k-1}(1+w_{k}^{j})^n )

所以设 ( f(x) = ( 1+x )^n ) ，求 k 次 ( f( w_{k}^{j} ) ) 就行。

对于这道题，为了凑一个二项式的形式，把 ( F[i] ) 看作斐波那契递推矩阵 A 的 ( A^{i}[0][0] ) ，就有

( ans = frac{1}{k}sumlimits_{j=0}^{k-1}sumlimits_{i=0}^{n}C_{n}^{i}A^{i}w_{k}^{j*i} )

有两个 i 次却没有 n-i 次，不能直接套。那个 ( w_{k}^{j*i} ) 的 ( w_{k}^{j} ) 与 i 无关，所以设 ( f(x) ) 的时候考虑把 ( w_{k}^{j} ) 作为 ( x ) 。

如果 ( f(x) = ( A+I*x )^n ) ，那么 ( f( w_{k}^{j} ) = sumlimits_{i=0}^{n}C_{n}^{i}A^{i}w_{k}^{j*(n-i)} )

想把 ( w_{k}^{j*(n-i)} ) 变成 ( w_{k}^{j*i} ) ，只要令 ( f(x) = x^{-n} ( A+I*x )^n ) ，然后求 ( f( w_{k}^{-j} ) ) 即可。

注意 n 是 long long 。

找原根是枚举 phi( mod ) 的质因子，然后看 ( g^{frac{phi(mod)}{pri}} ) 。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=2e4+5,K=35;
ll n;int k,mod,g,wn;
void upd(int &x){x>=mod?x-=mod:0;}
int pw(int x,int k)
{int ret=1;while(k){if(k&1)ret=(ll)ret*x%mod;x=(ll)x*x%mod;k>>=1;}return ret;}
void fnd_g()
{
  int pri[K],cnt=0,p=mod-1,d=p;////p=phi(mod)=mod-1!!
  for(int i=2;i*i<=d;i++)
    if(d%i==0)
      {pri[++cnt]=i;while(d%i==0)d/=i;}
  if(d>1)pri[++cnt]=d;
  for(g=2;;g++)
    {
      bool flag=1;
      for(int i=1;i<=cnt;i++)if(pw(g,p/pri[i])==1){flag=0;break;}
      if(flag)break;
    }
}
struct Mtr{
  int a[2][2];
  Mtr(){memset(a,0,sizeof a);}
  Mtr operator* (const Mtr &b)const
  {
    Mtr c;
    for(int i=0;i<=1;i++)
      for(int k=0;k<=1;k++)
    for(int j=0;j<=1;j++)
      c.a[i][j]=(c.a[i][j]+(ll)a[i][k]*b.a[k][j])%mod;
    return c;
  }
}A,I;
Mtr pw(Mtr x,ll k)//ll!!!
{Mtr ret=I;while(k){if(k&1)ret=ret*x;x=x*x;k>>=1;}return ret;}
int main()
{
  A.a[0][0]=A.a[0][1]=A.a[1][0]=1;
  I.a[0][0]=I.a[1][1]=1;
  int T;scanf("%d",&T);
  while(T--)
    {
      scanf("%lld%d%d",&n,&k,&mod);fnd_g();
      int ans=0,wn=pw(g,(mod-1)-(mod-1)/k),ml=(mod-1-n%(mod-1))%(mod-1);
      for(int i=0,w=1;i<k;i++,w=(ll)w*wn%mod)
    {
      Mtr t=A;t.a[0][0]+=w;t.a[1][1]+=w;upd(t.a[0][0]);upd(t.a[1][1]);
      t=pw(t,n);
      ans=(ans+(ll)t.a[0][0]*pw(w,ml))%mod;
    }
      ans=(ll)ans*pw(k,mod-2)%mod;
      printf("%d
",ans);
    }
  return 0;
}