循环卷积 -- 第K大C - FFT

虽然这题我并没有在其他题站上找到，但还是把贴在这了……

题意：给定数组a[N][M] 和b[N]，定义：

其中N是质数<2^18,M=2/3/4,求出第K 大的c[i]。

首先我们显然可以对b进行分类，将a拆成m个数组，那么我们就可以将问题转化为如下问题：

求x[i]=∑y[j]z[i*j mod n]

那么由于n是素数，那么我们求出n的原根g，则对于任何i∈[1,n)有g^p≡i (mod n)

那么我们就可以用离散对数解决问题，i=0和j=0可以单独处理。

令x'[p]表示x[i]，使得i为以p为底关于n的离散对数，对y和z也进行同样处理，那么有：

x'[p]=∑y'[q]z'[p+q mod n-1]

对于这样的问题我们显然可以将y'[]逆序为y''[],即y''[i]=y'[n-1-i]

则有x'[p]=∑y''[n-1-q]z'[p+q mod n-1]

所以我们只要求出A[i]=∑y''[n-1-j]z'[i+j]即可

#include"bits/stdc++.h"
using namespace std;
typedef long long ll;
#define gch getchar() 
template<class integer>
inline void read(integer&x){
	x=0;int f=1;char ch=gch;
	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=gch;}
	while(ch<='9'&&ch>='0'){x=x*10+ch-'0';ch=gch;}
	x*=f;
}

struct idouble{
	double r,i;
	idouble(double _r=0.0,double _i=0.0){r=_r,i=_i;}
	idouble operator +(const idouble&a)const{return idouble(r+a.r,i+a.i);}
	idouble operator -(const idouble&a)const{return idouble(r-a.r,i-a.i);} 
	idouble operator *(const idouble&a)const{return idouble(r*a.r-i*a.i,r*a.i+i*a.r);}
} ;

template<class iint>void rader(iint a[],int len){
	int i,j=len>>1,k;
	for(i=1;i<len-1;i++){
		if(i<j) swap(a[i],a[j]);
		for(k=len>>1;j>=k;j-=k,k>>=1);
		if(j<k) j+=k;
	}
}

const double PI = acos(-1.0);

void fft(idouble a[],int len,int dft){
	rader(a,len);int i,j,k,l;
	for(l=2;l<=len;l<<=1){
		idouble w(cos(-2*dft*PI/l),sin(-2*dft*PI/l));
		for(i=l>>1,j=0;j<len;j+=l){
			idouble wn(1.0,0.0);
			for(k=j;k<j+i;k++){
				idouble u=a[k],v=a[k+i]*wn;
				a[k]=u+v;a[k+i]=u-v;wn=wn*w;
			}
		}
	}
	if(dft==-1)for(i=0;i<len;i++)a[i].r/=len;
}

const int N=530005;
int n,m,g,K,a[N][4],b[N],logx[N],rev[N];

int kuai(int a,int t,int mod){
	int re=1;
	while(t){
		if(t%2)re=(ll)re*a%mod;
		a=(ll)a*a%mod;t>>=1;
	}
	return re;
}

int tmp[10],cnt;

void init(){
	int i,j,p;
	read(n),read(m),read(K);
	for(i=0;i<m;i++)
		for(j=0;j<n;j++)
			read(a[j][i]);
	for(i=0;i<n;i++)read(b[i]);
	for(i=2,p=n-1;i*i<=p;i++)
		if(p%i==0)for(tmp[++cnt]=i;p%i==0;p/=i);
	if(p>1)tmp[++cnt]=p;
	for(g=2;;g++){
		for(i=1;i<=cnt;i++)
			if(kuai(g,(n-1)/tmp[i],n)==1)
				break;
		if(i>cnt)break; 
	}
	for(i=0,j=1;i<n-1;i++){
		logx[j]=i;
		rev[i]=j;
		j=(ll)j*g%n;
	}
}

int ans[N],L;
idouble x[N],y[N];

void work(){
	for(L=1;L<=n*2;L<<=1);int i,j;
	for(i=1;i<n;i++)if(b[i]<m)ans[i]+=a[0][b[0]];
	for(i=0;i<n;i++)ans[0]+=a[i][b[0]];
	for(i=0;i<m;i++){
		for(j=1;j<n;j++){
			x[n-1-logx[j]].r=a[j][i]*1.0;
			y[logx[j]].r=b[j]==i;
		}
		fft(x,L,1);fft(y,L,1);
		for(j=0;j<L;j++)
			x[j]=x[j]*y[j];
		fft(x,L,-1);
		for(j=0;j<n-1;j++){
			x[j]=x[j]+x[j+n-1];
			ans[rev[j]]+=(int) (x[j].r+0.5);
		}
		for(j=0;j<L;j++)x[j]=y[j]=idouble(0.0,0.0);
	}
	sort (ans,ans+n);
	printf("%d\n",ans[n-K]);
}

int main(){
	init();work();
	return 0;
}