P13553 题解

首先要统计 $(i,j,k)$ 的数量。

枚举 $i$。如果 $a_i=j-i$ 或者 $a_i=k-i$ 是简单的——可以确定三元组中第二个数，那么第三个数只有 $O(1)$ 种可能，逐个判断即可。

否则，$a_i=k-j$。

令 $b_i=i-a_i$。分类讨论 $a_j$ 与 $a_k$ 的大小关系，可以得到：$b_j=b_k=i$ 或者 $b_j=i-a_i,b_k=i+a_i$。

用 bitset 直接维护每个 $b$ 的出现位置，求交即可，时空复杂度均为 $O(\frac{n^2}{w})$。

但是空间会爆，需要优化。

具体地，把 $b_i$ 按照大小分类，$>B$ 的用 bitset 维护，$\le B$ 的直接用 vector 维护，暴力枚举。

取 $B=\frac{n}{w}$，时间复杂度不变，空间复杂度降为 $O(n)$。为了减小常数，$B$ 可以开小一些。

然后就到了第二问，构造方案。

$n=20$ 和 $n=500$ 没啥好说的，直接用模拟退火之类的技巧搜索即可——实际上这个答案的界非常松，每次只保留更优解（卡在局部最优解的时候手动回退一下版本）就可以搜出来。

首先我们注意到，满足 $a_i$ 和 $j-i,k-i$ 匹配的东西是 $O(n)$ 的，满足 $b_j=b_k=i$ 的三元组也是 $O(n)$ 的（每个 $b$ 只会被算到一次）。

而题目要求我们构造的 $k$ 比它们能凑出来的三元组数量大一个数量级——考虑到所有构造方案都必定带一个小常数，这个量级接近 $O(n\sqrt n)$。

现在有 $b_j=i-a_i,b_k=i+a_i$。首先选一些根据 $b$ 划分的等价类，往每个等价类里面塞一大堆东西；然后对于所有的 $b_i,b_j(2\mid i+j,i<j)$，令 $a_\frac{i+j}{2}=\frac{j-i}{2}$（如果出现冲突，则任选一个值）。

这里每个 $b$ 平均要被覆盖根号次，考虑设置一个阈值 $B=O(\sqrt n)$，选择 $2x$ 与 $2Bx(0\le x<B)$ 作为 $b$ 的集合，这样每个 $b$ 都被覆盖了 $B+O(1)$ 次。

当然还有 $k-j=a_i$ 的限制。前面对 $a_i$ 进行了特殊的赋值，因此这个条件可以转化为 $j+a_j=k+a_k$。

现在考虑如下策略：每次选择一个起点 $p$，对于所有可能的 $b$，尝试在 $p+b$ 的位置放置 $p-b$（如果出现冲突或者越界就不放置），直到把所有位置放满。

把位置 $p$ 映射到 $([\frac{p}{B}],p\bmod B)$ 上，发现放置一组点相当于覆盖一个 L 形的区域。

一个简单的策略是，将平面划分成若干个 $B\times B$ 的块，每块再次划分成 $B$ 个 L 形区域——此时这一块的贡献约为 $\frac{B^3}{3}$。

粗略估计划分成的块数为 $\frac{n}{B^2}-1$，可以枚举 $B$ 求出一个优秀的块长——然后手动精细调整。

满分需要继续优化：在预处理部分，对于 $b_i,b_j(2\mid i+j,i<j)$，我们令 $a_\frac{i+j}{2}=\frac{j-i}{2}$——这一步不用做完，可以有一个后缀的 $a$ 没有被预处理。它们并不会被统计到太多次，因此没有太大的损失。

考虑这 $B$ 个 L 形区域，发现有 $B-1$ 个仍然有扩展空间——只要把划分成的块从 $B\times B$ 改为 $B\times (B+1)$（或者更大）即可。

具体的实现可以参考代码。

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const ll N=200007,B=500;
ll n,m,q,a[N],ans,id[N*3],cnt[N],bid[N],sum[N];
vector<int> vec[N];
bitset<N> bs[B];
mt19937 rng(20120712);
void mysort(ll& x,ll& y,ll& z){
	if (x>y) swap(x,y);
	if (y>z) swap(y,z);
	if (x>y) swap(x,y);
}
bool OK(ll x,ll y,ll z){
	if (z>n||z<=0) return 0;
	ll A[3]={abs(x-y),abs(y-z),abs(x-z)},B[3]={a[x],a[y],a[z]};
	mysort(A[0],A[1],A[2]);mysort(B[0],B[1],B[2]);
	return A[0]==B[0]&&A[1]==B[1]&&A[2]==B[2];
}
ll query(ll u,ll v,ll k){
	static bool ok[N];
	ll ans=0;
	if (cnt[u]<B&&cnt[v]<B){
		for (auto x:vec[v]) ok[x]=1;
		for (auto x:vec[u]) if (x+k<=n) ans+=ok[x+k];
		for (auto x:vec[v]) ok[x]=0;
		return ans;
	}
	if (cnt[u]<B){
		v=bid[v];
		for (auto x:vec[u]) if (x+k<=n) ans+=bs[v][x+k];
		return ans;
	}
	if (cnt[v]<B){
		u=bid[u];
		for (auto x:vec[v]) if (x-k>=0) ans+=bs[u][x-k];
		return ans;
	}
	return (bs[bid[u]]&(bs[bid[v]]>>k)).count();
}
ll count_triples(vector<int> H){
	n=H.size();
	for (int i=1;i<=n;++i) a[i]=H[i-1];
	
//	for (int i=1;i<=n;++i) for (int j=i+1;j<=n;++j) for (int k=j+1;k<=n;++k) if (OK(i,j,k)){
//		++sum[i];++sum[j];++sum[k];
//		++ans;
//	}
	
	for (int i=n;i;--i){
		//count
		int x=id[i+N],y,k=a[i];
		if (cnt[x]>=2) ans+=query(x,x,k);
		x=id[i-k+N];y=id[i+k+N];
		if (cnt[x]&&cnt[y]) ans+=query(x,y,k);
		//insert
		k=i-a[i]+N;
		if (!id[k]) id[k]=++m;
		k=id[k];
		vec[k].emplace_back(i);
		if ((++cnt[k])==B){
			bid[k]=++q;
			for (auto x:vec[k]) bs[q][x]=1;
		}
		if (cnt[k]>B) bs[bid[k]][i]=1;
//		cout<<ans<<endl;
	}
	for (int i=1;i<=n;++i){
		int j=i+a[i];
		if (j>n) continue;
		if (a[j]<a[i]) ans+=OK(i,j,j-a[j]);
		if (a[j]!=a[i]) ans+=OK(i,j,j+a[j]);
		if (a[j]!=2*a[i]&&a[j]*2!=a[i]) ans+=OK(i,j,i+a[j]);
	}
	return ans;
}
vector<int> construct_range(int M,int K){
	n=M;
	vector<int> vec,L,R;
	int B;
	if (n==20){
		vector<int> vec={3,1,2,1,4,3,6,5,6,7,2,3,4,1,2,1,8,1,6,5};
		return vec;
	}
	if (n==500){
		vector<int> vec={1,12,1,2,1,2,1,4,3,2,11,10,9,8,7,2,5,4,1,11,1,22,21,20,19,18,17,16,15,14,13,12,11,32,31,30,29,28,27,26,25,24,23,22,43,42,41,40,39,38,37,36,35,34,11,54,53,52,51,50,49,48,47,46,45,22,65,64,63,62,61,60,59,58,57,56,55,76,111,74,73,72,71,70,69,68,67,22,87,86,85,84,83,82,81,80,79,78,11,98,97,96,95,94,93,92,91,90,89,22,109,108,107,106,105,104,103,102,101,100,11,118,121,116,119,114,125,70,127,110,125,44,143,122,141,178,139,138,137,136,135,134,33,132,131,130,129,128,127,150,149,148,147,22,145,144,143,142,141,116,139,162,161,160,33,158,157,156,155,154,105,128,151,174,173,172,171,170,169,168,167,94,117,140,163,186,11,184,183,182,181,180,83,106,129,152,175,22,197,196,195,194,193,72,95,118,141,164,11,210,209,208,207,206,61,84,107,130,153,1,199,222,221,211,219,50,73,96,119,142,165,188,211,234,233,232,39,62,85,108,131,154,177,200,223,246,245,28,51,74,97,120,143,166,189,212,235,18,259,258,257,256,255,254,253,252,251,250,249,248,271,270,269,268,267,11,265,264,263,262,237,260,283,282,281,280,279,12,277,276,275,226,249,272,295,294,293,292,291,11,289,288,215,238,261,284,307,306,305,304,303,12,301,204,227,250,273,296,319,318,317,316,315,24,193,216,239,262,285,308,331,330,329,328,327,12,205,228,251,274,297,320,343,342,341,340,171,24,217,240,263,286,309,332,355,354,353,160,183,12,229,252,275,298,321,344,367,366,149,172,195,48,241,264,287,310,333,356,139,380,379,378,377,376,375,374,373,372,371,370,369,392,391,390,389,388,121,386,385,384,383,358,381,404,403,402,401,400,133,398,397,396,347,370,393,416,415,414,413,412,121,410,409,336,359,382,405,428,427,426,425,424,109,422,325,348,371,394,417,440,439,438,437,436,97,314,337,360,383,406,429,452,451,450,449,448,133,326,349,372,395,418,441,464,463,462,461,292,121,338,361,384,407,430,453,476,475,474,281,304,109,350,373,396,419,442,465,488,487,270,293,316,145,362,385,408,431,454,477,260};
		return vec;
	}
	if (n==5000) B=41;
	if (n==30000) B=100;
	if (n==100000) B=150;
	if (n==200000) B=258;
//	cin>>B;
	for (int i=0;i<B;++i) L.emplace_back(i);
	for (int i=1;i<B;++i) R.emplace_back(i*B);
	for (auto x:L) for (auto y:L) if (x!=y){
		if (!a[x+y]) a[x+y]=abs(y-x);
		else a[x+y]=min(a[x+y],(ll)abs(y-x));
	}
//	for (auto x:R) for (auto y:R) a[x+y]=y-x;
	for (auto x:L) for (auto y:R) a[x+y]=y-x;
//	if (n==5000) for (int i=2;i<B;i+=2) a[i]=B+1;
	int d=B*B+2800,lim=R.back()-3400;
	if (n==5000) d=B*B+3,lim=R.back()-1;
	if (n==30000) d=B*B+B,lim=R.back()-B-1;
	for (int i=lim+1;i<=n;++i) a[i]=0;
	vector<int> owo;
	for (int i=lim+1;i<=n;i+=d) owo.emplace_back(i);
	reverse(owo.begin(),owo.end());
	for (auto i:owo){
		for (int j=0;j<B;++j){
			int x=i+j*(B+1);
			for (auto d:L) if (x+d>0&&x+d<=n&&a[x+d]==0) a[x+d]=x-d;
			for (auto d:R) if (x+d>0&&x+d<=n&&a[x+d]==0) a[x+d]=x-d;
		}
	}
	for (int i=1;i<=n;++i) vec.emplace_back(min(max(a[i],1ll),n-1));
	return vec;
}