BZOJ4836 二元运算（分治FFT）

时间：2023-03-09 14:27:33
BZOJ4836 二元运算（分治FFT）

　　设A(n)为a中n的个数，B(n)为b中n的个数。如果只考虑加法显然是一个卷积，减法翻转一下也显然是一个卷积。

　　问题在于两者都有。容易想到分开处理。那么可以考虑分治。即对于值域区间[l,r]，分别计算A[l,mid]和B[mid+1,r]的贡献及A[mid+1,r]和B[l,mid]的贡献，然后再递归处理[l,mid]和[mid+1,r]。一定程度上类似于cdq分治。

　　注意结果可能爆int，用NTT的话不太方便。

#include<iostream>

#include<cstdio>

#include<cmath>

#include<cstdlib>

#include<cstring>

#include<algorithm>

using namespace std;

int read()

{

    int x=,f=;char c=getchar();

    while (c<''||c>'') {if (c=='-') f=-;c=getchar();}

    while (c>=''&&c<='') x=(x<<)+(x<<)+(c^),c=getchar();

    return x*f;

}

#define N 270000

const double PI=3.14159265358979324;

struct complex

{

    double x,y;

    complex operator +(const complex&a) const

    {

        return (complex){x+a.x,y+a.y};

    }

    complex operator -(const complex&a) const

    {

        return (complex){x-a.x,y-a.y};

    }

    complex operator *(const complex&a) const

    {

        return (complex){x*a.x-y*a.y,x*a.y+y*a.x};

    }

}c[N],d[N];

int T,n,m,q,a[N],b[N],r[N];

long long f[N];

void DFT(int n,complex *a,int p)

{

    for (int i=;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);

    for (int i=;i<=n;i<<=)

    {

        complex wn=(complex){cos(*PI/i),p*sin(*PI/i)};

        for (int j=;j<n;j+=i)

        {

            complex w=(complex){,};

            for (int k=j;k<j+(i>>);k++,w=w*wn)

            {

                complex x=a[k],y=w*a[k+(i>>)];

                a[k]=x+y,a[k+(i>>)]=x-y;

            }

        }

    }

}

void mul(int n,complex *a,complex *b)

{

    for (int i=;i<n;i++) r[i]=(r[i>>]>>)|(i&)*(n>>);

    for (int i=;i<n;i++) a[i].y=a[i].x-b[i].x,a[i].x=a[i].x+b[i].x;

    DFT(n,a,);

    for (int i=;i<n;i++) a[i]=a[i]*a[i];

    DFT(n,a,-);

    for (int i=;i<n;i++) a[i].x=a[i].x/n/;

}

void solve(int l,int r)

{

    if (l==r) {f[]+=1ll*a[l]*b[l];return;}

    int mid=l+r>>;

    solve(l,mid);

    solve(mid+,r);

    int t=;while (t<r-l+) t<<=;

    for (int i=;i<t;i++) c[i].x=c[i].y=d[i].x=d[i].y=;

    for (int i=l;i<=mid;i++) c[i-l].x=a[i];

    for (int i=mid+;i<=r;i++) d[i-mid-].x=b[i];

    mul(t,c,d);

    for (int i=l+mid+;i<=mid+r;i++) f[i]+=(long long)(c[i-l-mid-].x+0.5);

    for (int i=;i<t;i++) c[i].x=c[i].y=d[i].x=d[i].y=;

    for (int i=mid+;i<=r;i++) c[i-mid-].x=a[i];

    for (int i=l;i<=mid;i++) d[mid-i].x=b[i];

    mul(t,c,d);

    for (int i=;i<=r-l;i++) f[i]+=(long long)(c[i-].x+0.5);

}

int main()

{

#ifndef ONLINE_JUDGE

    freopen("bzoj4836.in","r",stdin);

    freopen("bzoj4836.out","w",stdout);

    const char LL[]="%I64d\n";

#else

    const char LL[]="%lld\n";

#endif

    T=read();

    while (T--)

    {

        n=read(),m=read(),q=read();

        memset(f,,sizeof(f));

        memset(a,,sizeof(a));memset(b,,sizeof(b));

        int x=,y;

        for (int i=;i<=n;i++) x=max(y=read(),x),a[y]++;

        for (int i=;i<=m;i++) x=max(y=read(),x),b[y]++;

        solve(,x);

        while (q--) printf(LL,f[read()]);

    }

    return ;

}

相关文章