大白二维计算几何模板

</pre><pre code_snippet_id="275812" snippet_file_name="blog_20140404_1_7192597" name="code" class="cpp">#include<iostream>#include<cstdio>#include<cstring>#include<cmath>#include<algorithm>#include<set>#include<map>#include<queue>#include<vector>#include<cstdlib>#include<stack>using namespace std;#define inf 0x3f3f3f3f#define eps 1e-7#define LL long long#define ULL unsigned long long#define MP make_pair#define pb push_back#define ls ( i << 1 )#define rs ( ls | 1 )#define md ( ( ll[i] + rr[i] ) >> 1 )#define mxn 5020#define PI acos( -1.0 )// 角度转换为弧度double torad( double deg ) {return deg / 180 * PI;}int dcmp( double x ) {  //OKif( fabs( x ) < eps )return 0;return x < 0 ? -1: 1;}struct point {double x, y;point( double x = 0, double y = 0 ): x( x ), y( y ) {}point operator + ( const point &b ) const {return point( x + b.x, y + b.y );}point operator - ( const point &b ) const {return point( x - b.x, y - b.y );}point operator * ( const double &k ) const {return point( x * k, y * k );}point operator / ( const double &k ) const {return point( x / k, y / k );}bool operator < ( const point &b ) const {if( dcmp( x - b.x ) != 0 )return x < b.x;return y < b.y;}bool operator == ( const point &b ) const {return dcmp( x - b.x ) == 0 && dcmp( y - b.y ) == 0;}double len() {return sqrt( x * x + y * y );}};struct Line {point v, p;double ang; //极角Line() {}Line( point v, point p ): v( v ), p( p ) {ang = atan2( v.y, v.x );}point get_point( double t ) {return v * t + p;}bool operator < ( const Line &L ) const {return ang < L.ang;}};// 圆struct circle {point c;double r;circle(){}circle( point c, double r ): c( c ), r( r ){}//圆心角为a的点point get_point( double a ) {return point( c.x + cos( a ) * r, c.y + sin( a ) * r );}};double cross( point A, point B ) {  // 叉乘 OKreturn A.x * B.y - A.y * B.x;}double dot( point A, point B ) { //OKreturn A.x * B.x + A.y * B.y;}// 极角double angle( point v ) {   //OKreturn atan2( v.y, v.x );}double angle( point A, point B ) {  //向量A和向量B的夹角return acos( dot( A, B ) / A.len() / B.len() );}double area( point A, point B, point C ) { // 三角形面积return fabs( cross( A - B, A - C ) ) / 2;}point rotate( point A, double rad ) {  // 向量A旋转rad弧度return point( A.x * cos( rad ) - A.y * sin( rad ), A.x * sin( rad ) + A.y * cos( rad ) );}point Normal( point A ) {// 单位法线double L = A.len();return point( - A.y / L, A.x / L );}// 直线的交点,P 和 Q是直线上的点， v和w是方向向量point GetLineIntersection( point P, point v, point Q, point w ) {point u = P - Q;double t = cross( w, u ) / cross( v, w );return P + v * t;}// 点到直线的距离double DisToLine( point P, point A, point B ) {  //OKpoint v1 = B - A, v2 = P - A;return fabs( cross( v1, v2 ) / v1.len() );}// 点到线段的距离double DisToSeg( point P, point A, point B ) {if( A == B )return ( P - A ).len();point v1 = B - A, v2 = P - A, v3 = P - B;if( dcmp( dot( v1, v2 ) ) < 0 )return v2.len();if( dcmp( dot( v1, v3 ) ) > 0 )return v3.len();return fabs( cross( v1, v2 ) / v1.len() );}// 点在直线的投影point GetLinePro( point P, point A, point B ) {point v = B - A;return A + v * ( dot( v, P - A ) / dot( v, v ) );}// 线段规范相交bool SegProInters( point a1, point a2, point b1, point b2 ) {double c1 = cross( a2 - a1, b1 - a1 ), c2 = cross( a2 - a1, b2 - a1 ),c3 = cross( b2 - b1, a1 - b1 ), c4 = cross( b2 - b1, a2 - b1 );return dcmp( c1 ) * dcmp( c2 ) < 0 && dcmp( c3 ) * dcmp( c4 ) < 0;}// 点在线段上bool Ons( point P, point a1, point a2 ) {return dcmp( cross( a1 - P, a2 - P ) ) == 0 && dcmp( dot( a1 - P, a2 - P ) ) < 0;}//多边形有向面积double PolygonArea( point *p, int n ) {double area = 0;for( int i = 1; i < n - 1; ++i ) area += cross( p[i] - p[0], p[i+1] - p[0] );return area / 2;}int GetLineCirIns( Line L, circle C, double &t1, double &t2, vector<point> &sol ) {double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;double e = a * a + c * c, f = 2 * ( a * b + c * d ), g = b * b + d * d - C.r * C.r;double delta = f * f - 4 * e * g;if( dcmp( delta ) < 0 )return 0;if( dcmp( delta ) == 0 ) {t1 = t2 = - f / ( 2 * e );sol.pb( L.get_point( t1 ) );return 1;}t1 = ( - f - sqrt( delta ) ) / ( 2 * e );sol.push_back( L.get_point( t1 ) );t2 = ( - f + sqrt( delta ) ) / ( 2 * e );sol.push_back( L.get_point( t2 ) );return 2;}// 两圆相交int getCCins( circle c1, circle c2, vector<point> &sol ) {double d = ( c1.c - c2.c ) .len();if( dcmp( d ) == 0 ) {if( dcmp( c1.r - c2.r ) == 0 )return -1;return 0;}if( dcmp( c1.r + c2.r - d ) < 0 )return 0;if( dcmp( fabs( c1.r - c2.r ) - d ) > 0 )return 0;double a = angle( c2.c - c1.c );double da = acos( ( c1.r * c1.r + d * d - c2.r * c2.r ) / ( 2 * c2.r * d ) );point p1 = c1.get_point( a - da ), p2 = c1.get_point( a + da );sol.push_back( p1 );if( p1 == p2 )return 1;sol.push_back( p2 );return 2;}// 过点p到圆c的切线，v[i]是第i条切线。返回切线条数int getTangents( point p, circle c, point *v ) {point u = c.c - p;double dist = u.len();if( dist < c.r )return 0;if( dcmp( dist - c.r ) == 0 ) {v[0] = rotate( u, PI / 2 );return 1;}double ang = asin( c.r / dist );v[0] = rotate( u, - ang );v[1] = rotate( u, + ang );return 2;}// 求两圆的切线，返回切线条数，-1表示无穷条，// a[i] 和 b[i]分别是第i条切线在圆A和圆B上的切点int getTangents( circle A, circle B, point *a, point *b ) {int cnt = 0;if( A.r < B.r )swap( a, b ), swap( A, B );double d2 = ( A.c.x - B.c.x ) * ( A.c.x - B.c.x ) + ( A.c.y - B.c.y ) * ( A.c.y - B.c.y );double rdiff = A.r - B.r;double rsum = A.r + B.r;if( dcmp( rdiff * rdiff - d2 ) > 0 )return 0;double base = atan2( B.c.y - A.c.y, B.c.x - A.c.x );if( dcmp( d2 ) == 0 && dcmp( A.r - B.r ) == 0 )return -1;if( dcmp( d2 - rdiff * rdiff ) == 0 ) {a[cnt] = A.get_point( base );b[cnt] = B.get_point( base );cnt++;return 1;}double ang = acos( ( A.r - B.r ) / sqrt( d2 ) );a[cnt] = A.get_point( base + ang );b[cnt++] = B.get_point( base + ang );a[cnt] = A.get_point( base - ang );b[cnt++] = B.get_point( base - ang );if( dcmp( d2 - rsum * rsum ) == 0 ) {a[cnt] = A.get_point( base );b[cnt++] = B.get_point( PI + base );}else {if( dcmp( d2 - rsum * rsum ) > 0 ) {double ang = acos( ( A.r + B.r ) / sqrt( d2 ) );a[cnt] = A.get_point( base + ang );b[cnt++] = B.get_point( PI + base + ang );a[cnt] = A.get_point( base - ang );b[cnt++] = B.get_point( PI + base - ang );}}return cnt;}// 点在多边形内判定, -1 在边界，0,不在，1在// 点在凸多边形内的判定, 判断点是不是在所有边的左边int inPoly( point p, point * poly, int n ) {int wn = 0;for( int i = 0; i < n; ++i ) {if( Ons( p, poly[i], poly[(i+1)%n] ) )return -1;int k = dcmp( cross( poly[(i+1)%n] - poly[i], p - poly[i] ) );int d1 = dcmp( poly[i].y - p.y );int d2 = dcmp( poly[(i+1)%n].y - p.y );if( k > 0 && d1 <= 0 && d2 > 0 )wn++;if( k < 0 && d2 <= 0 && d1 > 0 )wn--;}if( wn )return 1;return 0;}// 凸包int tb( point *p, int n, point *ch ) {   //OKsort( p, p + n );int m = 0;for( int i = 0; i < n; ++i ) {while( m > 1 && cross( ch[m-1] - ch[m-2], p[i] - ch[m-2] ) <= 0 )m--;ch[m++] = p[i];}int k = m;for( int i = n - 2; i >= 0; --i ) {while( m > k && cross( ch[m-1] - ch[m-2], p[i] - ch[m-2] ) <= 0 )m--;ch[m++] = p[i];}if( n > 1 )m--;return m;}// 点在有向线段左边bool OnLeft( Line L, point p ) {   //OKreturn dcmp( cross( L.v, p - L.p ) ) > 0;}// 二直线交点，假设交点唯一存在point GetIntersection( Line a, Line b ) {point u = a.p - b.p;double t = cross( b.v, u ) / cross( a.v, b.v );return a.p + a.v * t;}// 半平面交int Halfins( Line *L, int n, point *poly ) { //OKsort( L, L + n );int first, last;point *p = new point[n];Line *q = new Line[n];q[first=last=0] = L[0];for( int i = 1; i < n; ++i ) {while( first < last && !OnLeft( L[i], p[last-1] ) ) last --;while( first < last && !OnLeft( L[i], p[first] ) )first++;q[++last] = L[i];if( fabs( cross( q[last].v, q[last-1].v ) ) < eps ) {last--;if( OnLeft( q[last], L[i].p ) )q[last] = L[i];}if( first < last )p[last-1] = GetIntersection( q[last-1], q[last] );}while( first < last && !OnLeft( q[first], p[last-1] ) )last--;if( last - first <= 1 )return 0;p[last] = GetIntersection( q[last], q[first] );int m = 0;for( int i = first; i <= last; ++i )poly[m++] = p[i];return m;}//*/