当前位置: 代码迷 >> 综合 >> 【NTT】【二维卷积】CodeChef BuyingLand
  详细解决方案

【NTT】【二维卷积】CodeChef BuyingLand

热度:61   发布时间:2023-09-27 04:07:20.0

分析:

非常脑残的二维卷积题

把式子平方暴力拆开,然后发现得到的6个式子居然可以算。
3个直接O(1)求得。
2个需要二维前缀和。
剩下一个用二维卷积。

#include<cstdio>
#include<cstring>
#include<algorithm>
#define SF scanf
#define PF printf
#define MAXN 1410
#define MAXM 4000010
#define MOD 998244353
using namespace std;
typedef long long ll;
int w,h,n,m;
ll a[MAXN][MAXN],g[MAXN][MAXN];
ll pre[MAXN][MAXN];
ll G[MAXM],F[MAXM],res[MAXM];
ll ans[MAXN][MAXN];
pair<ll,int> pri[MAXM];
ll fsp(ll x,int y){
    ll res=1;while(y){
    if(y&1)res=res*x%MOD;x=x*x%MOD;y>>=1;}return res;
}
void NTT(ll A[],int N,int flag){
    for(int i=1,j=0;i<N;i++){
    for(int d=N;j^=d>>=1,~j&d;);if(i<j)swap(A[i],A[j]);}for(int i=1;i<N;i<<=1){
    ll wn=fsp(3,(MOD-1)/(i<<1));if(flag)wn=fsp(wn,MOD-2);for(int j=0;j<N;j+=(i<<1)){
    ll w=1;for(int k=0;k<i;k++,w=w*wn%MOD){
    ll x=A[j+k],y=A[i+j+k]*w%MOD;A[j+k]=(x+y)%MOD;A[i+j+k]=(x-y+MOD)%MOD;}}}if(flag){
    ll invN=fsp(N,MOD-2);for(int i=0;i<N;i++)A[i]=A[i]*invN%MOD;}
}
void Mul(ll A[],ll B[],ll res[],int len){
    int p=1;while(p<=len)p<<=1;NTT(A,p,0);NTT(B,p,0);for(int i=0;i<p;i++)res[i]=A[i]*B[i]%MOD;NTT(res,p,1);for(int i=0;i<len;i++)if(res[i]>100000000ll)res[i]-=MOD;
}
int main(){
    SF("%d%d",&n,&m);for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)SF("%lld",&a[i][j]);SF("%d%d",&w,&h);for(int i=1;i<=w;i++)for(int j=1;j<=h;j++)SF("%lld",&g[i][j]);int posx,posy;SF("%d%d",&posx,&posy);ll tot2=0,tot=0;for(int i=1;i<=w;i++)for(int j=1;j<=h;j++)tot2+=g[i][j]*g[i][j],tot+=g[i][j];for(int i=1;i<=n-w+1;i++)for(int j=1;j<=m-h+1;j++)ans[i][j]=a[i+posx-1][j+posy-1]*a[i+posx-1][j+posy-1]*w*h+tot2+2ll*a[i+posx-1][j+posy-1]*tot;for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)pre[i][j]=a[i][j]*a[i][j]+pre[i-1][j]+pre[i][j-1]-pre[i-1][j-1];for(int i=1;i<=n-w+1;i++)for(int j=1;j<=m-h+1;j++)ans[i][j]+=(pre[i+w-1][j+h-1]-pre[i-1][j+h-1]-pre[i+w-1][j-1]+pre[i-1][j-1]);for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)pre[i][j]=a[i][j]+pre[i-1][j]+pre[i][j-1]-pre[i-1][j-1];for(int i=1;i<=n-w+1;i++)for(int j=1;j<=m-h+1;j++)ans[i][j]-=2ll*a[i+posx-1][j+posy-1]*(pre[i+w-1][j+h-1]-pre[i-1][j+h-1]-pre[i+w-1][j-1]+pre[i-1][j-1]);for(int i=1;i<=w;i++)reverse(g[i]+1,g[i]+1+h);for(int i=1;i<w-i+1;i++)swap(g[i],g[w-i+1]);for(int i=1;i<=n;i++)for(int j=1;j<=m;j++)G[(i-1)*m+j-1]=g[i][j],F[(i-1)*m+j-1]=a[i][j];Mul(G,F,res,n*m*2);int cnt=0;
// for(int i=1;i<=n;i++)
// for(int j=1;j<=m;j++)
// PF("[%d %d %lld]\n",i,j,res[(i-1)*m+j-1]);for(int i=1;i<=n-w+1;i++)for(int j=1;j<=m-h+1;j++)pri[cnt++]=make_pair(ans[i][j]-res[(i-2+w)*m+j-2+h]*2ll,(i-1)*m+j-1);sort(pri,pri+cnt);int k;SF("%d",&k);for(int i=0;i<k;i++){
    int idx=pri[i].second/m+1;int idy=pri[i].second%m+1;PF("%d %d %lld\n",idx,idy,pri[i].first);}
}