[2013-2014 Samara SAU ACM ICPC Quarterfinal Qualification Contest pK][BIT][sqrt method] Three Contests

link
You can view the official tutorial here
This is a simple but interesting data structure problem.
Let's talk about some important ideas first. The same discussion can also be read in the official tutorial.

The problem asks for the number of pairs of points(i,j) such that x_i < x_j or y_i<y_j or k_i<k_j.

Since directly calculating the number would be difficult,we can restate the problem first.

The first thing we do is transfer the form of the answer to:
C(n,2) - the number of pairs of points(i,j) such that x_i<x_j and y_i<y_j and z_i<z_j.
we'll see the "and" condition is more convenient to deal with than the "or" condition in this problem.

The other important idea is that three dimensions will be too difficult to deal with, we can make an offline algorithm with z increasing so we can make the problem to 2D.
Assume we are dealing with point i. Note that the points we added in before must have smaller z. So the thing we want to know is the number of points in the rectangle with border points (0,0) and (x_i,y_i).

This query is a classic problem which can be done with 2D Segment Tree. The official tutorial says it is too ugly to implement. Nevertheless , I found the method not only ugly but also unacceptable under the memory constraint. I do believe it's possible to get AC with a 2D Segment Tree, but it surely requires some techniques in memory optimization. (I'm really poor of it orz)

The interesting part comes here. We can use square root method to make time  and memory acceptable. The official tutorial makes it clear enough. I would explain more clearly and completely.

Now , let's develop the intuition behind the square method. Since 2D Segment Tree takes too much memory , the first thing comes to our minds is the necessity of saving the whole graph. Consider a value K . If we only save the sum of rectangle with border points (0,0) and ( (K , 2K , 3K .... ) , (K, 2K , 3K .... ) ) , the number of points on the graph would become n/K * n/K . By maintaining the graph with 2D Fenwick Tree (BIT), we can quickly answer all the queries where both x and y are divisible by K. It takes O(n/K log (n/K) log (n/K)) both in time and space complexity.

What about the points whose x or y is not divisible by K ? We'll just do the simple numeration.The points we can't count with 2D BIT are like
( ( x , x-1 , x-2 ...... x-(x mod K) +1) , (some y < y_i ) ) and y alike. Note that there will be at most K-1 x and K-1 y to enumerate , and the constraint of all the points have different x and y makes it even more convenient to count. We can simply save the corresponding x to every y and y to every x. The enumeration takes O(K) in time complexity and O(n) in space complexity.

Note that you can take K=sqrt(200000) or any other number as long as the memory allows. Basically , more space being used , less time this algorithm takes.

#include<cstdio>
#include<cmath>
#include<algorithm>
#define LL long long
using namespace std;
const int maxn=200000+5;
const int K = 450;
struct Point
{
    int x,y,z;
    bool operator < (const Point& rhs) const {return z<rhs.z;}
}p[maxn];
int n,c[maxn/K+5][maxn/K+5];
inline int lowbit(int x) {return x&(-x);}
int query2D(int qx,int qy)
{
    int ret=0;
    while(qy)
    {
        ret+=c[qx][qy];
        qy-=lowbit(qy);
    }
    return ret;
}
int query1D(int qx,int qy)
{
    int ret=0;
    while(qx)
    {
        ret+=query2D(qx,qy);
        qx-=lowbit(qx);
    }
    return ret;
}
void modify2D(int px,int py,int v)
{
    while(py<=(n/K))
    {
        c[px][py]+=v;
        py+=lowbit(py);
    }
}
void modify1D(int px,int py,int v)
{
    while(px<=(n/K))
    {
        modify2D(px,py,v);
        px+=lowbit(px);
    }
}
int X[maxn],Y[maxn];
int main()
{
    scanf("%d",&n);
    for(int i=0;i<n;i++) scanf("%d%d%d",&p[i].x,&p[i].y,&p[i].z);
    for(int i=0;i<=n;i++) X[i]=n+1,Y[i]=n+1;
    sort(p,p+n);
    LL ans=(LL)n*(n-1)/2LL;
    for(int i=0;i<n;i++)
    {
        ans-=query1D(p[i].x/K,p[i].y/K);
        for(int j=p[i].x;j>p[i].x-(p[i].x%K);j--) if(Y[j]<p[i].y) ans--;
        for(int j=p[i].y;j>p[i].y-(p[i].y%K);j--) if(X[j]<=p[i].x-(p[i].x%K)) ans--;
        X[p[i].y]=p[i].x,Y[p[i].x]=p[i].y;
        modify1D((p[i].x+K-1)/K,(p[i].y+K-1)/K,1);
    }
    printf("%I64d\n",ans);
}