把ABCD看作数轴上的四个点, 满足\(1\le A<B<C<D\le N\). 并记有cnt[i]个点的坐标为i.
枚举CD的距离lenCD, 则应该满足\(1\le lenCD\)且\(lenCD \times 9 < N\).

假设已知D的坐标\(x_D\), 则应该有\(lenCD \times 9 + 1 < x_D \le N\). 那么C的坐标\(x_C = x_D - lenCD\). 只需要知道有多少对满足条件的A和B, 不妨记为s, 则以\(x_D\)为D的方案数 (在lenCD固定的情况下) 为\(s \times cnt[x_C]\). 同理以\(x_C\)为C的方案数为\(s \times cnt[x_D]\).

若我们从小到大依次枚举\(x_D\), 则对于当前的\(x_D\)来说s只比\(x_D\)时的s多了\(cnt[x_D-9\times lenCD-1] \times cnt[x_D-7\times lenCD-1]\), 这恰好对应了唯一多出来的一部分可能的二元点对A和B. 至此可以在\(O(C^2/9)\)内求出来所有值\(1\le x\le N\)作为C点坐标和D点坐标的方案数.

相似的思想, 在枚举\(lenCD\)后去枚举A的坐标, 即可找出所有点作为A和B的方案总数.

我们可以把这4个数看作是个数轴上的点，根据题目给的条件，可知AB=2*CD，BC＞6*CD，AD>9*CD

那么，我们只需要确定D，就可以确定C点，然后再找AB。我们也可以通过找C来确定ABD。
100分代码

#include <cstdio>
int n,m,x[40001],vis[15001],a[15001],b[15001],c[15001],d[15001];
int main() {
    //freopen("magic.in","r",stdin);
    //freopen("magic.out","w",stdout);
    scanf("%d%d",&n,&m);
    for (int i = 1;i <= m;i++) {
        scanf("%d",&x[i]);
        vis[x[i]]++;
    }
    for (int i = 1;i <= n/9;i++) {
        int sum = 0;
        for (int j = i*9+2;j <= n;j++) {
            sum += vis[j-(9*i+1)]*vis[j-(9*i+1)+2*i];
            d[j] += sum*vis[j-i];
            c[j-i]+=sum*vis[j];
        }
        sum = 0;
        for (int j = n-(9*i+1);j >= 1;j--) {
            sum += vis[j+(9*i+1)]*vis[j+(9*i+1)-i];
            a[j] += sum*vis[j+2*i];
            b[j+2*i] += sum*vis[j];

        }
    }
    for (int i = 1;i <= m;i++) printf("%d %d %d %d\n",a[x[i]],b[x[i]],c[x[i]],d[x[i]]);
    return 0;
}

首先可以发现每个x都小于n，而n最大值只是15000，所以可以开一个桶来存每个魔法值出现的次数

回忆一下3个约束条件

xa<xb<xc<xd ①

xb−xa=2(xd−xc) ②

xb−xa<(xc−xb)/3 ③

现在魔改一下这三个式子

设t=xd−xc
所以②可化为xb−xa=2t ④

将④代入③

2t<(xc−xb)/3
移项一下，就变成

6t<xc−xb ⑤

再魔改一下

设6t+k=xc−xb（就是把差的部分补上去）

于是可以画出来一个图
https://images2018.cnblogs.com/blog/1113423/201808/1113423-20180801210708340-1039644464.png

显然，A的最小值为1，D的最大值为n
由图可得AD=9t+k
所以我们可以尝试着枚举t，用t来表示各个魔法值的值

由上易得t的范围为1<=t<=(n−1)/9
在代码中为了避免除法写成t∗9<n
再枚举D，因为我们已经枚举出了t，所以C的值是可以直接算出来的

C=D−t
又因为使A,B,C,D满足条件的k的最小值为1，所以对于当前的C和D，最大的A和B为A=D−9t−1,B=D−7t−1
那么如果A和B更小怎么办？

观察到在其他条件不变的情况下，只要C和B满足Xc−Xb>6t，那么这个魔法阵就一定成立，所以当(a1<a2,b1<b2)时，只要a2和b2能够和C,D组成魔法阵，a1,b1也一定能和C,D组成魔法阵，所以可以使用前缀和优化

然后又由乘法原理可得，当前魔法值作为D物品的个数为SumD=SumA∗SumB∗SumC
所以我们利用前缀和优化SumA∗SumB
C的情况可以顺便在算D的时候算出来

那么还有一个问题是，我们枚举的D的范围是多少？

因为要统计前缀和，所以一定是要顺推下去的，由上面那张图我们可以知道，D的最大值为n，最小值则为当k=1且A=1的时候，所以D的最小值为9∗t+2，再小是无法组成魔法阵的

同理可以枚举A

但是这个的情况又和枚举D的情况有一点不同

在其他条件不变的情况下，只要C和B满足Xc−Xb>6t，那么这个魔法阵就一定成立，所以当(c1<c2,d1<d2)时，只要c1和d1能够和A,B组成魔法阵，c2,d2也一定能和A,B组成魔法阵，所以可以使用后缀和优化

因为需要统计后缀和，所以需要逆推

枚举的范围：A的最大值为(n−t∗9−1)(因为当k=1，D=n的时候A才最大)，A的最小值则为1

所以就可以算出每个魔法值作为A,B,C,D物品的次数了，输出时直接输出当前魔法物品的魔法值的次数就可以了

#include<cstdio>

#include<iostream>

#include<cmath>

#include<cstring>

#include<string>

#include<algorithm>

#include<functional>

#include<queue>

#include<stack>

#include<set>

#include<vector>

#include<map>

using namespace std;

int a[15005],b[15005],c[15005],d[15005],w[15005],h[40005],n,m,x,y;

//abcd表示某个点作为abcd物品出现的次数；w表示数轴上每个点出现的次数；h表示每个物品的魔法值

//n表示最大魔法值，m表示物品数量

int main()

{

cin>>n>>m;

if(n<11)
{
for(int i=1;i<=m;i++) printf("0 0 0 0\n");
return 0;
}
for(int i=1;i<=m;i++)
{
scanf("%d",&h[i]);
w[h[i]]++;
}//把这些点标记在数轴上
for(int i=1;i<=n/9;i++)
{
//若数轴上有一个魔法阵：ABCD，其中有AB=2*CD，BC>6*CD

//所以只需枚举CD的长度就可以了

x=1+9*i,y=0;//x为AD最短长度

for(int j=2+9*i;j<=n;j++)

{

//因为数轴是从1开始的，所以从1+x开始枚举

//枚举D点即j，则C点为j-i，A点为j-x，B点为j-x+2*i

//CD的个数取决于AB有多少组，所以我们用y表示AB的组数

y+=w[j-x]*w[j-x+i+i];//y为AB的对数

//D点是不定的。但是D点变化时，之前合格的AB两点仍然合格，所以要累加

d[j]+=y*w[j-i];//有几组AB，就有几个C点，就有几个D点

c[j-i]+=y*w[j];//有几组AB，就有几个D点，就有几个C点

}

//注意，魔法值可能重复，所以在加的时候，注意不要直接加。

//同理，枚举CD两点，确定AB的个数

x=8*i+1,y=0;

for(int j=n-9*i-1;j>=1;j--)

{

y+=w[j+x]*w[j+x+i];

a[j]+=y*w[j+i+i];

b[j+i+i]+=y*w[j];

}

}

for(int i=1;i<=m;i++)//输出每个物品对应的魔法值的作为abcd物品的次数

cout<<a[h[i]]<<' '<<b[h[i]]<<' '<<c[h[i]]<<' '<<d[h[i]]<<endl;

}

60分代码，应该是最好理解的60分了。

#include<stdio.h>
#include<stdlib.h>
#include<iostream>
#include<string.h>
#include<algorithm>
using namespace std;
struct thing
{
    int index;
    int value;
    int A=0;
    int B=0;
    int C=0;
    int D=0;
}t[40000];
bool cmp_value(thing x,thing y)
{
    return x.value<y.value;
}
bool cmp_index(thing x,thing y)
{
    return x.index<y.index;
}
int main()
{
    int n,m;
    double temp;
    cin>>n>>m;
    for(int i=0;i<m;i++)
    {
        cin>>t[i].value;
        t[i].index=i;
    }
    sort(t,t+m,cmp_value);
    for(int a=0;a<m;a++)
    {
        for(int b=a+1;b<m;b++)
        {
            for(int c=b+1;c<m;c++)
            {
                for(int d=c+1;d<m;d++)
                {
                    temp=(t[c].value-t[b].value)/3.0;
                    if(t[a].value<t[b].value  && t[b].value<t[c].value && t[c].value<t[d].value && t[b].value-t[a].value==2*(t[d].value-t[c].value) && t[b].value-t[a].value<temp)
                    {
                        t[a].A++;
                        t[b].B++;
                        t[c].C++;
                        t[d].D++;
                    }
                }
            }
        }
    }
    sort(t,t+m,cmp_index);
    for(int i=0;i<m;i++)
    {
        cout<<t[i].A<<' '<<t[i].B<<' '<<t[i].C<<' '<<t[i].D<<endl;
    }
    return 0;
}