题解:P10666 COUNTARI

· · 题解

题目传送门

神秘做法。

题目分析

c_i 表示 i 的出现次数。

首先我们先不算三个数相等的情况,最后再加上去。这样方便后面处理。

先不考虑位置顺序的限制,发现是好求的。您都开到这道题了那一定会的对吧!。记 s 表示 c 自己卷起来的结果,那这一部分的答案就是 \sum c_i(s_i-c_i)。但发现其实每个都算了两次,不过问题不大。

考虑减掉多算的部分,发现只有两种情况:x-d,x+d,xx,x-d,x+d,其中 d \neq 0。第二种其实和第一种是类似的,只要翻转一下重新求即可。下面考虑求第一种情况。

记现在处理到 p,之前每个数的出现个数为 f_i,那么我们有两种求法:

发现上面两种方法都是很慢的。将两种方法结合一下,类似于分块,我们可以快速将前几个块末求出的 h 与这块内的散点合并。具体地,记前面几个整块末求出的个数为 f,自己卷积的结果为 h,这块内每个数的出现次数为 g,则答案为 h_{2a_p}-f_{a_p}^2+\sum_{i>R} 2f_{2a_p-a_i}+g_{2a_p-a_i}。(乘 2 是因为和其它的统一一下也算两次)。

设块长为 B。则卷积部分复杂度为 O(\frac{n}{B}n \log n),枚举计算部分为 O(nB)。于是 B\sqrt{n\log n} 时最小,但由于卷积常数很大,可以再取大些。

于是就做完啦。

Code

// by wangyizhi(571247)
#include<bits/stdc++.h>
//#include<bits/extc++.h>
using namespace std;
using ll=long long;
using ld=long double;
//#define int ll
using pii=pair<int,int>;
//bool Mst;
namespace wyzmath
{
    /*complex*/
    template<typename _Tp=double>
    class complex
    {
    public:
        typedef _Tp vtype;
        vtype x,y;
        complex(vtype _x=0,vtype _y=0):x(_x),y(_y){}
        vtype length(){return std::sqrt(x*x+y*y);}
    }; // class complex
    /*about complex*/
    template<typename _Tp=double>
    complex<_Tp> operator+(complex<_Tp> a,complex<_Tp> b)
    {return complex<_Tp>(a.x+b.x,a.y+b.y);}
    template<typename _Tp=double>
    complex<_Tp> operator-(complex<_Tp> a,complex<_Tp> b)
    {return complex<_Tp>(a.x-b.x,a.y-b.y);}
    template<typename _Tp=double>
    complex<_Tp> operator*(complex<_Tp> a,complex<_Tp> b)
    {return complex<_Tp>(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
    template<typename _Tp=double>
    complex<_Tp> operator+=(complex<_Tp> &a,complex<_Tp> b){return a=a+b;}
    template<typename _Tp=double>
    complex<_Tp> operator-=(complex<_Tp> &a,complex<_Tp> b){return a=a-b;}
    template<typename _Tp=double>
    complex<_Tp> operator*=(complex<_Tp> &a,complex<_Tp> b){return a=a*b;}
}
using namespace std;
using ld=long double;
using compld=wyzmath::complex<ld>;
const compld I(1,0);
const int N=1<<21,V=3e4,VV=6e4,B=5000;
const ld PI=3.141592653589793;
int rev[N];
ll a[N],c[N],s[N],d[N],e[N],p[N];
void FFT(compld *a,int n,int op)
{
    for(int i=0;i<n;i++) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
    for(int mid=1;mid<n;mid<<=1)
    {
        compld omega(cos(PI/mid),op*sin(PI/mid));
        for(int len=mid<<1,i=0;i<n;i+=len)
        {
            compld o(1,0),x,y;
            for(int j=0;j<mid;j++,o=o*omega)
                x=a[i+j],y=o*a[i+j+mid],a[i+j]=x+y,a[i+j+mid]=x-y;
        }
    }
}
void conv(ll *a,ll *b,ll *c,int n)
{
    static compld tmpa[N],tmpb[N];
    int sz=1,l=0;
    while(sz<n) sz<<=1,l++;
    for(int i=1;i<sz;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    for(int i=0;i<sz;i++) tmpa[i]=a[i],tmpb[i]=b[i];
    FFT(tmpa,sz,1),FFT(tmpb,sz,1);
    for(int i=0;i<sz;i++) tmpa[i]=tmpa[i]*tmpb[i];
    FFT(tmpa,sz,-1);
    for(int i=0;i<sz;i++) c[i]=(ll)(tmpa[i].x/sz+0.5);
    for(int i=n;i<sz;i++) c[i]=0;
}
//bool Med;
signed main()
{
//  cerr<<"Memory Size: "<<abs((&Med)-(&Mst))/1024.0/1024<<" MB\n";
    ios::sync_with_stdio(0);cin.tie(0),cout.tie(0);
//  freopen("in.in","r",stdin);
//  freopen("out.out","w",stdout);
    int n;ll ans=0;
    cin>>n;
    for(int i=1;i<=n;i++) cin>>a[i],a[i]--,c[a[i]]++;
    conv(c,c,s,VV);
    for(int i=0;i<V;i++) s[i<<1]-=c[i]*c[i];
    for(int i=0;i<V;i++) ans+=s[i<<1]*c[i];
//  cout<<ans<<"\n";
    auto sol=[&]()->void
    {
        for(int i=0;i<VV;i++) s[i]=0;
        for(int i=0;i<V;i++) d[i]=e[i]=0;
        for(int i=1,le=0;i<=n;i++)
        {
            ans-=s[a[i]*2]-d[a[i]]*d[a[i]];
            for(int j=le+1;j<=i;j++) if(a[i]!=a[j]&&a[i]*2>=a[j])
                ans-=d[a[i]*2-a[j]]*2+e[a[i]*2-a[j]];
            e[a[i]]++;
            if(i%B==0&&i!=n)
            {
                for(int j=0;j<V;j++) d[j]+=e[j],e[j]=0;
                conv(d,d,s,VV);
                le=i;
            }
        }
    };
    sol();reverse(a+1,a+n+1);sol();ans/=2;
    for(int i=0;i<V;i++) ans+=c[i]*(c[i]-1)*(c[i]-2)/6;
    cout<<ans<<"\n";
    return 0;
}