[BZOJ4445]-[Scoi2015]小凸想跑步-计算几何

说在前面

由于me的推导方向错误
导致这么一道简单的半平面交题,被me写成了模板大全
“谨”以此文,来记录自己的智障


题目

BZOJ4445传送门
非权限题,看题可进传送门


解法

因为要求三角形「P,0,1」的面积最小,也就是小于等于所有的三角形,那么不妨对每一条边考虑
把边0,1设为A,把凸包上的每一条边都当作底边B。底边确定之后,如果可以找到一条线,使「这条线上的每一个点」与A和B形成的三角形面积都相等,那么以这条线为界,A所在的一侧就是可行区域。对所有可行区域取个并集即可

然后怎么求出这条线呢,正常的做法是先设出点的坐标(x,y),然后用叉积写出面积,不等式化一下,就得到形如 Ax+By+C0 A x + B y + C ≤ 0 的式子,然后就可以上半平面交了…

然后me的做法是:
这里写图片描述
要找出一条线,让这条线上的点与A和B构成的三角形面积相等(这里先不考虑A,B平行的情况),那么A,B的延长线交点P肯定在这个线上(然而P的坐标极限可以达到1e18级别,强行掉了一波精度)。那么已经确定了一个点,如果能再确定斜率k或者另一个点Q,那么这条直线就出来了
这时注意到 AhA=BhB A ⋅ h A = B ⋅ h B ,并且还有公共边PQ,所以可以得到 sinαsinβ=hAhB=BA sin ⁡ α sin ⁡ β = h A h B = B A ,而且 α+β α + β 的角度可以算(就是用PA的幅角减去PB的幅角),那么现在只要解出这个方程,就可以获得 α α β β 的角度,从而可以通过将PA旋转,来获得PQ。(获得幅角使用了atan2,旋转的时候需要用到sin和cos,掉精度*3)

me发现me貌似并不会解这个方程,不过这并不能阻挡me前进的步伐。注意到,当 α+β α + β 不是 180 180 ∘ 的时候,较大的角的 sin s i n 值始终比较小的角要大一些,且随着角度的变大,这个比值也会不断加大(可以画两个 sin s i n 图像感受一下,较小的角降落的快一些),也就是说, sinαsinβ sin ⁡ α sin ⁡ β 具有单调性,所以 α α 的角度可以使用实数二分得出(精度已经掉没了)

那么这条线就被成功的找了出来!

别急,万一A,B平行怎么办?显然,当A,B平行的时候,PQ也平行A和B,并且与A和B的距离的比值就是 hA h A hB h B 。那么只需要算出AB的距离,然后再把A或者B平移就可以了

然后me就成功解决掉了此题

其实思维过程很煞笔,关键是写起来恶心
比如说,因为半平面固定保留直线的左半部分,所以还要注意线的方向问题,需要分类讨论

模拟考,考场上用了4个小时干死这题,写完发现死活过不去自己造的大样例,输出答案都是1e10级别的数。me猜测这是由于精度问题,于是输出的时候加了特判(ans>1||ans<=0)?0:ans。最后得到了60分(半平面交忘了判断直线平行)

最后的最后,在ZJC的帮助下,强行用long double踩过= =


下面是大常数,大长度,还自带掉精度的代码

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;

# define double long double

int N ;
const double PI = acos ( -1 ) , eps = 1e-6 ;
struct Vector{
    double x , y ;
    double len(){ return sqrt( x * x + y * y ) ; }
} p[100005] ;
typedef Vector Point ;
struct Line{
    Point p1 , p2 ;
    Vector v ;
}L[100005] , pl[100005] ;
double Lang[100005] ;
Vector operator + ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x + B.x , A.y + B.y } ; }
Vector operator - ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x - B.x , A.y - B.y } ; }
Vector operator * ( const Vector &A , const double &p ){ return ( Vector ){ A.x * p , A.y * p } ; }
Vector operator / ( const Vector &A , const double &p ){ return ( Vector ){ A.x / p , A.y / p } ; }
double cross( const Vector &A , const Vector &B ){ return A.x * B.y - A.y * B.x ; }
double dot( const Vector &A , const Vector &B ){ return A.x * B.x + A.y * B.y ; }
double fix( double x ){ return x < 0 ? x + 2*PI : x ; }
int dcmp( double x ){
    if( x < eps && x > -eps ) return 0 ;
    return x > 0 ? 1 : -1 ;
}

Point glt( const Point &u , const Vector &P , const Point &v , const Vector &Q ){
    double t1 = cross( Q , v - u ) / cross( Q , P ) ; 
    return u + P * t1 ;
}
Point glt( const Line &A , const Line &B ){ return glt( A.p1 , A.v , B.p1 , B.v ) ; }

Vector Rotate( Vector v , double B ){
    double cosB = cos( B ) , sinB = sin( B ) ;
    return ( Vector ){ v.x * cosB - v.y * sinB , v.y * cosB + v.x * sinB } ;
}

Line Rotate( Line x , double ang ){
    Vector v = Rotate( x.v , ang ) ;
    return ( Line ){ x.p1 , x.p1 + v , v } ;
}

bool Onleft( const Line &u , const Point &P ){ return dcmp( cross( u.p1 - P , u.p2 - P ) ) <= 0 ; }

double getAng( double x , double ang ){
    double lf = 0 , rg = ang , mid ;
    while( rg - lf > eps ){
        mid = ( lf + rg ) / 2 ;
        if( sin( mid ) / sin( ang - mid ) > x ) rg = mid ;
        else lf = mid ;
    } return ( lf + rg ) / 2 ;
}

void preWork(){
    Line o = L[1] , tmp ;
    double Oang = atan2( o.v.y , o.v.x ) , olen = o.v.len() ;
    pl[1] = o ;
    for( int i = 2 ; i <= N ; i ++ ){
        if( dcmp( cross( o.v , L[i].v ) ) == 0 ){
            Line chuizhi = Rotate( L[i] , PI / 2 ) ;
            Point pp = glt( chuizhi , o ) ;
            double dist = ( pp - L[i].p1 ).len() , B = L[i].v.len() , cc = ( dist * olen / ( B + olen ) ) ;
            Vector fa = chuizhi.v / chuizhi.v.len() ;
            pl[i] = ( Line ){ L[i].p1 + fa * cc , L[i].p2 + fa * cc , L[i].v } ;
        } else {
            Point pp = glt( o , L[i] ) ;
            if( dot( pp - o.p2 , o.p1 - o.p2 ) <= 0 ){
                double ang = PI - fix( ( atan2( L[i].v.y , L[i].v.x ) - Oang ) ) ;
                double rot = getAng( L[i].v.len() / olen , ang ) ;
                tmp = ( Line ){ pp , o.p1 , o.p1 - pp } ;
                pl[i] = Rotate( tmp , 2 * PI - rot ) ;
            } else{
                double ang = PI - fix( Oang - atan2( L[i].v.y , L[i].v.x ) ) ;
                double rot = getAng( L[i].v.len() / olen , ang ) ;
                tmp = ( Line ){ pp , o.p2 , o.p2 - pp } ;
                tmp = Rotate( tmp , rot ) ;
                pl[i] = ( Line ){ tmp.p2 , tmp.p1 , tmp.p1 - tmp.p2 } ;
            }
        }
    }
}

double getSiz( Point *cvx , int siz ){
    double rt = 0 ;
    for( int i = 2 ; i <= siz ; i ++ )
        rt += cross( cvx[i-1] , cvx[i] ) ;
    rt += cross( cvx[siz] , cvx[1] ) ;
    rt /= 2 ;
    return rt > 0 ? rt : -rt ;
}

bool cmp( const int &A , const int &B ){ return Lang[A] < Lang[B] ; }

Point cvx[100005] ;
int id[100005] , que[100005] , fr , ba , tot ;
void solve(){
    for( int i = 1 ; i <= N ; i ++ ){
        id[i] = i ; Lang[i] = fix( atan2( pl[i].v.y , pl[i].v.x ) ) ;
    } sort( id + 1 , id + N + 1 , cmp ) ;
    fr = 1 , ba = 0 ;
    for( int iii = 1 ; iii <= N ; iii ++ ){
        int i = id[iii] ;
        while( ba > fr && Onleft( pl[i] , cvx[ba-1] ) ) ba -- ;
        while( ba > fr && Onleft( pl[i] , cvx[fr] ) ) fr ++ ;
        que[++ba] = i ;
        if( ba > fr && dcmp( cross( pl[ que[ba] ].v , pl[ que[ba-1] ].v ) ) == 0 ){
            ba -- ;
            if( Onleft( pl[i] , cvx[ba] ) ) que[ba] = i ;
        } if( ba > fr ) cvx[ba-1] = glt( pl[ que[ba] ] , pl[ que[ba-1] ] ) ;
    }
    while( ba > fr && Onleft( pl[ que[fr] ] , cvx[ba-1] ) ) ba -- ;
    cvx[ba] = glt( pl[ que[ba] ] , pl[ que[fr] ] ) ;

    for( int i = fr ; i <= ba ; i ++ ) cvx[++tot] = cvx[i] ;
    double ans = getSiz( cvx , tot ) / getSiz( p , N ) ;
    printf( "%.4f" , ( float ) ans ) ;
}

int main(){
    scanf( "%d" , &N ) ;
    for( int i = 1 ; i <= N ; i ++ )  {
        int x, y ;
        scanf ( "%d%d", & x, & y ) ;
        p [i] = ( Point )  { ( double ) x, ( double ) y } ;
    }
    for( int i = 1 ; i < N ; i ++ )
    L[i] = ( Line ){ p[i] , p[i+1] , p[i+1] - p[i] } ;
    L[N] = ( Line ){ p[N] , p[1] , p[1] - p[N] } ;
    preWork() ; solve() ;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值