2021杭电多校第五场1002（单位根反演）
作者：互联网
2021杭电多校第五场1002 Problem - 7013 (hdu.edu.cn)
题意：
给一个长度为 \(L\) 的字符串，包含前 \(k(k>=2)\) 个小写字母，可以得到不同的字符串有 \(k^L\) 种
对于每一对 \((i,j),(0\le i,j)\) ,找出包含 \(p\) 个 \('a'\) , \(q\) 个 \('b'\)，满足 \(q\equiv i(mod\ n),p\equiv j(mod\ n)\)的字符串个数
\(2\le k\le26,1\le L\le 10^{18},1\le n\le 500\)，对 \(P=1e9+9\) 取模，同时保证\(n|P-1\)
思路：
枚举 \('a'\) 和 \('b'\) 出现的次数可得
\(ans[i][j]=\sum\limits_{x=0}^L\sum\limits_{y=0}^{L-x}[n|x-i][x|y-j]\left(\matrix{L\\x}\right)\left(\matrix{L-x\\y}\right)(k-2)^{L-x-y}\)
对于\([n|k]\),由单位根反演的式子可得\(\frac{1}{n}\sum\limits_{i=0}^{n-1}{\omega}^{ki}\)
所以\(ans[i][j]=\sum\limits_{x=0}^L\sum\limits_{y=0}^{L-x}\frac{1}{n}\sum\limits_{p=0}^{n-1}{\omega}_n^{p\times(x-i)}\frac{1}{n}\sum\limits_{q=0}^{n-1}{\omega}_n^{q\times(y-j)}\left(\matrix{L\\x}\right)\left(\matrix{L-x\\y}\right)(k-2)^{L-x-y}\)
更换求和顺序有\(\frac{1}{n^2}\sum\limits_{p=0}^{n-1}\sum\limits_{q=0}^{n-1}\sum\limits_{x=0}^L\sum\limits_{y=0}^{L-x}{\omega}_n^{px}{\omega}_n^{qy}\left(\matrix{L\\x}\right)\left(\matrix{L-x\\y}\right)(k-2)^{L-x-y}{\omega}_n^{-pi}{\omega}_n^{-qj}\)
考虑m项式的n次展开有
\((x_1+x_2+\cdots+x_m)^n=\sum\limits_{a_1+a_2+\cdots+a_m=n}\left(\matrix{n\\a_1,\cdots,a_n}\right)\prod\limits_{k=1}^{m}x_k^{a_k}\)
其中\(\left(\matrix{n\\a_1,\cdots,a_m}\right)\)代表从n个小球中依次选取\(a_1,a_2,\cdots,a_m\)个小球的方案数化简后变为
\(\left(\matrix{n\\a_1,\cdots,a_m}\right)=\frac{n!}{a_1!a_2!\cdots a_m!}\)
我们返回来看\(\sum\limits_{x=0}^L\sum\limits_{y=0}^{L-x}{\omega}_n^{px}{\omega}_n^{qy}\left(\matrix{L\\x}\right)\left(\matrix{L-x\\y}\right)(k-2)^{L-x-y}\)这一段，令\(a={\omega}_n^p,b={\omega}_n^q,c=k-2\)
则其可以进一步化简为
\(\sum\limits_{x+y+z=L}\frac{L!}{x!(L-x)!}\frac{(L-x)!}{y!(L-x-y)!}a^xb^yc^{L-x-y}=\sum\limits_{x+y+z=L}\frac{L!}{x!y!z!}\)正好为三项式展开的形式
那么这部分就可以化为 \(({\omega}_n^p+{\omega}_n^q+k-2)^L\)
我们计 \(A[i][p]={\omega}_n^{-pi},B[p][q]=\frac{1}{n^2}({\omega}_n^p+{\omega}_n^q+k-2)^L,C[q][j]={\omega}_n^{-qj}\)
那么\(ans=A\times B\times C\)
莫名奇妙被卡常了，拿std也跑了5s，时限5.5s，故下面分别贴std和没通过的代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<ll,ll> pii;
#define rep(i,x,y) for(int i=x;i<y;i++)
#define rept(i,x,y) for(int i=x;i<=y;i++)
#define all(x) x.begin(),x.end()
#define fi first
#define se second
#define mes(a,b) memset(a,b,sizeof a)
#define mp make_pair
#define pb push_back
#define dd(x) cout<<#x<<"="<<x<<" "
#define de(x) cout<<#x<<"="<<x<<"\n"
const int inf=0x3f3f3f3f;
const int maxn=500;
const int mod=1e9+9;
ll x[maxn],y[maxn];
int p[maxn],q[maxn];
class matrix
{
public:
ll arrcy[maxn][maxn];//?????С???±?0??row-1,0??column-1
int row,column;//row??У?column???
matrix()
{
memset(arrcy,0,sizeof arrcy);
column=row=0;
}
friend matrix operator *(matrix s1,matrix s2)
{
int i,j;
matrix s3;
for (i=0;i<s1.row;i++)
{
for (j=0;j<s2.column;j++)
{
for (int k=0;k<s1.column;k++)
{
s3.arrcy[i][j]+=s1.arrcy[i][k]*s2.arrcy[k][j];
s3.arrcy[i][j]%=mod;
}
}
}
s3.row=s1.row;
s3.column=s2.column;
return s3;
}
void show()
{
for(int i=0;i<row;i++)
{
for (int j=0;j<column;j++)
cout<<arrcy[i][j]<<" ";
cout<<"\n";
}
}
}mat1,mat2,mat3;
matrix mul(matrix &s1,matrix &s2)
{
int i,j;
matrix s3;
for (i=0;i<s1.row;i++)
{
for (j=0;j<s2.column;j++)
{
for (int k=0;k<s1.column;k++)
{
s3.arrcy[i][j]+=s1.arrcy[i][k]*s2.arrcy[k][j];
s3.arrcy[i][j]%=mod;
}
}
}
s3.row=s1.row;
s3.column=s2.column;
return s3;
}
ll qpow(ll a,ll b)
{
ll ans=1;
for(;b;b>>=1,a=a*a%mod)
if(b&1)
ans=ans*a%mod;
return ans;
}
/*
matrix quick_pow(matrix s1,long long n)
{
matrix mul=s1,ans;
ans.row=ans.column=s1.row;
memset(ans.arrcy,0,sizeof ans.arrcy);
for(int i=0;i<ans.row;i++)
ans.arrcy[i][i]=1;
while(n)
{
if(n&1) ans=ans*mul;
mul=mul*mul;
n/=2;
}
return ans;
}
*/
int g=13;
void solve()
{
int k,n;
ll l;
cin>>k>>l>>n;
k-=2;
x[0]=1;
g=qpow(13,(mod-1)/n);
rept(i,1,n) x[i]=x[i-1]*g%mod;
mat1.row=mat1.column=n;
rep(i,0,n)
rep(j,0,n)
mat1.arrcy[i][j]=qpow(g,i*j);
// mat1.show();
mat3=mat1;
rep(i,0,n) rep(j,i+1,n) swap(mat3.arrcy[i][j],mat3.arrcy[j][i]);
mat2.row=mat2.column=n;
rep(i,0,n) rep(j,0,n) mat2.arrcy[i][j]=qpow(x[i]+x[j]+k,l);
//mat1.show();
//mat2.show();
//mat3.show();
// mat1=mat1*mat2;
// mat1=mat1*mat3;
mat1=mul(mat1,mat2);
mat1=mul(mat1,mat3);
p[0]=q[0]=0;
rep(i,1,n) p[i]=q[i]=n-i;
ll nn=qpow(n*n,mod-2);
rep(i,0,n)
{
rep(j,0,n)
{
cout<<mat1.arrcy[p[i]][q[j]]*nn%mod;
if(j==n-1) cout<<"\n";
else cout<<" ";
}
}
return ;
}
int main(){
ios::sync_with_stdio(false);
cin.tie(0);cout.tie(0);
int T;
cin>>T;
while(T--)
solve();
return 0;
}
#include<bits/stdc++.h>
using namespace std;
const int mod=1e9+9;
const int maxn=1e3+50;
struct tnode{
int n,m;
int a[maxn][maxn];
}A,B,C,D;
int p[maxn],q[maxn];
int ksm(int a,long long b)
{
int res=1;
while(b)
{
if(b&1)res=1ll*res*a%mod;
a=1ll*a*a%mod;
b>>=1;
}
return res;
}
int g=13;
int w[maxn];
int n,k;
long long L;
void mul(tnode &a,tnode &b,tnode &c)
{
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
c.a[i][j]=0;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
{
for(int k=0;k<n;k++)
c.a[i][j]=(c.a[i][j]+1ll*a.a[i][k]*b.a[k][j]%mod)%mod;
}
}
void show(tnode &a)
{
for(int i=0;i<n;i++)
{
for(int j=0;j<n;j++)
cout<<a.a[i][j]<<" ";
cout<<"\n";
}
}
int main()
{
int t;
w[0]=1;
scanf("%d",&t);
while(t--)
{
scanf("%d%lld%d",&k,&L,&n);
int g=13;
g=ksm(g,(mod-1)/n);
for(int i=1,j=g;i<=1000;i++)
w[i]=j,j=1ll*j*g%mod;
int nn=ksm(n*n,mod-2);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
// A.a[i][j]=ksm(g,i*j);
{
if(j==0)A.a[i][j]=1;
else A.a[i][j]=1ll*A.a[i][j-1]*w[i]%mod;
}
// show(A);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
B.a[i][j]=ksm(w[j]+w[i]+k-2,L);
// show(B);
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
C.a[i][j]=A.a[j][i];
// show(C);
mul(A,B,D);
mul(D,C,A);
p[0]=q[0]=0;
for(int i=1;i<n;i++)
p[i]=q[i]=n-i;
for(int i=0;i<n;i++)
{
for(int j=0;j<n-1;j++)
printf("%d ",1ll*nn*A.a[p[i]][q[j]]%mod);
printf("%d\n",1ll*nn*A.a[p[i]][q[n-1]]%mod);
}
}
return 0;
}
标签：杭电多校,matrix,limits,mat1,int,第五场,sum,单位根,omega 来源： https://www.cnblogs.com/1427314831a/p/15195352.html