首先变换为 $p_i/p_j=p_l/p_k$，设 $q_{p_i}=i$，也就是数 $i$ 的下标，考虑计算数对 $(a,b,c,d)$ 满嘴 $q_{ac} < q_{bc} < q_{bd} < q_{ad}$ 并且 $\gcd (a,b) =1$ 的个数。
考虑把 $\gcd(a,b)=1$ 去掉，考虑莫比乌斯反演，考虑全部 $a,b$ 的倍数 $x$，此时只需要把所有的下标都除以 $x$，乘上 $\mu(x)$ 即可。此时我们就转化成了若干个去掉 $\gcd(a,b)=1$ 的原问题。
考虑式子 $q_{ac} < q_{bc} < q_{bd} < q_{ad}$ 如何计算，考虑到这是乘积的形式，因此可以考虑根号分治：
当 $\max(a,b) \leq B$ 时，我们先枚举 $b$，考虑序列 $v_i=(q_{bi},i)$，相当于
std::pair，此时我们对序列 $v$ 排序。此时我们再枚举 $v$ 中的元素 $(q_{bd},d)$，再枚举 $a$，此时我们需要满足 $a \leq \min(B,n/d)$，只需要统计此时是否满足 $q_{bd} < q_{ad}$，并且在之前存在多少个 $d_0$，满足 q_{a d_0 }<q_{d_0 a}（这里 $\text{Latex}$ 渲染有问题，只能直接打出来了）。
由于我们对 $q$ 执行了排序操作，所以 $q_{bd'}< q_{bd}$ 是天然就满足的了。
时间复杂度瓶颈为枚举 $(a,b,d)$ 的复杂度，为 $O(nB)$。具体见代码。
当 $\max(a,b) > B$ 时，我们有 $\max(c,d) < \frac{n}{B}$，我们可以暴力枚举 $(c,d)$，并且对于全部的 $i$，计算点对 $(q_{ci},q_{di})$。
限制形如 $q_{ac} < q_{bc} < q_{bd} < q_{ad}$，相当于是一个二维数点问题，统计在一个矩形中存在多少个点。
有一个细节就是我们不能让 $a,b$ 同时 $\leq B$，否则会与上一种情况算重，此时可以考虑分离查询点与插入点，插入点正常插入，查询的时候，设我们查询的点为 $(x,y)$ 我们需要知道在平面内有多少个点 $(x',y')$ 满足 $x' < x \wedge y < y'$ 即可。此时我们只对 $b>B$ 的情况查询。
这一部分的时间复杂度为 $O\left(\frac{n^2\log n}{B}\right)$。
平衡一下可有 $O\left(n\sqrt{n \log n}\right)$ 的时间复杂度。
总复杂度为 $O\left(\sum\limits_{x=1}^n (n/x)^{1.5} \log n\right)$，采用积分近似：
$$\sum\limits_{x=1}^n (n/x)^{1.5} \sqrt{\log n} \approx \int_{1}^n x^{-1.5} \mathrm{d}x n\sqrt{n\log n} = \frac{1}{-0.5}(n^{-0.5}-1) n\sqrt{n\log n} = O(n\sqrt{n\log n})$$
所以时间复杂度为 $O(n\sqrt{n\log n})$。