zoukankan      html  css  js  c++  java
  • 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=0,f=1;char c=getchar();
        while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
        while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),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=0;i<n;i++) if (i<r[i]) swap(a[i],a[r[i]]);
        for (int i=2;i<=n;i<<=1)
        {
            complex wn=(complex){cos(2*PI/i),p*sin(2*PI/i)};
            for (int j=0;j<n;j+=i)
            {
                complex w=(complex){1,0};
                for (int k=j;k<j+(i>>1);k++,w=w*wn)
                {
                    complex x=a[k],y=w*a[k+(i>>1)];
                    a[k]=x+y,a[k+(i>>1)]=x-y;
                }
            }
        }
    }
    void mul(int n,complex *a,complex *b)
    {
        for (int i=0;i<n;i++) r[i]=(r[i>>1]>>1)|(i&1)*(n>>1);
        for (int i=0;i<n;i++) a[i].y=a[i].x-b[i].x,a[i].x=a[i].x+b[i].x;
        DFT(n,a,1);
        for (int i=0;i<n;i++) a[i]=a[i]*a[i];
        DFT(n,a,-1);
        for (int i=0;i<n;i++) a[i].x=a[i].x/n/4;
    }
    void solve(int l,int r)
    {
        if (l==r) {f[0]+=1ll*a[l]*b[l];return;}
        int mid=l+r>>1;
        solve(l,mid);
        solve(mid+1,r);
        int t=1;while (t<r-l+1) t<<=1;
        for (int i=0;i<t;i++) c[i].x=c[i].y=d[i].x=d[i].y=0;
        for (int i=l;i<=mid;i++) c[i-l].x=a[i];
        for (int i=mid+1;i<=r;i++) d[i-mid-1].x=b[i];
        mul(t,c,d);
        for (int i=l+mid+1;i<=mid+r;i++) f[i]+=(long long)(c[i-l-mid-1].x+0.5);
        for (int i=0;i<t;i++) c[i].x=c[i].y=d[i].x=d[i].y=0;
        for (int i=mid+1;i<=r;i++) c[i-mid-1].x=a[i];
        for (int i=l;i<=mid;i++) d[mid-i].x=b[i];
        mul(t,c,d);
        for (int i=1;i<=r-l;i++) f[i]+=(long long)(c[i-1].x+0.5);
    }
    int main()
    {
    #ifndef ONLINE_JUDGE
        freopen("bzoj4836.in","r",stdin);
        freopen("bzoj4836.out","w",stdout);
        const char LL[]="%I64d
    ";
    #else
        const char LL[]="%lld
    ";
    #endif
        T=read();
        while (T--)
        {
            n=read(),m=read(),q=read();
            memset(f,0,sizeof(f));
            memset(a,0,sizeof(a));memset(b,0,sizeof(b));
            int x=0,y;
            for (int i=1;i<=n;i++) x=max(y=read(),x),a[y]++;
            for (int i=1;i<=m;i++) x=max(y=read(),x),b[y]++;
            solve(0,x);
            while (q--) printf(LL,f[read()]);
        }
        return 0;
    }
  • 相关阅读:
    遍历当前窗口名字
    Delphi获取其它进程窗口句柄的3种方法
    Delphi判断一个字符是否为汉字的最佳方法
    Delphi拷贝目录(含子目录)的方法
    贴一份用delphi修改注册表改网卡MAC地址的代码
    delphi的TFileStream 内存流
    Delphi的ListView自动排序
    Delphi中上指定进程(进程名)
    Delphi-IP地址的隐藏
    C# 简单线程实例
  • 原文地址:https://www.cnblogs.com/Gloid/p/9445977.html
Copyright © 2011-2022 走看看