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;
}
上一篇:boost::multiprecision模块将 std::numeric_limits 用作 multiprecision.qbk 上的多精度文档片段的示例


下一篇:狄利克雷卷积重要公式及定义