hdu6158 The Designer

链接

点击跳转

题解

这都$21$世纪了，还不写$spj$？

做法很容易，直接反演圆，然后所求的圆就变成夹在两个平行线之间的一些圆，枚举这些圆，反演回去求面积，面积太小就可以停止计算了（这里精度要高一些，实测1e-12可过，1e-8过不了）

代码

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#define iinf 0x3f3f3f3f
#define linf (1ll<<60)
#define eps 1e-8
#define maxn 1000010
#define maxe 1000010
#define cl(x) memset(x,0,sizeof(x))
#define rep(_,__) for(_=1;_<=(__);_++)
#define em(x) emplace(x)
#define emb(x) emplace_back(x)
#define emf(x) emplace_front(x)
#define fi first
#define se second
#define de(x) cerr<<#x<<" = "<<x<<endl
#define sqr(x) ((x)*(x))
using namespace std;
using namespace __gnu_pbds;
typedef long long ll;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
const double pi = acos(-1);
ll read(ll x=0)
{
    ll c, f(1);
    for(c=getchar();!isdigit(c);c=getchar())if(c=='-')f=-f;
    for(;isdigit(c);c=getchar())x=x*10+c-0x30;
    return f*x;
}
int sign(double x)
{
    if(fabs(x)<eps)return 0;
    return x<0?-1:1;
}
struct point
{
    double x, y;
    point(double x, double y):x(x),y(y){}
    point(){}
}_;
struct vec
{
    point o; double x, y;
    vec(point o, double x, double y):o(o),x(x),y(y){}
    vec(){};
    vec operator+(vec v){return vec(_, x+v.x, y+v.y);} //向量和
    vec operator-(vec v){return vec(_, x-v.x, y-v.y);} //向量差
    double operator^(vec v){return x*v.y - v.x*y;} //向量内积
    double operator*(vec v){return x*v.x + y*v.y;} //向量外积
    vec operator*(double t){return vec(_, x*t, y*t);} //向量数乘
    double abs2(){return sqr(x)+sqr(y);}   //模的平方
    double abs(){return sqrt(abs2());}         //模
    vec rot(double th)   //逆时针旋转，精度很差
    {return vec(o, -sin(th)*y+cos(th)*x, cos(th)*y+sin(th)*x );}
    vec rot(double sinth, double costh)  //逆时针旋转，精度取决于sinth和costh
    {return vec(o, -sinth*y+costh*x, costh*y+sinth*x );}
    vec trunc(double l){return (*this)*(l/abs());}    //把向量v的长度改为l，方向不变
};
point operator+(point p, vec v){return point(p.x+v.x, p.y+v.y);} //点沿着向量方向平移
vec operator-(point p1, point p2){return vec(p2,p1.x-p2.x,p1.y-p2.y);}//用两点构造向量
bool on(point p, vec v){return sign((p-v.o)^v)==0;}   //点在直线上
point projection(point p, vec v){return v.o + v*((p-v.o)*v/v.abs2());}  //求投影点
struct circle
{
    point o; double r;
    circle(point o, double r):o(o),r(r){}
    circle(){};
    int judge(point p){return sign(sqr(r)-(o-p).abs2());}//点和圆的位置关系,1内,0上,-1外
    int judge(vec v){return judge(projection(o,v));}//直线与圆的位置关系,1交,0切,-1离
    double area(){return pi*sqr(r);}
    tuple<point,point> tangent_points(point p)    //过圆外一点求圆的切点
    {
        auto d = (o-p).abs(), x = sqr(r)/d, h = sqrt(sqr(r)-sqr(x));
        auto v = p-o;
        auto pp = o+v.trunc(x);
        return tuple<point,point>(pp+v.trunc(h).rot(1,0), pp+v.trunc(h).rot(-1,0));
    }
    tuple<int,point,point> cross_points(vec v)  //返回(交点个数,第一个点，第二个点)
    {
        auto type=judge(v);
        if(type==-1)return make_tuple(0,_,_);
        auto p=projection(o,v);
        if(type==0)return make_tuple(1,p,_);
        auto op=p-o;
        auto d=op.abs(), x=sqrt(sqr(r)-op.abs2());
        return make_tuple(2,p+v.trunc(x),p+v.trunc(-x));
    }
};
struct Inversion
{
    circle O;
    point inv(point p)//点的反演（不过反演圆圆心）
    {
        auto v=p-O.o;
        v=v*(sqr(O.r)/v.abs2());
        return O.o+v;
    }
    vec inv_on(vec l){return l;}  //直线反演（过反演圆圆心）
    circle inv_off(vec v)         //直线反演（不过反演圆圆心）
    {
        auto p=inv(projection(O.o,v));
        auto op=(p-O.o)*0.5;
        auto cc=circle(O.o+op,op.abs());
        return circle(O.o+op,op.abs());
    }
    vec inv_on(circle c)      //圆反演（过反演圆圆心）
    {
        auto p = inv( O.o + (c.o-O.o)*2 );
        auto v = p-O.o;
        return vec(O.o+v,v.y,-v.x);
    }
    circle inv_off(circle c)  //圆反演（不过反演圆圆心）
    {
        auto oc=c.o-O.o, v=oc*(1.0/oc.abs());
        auto a=c.o+v*(-c.r), b=c.o+v*c.r, aa=inv(a), bb=inv(b);
        return circle( point( (aa.x+bb.x)/2, (aa.y+bb.y)/2 ), (aa-bb).abs()/2 );
    }
}Inv;
int main()
{
    ll T=read(), N;
    double ans;
    Inv.O = circle(_,500);
    while(T--)
    {
        circle c1, c2;
        c1.o.x=c1.r=read(), c1.o.y=0;
        c2.o.x=c2.r=read(), c2.o.y=0;
        if(c1.r>c2.r)swap(c1,c2);
        auto l1 = Inv.inv_on(c1), l2 = Inv.inv_on(c2);
        auto r = (l1.o.x-l2.o.x)/2;
        auto now = circle( point(l2.o.x+r,0) , r );
        ans=0;
        N=read();
        while(N--)
        {
            auto c = Inv.inv_off(now);
            // de(now.o.x), de(now.o.y), de(now.r);
            // de(c.o.x), de(c.o.y), de(c.r);
            if(c.area()<1e-12)break;
            ans += c.area();
            if(sign(now.o.y)==1)
            {
                now.o.y*=-1;
            }
            else
            {
                now.o.y*=-1;
                now.o.y+=2*r;
            }
        }
        printf("%.5lf\n",ans);
    }
    return 0;
}