标签：杭电多校 matrix limits mat1 int 第五场 sum 单位根 omega
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
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。