【BZOJ3616】War，KD树+bitset压位

思路：
一开始想到顺序无关、轮数独立啥的，想用DP来做，但发现怎么定义状态都有后效性，想套容斥发现也没有什么用，后来聪哥提供思路(bao chu zheng jie)，对于一些概率期望的题目，除DP外还可以考虑每一个元素的贡献
沿着这种思路想一下，求解就可以考虑每个阵营的贡献，即阵营i的所有炮塔一轮中都不会受攻击的概率$P_i$的$m$次方，累和起来答案就是$\displaystyle \sum^k_{i=1} P_i^m$
$P_i=\frac {n-s_i} n$
$s_i$指能攻击到种类i炮塔的炮塔数，所以问题的关键就是求$s_i$了
“二维平面中求与某点距离在某一范围内的元素个数”可以用kd树，所以枚举每个炮塔，把能打到的其他炮塔标记一下就可以了
但是对于某一炮塔，要存每个类型是否已经打过它了，暴力搞是$O(nk)$的内存，然后就是第一次用bitset压位来搞这种题……
（感觉bitset处理这种记录状态的题目好劲啊，以后还是多学习一个的）
如果打标记+存节点本身的状态，好像要用两个35000×35000的bitset？爆M了啊
像标记永久化线段树那样搞倒是只用一个bitset，但节点是2n的，所以并没有什么卵用:-)
mrazer提醒了我，每次kd树的修改最多遍历$\sqrt n$个节点，也就是说除去标记，总共修改的状态数不超过$n\sqrt n$个，所以我们可以把这些修改暴力存下来，至于标记，可以让它永久化，操作全部结束后再放下来
最后统计某一类型的答案时把该类型的bitset都or起来，再处理一下暴力存的点就可以了
不用long double ，double就够了
写的挺快，改了好久，拍了好久，最后还因为没有测试极限数据而没有发现建树写挫了，没能1A:-(
问了问大佬，调用n位的bitset的count()，reset()，set()什么的时间复杂度大概是$O(\frac n {32})$，修改某一位是$O(1)$
时间复杂度应该是$O(n\sqrt n+\frac {n^2} {32}+\frac {nk} {32})$
代码：

#include<cstdio>
#include<bitset>
#include<algorithm>
#include<iostream>
#include<vector>
#define M 35003
#define ls tr[rt].l
#define rs tr[rt].r
#define abs(x) ((x)>0?(x):-(x))
using namespace std;
int n,m,k,root,cnt;
bitset<M>bit[M],tmp;
vector<int>e[M],belong[M];
int data[2][M],R[M],A[M],P[M],id[M];
bool key;
struct kdtree
{
    int l,r,mi[2],mx[2],dat[2],ID;
}tr[M];
bool cmp(int x,int y)
{
    return data[key][x]<data[key][y];
}
inline int dis1(int xa,int ya,int xb,int yb)
{
    return abs(xa-xb)+abs(ya-yb);
}
inline int dis2(int xa,int ya,int xb,int yb)
{
    return (xa-xb)*(xa-xb)+(ya-yb)*(ya-yb);
}
inline int suan1(int rt,int x,int y,bool T)
{
    if (T)//求最小的曼哈顿距离 
    {
        int mi=0;
        if (x-tr[rt].mx[0]>0) mi+=x-tr[rt].mx[0];
        else if (tr[rt].mi[0]-x>0) mi+=tr[rt].mi[0]-x;
        if (y-tr[rt].mx[1]>0) mi+=y-tr[rt].mx[1];
        else if (tr[rt].mi[1]-y>0) mi+=tr[rt].mi[1]-y;
        return mi;

    }
    else//求最大的曼哈顿距离
        return max(abs(y-tr[rt].mx[1]),abs(tr[rt].mi[1]-y))+max(abs(x-tr[rt].mx[0]),abs(tr[rt].mi[0]-x));
}
inline int suan2(int rt,int x,int y,bool T)
{
    int X=0,Y=0;
    if (T)//求最小的欧氏距离
    {
        if (x-tr[rt].mx[0]>0) X=x-tr[rt].mx[0];
        else if (tr[rt].mi[0]-x>0) X=tr[rt].mi[0]-x;
        if (y-tr[rt].mx[1]>0) Y=y-tr[rt].mx[1];
        else if (tr[rt].mi[1]-y>0) Y=tr[rt].mi[1]-y;
        return X*X+Y*Y;
    }
    else
    {
        X=max(abs(x-tr[rt].mx[0]),abs(tr[rt].mi[0]-x));
        Y=max(abs(y-tr[rt].mx[1]),abs(tr[rt].mi[1]-y));
        return X*X+Y*Y;
    }
}
int build(int l,int r,bool T)
{
    int mid=(l+r)>>1,rt=++cnt;
    key=T; 
    nth_element(id+l,id+mid,id+r,cmp);
    tr[rt].ID=id[mid];
    tr[rt].dat[0]=tr[rt].mi[0]=tr[rt].mx[0]=data[0][id[mid]];
    tr[rt].dat[1]=tr[rt].mi[1]=tr[rt].mx[1]=data[1][id[mid]];
    belong[P[id[mid]]].push_back(rt);
    if (l<mid) ls=build(l,mid-1,T^1);
    if (mid<r) rs=build(mid+1,r,T^1);
    for (int i=0;i<2;++i)
    {
        if (ls)
            tr[rt].mi[i]=min(tr[rt].mi[i],tr[ls].mi[i]),
            tr[rt].mx[i]=max(tr[rt].mx[i],tr[ls].mx[i]);
        if (rs)
            tr[rt].mi[i]=min(tr[rt].mi[i],tr[rs].mi[i]),
            tr[rt].mx[i]=max(tr[rt].mx[i],tr[rs].mx[i]);
    }
    return rt;
}
inline void update(int rt,int id)
{
    int mx1=suan1(rt,data[0][id],data[1][id],0),mx2=suan2(rt,data[0][id],data[1][id],0),
        mi1=suan1(rt,data[0][id],data[1][id],1),mi2=suan2(rt,data[0][id],data[1][id],1);
    if (mx1<=A[id]||mx2<=R[id]*R[id]) return void(bit[rt].set(id));
    if (mi1>A[id]&&mi2>R[id]*R[id]) return;
    int d1=dis1(tr[rt].dat[0],tr[rt].dat[1],data[0][id],data[1][id]),
        d2=dis2(tr[rt].dat[0],tr[rt].dat[1],data[0][id],data[1][id]);
    if (d2<=R[id]*R[id]||d1<=A[id]) e[rt].push_back(id);
    if (ls) update(ls,id);
    if (rs) update(rs,id);
}
void down(int rt)
{
    if (ls) bit[ls]|=bit[rt],down(ls);
    if (rs) bit[rs]|=bit[rt],down(rs);
}
inline double qr(double x,int y)
{
    double t=1.0;
    for (;y;y>>=1,x=x*x)
        if (y&1)
            t=t*x;
    return t;
}
main()
{
    scanf("%d%d%d",&n,&m,&k);
    for (int i=1;i<=n;++i)
        scanf("%d%d%d%d%d",data[0]+i,data[1]+i,R+i,A+i,P+i),
        id[i]=i;
    root=build(1,n,0);
    for (int i=1;i<=n;++i) update(root,i);
    down(root);
    double ans=0;
    for (int i=1;i<=k;++i)
    {
        tmp.reset();
        for (int v,j=0;j<belong[i].size();++j)
        {
            v=belong[i][j];
            tmp|=bit[v];
            for (int p=0;p<e[v].size();++p) tmp.set(e[v][p]);
        }
        for (int j=0;j<belong[i].size();++j)
            tmp.reset(tr[belong[i][j]].ID);
        ans+=qr(1.0*(n-tmp.count())/n,m);
    }
    printf("%.5lf\n",ans);
}