bzoj 2244: [SDOI2011]拦截导弹 （CDQ分治+DP）

题目描述

传送门

题目大意：某国为了防御敌国的导弹袭击，发展出一种导弹拦截系统。但是这种导弹拦截系统有一个缺陷：虽然它的第一发炮弹能够到达任意的高度、并且能够拦截任意速度的导弹，但是以后每一发炮弹都不能高于前一发的高度，其拦截的导弹的飞行速度也不能大于前一发。某天，雷达捕捉到敌国的导弹来袭。由于该系统还在试用阶段，所以只有一套系统，因此有可能不能拦截所有的导弹。
在不能拦截所有的导弹的情况下，我们当然要选择使国家损失最小、也就是拦截导弹的数量最多的方案。但是拦截导弹数量的最多的方案有可能有多个，如果有多个最优方案，那么我们会随机选取一个作为最终的拦截导弹行动蓝图。
我方间谍已经获取了所有敌军导弹的高度和速度，你的任务是计算出在执行上述决策时，每枚导弹被拦截掉的概率。

题解

$f[i]$ 表示以i结尾的最长不升子序列的长度
$g[i]$表示以i结尾的最长不升子序列的方案数
$f1[i]$表示以i开头的最长不升子序列的长度
$g1[i]$表示以i开头的最长不升子序列的方案数
那么第一问的答案就是$ans=max\{f[i]+f1[i]-1\}$
第二问的答案，如果$f[i]+f1[i]-1=ans$,那么$p[i]={g[i]*g1[i]\over sum}$
其中$sum={{\sum_{i=1}^n g[i]*g1[i]*[f[i]+f1[i]-1=ans]}\over ans}$
问什么要除以ans,因为每种方案都会被重复计算呀。
剩下的就是CDQ解决三维偏序的问题了，这道题f,g要一起计算和维护，而且还有注意某一维的值相同的问题。
这道题的精度有些差，建议跟g有关的计算都用double。
话说我稀里糊涂的成了code vs 第一个没有被卡精度的人。。。。。

代码

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 100003
#define LL long long 
#define eps 1e-7
using namespace std;
struct data{
    int x,y,id;
    bool pd;
}a[N];
int b[N],c[N],n,f[N],f1[N],tr[N],ans[N],cntx,cnty;
double sum1[N],g[N],g1[N],sum[N];
struct point{
    int x; double y;
};
int cmp(data a,data b)
{
    return a.x<b.x||a.x==b.x&&a.y<b.y||a.x==b.x&&a.y==b.y&&a.pd>b.pd;
}
int cmp1(data a,data b)
{
    return a.id<b.id;
}
int lowbit(int x)
{
    return x&(-x);
}
point query(int x)
{
    point ans; ans.x=ans.y=0;
    for (int i=x;i>=1;i-=lowbit(i)) {
        if (tr[i]>ans.x) ans.x=tr[i],ans.y=sum[i];
        else if (tr[i]==ans.x) ans.y+=sum[i];
    }
    ans.x++;
    return ans;
}
void change(int x,int v,double t)
{
    for (int i=x;i<=cnty;i+=lowbit(i)){
        if (tr[i]<v) tr[i]=v,sum[i]=t;
        else if (tr[i]==v) sum[i]+=t;
    }
}
void clear(int x)
{
    for (int i=x;i<=cnty;i+=lowbit(i)) tr[i]=sum[i]=0;
}
void divide(int l,int r)
{
    if (l==r) return;
    int mid=(l+r)/2;
    divide(l,mid);
    for (int i=l;i<=mid;i++) a[i].pd=1;
    for (int i=mid+1;i<=r;i++) a[i].pd=0;
    sort(a+l,a+r+1,cmp);
    int j=l;
    for (int i=l;i<=r;i=j+1) {
        j=i;
        while (j<r&&a[i].x==a[j+1].x) j++;
        int lx=0;
        for (int k=i;k<=j;k=lx+1) {
           lx=k;
           while (lx<j&&a[k].y==a[lx+1].y) lx++;
           for (int t=k;t<=lx;t++) 
            if (a[t].pd) change(a[t].y,f[a[t].id],g[a[t].id]);
           for (int t=k;t<=lx;t++)
            if (!a[t].pd) {
                point ans=query(a[t].y); int pos=a[t].id;
                if (ans.x>f[pos]) f[pos]=ans.x,g[pos]=ans.y;
                else if (ans.x==f[pos]) g[pos]+=ans.y;
            }
        }
    }
    for (int i=l;i<=r;i++)
     if (a[i].pd) clear(a[i].y);
    sort(a+l,a+r+1,cmp1);
    divide(mid+1,r);
}
void divide1(int l,int r)
{
    if (l==r) return;
    int mid=(l+r)/2;
    divide1(mid+1,r);
    for (int i=l;i<=mid;i++) a[i].pd=0;
    for (int i=mid+1;i<=r;i++) a[i].pd=1;
    sort(a+l,a+r+1,cmp);
    int j=l;
    for (int i=l;i<=r;i=j+1) {
        j=i;
        while (j<r&&a[i].x==a[j+1].x) j++;
        int lx=0;
        for (int k=i;k<=j;k=lx+1) {
           lx=k;
           while (lx<j&&a[k].y==a[lx+1].y) lx++;
           for (int t=k;t<=lx;t++) 
            if (a[t].pd) change(a[t].y,f1[a[t].id],g1[a[t].id]);
           for (int t=k;t<=lx;t++)
            if (!a[t].pd) {
                point ans=query(a[t].y); int pos=a[t].id;
                if (ans.x>f1[pos]) f1[pos]=ans.x,g1[pos]=ans.y;
                else if (ans.x==f1[pos]) g1[pos]+=ans.y;
            }
        }
    }
    for (int i=l;i<=r;i++)
     if (a[i].pd) clear(a[i].y);
    sort(a+l,a+r+1,cmp1);
    divide1(l,mid);
}
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    scanf("%d",&n);
    for (int i=1;i<=n;i++) {
      scanf("%d%d",&a[i].x,&a[i].y);
      b[i]=a[i].x; c[i]=a[i].y;
    }
    sort(b+1,b+n+1); sort(c+1,c+n+1);
    cntx=unique(b+1,b+n+1)-b-1;
    cnty=unique(c+1,c+n+1)-c-1;
    for (int i=1;i<=n;i++) {
        a[i].x=cntx-(lower_bound(b+1,b+cntx+1,a[i].x)-b)+1;
        a[i].y=cnty-(lower_bound(c+1,c+cnty+1,a[i].y)-c)+1;
        a[i].id=i;
        //cout<<a[i].x<<" "<<a[i].y<<endl;
    }
    for (int i=1;i<=n;i++) f[i]=1,g[i]=1;
    divide(1,n);
    //for (int i=1;i<=n;i++) cout<<f[i]<<" "<<g[i]<<endl;
    for (int i=1;i<=n;i++) f1[i]=1,g1[i]=1;
    sort(a+1,a+n+1,cmp1);
    memset(tr,0,sizeof(tr));
    memset(sum,0,sizeof(sum));
    for (int i=1;i<=n;i++) 
     a[i].x=cntx-a[i].x+1,a[i].y=cnty-a[i].y+1;
    divide1(1,n);
    //for (int i=1;i<=n;i++) cout<<f1[i]<<" "<<g1[i]<<endl;
    double tot=0; LL mx=0;
    for (int i=1;i<=n;i++) {
        ans[i]=f[i]+f1[i]-1;
        sum1[i]=g[i]*g1[i];
        if (ans[i]>mx) mx=ans[i],tot=sum1[i];
        else if (ans[i]==mx) tot+=sum1[i];
    }
    tot/=(double)mx;
    printf("%d\n",mx);
    for (int i=1;i<=n;i++) {
        if (ans[i]!=mx) printf("0.00000 ");
        else printf("%.5lf ",fabs(sum1[i]/tot));
    }
    printf("\n");
}

补充

这道题其实我刚开始想的DP转移不是这样的，用的是期望DP，也可以做，但是写起来比较麻烦，而且因为除法的问题可能精度比较差，所以没有改成CDQ分治。不过也是一种思路。
f[i]表示以i结尾的最长不升子序列的长度，g[i]表示的是方案数。
那么设$tot=\sum_{i=1}^n g[i] *[f[i]==ans]$
那么我们到这统计概率即可。如果$f[i]==ans，p[i]={g[i]\over tot}$
然后用$p[i]$去更新$f[j]+1==f[i]，p[j]+={g[j]*p[i]\over g[i]}$

暴力DP

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 100003
using namespace std;
int n,a[N],b[N],f[N],g[N];
double p[N];
int main()
{
    freopen("a.in","r",stdin);
    freopen("my.out","w",stdout);
    scanf("%d",&n);
    for (int i=1;i<=n;i++) scanf("%d%d",&a[i],&b[i]);
    f[1]=1; g[1]=1;
    for (int i=2;i<=n;i++){
        int t=0;
        for (int j=1;j<i;j++)
         if (a[j]>=a[i]&&b[j]>=b[i]) t=max(t,f[j]);
        f[i]=t+1;
        if (f[i]==1) g[i]=1;
        for (int j=1;j<i;j++)
         if (a[j]>=a[i]&&b[j]>=b[i]&&f[j]+1==f[i]) g[i]+=g[j];
    }
    int ans=0;
    for (int i=1;i<=n;i++) ans=max(ans,f[i]);
    printf("%d\n",ans);
    int t=0;
    //for (int i=1;i<=n;i++) cout<<g[i]<<" "; cout<<endl;
    for (int i=1;i<=n;i++)
     if (ans==f[i]) t+=g[i];
    //cout<<t<<endl;
    for (int i=1;i<=n;i++)
     if (ans==f[i]) p[i]=(double)g[i]*1.0/t;
    for (int i=n;i>=2;i--) 
     if (p[i]){
        for (int j=i-1;j>=1;j--)
         if (a[i]<=a[j]&&b[i]<=b[j]&&f[i]==f[j]+1)
          p[j]+=double(1.0*g[j]/g[i])*p[i];
    }
    for (int i=1;i<=n;i++) printf("%.5lf ",p[i]);
    printf("\n");
}