C++博客-From A Start,As An Acmer

CF 159 A [没啥，只是想记录一下自己想的贪心是错的]

BUPT-[aswmtjdsj] @ Penalty — Tue, 27 Mar 2012 05:12:00 GMT

orz。如此一个水题让我想了很久用了无数STL之后已然WA到死。直到想明白自己哪里错了。
哎，自己的贪心果断不靠谱啊。但是为啥他们那些如此暴力的做法都能过呢- -构造卡N^3的数据太难了是么。。

题目大意：
n(int,[1,1000]), d(int,[1,1000])
n条记录，每条记录有a,b(string,len[1,20]),t(int,[0,10000])。按非降序列给出。
a,b,t表示在t时刻，a给b发了一条信息。
约束条件，如果a给b在t1时刻发了信息，且b给a在t2时刻发了信息，且0 < t2 - t1 <= d。那么a和b就是朋友关系。
(不存在自己给自己发消息的情况）

让你给出朋友关系的条数，和每条朋友关系，spj。

我的做法：
错误做法，N*logN的贪心：
逐条读入，用string存unique的string，map存映射。用int[][]存发送时间，abt则rel[index(a)][index(b)] = max(t, rel[index(a)][index(b)]。并与rel[index(b)][index(a)]比较，如果满足约束条件则将(a,b)插入一个set中。
错误关键点：在线更新rel[][]是不对的，虽说两个人的朋友关系与否只和最近一次的两人互通消息有关，但是这中间存在一个反例（由于题目条件的特殊性）：
a b 1
a b 2
b a 2
这组例子a、b是朋友关系。但是由于rel[index(a)][index(b)] = 2 其矩阵对称值也为2，所以不满足约束关系，所以就判错了。
所以一个简单的正确做法为离线做。
正确做法，N^2*logN的枚举：
对于每条记录枚举之前所有记录，看满足约束条件否，满足则扔set。。
（ps，有人用vector代替set都不超时。。。N^3啊，泪目）
（贪心什么的，一证明，二想反例，一定要慎重啊。自己还是太弱了。

太伤了。。。

BUPT-[aswmtjdsj] @ Penalty 2012-03-27 13:12 发表评论

TC 538 DIV II 300+500

BUPT-[aswmtjdsj] @ Penalty — Wed, 21 Mar 2012 02:32:00 GMT

最近TC一路猛跌，现在光荣的在div2厮混。昨晚打了一场，全是ws题，可惜自己手速不够快，没抢到分也没cha够人。

300：
给你一个只包含LR?的序列，表示你每一步走的方向，？表示L、R都有可能，问你最远你能走多远（过程中）。

标准做法：贪心，枚举所有？为L或R。两遍扫。大略的证明，如果最远点在左边，那么过程中向右绝对没好处；反之亦然。
for i = 0 ,a = 0, i < len , i++
if s[i] = 'L' , a++
else a--
ans = max(abs(a),ans)

for i = 0 ,a = 0, i < len , i++
if s[i] = 'R' , a++
else a--
ans = max(abs(a),ans)

我的做法：由于本人经研究发现，自己贪心最弱，所以遇题从不敢写贪心，于是剪枝爆搜了一把。
O(2^(len/2))的搜索，亏着len <= 50
将原串中连续的L、R、？压缩成一个pair。LLL=pair(L,3)，塞到一个vector里面。
dfs之，遇到？pair的时候分化。
由于最坏情况为L?R?L?R?......所以最坏有25个离散的?，离散的LR不影响，因为不产生分支。
所以为O(2^25)的搜索。tc服务器比较快，俨然过了。

500：
没有发现水题本质哇，就是一个猜结论题。
你一开始在(0,0)，给你一个xy坐标序列，一个奇偶预测p。问你是否能找到一个走法序列使得你的步数满足这个预测。
举例：序列为(1,1),(1,0),p=1
则(0,0)-偶->(1,1)-奇->(1,0) ：奇，所以存在这样的序列。

结论：
咱们画个图表示一下：上层点代表偶点，下层点代表奇点。
奇解的情况：
________
\
--------------
偶解的情况：

________ __________
\ /
--------------

将序列的每个元素自加，奇就为奇点，偶就为偶点。
因为初始点为偶点，那么序列中存在的偶点有两个选择，要么缀在这个偶点之后，要么缀在之后的某个奇点之后。
假如存在奇点，那么也有两个选择，要么缀在偶点之后，要么插入偶点之间。
如果奇缀偶后为最终结果，那么该解为奇解。
如果奇插偶中为最终结果，那么该解为偶解。
那么显然，当序列中不存在偶点一定得不到偶解；不存在奇点则一定得不到奇解。

竟然想了半个小时，泪目。

BUPT-[aswmtjdsj] @ Penalty 2012-03-21 10:32 发表评论

sgu 108[卡内存什么的]

BUPT-[aswmtjdsj] @ Penalty — Fri, 16 Mar 2012 12:35:00 GMT

题目大意：
给一个n，另d(n)= n + sum{digit_i}(for i = 1 to length(n), digit_i is n's no.i digit)
若y = d(n)，则n就是y的generator。
若某个数没有generator，那么这个数就叫self number。
要求：给你一个N[1,10^7]，一个k[1,5000]
k个数，s1...sk。
让你求出[1,N]中有多少个self numbers。并输出k个a[si]，即第si个self number。保证si不大于N中self number的数量。

做法，本来很水的一题，用类似筛法的思想筛一遍就搞定了。
但是因为当年的SGU很穷，所以没钱买内存，于是乎内存只给了4M。必须用点ws的方法了。

10^7的bool数组大小约为9.5M，10^6的int数组约为4M内存（10^7中一共有977xxx个self numbers），省内存就要从两方面入手。

节省int数组的开销，即不存储这个数组，而是把输入的si排序，当扫描到对应的si的时候，直接标记下来。注意wa点：si不一定有序，且有可能重复。
节省bool数组的开销：
法1：位压缩，bitset或者手写位运算均可。
法2：滚动数组。由于d(n) <= n + 9*7在10^7的范围内。所以只需开一个64大小的滚动数组即可。
实现：
if(vis[i%64])//if(vis[i&63])
{
//blabla
}
vis[i%64] = true;//这时候i已经扫描过去了，所以对于下一个等于i%64的i，默认是true的
vis[d(i)%64] = false;//筛掉d(i)

虽然我莫名其妙的G++WA on testcase9，C++AC了吧。。

各种小技巧还是挺有意思的。

BUPT-[aswmtjdsj] @ Penalty 2012-03-16 20:35 发表评论

三分法 for 单峰函数[数值计算][题目sgu 114, zoj 3421]

BUPT-[aswmtjdsj] @ Penalty — Mon, 20 Feb 2012 05:22:00 GMT

（好多地方其实画图比较好，奈何linux的各种画图工具用着不舒服。。又不愿意用opengl什么的画图，那么各位看官就锻炼一下自己的yy能力吧～）
所谓单调函数二分法，单峰函数三分法。

单调函数二分法，函数在给定定义域内单调，可以根据单调性来二分逼近答案。
单峰函数三分法，函数在给定定义域内有不超过一个的峰值，可以根据其单峰特性来三分逼近答案。

只说三分：
迭代到当前轮：
迭代的当前区间[l,r]，区间长度为d=r-l，那么三分点为m1=l+d，m2=r-d。4点三等分区间，这样保证每次迭代后，新区间变为原区间的2/3;如果用1/2,3/4分点的话，迭代后有可能变为原来的3/4,这样比较慢。log1.25(n)毕竟比log1.3333(n)要大啊～
设目标函数是f(x)，比较f(m1) , f(m2)，如果f(m1)比f(m2)更趋近峰值，那么就让r=m2,即右端点左移；反之左端点右移。仔细想想就是，如果左边更靠近峰值的话，那么肯定让右侧向左侧移动；反之亦然。
伪码：
while(l < r)//dblcmp is preferred
{
d = (r - l) / 3;
m1 = l +d;
m2 = m1 +d;
if(f(m1) better than f(m2))
r = m2;
else
l = m1;
}
l or r is the solution

下面谈到具体问题：
sgu 114, 题目大意就不说了，自己看去，抽象出来的方程就是。
f(x) = sum{|x-y[i]| * p[i]},i = 1 to n
ans = min(f(x)) 其中y[i]、p[i]给定
这题做法很多，求中位数、分段二分、三分皆可做；
说一下思路吧：
仔细观察发现，这是一个分段单调的函数。
以ｘ为自变量会发现，在每一个分段内，f(x) = f(floor(x)) + (sum{p[i]}(小于x的每个ｉ} - sum{p[i]}(大于x的每个i)) * (x - floor(x))
其中含有sum的那部分就是每个分段的不同点：斜率或者说单调性；可以发现这个单调性是随着x的增长（分段），不断下降的，从单调降逐渐过度到单调升（如果存在这两个过程的话）。
所以此函数单峰。。其他的就不用多说了。

zoj 3421,10年成都现场赛F题，印象深刻啊，想当年不知道单峰函数可以三分，不然没准就可以拿到silver了。。

一堆a>=0的二次函数Si(x)，f(x) = max(Si(x))，求在x belongs to [0,1000]内的min(f(x))
稍微画一下图便可以看出来，f(x)的图像是一个单峰的函数，根据所有的Si(x)a>=0也可以看出来，f(x)的函数值越靠近两侧应该越大，靠近中间应该越低，所以出现了谷（反着的峰哈。。）。。。然后三分之。。。
精度要求挺高，输出1e-4,而eps到1e-8竟然还不够，需要1e-10。。。
另外，输出的是f(ans)，不是ans。。。坑死我了，10多炮wa啊。。。
代码什么的就不attach了。
以上。。三分什么的，也是折半思想一种典型应用。欢迎吐槽+拍砖。

BUPT-[aswmtjdsj] @ Penalty 2012-02-20 13:22 发表评论

一份题目包中的内容[for 出题者]

BUPT-[aswmtjdsj] @ Penalty — Sat, 04 Feb 2012 17:25:00 GMT

好久没写东西了，随便灌点水充充门面。

最近总出题，作为一道比较正规的题目，它的题目包中应包含以下部分，我认为：

[] means '1 or more'

* means '0 is okay'

| means 'or'

*[Problem Description(zh_cn)]

Problem Description(en_us)

Additional Parts(img or others)

Data Generator(py or cpp or whatever)

Data(data.in and data.out)

*[Data Checker(py or cpp)]

*[Standard Program(brute-force)]//maybe the same as the one just above

*[Standard Program(AC with no STL)]

Standard Program(AC with total STL)

*[Standard Program(others')](For Check)

[Solution Report(zh_cn | en_us)]

BUPT-[aswmtjdsj] @ Penalty 2012-02-05 01:25 发表评论

HDU 4087 & BOJ 244 A Letter to Programmers 【三维仿射变幻的齐次矩阵乘+快速幂】

BUPT-[aswmtjdsj] @ Penalty — Mon, 07 Nov 2011 11:55:00 GMT

题意很裸，亮点在于将齐次仿射变换矩阵和快速幂结合起来。

题目读入延续了pku出题的ws风格，还得用一个dfs处理读入。
题目输出也很蛋疼，标准写法竟然会导致-0.0出现，必须手动加一个eps来修改偏差。没有spj的二货。

附上仿射变换矩阵：

平移

translate tx ty tz

1 0 0 tx

0 1 0 ty

0 0 1 tz

0 0 0 1

缩放

scale kx ky kz

kx 0 0 0

0 ky 0 0

0 0 kz 0

0 0 0 1

绕任意轴（过原点）旋转（注意要把轴向量归一化，不然会在“点在轴上”这个情况下出问题）

rotate x y z d

(1-cos(d))*x*x+cos(d) (1-cos(d))*x*y-sin(d)*z (1-cos(d))*x*z+sin(d)*y 0

(1-cos(d))*y*x+sin(d)*z (1-cos(d))*y*y+cos(d) (1-cos(d))*y*z-sin(d)*x 0

(1-cos(d))*z*x-sin(d)*y (1-cos(d))*z*y+sin(d)*x (1-cos(d))*z*z+cos(d) 0

0 0 0 1

附代码，一开始装b用stack模拟读入过程，re啊wa啊各种sb。最后被迫用了dfs。唉，naive。

#include
#include
#include
#include
#include
#include
using namespace std;
const double eps = 1e-6;
const double pi = acos(-1.0);
const int maxn = 1005;
int sgn(double x){return x < -eps?-1:x > eps;}
#define rep(i,j,k) for(int (i) = (j);(i) < (k);(i) ++)
struct matrix
{
    int r,c;
    double v[5][5];
    matrix(){}
    matrix(int a,int b):r(a),c(b)
    {
        rep(i,0,r)
            fill(v[i],v[i] + c,0.0);
        if(r == c)
        {
            rep(i,0,r)
                v[i][i] = 1.0;
        }
    }
    matrix init()
    {
        return matrix(r,c);
    }
    matrix operator * (const matrix p)
    {
        matrix ans(r,p.c);
        rep(i,0,r)
            rep(j,0,ans.c)
            {
                double tmp = 0.0;
                rep(k,0,c)
                    tmp += v[i][k] * p.v[k][j];
                ans.v[i][j] = tmp;
            }
        return ans;
    }
};
matrix power(matrix a,int n)
{
    matrix e(a.c,a.c);//square matrix of a.c * a.c
    if(n == 0)
        return e;
    matrix ans = power(a,n/2);
    ans = ans * ans;
    return (n%2)?ans * a : ans;
}
struct point
{
    double x,y,z;
    point(){}
    point(double a,double b,double c):x(a),y(b),z(c){}
    void in()
    {
        scanf("%lf %lf %lf",&x,&y,&z);
    }
    void print()
    {
        printf("%.2lf %.2lf %.2lf\n",x + eps,y + eps,z + eps);
    }
}p[maxn];
int n;
point fuck(point now,matrix a)
{
    matrix o(4,1);
    o.v[0][0] = now.x;
    o.v[1][0] = now.y;
    o.v[2][0] = now.z;
    o.v[3][0] = 1.0;
    matrix ans = a * o;
    return point(ans.v[0][0],ans.v[1][0],ans.v[2][0]);
}
char op[20];
matrix dfs(int times)
{
    matrix ans(4,4);
    double v[5];
    int what;
    while(scanf("%s",op) == 1)
    {
        matrix temp(4,4);
        if(op[0] == 'e')
            break;
        else if(op[0] == 'r')
        {
            if(op[1] == 'e')
            {
                scanf("%d",&what);
                temp = dfs(what);
            }
            else
            {
                rep(i,0,4)
                    scanf("%lf",&v[i]);
                double l = 0.0;
                rep(i,0,3)
                    l += v[i] * v[i];
                l = sqrt(l);
                v[3] = v[3] / 180.0 * pi;
                double ct = cos(v[3]),st = sin(v[3]);
                rep(i,0,3)
                    v[i] /= l;
                temp.v[0][0] = (1.0 - ct) * v[0] * v[0] + ct;
                temp.v[0][1] = (1.0 - ct) * v[0] * v[1] - st * v[2];
                temp.v[0][2] = (1.0 - ct) * v[0] * v[2] + st * v[1];
                temp.v[1][0] = (1.0 - ct) * v[1] * v[0] + st * v[2];
                temp.v[1][1] = (1.0 - ct) * v[1] * v[1] + ct;
                temp.v[1][2] = (1.0 - ct) * v[1] * v[2] - st * v[0];
                temp.v[2][0] = (1.0 - ct) * v[2] * v[0] - st * v[1];
                temp.v[2][1] = (1.0 - ct) * v[2] * v[1] + st * v[0];
                temp.v[2][2] = (1.0 - ct) * v[2] * v[2] + ct;
            }
        }
        else if(op[0] == 't')
        {
            rep(i,0,3)
            {
                scanf("%lf",&v[i]);
                temp.v[i][3] = v[i];
            }
        }
        else if(op[0] == 's')
        {
            rep(i,0,3)
            {
                scanf("%lf",&v[i]);
                temp.v[i][i] = v[i];
            }
        }
        ans = temp * ans;
    }
    return power(ans,times);
}
int main()
{
    while(scanf("%d",&n) == 1 && n)
    {
        matrix temp(4,4);
        temp = dfs(1) * temp;
        for(int i = 0;i < n;i++)
        {
            p[i].in();
            p[i] = fuck(p[i],temp);
            p[i].print();
        }
        puts("");
    }
}

BUPT-[aswmtjdsj] @ Penalty 2011-11-07 19:55 发表评论

HDU 2297 Run 【半平面交 or 维护单调性】

BUPT-[aswmtjdsj] @ Penalty — Sat, 05 Nov 2011 13:32:00 GMT

题目大意：
给n个人参加赛跑，给出某一时刻每个人的位置和速度（保证单向，速度恒定，且互不相同，且不为0），问，假如无限跑下去，过程中有可能有多少个不同的人领先过？（如果某一时刻两人同时领先，则都不算）

做法：
画成s-t图的话，很容易得出半平面交的做法，即以y轴和无穷远向这n条射线包围，求得围成的多边形内核上有多少条线段的问题（去掉y轴和无穷远）。
半平面交N^logN的我不会写，于是就用了himdd大神的模板，算法来源zzy大神的论文。
由于最后求得的点集可能有重点问题，需要先sort再unique一下去重。（表示这份代码用到了complex类，unique，distance等stl。。。太orz了）
代码：

#include
#include
#include
#include
#include
#include
using namespace std;
const double inf=1e-10;
const int Max=60010;
typedef complex<double> Point;
inline int dbcmp(double tp){return tp<-inf?-1:tp>inf;}
double operator^(Point a,Point b)
{return imag(conj(a)*b);}
struct Line{
    Point a,b;
    double angle;
    Line(){}
    Line(Point a,Point b):a(a),b(b){angle=arg(a-b);}
};
bool operator<(Line p,Line l){
    return dbcmp(p.angle-l.angle)<0||
        (dbcmp(p.angle-l.angle)==0&&dbcmp(p.b-p.a^l.b-p.a)<0);
}

Point operator*(const Line& u,const Line&v)
{
    Point tp=u.a;
    double fz=u.a-v.a^v.b-v.a;
    double fm=u.a-u.b^v.b-v.a;
    return tp+(u.b-u.a)*fz/fm;
}
bool onleft(Point p,Line u)
{
    return dbcmp(p-u.a^u.b-u.a)<=0;
}
bool uniquecmp(Line u,Line v)
{
    int a = dbcmp(u.angle-v.angle)==0;
    if(a)
        puts("cao");
    return a;
}
bool iter(Point p,Point q)
{
    return dbcmp(real(p) - real(q)) < 0 || (dbcmp(real(p) - real(q)) == 0 && dbcmp(imag(p) - imag(q)) < 0);
}
Point p[Max];
Line l[Max];
int e[Max];//line的下标
int get_kernel(Line *l,int lsz,Point *p )
{
    int eb=0,ee=2,pb=0,pe=1;
    sort(l,l+lsz);
    p[lsz]=p[0];
    lsz=distance(l,unique(l,l+lsz,uniquecmp));
    e[0]=0;e[1]=1;
    p[0]=l[0]*l[1];
    for(int i=2;i        while(pb!=pe&&!onleft(p[pe-1],l[i]))pe--,ee--;
        while(pb!=pe&&!onleft(p[pb],l[i]))pb++,eb++;
        p[pe++]=l[i]*l[e[ee-1]];
        e[ee++]=i;
    }
    while(pb!=pe&&!onleft(p[pe-1],l[e[eb]]))pe--,ee--;
    //这个循环好像是没有用的。。
    while(pb!=pe&&!onleft(p[pb],l[e[ee-1]]))pb++,eb++;
    p[pe++]=l[e[pb]]*l[e[ee-1]];
    int psz=pe-pb;
    for(int i=0;pb!=pe;i++)p[i]=p[pb++];
    return psz;//点的个数
}
bool uniqueccc(Point p,Point q)
{
    return  dbcmp(real(p)-real(q)) == 0 && dbcmp(imag(p)-imag(q))==0;
}
void gao()
{
    int n;int lsz=0;
    scanf("%d",&n);
    double x,y;
    for(int i=0;i    {
        scanf("%lf%lf",&x,&y);
        l[lsz++]=Line(Point(0,x),Point(1,x+y));
    }
    l[lsz++]=Line(Point(0,0),Point(1,0));
    l[lsz++]=Line(Point(1e22,0),Point(1e22,1));
    l[lsz++]=Line(Point(1,1e22),Point(0,1e22));
    l[lsz++]=Line(Point(0,1e22),Point(0,0));
    int psz=get_kernel(l,lsz,p);
    sort(p,p + psz,iter);
    psz = distance(p,unique(p,p + psz,uniqueccc));
    printf("%d\n",psz-2);
}
int main()
{
    int t;
    scanf("%d",&t);
    for(int i = 0;i < t;i++)
        gao();
}

另一个做法：
学自各种大牛的单调栈。
对于所有元素按照p进行从大到小的排序。因为初始时刻p最大的一定是解集元素之一，
所以虚拟一个元素（p与p最大值相同，速度为0），把这个虚拟元素和第一元素压入栈中，作为初始化。
然后对于之后的第i个元素，先判断它是否比所有已出现过的元素v都大，如果不大于，则直接跳过（因为p小，v同，肯定追不上），再之后就比较该元素追上栈顶下一个元素和栈顶元素追上栈顶下一个元素的时间，如果该元素时间短，就弹栈，直到出现>=为止，把该元素入栈。
最后栈中元素减一即为答案。
算法的原理在于，第一个元素肯定可以领先，然后就比较谁能更早的超过曾经领先过的元素，谁也就曾经“真的”领先过。
栈中元素保证了p单调非增，v单调增。维护出来的依旧是一个凸集，和半平面交效果相同。

有一个问题，很值得思考，这个做法对于所有约束条件符号完全相同的线性规划问题，是有通解性的吧？ax + by + c < 0，保证a，b，< 分别同号的情况下。
“维护单调性”。orz这个神想法。
附代码：

#include
#include
#include
#include
using namespace std;
#include
#define maxn 50005
typedef long long ll;
int t,n;
struct node
{
    int p,v;
    node(){}
    node(int _p,int _v):p(_p),v(_v){}
}nn[maxn];
int ss[maxn];
bool operator < (node a,node b){return a.p > b.p || (a.p == b.p && a.v < b.v);}
bool ok(node a,node b,node c)
{
    return (ll) (b.p - c.p) * (ll) (a.v - b.v) <= (ll)(b.p - a.p) * (ll)(c.v-b.v);
}
void gao()
{
    scanf("%d",&n);
    int a,b;
    for(int i = 0;i < n;i++)
    {
        scanf("%d %d",&a,&b);
        nn[i] = node(a,b);
    }
    sort(nn,nn + n);
    int top = 1;
    nn[n] = node(nn[0].p,0);
    ss[0] = n;
    ss[1] = 0;
    int fuck = nn[0].v;
    for(int i = 1;i < n;i++)
    {
        if(nn[i].v < fuck)
            continue;
        fuck = nn[i].v;
        while(top > 0 && ok(nn[ss[top]],nn[ss[top-1]],nn[i]))
            top--;
        ss[++top]=i;
    }
    printf("%d\n",top);
}
int main()
{
    scanf("%d",&t);
    for(int i = 0;i < t;i++)
        gao();
}

BUPT-[aswmtjdsj] @ Penalty 2011-11-05 21:32 发表评论

UVA 12303 Composite Transformations ［三维仿射矩阵］［更正：旋转轴归一化和平面三点式问题］

BUPT-[aswmtjdsj] @ Penalty — Thu, 03 Nov 2011 17:40:00 GMT

摘要: 三维仿射变换什么的，纠结了很多天想把绕任意轴旋转的齐次矩阵手推出来，最后，结果连sample都过不了，只能等回头再调了，目前先用着模板。所谓仿射变换，咱们通俗的讲就是通过补维把坐标转化为齐次矩阵的形式进行矩阵乘法，这样就能找到平移、拉伸、旋转的通式解。此处要用到的有三个变换矩阵：平移：(偏移量(dx,dy,dz))1 0 0 dx0 1 0 dy0 0 1 dz0 0 0 1放缩：(比例系数(sx... 阅读全文

BUPT-[aswmtjdsj] @ Penalty 2011-11-04 01:40 发表评论

ZJU 3551 Bloodsucker 【简单概率dp】

BUPT-[aswmtjdsj] @ Penalty — Sat, 29 Oct 2011 10:42:00 GMT

今天浙大月赛dp+数学专场。。俩小盆友不在，哥这个从来不写dp和数学的人表示很蛋疼。
很少写dp，很少写概率。哪怕写了一道水题也很高兴哇～
写一下吧～

题目大意：
说你手里有n个生物，n-1个人，1个吸血鬼。每一天，你从这堆里面随机挑出两头来，如果是不同种，那么人就有p的概率会被同化为吸血鬼。
求一个期望，问给定n和p的情况下，所有人变吸血鬼的期望天数D。

做法：
首先是要看出最基本的式子，就是在n个生物，n-x个人，x个吸血鬼的情况下抽出两头不同生物的概率：
假设分先后抽，那么先抽吸血鬼x / n , 后抽人 (n - x) / (n - 1) , 概率为 x * ( n - x) / (n * (n - 1))；先抽人 (n - x) / n , 后抽吸血鬼 x / (n - 1) , 概率为 x * (n - x) / (n * ( n - 1))。
所以th[x]（表示在n个生物，n-x个人，x个吸血鬼的情况下抽出两头不同生物的概率） = 2 * x * (n - x) / (n * (n - 1))。
然后就可以得出在某轮已有x个吸血鬼的前提下，吸血鬼数量增加1的概率为f[x] = th[x] * p 。而不变的概率则为g[x] = 1.0 - f[x]。
则可以建立如下概率转移图。Sx表示当前有x个吸血鬼。第k层的Sx的表示，到了第k-1天，有x个吸血鬼（的概率转移）。

所以可以看出Sx只和Sx-1和Sx本身有关。假设DP[x]表示由初始状态变成有x个吸血鬼所需天数的期望。
DP[x] = sigma(i = 1 to infinity)[f[x-1] * i * g[x-1] ^ (i - 1)]
而最终的期望即为所有期望之和ans = sigma(i=2 to n) DP[i]。
题目的关键就是DP[x]的求解，仔细观察这是一个等差乘等比的无穷级数求和，根据高中所学的错位相减法得到公式：
设Sn=DP[x],a=f[x-1],p=g[x-1]。则Sn-p*Sn=a+a*p+a*p^2....+a*p^n-1-n*a*p^n
Sn=[a*(1-p^n) - n*(1-p)*p^n]/(1-p)^2
lim(n->inf)Sn=a/(1-p)^2=1/a
所以最后的ans=sigma(1/a[i])=sigma(1.0/(2.0*p*i*(n-i)/(n*(n-1))))。
因为n在10^5级别，所以要用double乘，不然int会爆精度。

我也是可以做数学题的哈。

BUPT-[aswmtjdsj] @ Penalty 2011-10-29 18:42 发表评论

UVALIVE 5685 Difficult Routes ［最短路问题变形，坑爹题］

BUPT-[aswmtjdsj] @ Penalty — Fri, 28 Oct 2011 17:01:00 GMT

题目大意：
给你一个三维地图，里面N个点，给出对应的xyz坐标；给你M条边，连接对应编号的点。
每条边有两个值。
cost：两点的欧式距离。
d：两点的高度差（z轴差值的绝对值）除以两点在XY平面的投影的欧式距离，再乘以100，最后下取整，如果a点比b点高或者和b点持平，则d值为0。=>所以边有向。
给你一个起始点编号一个终点编号，再给你一个目标值d。
让你找一条最短路径（cost和最小），要求起点终点为对应的点，且这条路径上的最大d值为给定d值。起终点可能相同。

做法：
就是最短路嘛，特点是要求必须经过满足约束的边，那么就在松弛条件上面做手脚。
把dis数组开成两维的：dis[0][p]表示到该点还没满足条件时候的最短路，dis[1][p]表示到这点已经满足条件时候的最短路。
由于稀疏图，那么就开始spfa了。
从dis[0][s]=0开始，把其他所有值初始化置为－1或者inf。
松弛的时候这样，对于dis[0]，该咋松弛就咋松弛。
对于邻接表中当前边的d值等于给定d值，那么就用当前点的min(dis[0][now],dis[1][now])去更新下一点的dis[1][p]；
对于不等的情况，如果当前点的dis[1][now]不等于－1或者inf，那么就用dis[1][now]去更新dis[1][p]。
然后就完了。答案就是dis[1][t]。inf或者ans。

坑爹之处在于题目中N和M的数据范围少给了个0，RE到死哇。。一开始还tmd TLE。。。擦。。。

话说某些小盆友还用了spfa求最小环的想法，我觉得没必要的说。把满足条件后的点和不满足条件的点当成不同的点，那么根本就不存在环的问题了。。
对于起终点相同的情况，最短路径实际上的松弛方式就是i->j~k->i'。。。看嘛，哪里有环了。

附代码：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#include <queue>
using namespace std;
int n,m,s,t,d;
#define maxn 110005
#define maxm 110005
const double eps = 1e-8;
int sgn(double x){return fabs(x) < eps ? 0 : x < -eps ? -1 : 1;}
struct edge
{
    int p,next,d;
    double c;
    edge(){}
    edge(int a,int b,int _c,double _d):p(a),next(b),d(_c),c(_d){}
}e[maxm];
int tot;
struct point
{
    double x,y,z;
    point(){}
    point(double a,double b,double c):x(a),y(b),z(c){}
    point operator - (const point p)
    {
        return point(x - p.x,y - p.y,z - p.z);
    }
    double norm()
    {
        return sqrt(x * x + y * y + z * z);
    }
    double norm2()
    {
        return sqrt(x * x + y * y);
    }
}p[maxn];
int st[maxn];
void init()
{
    tot = 0;
    fill(st,st + n + 1,-1);
}
void add(int p,int q,int d,double c)
{
    e[tot] = edge(q,st[p],d,c);
    st[p] = tot++;
}
double dis[2][maxn];
bool in[maxn];
void bfs()
{
    for(int i = 0;i < 2;i++)
        for(int j = 1;j <= n;j++)
            dis[i][j] = -1.0;
    fill(in,in + n + 1,0);
    queue <int> Q;
    Q.push(s);
    in[s] = true;
    dis[0][s] = 0.0;
    while(!Q.empty())
    {
        int now = Q.front();
        Q.pop();
        in[now] = false;
        for(int k = st[now];~k;k = e[k].next)
        {
            bool mark = false;
            int p = e[k].p;
            int dd = e[k].d;
            double cost = e[k].c;
            if(dd > d)
                continue;
            //cout << p << ' ' << dd << ' ' << cost << endl;
            if(sgn(dis[0][p]+1) == 0 || sgn(dis[0][p] - dis[0][now] - cost) > 0)
                dis[0][p] = dis[0][now] + cost,mark = true;
            if(dd == d)
            {
                double tmp ;
                if(sgn(dis[1][now] + 1.0) != 0)
                    tmp = min(dis[0][now],dis[1][now]);
                else
                    tmp = dis[0][now];
                if(sgn(dis[1][p] + 1.0) == 0 || sgn(dis[1][p] - tmp - cost) > 0)
                    dis[1][p] = tmp + cost,mark=true;
            }
            else if(sgn(dis[1][now]+1) != 0)
            {
                double tmp = dis[1][now];
                if(sgn(dis[1][p] + 1.0) == 0 || sgn(dis[1][p] - tmp - cost) > 0)
                    dis[1][p] = tmp + cost,mark=true;
            }
            if(mark && !in[p])
            {
                Q.push(p);
                in[p] = true;
            }
            //if(Q.size() > 100000)
            //    cout << "fuck" << endl;
        }
    }
}
void gao()
{
    double x,y,z;
    for(int i = 1;i <= n;i++)
    {
        scanf("%lf %lf %lf",&x,&y,&z);
        p[i] = point(x,y,z);
    }
    init();
    for(int i = 0;i < m;i++)
    {
        int a,b;
        scanf("%d %d",&a,&b);
        int fuck = (int)floor(100 * fabs((p[b].z - p[a].z)) / (p[a] - p[b]).norm2());
        double l = (p[a] - p[b]).norm();
        if(sgn(p[a].z - p[b].z) < 0)
        {
            add(a,b,fuck,l);
            add(b,a,0,l);
        }
        else
        {
            add(b,a,fuck,l);
            add(a,b,0,l);
        }
    }
    scanf("%d %d %d",&s,&t,&d);
    bfs();
    double ans = dis[1][t];
    //for(int i = 1;i <= n;i++)
    //    printf("%.3lf %.3lf\n",dis[0][i],dis[1][i]);
    if(sgn(ans + 1.0) == 0)
        puts("None");
    else
        printf("%.3lf\n",ans);
}
int main()
{
    while(scanf("%d %d",&n,&m) == 2 && (n || m))
        gao();
}

BUPT-[aswmtjdsj] @ Penalty 2011-10-29 01:01 发表评论