【Luogu4191】[CTSC2010] 性能优化

题目链接

题意简述

求循环卷积意义下的 (A(x)*B(x)^C)。
模数为 n+1 ，长度为 n。

Sol

板子题。

循环卷积可直接把点值快速幂来解决。
所以问题就是要快速 (DFT)，由于长度是 n 不一定是NTT模数，我们要解决任意长度的 (DFT)

这道题保证了 (n) 质因数分解之后的质因子最大不超过 10 。
我们可以模仿朴素 (FFT) 对点值分组分别计算然后合并的方法。

每次分成 (p) 组然后合并点值即可。根据如下式子:
(F(x)=sum a_ix^i) ，(F_r(x)=sum a_{ip+r}x^i)
(F(x)=sum x^rF_r(x^p))

[F(w_n^{an+b})=sum_{r=0} w_{np}^{(an+b)r}F_r(w_n^b) ]

写的时候可以类似的用 dfs 预处理出每一个数最后到达的位置

code:

#include<bits/stdc++.h>
#define Set(a,b) memset(a,b,sizeof(a))
using namespace std;
const int N=1e6+10;
int mod,g;
template <typename T> inline void init(T&x){
    x=0;char ch=getchar();bool t=0;
    for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
    for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
    if(t) x=-x;return;
}
typedef long long ll;
template <typename T>inline void Inc(T&x,int y){x+=y;if(x>=mod) x-=mod;return;}
template <typename T>inline void Dec(T&x,int y){x-=y;if(x <  0) x+=mod;return;}
template <typename T>inline int fpow(int x,T k){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod) if(k&1) ret=(ll)ret*x%mod;return ret;}
int Sum(int x,int y){x+=y;if(x>=mod) return x-mod;return x;}
int Dif(int x,int y){x-=y;if(x < 0 ) return x+mod;return x;}
int n,C;int pri[N],cur=0;
inline int Getroot(int p){
    int x=p-1;cur=0;
    for(int i=2;i*i<=x;++i) while(x%i==0) pri[++cur]=i,x/=i;
    if(x>1) pri[++cur]=x;int g;
    for(g=2;;++g){bool fl=1;
        for(int j=1;j<=cur;++j) if(pri[j]!=pri[j+1]&&fpow(g,p/pri[j])==1) {fl=0;break;}
        if(fl) break;
    }return g;
}
int rader[N],Po[N],IP[N],A[N],B[N];
int dfs(int s,int p,int now,int blk){
    if(now==cur+1)return s+p;
    int nxt=blk/pri[now];
    return dfs(s+nxt*(p%pri[now]),(p-p%pri[now])/pri[now],now+1,nxt);
}
inline void NTT(int*A,int n,int f){
    static int tmp[N];
    for(int i=0;i<n;++i) tmp[rader[i]]=A[i];
    for(int i=0;i<n;++i) A[i]=tmp[i],tmp[i]=0;
    for(int i=1,now=cur;i<n;i*=pri[now],--now){// 模拟 FFT
        for(int t=i*pri[now],j=0;j<n;j+=t)
            for(int k=0;k<t;k+=i)
                for(int l=0;l<i;++l)
                    for(int o=0;o<pri[now];++o){
                        if(~f) Inc(tmp[j+k+l],(ll)Po[n/t*(k+l)*o%n]*A[j+i*o+l]%mod);
                        else   Inc(tmp[j+k+l],(ll)IP[n/t*(k+l)*o%n]*A[j+i*o+l]%mod);
                    }
        for(int j=0;j<n;++j) A[j]=tmp[j],tmp[j]=0;
    }
    if(f==-1) for(int i=0,inv=fpow(n,mod-2);i<n;++i) A[i]=(ll)A[i]*inv%mod;
    return;
}
int main()
{
    init(n);mod=n+1;init(C);g=Getroot(mod);
    Po[0]=IP[0]=1,Po[1]=g,IP[1]=fpow(g,mod-2);
    for(int i=0;i<n;++i) init(A[i]);
    for(int i=0;i<n;++i) init(B[i]);
    for(int i=2;i<n;++i) Po[i]=(ll)Po[i-1]*Po[1]%mod,IP[i]=(ll)IP[i-1]*IP[1]%mod;
    for(int i=0;i<n;++i) rader[i]=dfs(0,i,1,n);
    NTT(A,n,1);NTT(B,n,1);
    for(int i=0;i<n;++i) A[i]=(ll)A[i]*fpow(B[i],C)%mod;
    NTT(A,n,-1);
    for(int i=0;i<n;++i) printf("%d
",A[i]);
    return 0;
}

本题还有一种常数巨大的完全过不了的算法。(用上合并DFT的科技说不定能过)
还是快速求解任意长度 (DFT)。
有一种叫做(Bluestein’s; Algorithm) 的算法。

考虑我们要求解:
(A(w_n^k)=sum_{i=0}^{n-1}a_iw_n^{ki})
把 (ki) 以一种能够构成卷积的方式代换，为了防止出现单位根不存在二次剩余的情况这里选择用组合数替换: (ki={k+ichoose 2}-{kchoose 2}-{ichoose 2})
所以要求的就是:
(A(w_n^k)=sum_{i=0}^{n-1}a_iw_n^{{k+ichoose 2}-{kchoose 2}-{ichoose 2}})
(A(w_n^k)=w_n^{-{kchoose 2}} sum_{i=0}^{n-1}a_iw_n^{{k+ichoose 2}-{ichoose 2}})
(A(w_n^k)=w_n^{-{kchoose 2}} sum_{i=0}^{n-1}a_iw_n^{-{ichoose 2}} w_n^{{k+ichoose 2}})
后面那个东西翻转一下就是一个卷积的形式了，所以我们可以用 (FFT) 等多项式卷积算法来计算任意长度(DFT) !

但是这道题里显然不能朴素 (NTT)，因为模数并不是 (NTT) 模数。
如果用 (FFT) 代替，存在精度误差过不了，那么只能 (MTT)。
这样一算下来，我们总共需要 (DFT)三次。((IDFT) 和 (DFT)没有区别)
每次 (DFT) 里我们要用上一个 (MTT) ，而(MTT)每次要做 (7) 次 (DFT)。

也就是说我们总共做了 (3*7=21) 次 (DFT) 。
也就是说 复杂度 (O(nlog^2n))...这个东西比直接倍增算还慢吧...

(哪位大佬有更加优秀的非合并 (DFT) 的做法来教教我啊)

80'code:

#include<bits/stdc++.h>
#define Set(a,b) memset(a,b,sizeof(a))
using namespace std;
const int N=5e5+10;
const int MAXN=2097152;
int mod,g;
template <typename T> inline void init(T&x){
	x=0;char ch=getchar();bool t=0;
	for(;ch>'9'||ch<'0';ch=getchar()) if(ch=='-') t=1;
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<1)+(x<<3)+(ch-48);
	if(t) x=-x;return;
}
typedef long long ll;
template <typename T>inline void Inc(T&x,int y){x+=y;if(x>=mod) x-=mod;return;}
template <typename T>inline void Dec(T&x,int y){x-=y;if(x <  0) x+=mod;return;}
template <typename T>inline int fpow(int x,T k){int ret=1;for(;k;k>>=1,x=(ll)x*x%mod) if(k&1) ret=(ll)ret*x%mod;return ret;}
int Sum(int x,int y){x+=y;if(x>=mod) return x-mod;return x;}
int Dif(int x,int y){x-=y;if(x < 0 ) return x+mod;return x;}
int n,C;int pri[N],cur=0;
inline int Getroot(int p){
	int x=p-1;cur=0;
	for(int i=2;i*i<=x;++i) while(x%i==0) pri[++cur]=i,x/=i;
	if(x>1) pri[++cur]=x;int g;
	for(g=2;;++g){bool fl=1;
		for(int j=1;j<=cur;++j) if(pri[j]!=pri[j+1]&&fpow(g,p/pri[j])==1) {fl=0;break;}
		if(fl) break;
	}return g;
}
int rader[MAXN],Po[N],IP[N],A[N],B[N];
inline int Init(int n){
	int len=1,up=-1;for(;len<=n;len<<=1,++up);
	for(int i=0;i<len;++i) rader[i]=(rader[i>>1]>>1)|((i&1)<<up);
	return len;
}
typedef double db;
namespace MTT{
	const db PI=acos(-1);
	struct Complex{
		db x,y;Complex(db _x=0.0,db _y=0.0){x=_x,y=_y;}
		inline Complex operator +(const Complex B){return Complex(x+B.x,y+B.y);}
		inline Complex operator -(const Complex B){return Complex(x-B.x,y-B.y);}
		inline Complex operator *(const Complex B){return Complex(x*B.x-y*B.y,x*B.y+y*B.x);}
	}w[MAXN];
	inline void Calc(int n){for(int i=1;i<n;i<<=1) for(int j=0;j<i;++j) w[n/i*j]=Complex(cos(PI/i*j),sin(PI/i*j));return;}
	inline void FFT(Complex*A,int n,int f){
		for(int i=0;i<n;++i) if(rader[i]>i) swap(A[rader[i]],A[i]);
		for(int i=1;i<n;i<<=1)
			for(int j=0,p=i<<1;j<n;j+=p)
				for(int k=0;k<i;++k){
					Complex X=A[j|k],Y=A[j|k|i]* ((~f)? w[n/i*k]:Complex(w[n/i*k].x,-w[n/i*k].y));
					A[j|k]=X+Y,A[j|k|i]=X-Y;
				}
		if(!~f) for(int i=0;i<n;++i) A[i].x/=(db)n;return;
	}
	inline void Mul(int*A,int*B,int*C,int len){
		static Complex A1[MAXN],A2[MAXN],B1[MAXN],B2[MAXN];
		int MO=sqrt(mod);
		for(int i=0;i<len;++i) A1[i]=Complex(A[i]/MO,0.0),B1[i]=Complex(A[i]%MO,0.0),A2[i]=Complex(B[i]/MO,0.0),B2[i]=Complex(B[i]%MO,0.0);
		FFT(A1,len,1),FFT(A2,len,1),FFT(B1,len,1),FFT(B2,len,1);
		for(int i=0;i<len;++i) {Complex X;
			X=A1[i]*A2[i],A2[i]=A2[i]*B1[i];
			B1[i]=B1[i]*B2[i];B2[i]=B2[i]*A1[i];
			A1[i]=X,A2[i]=A2[i]+B2[i];
		}int MOD=MO*MO%mod;
		FFT(A1,len,-1),FFT(B1,len,-1),FFT(A2,len,-1);
		for(int i=0;i<len;++i) {
			int X=(ll)(A1[i].x+0.5)%mod,Y=(ll)(B1[i].x+0.5)%mod,Z=(ll)(A2[i].x+0.5)%mod;
			int ans=(ll)MOD*X%mod;Inc(ans,(ll)MO*Z%mod);Inc(ans,Y);
			C[i]=ans;
		}return;
	}
}using MTT::Calc;
inline int Co(int x){return (ll)x*(x-1)/2%n;}
inline void DFT(int*A,int n,int len,int f){
	int m=2*n-1;static int F[MAXN],G[MAXN];
	if(~f) {
		for(int i=0;i<n;++i) F[i]=(ll)A[i]*IP[Co(i)]%mod;for(int i=n;i<len;++i) F[i]=0;
		for(int i=0;i<m;++i) G[i]=Po[Co(i)];for(int i=m;i<len;++i) G[i]=0;
	}
	else   {
		for(int i=0;i<n;++i) F[i]=(ll)A[i]*Po[Co(i)]%mod;for(int i=n;i<len;++i) F[i]=0;
		for(int i=0;i<m;++i) G[i]=IP[Co(i)];for(int i=m;i<len;++i) G[i]=0;
	}reverse(F,F+n);MTT::Mul(F,G,F,len);
	for(int k=0,i=n-1;i<m;++i,++k) {
		if(~f) A[k]=(ll)F[i]*IP[Co(k)]%mod;
		else   A[k]=(ll)F[i]*Po[Co(k)]%mod;
	}
	if(!~f) for(int i=0,inv=fpow(n,mod-2);i<n;++i) A[i]=(ll)A[i]*inv%mod;
	return;
}
int main()
{
	init(n);mod=n+1;init(C);g=Getroot(mod);
	Po[0]=IP[0]=1,Po[1]=g,IP[1]=fpow(g,mod-2);
	for(int i=0;i<n;++i) init(A[i]),A[i]%=mod;
	for(int i=0;i<n;++i) init(B[i]),B[i]%=mod;
	for(int i=2;i<n;++i) Po[i]=(ll)Po[i-1]*Po[1]%mod,IP[i]=(ll)IP[i-1]*IP[1]%mod;
	int len=Init(3*n-3);Calc(len);
	DFT(A,n,len,1);DFT(B,n,len,1);
	for(int i=0;i<n;++i) A[i]=(ll)A[i]*fpow(B[i],C)%mod;
	DFT(A,n,len,-1);
	for(int i=0;i<n;++i) printf("%d
",A[i]);
	return 0;
}