2012-A3-0025

从 Trac 迁移的文章

这是从旧校内 Wiki 迁移的文章，可能存在一些样式问题，您可以向 memset0 反馈。

原文章内容如下：

代码写得比较长，还剩了好些调试信息，懒得删了
{{{
#include<cstdio>
#include<cstring>
#include<cmath>
#include<cassert>
#include<cstdlib>
#include<iostream>
using namespace std;
#define InCase(x) if(__builtin_expect((x),0))
typedef double real;

const real
    g   =   10,
    inf =   1e64,
    eps =   1e-8;

static real
    ans;

#define sqr(x) ((x)*(x))

inline void check(real tmp){
    InCase(tmp<0)
        throw "At check(real): tmp < 0!";
    if(tmp<ans)
        ans= tmp;
}

inline real vat(real y){
// return the speed at y
    const static real k = g*-2.0;
    InCase(y>0){
        throw "At vat(real): Sqrt a negetive real";
    }
    return sqrt(k*y);
}
inline real walk(real xa,real ya,real xb,real yb){
//return the time walk between (xa,yb) and (xb,yb)
    const static real ka = sqrt(2.0/g);

    InCase(ya>0 || yb > 0){
        return inf;
        throw "At walk(real, ...): y higher than 0";
    }
    InCase(ya>eps||yb>eps)return inf;
    if(fabs(ya - yb)<eps){
        if(ya<-eps)
            return fabs(xa - xb)/vat(yb);
        else return inf;
    }
    //printf("%.8lf\n",ka*hypot(xa-xb,ya-yb)/(sqrt(-ya)+sqrt(-yb)));
    return ka*hypot(xa-xb,ya-yb)/(sqrt(-ya)+sqrt(-yb));
}
pair<real,real> cross(real r0,real r1,real xc,real yc){
// a circle at (0,0) radical = r0
// a circle at (xc,yc) radical = r1 
// return the cross point choosing the one with small y
// when xc == 0 ,return the one with positive x
    //x^2 + y ^2 =r0^2
    //(x-xc)^2 + (y-yc)^2 =r1^2
    //2xc*x + 2yc*y = r0^2 -r1^2+xc^2 + yc^2

    static real x,y;
    InCase(fabs(xc)+fabs(yc)<eps){
        throw "At cross(real,...): circle has same middle";
    }
    if(fabs(xc)>eps){
        real    oc = sqr(xc)+sqr(yc), t = oc + sqr(r0)-sqr(r1);
        real    A = oc * 4, B = -4 * t * yc, C = sqr(t)-sqr(r0*2*xc);
        y= (-B-sqrt(sqr(B)-4*A*C))/(2*A);
        x=(t*0.5-y*yc)/xc;
    }else{
        y =(sqr(r0)+sqr(xc)+sqr(yc)-sqr(r1))/(yc*2);
        x = sqrt(sqr(r0)-sqr(y));
    }
    return make_pair(x,y);
}

real walk2(real xd,real yd,real a,real b){
// time from (0,0) to (xd,yd) through pipe length a then lenght b.
    InCase(yd>0){
        return inf;
        throw "At walk2(real, ...): y higher than 0";
    }
    InCase(fabs(xd)+fabs(yd)<eps){
        assert(fabs(a-b)<eps);
        return walk(0,0,0,-a)*2;
    }
    static pair<real,real> p;
    InCase(fabs(hypot(xd,yd)-(a+b))<eps)
        return walk(0,0,xd,yd);
    p = cross(a,b,xd,yd);
    return (p.second< -eps)?
        walk(0,0,p.first,p.second)+walk(p.first,p.second,xd,yd)
    :   inf;
}

class FindMin{
    real x,y,a,b,c;
    const static real
        g1  =   0.6180339887498948482045868343656381177203091798057628,
        g2  =   0.3819660112501051517954131656343618822796908201942372;
    real f(real ang){
        static real tx,ty;
        tx =    x + c * cos(ang);
        ty =    y + c * sin(ang);
        return walk2(tx,ty,a,b)+walk(tx,ty,x,y);
    }
public:
    FindMin(real _x,real _y,real _a,real _b ,real _c){
        x = _x;y =  _y;a =  _a;b =  _b;c = _c;
    }
    real Min(real l,real r){
        static real ml,mr,fml,fmr;
        if(l>r+eps){
            return inf;
            printf("%d %d %d %d %d\n",(int)(x),(int)(y),(int)(a+eps),(int)(b+eps),(int)(c+eps));
            throw "at Min(real,real): bound error";
        }
        if(r-l>100){
            return min(Min(l,(l+r)*0.5+eps),Min((l+r)*0.5-eps,r));
        }
        fml = f(ml = l * g1 + r * g2);
        fmr = f(mr = r * g1 + l * g2);
        for(size_t i=0;i<96;++i)
            if(fml>fmr){
                l = ml;
                fml = fmr;
                ml  = mr;
                fmr = f(mr = r *g1 + l * g2);
            }else{
                r = mr;
                fmr = fml;
                mr = ml;
                fml = f(ml = l * g1 + r * g2);
            }
        /*
        for(size_t i=0;i<200;++i){
            ml = l*0.666666666666666666666666666666667 + r * 0.33333333333333333333333333333;
            mr = r*0.666666666666666666666666666666667 + l * 0.33333333333333333333333333333;
            if(f(ml)>f(mr))l=ml;
                else r= mr;
        }*/
        return f((l+r)*0.5);
    }
};

real walk3(real xd,real yd,real a ,real b, real c){
// time from (0,0) to (xd,yd) through pipe length a then lenght b ,length c.

    static pair<real , real> p ;
    static real l,r,Jl,Jr,tm;
    static real ri,ro,rd;

    ri = fabs(a-b);
    ro = a+b;
    rd = hypot(xd,yd);

    InCase( ri  > rd + c + eps || fabs(rd -c)> ro +eps )
        return inf;

    InCase(fabs(rd -(ro+c))<eps)
        return walk(0,0,xd,yd);

    if(xd>0)xd = -xd;//make it in 3rd Quadrant
    FindMin app(xd,yd,a,b,c);
    InCase(fabs(xd)+fabs(yd)<eps){
        // handle (0,0)to(0,0)
        l   =   -M_PI;
        r       =   -acos((sqr(a)+sqr(c)-sqr(b))/(2*a*c));
        return app.Min(l,r);
    }else{
        l = (yd<-eps)?atan2(yd,xd):-M_PI;
        r = l+ M_PI;
        if( ro < rd + c){
            p = cross (ro,c,xd,yd);
            l = max(l,atan2(p.second-yd,p.first-xd));
        }
        if( ri > fabs(rd - c)){
            p = cross (ri,c,xd,yd);
            r = min(r,atan2(p.second-yd,p.first-xd));
        }
        tm = atan2(yd,xd);
        Jr = tm*2-l; 
        Jl = Jr - (r-l);
        InCase(c > -yd){
            r = min(r , atan2(-yd, sqrt(sqr(c)-sqr(yd))));
            Jl = max(Jl , atan2(-yd, -sqrt(sqr(c)-sqr(yd)))-2*M_PI);
        }
        return min(app.Min(l,r),app.Min(Jl,Jr));
    }
    //if(l>r+eps)return inf;
    //return min(FindMin(xd,yd,a,b,c).Min(0,M_PI),FindMin(xd,yd,a,b,c).Min(-M_PI,0));
}
void test(){
    printf("%.12lf\n",walk(0,0,0,-5));
    printf("%.12lf\n",walk(0,0,5,-5));
    printf("%.12lf\n",walk(0,-5,10,-5));
    printf("%.12lf\n",walk(9,-5,9,-20));

    printf("%.12lf\n",walk2(5,-5,5,5));
    printf("%.12lf\n",walk2(-3,-4,1,4));
}
int main(){
    //test();return 0;
    for(int l[3],x,y;scanf("%d%d%d%d%d",&x,&y,l,l+1,l+2)==5;
      (ans<0.5*inf)?printf("%.15lf\n",ans):puts("Impossible!")){
        assert(abs(x)<=1000);
        assert(abs(y)<=1000);
        try{
            for(int i=0;i<3;++i){
                InCase(l[i]<1||l[i]>1000)throw "error input";
                assert(l[i]>=1 && l[i]<=1000);
            }
            ans = inf;
            if(y>0)continue; 
            real dis = hypot((real)x,(real)y);
            for(int i=0;i<3;++i)
                if(fabs(dis - l[i])<eps) 
                    check(walk(0,0,x,y));
            for(int i=0;i<3;++i)
                for(int j=0;j<3;++j)
                if(i!=j){
                    if(l[i]+l[j]+eps>dis && fabs(l[i]-l[j])<eps+dis)
                        check(walk2(x,y,l[i],l[j]));
                    check(walk3(x,y,l[i],l[j],l[3-i-j]));
                    //printf("%d %d %d %d %d:%.15lf\n",x,y,l[i],l[j],l[3-i-j],walk3(x,y,l[i],l[j],l[3-i-j]));
                }
        }
        catch(const char* error){
            puts("_____________________________________");
            puts(error);
            puts("at case:");
            printf("%d %d %d %d %d\n",x,y,l[0],l[1],l[2]);
            puts("_____________________________________");
        }
        catch(...){
            exit(-1);
        }
    }
}

}}}

代码写得比较长，还剩了好些调试信息，懒得删了

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cassert>
#include<cstdlib>
#include<iostream>
using namespace std;
#define InCase(x) if(__builtin_expect((x),0))
typedef double real;
const real
    g   =   10,
    inf =   1e64,
    eps =   1e-8;
static real
    ans;
#define sqr(x) ((x)*(x))
inline void check(real tmp){
    InCase(tmp<0)
        throw "At check(real): tmp < 0!";
    if(tmp<ans)
        ans= tmp;
}
inline real vat(real y){
// return the speed at y
    const static real k = g*-2.0;
    InCase(y>0){
        throw "At vat(real): Sqrt a negetive real";
    }
    return sqrt(k*y);
}
inline real walk(real xa,real ya,real xb,real yb){
//return the time walk between (xa,yb) and (xb,yb)
    const static real ka = sqrt(2.0/g);
    InCase(ya>0 || yb > 0){
        return inf;
        throw "At walk(real, ...): y higher than 0";
    }
    InCase(ya>eps||yb>eps)return inf;
    if(fabs(ya - yb)<eps){
        if(ya<-eps)
            return fabs(xa - xb)/vat(yb);
        else return inf;
    }
    //printf("%.8lf\n",ka*hypot(xa-xb,ya-yb)/(sqrt(-ya)+sqrt(-yb)));
    return ka*hypot(xa-xb,ya-yb)/(sqrt(-ya)+sqrt(-yb));
}
pair<real,real> cross(real r0,real r1,real xc,real yc){
// a circle at (0,0) radical = r0
// a circle at (xc,yc) radical = r1 
// return the cross point choosing the one with small y
// when xc == 0 ,return the one with positive x
    //x^2 + y ^2 =r0^2
    //(x-xc)^2 + (y-yc)^2 =r1^2
    //2xc*x + 2yc*y = r0^2 -r1^2+xc^2 + yc^2
    static real x,y;
    InCase(fabs(xc)+fabs(yc)<eps){
        throw "At cross(real,...): circle has same middle";
    }
    if(fabs(xc)>eps){
        real    oc = sqr(xc)+sqr(yc), t = oc + sqr(r0)-sqr(r1);
        real    A = oc * 4, B = -4 * t * yc, C = sqr(t)-sqr(r0*2*xc);
        y= (-B-sqrt(sqr(B)-4*A*C))/(2*A);
        x=(t*0.5-y*yc)/xc;
    }else{
        y =(sqr(r0)+sqr(xc)+sqr(yc)-sqr(r1))/(yc*2);
        x = sqrt(sqr(r0)-sqr(y));
    }
    return make_pair(x,y);
}
real walk2(real xd,real yd,real a,real b){
// time from (0,0) to (xd,yd) through pipe length a then lenght b.
    InCase(yd>0){
        return inf;
        throw "At walk2(real, ...): y higher than 0";
    }
    InCase(fabs(xd)+fabs(yd)<eps){
        assert(fabs(a-b)<eps);
        return walk(0,0,0,-a)*2;
    }
    static pair<real,real> p;
    InCase(fabs(hypot(xd,yd)-(a+b))<eps)
        return walk(0,0,xd,yd);
    p = cross(a,b,xd,yd);
    return (p.second< -eps)?
        walk(0,0,p.first,p.second)+walk(p.first,p.second,xd,yd)
    :   inf;
}
class FindMin{
    real x,y,a,b,c;
    const static real
        g1  =   0.6180339887498948482045868343656381177203091798057628,
        g2  =   0.3819660112501051517954131656343618822796908201942372;
    real f(real ang){
        static real tx,ty;
        tx =    x + c * cos(ang);
        ty =    y + c * sin(ang);
        return walk2(tx,ty,a,b)+walk(tx,ty,x,y);
    }
public:
    FindMin(real _x,real _y,real _a,real _b ,real _c){
        x = _x;y =  _y;a =  _a;b =  _b;c = _c;
    }
    real Min(real l,real r){
        static real ml,mr,fml,fmr;
        if(l>r+eps){
            return inf;
            printf("%d %d %d %d %d\n",(int)(x),(int)(y),(int)(a+eps),(int)(b+eps),(int)(c+eps));
            throw "at Min(real,real): bound error";
        }
        if(r-l>100){
            return min(Min(l,(l+r)*0.5+eps),Min((l+r)*0.5-eps,r));
        }
        fml = f(ml = l * g1 + r * g2);
        fmr = f(mr = r * g1 + l * g2);
        for(size_t i=0;i<96;++i)
            if(fml>fmr){
                l = ml;
                fml = fmr;
                ml  = mr;
                fmr = f(mr = r *g1 + l * g2);
            }else{
                r = mr;
                fmr = fml;
                mr = ml;
                fml = f(ml = l * g1 + r * g2);
            }
        /*
        for(size_t i=0;i<200;++i){
            ml = l*0.666666666666666666666666666666667 + r * 0.33333333333333333333333333333;
            mr = r*0.666666666666666666666666666666667 + l * 0.33333333333333333333333333333;
            if(f(ml)>f(mr))l=ml;
                else r= mr;
        }*/
        return f((l+r)*0.5);
    }
};
real walk3(real xd,real yd,real a ,real b, real c){
// time from (0,0) to (xd,yd) through pipe length a then lenght b ,length c.
    static pair<real , real> p ;
    static real l,r,Jl,Jr,tm;
    static real ri,ro,rd;
    ri = fabs(a-b);
    ro = a+b;
    rd = hypot(xd,yd);
    InCase( ri  > rd + c + eps || fabs(rd -c)> ro +eps )
        return inf;
    InCase(fabs(rd -(ro+c))<eps)
        return walk(0,0,xd,yd);
    if(xd>0)xd = -xd;//make it in 3rd Quadrant
    FindMin app(xd,yd,a,b,c);
    InCase(fabs(xd)+fabs(yd)<eps){
        // handle (0,0)to(0,0)
        l   =   -M_PI;
        r       =   -acos((sqr(a)+sqr(c)-sqr(b))/(2*a*c));
        return app.Min(l,r);
    }else{
        l = (yd<-eps)?atan2(yd,xd):-M_PI;
        r = l+ M_PI;
        if( ro < rd + c){
            p = cross (ro,c,xd,yd);
            l = max(l,atan2(p.second-yd,p.first-xd));
        }
        if( ri > fabs(rd - c)){
            p = cross (ri,c,xd,yd);
            r = min(r,atan2(p.second-yd,p.first-xd));
        }
        tm = atan2(yd,xd);
        Jr = tm*2-l; 
        Jl = Jr - (r-l);
        InCase(c > -yd){
            r = min(r , atan2(-yd, sqrt(sqr(c)-sqr(yd))));
            Jl = max(Jl , atan2(-yd, -sqrt(sqr(c)-sqr(yd)))-2*M_PI);
        }
        return min(app.Min(l,r),app.Min(Jl,Jr));
    }
    //if(l>r+eps)return inf;
    //return min(FindMin(xd,yd,a,b,c).Min(0,M_PI),FindMin(xd,yd,a,b,c).Min(-M_PI,0));
}
void test(){
    printf("%.12lf\n",walk(0,0,0,-5));
    printf("%.12lf\n",walk(0,0,5,-5));
    printf("%.12lf\n",walk(0,-5,10,-5));
    printf("%.12lf\n",walk(9,-5,9,-20));
    printf("%.12lf\n",walk2(5,-5,5,5));
    printf("%.12lf\n",walk2(-3,-4,1,4));
}
int main(){
    //test();return 0;
    for(int l[3],x,y;scanf("%d%d%d%d%d",&x,&y,l,l+1,l+2)==5;
      (ans<0.5*inf)?printf("%.15lf\n",ans):puts("Impossible!")){
        assert(abs(x)<=1000);
        assert(abs(y)<=1000);
        try{
            for(int i=0;i<3;++i){
                InCase(l[i]<1||l[i]>1000)throw "error input";
                assert(l[i]>=1 && l[i]<=1000);
            }
            ans = inf;
            if(y>0)continue; 
            real dis = hypot((real)x,(real)y);
            for(int i=0;i<3;++i)
                if(fabs(dis - l[i])<eps) 
                    check(walk(0,0,x,y));
            for(int i=0;i<3;++i)
                for(int j=0;j<3;++j)
                if(i!=j){
                    if(l[i]+l[j]+eps>dis && fabs(l[i]-l[j])<eps+dis)
                        check(walk2(x,y,l[i],l[j]));
                    check(walk3(x,y,l[i],l[j],l[3-i-j]));
                    //printf("%d %d %d %d %d:%.15lf\n",x,y,l[i],l[j],l[3-i-j],walk3(x,y,l[i],l[j],l[3-i-j]));
                }
        }
        catch(const char* error){
            puts("_____________________________________");
            puts(error);
            puts("at case:");
            printf("%d %d %d %d %d\n",x,y,l[0],l[1],l[2]);
            puts("_____________________________________");
        }
        catch(...){
            exit(-1);
        }
    }
}