罗书的计算几何：
Lifting the Stone
Lifting the Stone：
求多边形的重心。
- 质量集中在顶点上。n个顶点坐标为(xi, yi), 质量为mi，则重心坐标为：
X = ∑( xi×mi ) / ∑mi
Y = ∑( yi×mi ) / ∑mi
特殊地，若每个点的质量相同，则
X = ∑xi / n
Y = ∑yi / n - 质量分布均匀。这个题就是如此，特殊的，质量均匀的三角形重心：
X = ( x0 + x1 + x2 ) / 3
Y = ( y0 + y1 + y2 ) / 3
解析：
- 以一个点为顶点，做三角形，然后求出每个三角形的重心与质量。
- 然后用每个三角形再重新构造一个多边形。
- 这个多边形就是第一种情况，就可以求解了。
(由于均匀分布，所以质量等同于面积)
用到的知识点：求多边形的重心，质量均匀用面积来替代，所以我们通过每两个向量把它们分成多个三角形，然后求出每个三角形的重心，然后开始求当下多边形的重心，来替代多边形的重心。
#include <bits/stdc++.h>
using namespace std;
const int N = 1e6 + 10;
struct point{
double x, y;
}pnt[N];
point bcenter(int n)
{
point p, s;
double tp, area = 0, tpx = 0, tpy = 0;
p.x = pnt[0].x, p.y = pnt[0].y;
for(int i = 1; i <= n; i ++ )
{
s.x = pnt[i == n ? 0 : i].x;
s.y = pnt[i == n ? 0 : i].y;
tp = (s.y * p.x - s.x * p.y);//叉积，求三角形的面积
area += tp / 2.0;//除以二是因为叉积求出来的是四边形的面积
tpx += (p.x + s.x) * tp;//上述公式求和
tpy += (p.y + s.y) * tp;//同理
p.x = s.x;//两个向量来求面积
p.y = s.y;
}
s.x = tpx / (6 * area);//除以6是因为重心要除以3，然后加了两边的.x所以除以2
s.y = tpy / (6 * area);
return s;
}
int main()
{
ios::sync_with_stdio(false);
int t, n;
point z;
cin >> t;
while(t -- )
{
cin >> n;
for(int i = 0; i < n; i ++ )
{
cin >> pnt[i].x >> pnt[i].y;
}
z = bcenter(n);
printf("%.2lf %.2lf\n", z.x, z.y);
}
return 0;
}
Surround the Trees:
Surround the Trees
分析：这道题是通过一个绳子围城一颗树，这就是典型的凸包问题。
题解：先通过存入两个点，然后加入第三个点，然后进行check，如果当前i点是逆时针的点，那么就删除掉len - 1的点，然后进行下一次check，通过凸包的性质，进行顺时针把点加入即可，正着来一遍反着来一遍。
#include <bits/stdc++.h>
using namespace std;
const int N = 110;
struct node{
double x, y;
}edge[N], ans[N];
bool cmp(node a, node b)
{
if(a.x != b.x) return a.x < b.x;
return a.y < b.y;
}
double mul(node a, node b, node c)
{
return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
}
double dis(node a, node b)
{
return sqrt((a.x*1.0-b.x)*(a.x*1.0-b.x) + (a.y*1.0-b.y)*(a.y*1.0-b.y));
}
int solve(int n)
{
int len = 0;
for(int i = 0; i < n; i ++ )
{
//判断凸包问题，我们所要的是逆时针旋转的每一个点，而顺时针的点不需要已经包含在内了。
while(len > 1 && mul(ans[len - 2], ans[len - 1], edge[i]) <= 0) len --;
ans[len ++ ] = edge[i];
}
int k = len;
for(int i = n - 1; i >= 0; i -- )
{
while(len > k && mul(ans[len - 2], ans[len - 1], edge[i]) <= 0) len --;
ans[len ++ ] = edge[i];
}
return len - 1;
}
int main()
{
int t, n;
while(~scanf("%d", &n) && n)
{
for(int i = 0; i < n; i ++ ) scanf("%lf%lf", &edge[i].x, &edge[i].y);
sort(edge, edge + n, cmp);
int len = solve(n);
double ans1 = 0;
for(int i = 0; i < len; i ++ )
{
ans1 += dis(ans[i], ans[i + 1]);
}
if(len == 2)
ans1 /= 2;
printf("%.2lf\n", ans1);
}
}
Segment set:
Segment set
题意：求某一条线段与其他已知线段相交的线段个数。
这道题用到了线段并查集以及判断线段相交的性质
- 并查集就用序号代替
- 判断线段相交需要用到叉积，用当前线段与另一条线段的两点进行叉积，之后的乘积如果为负值，则说明两条线段相交。
#include<cstdio>
#include<cmath>
#include<cstring>
#define lk (k<<1)
#define rk (k<<|1)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const double pi=acos(-1);
const double eps=1e-8;
const int inf=0x3f3f3f3f;
const int N=1e3+10;
int fa[N],cnt,num[N];
int sgn(double x)
{
if(fabs(x)<eps) return 0;
else return x<0?-1:1;
}
int find(int k)
{
if(fa[k]!=k) fa[k]=find(fa[k]);
return fa[k];
}
struct point
{
double x,y;
point(double xx=0,double yy=0):x(xx),y(yy){}
point operator + (point a){return point(x+a.x,y+a.y);}
point operator - (point a){return point(x-a.x,y-a.y);}
point operator * (double k){return point(x*k,y*k);}
point operator / (double k){return point(x/k,y/k);}
};
double cross(point a,point b)
{
return a.x*b.y-a.y*b.x;
}
struct line
{
point p1,p2;
line(){}
line(point p11,point p22):p1(p11),p2(p22){}
};
bool cross_seg(point a,point b,point c,point d)
{
//这个应用叉积的正负来确定，两个线段是否相交
//把当前线段定下来，然后去另一个线段的两点，
//然后计算两点与当前线段的叉积的乘积为负即相交
double c1=cross(c-a,b-a),c2=cross(d-a,b-a);
double d1=cross(d-c,b-c),d2=cross(d-c,a-c);
if(sgn(c1)*sgn(c2)<=0&&sgn(d1)*sgn(d2)<=0) return 1;
return 0;
}
line l[N];
int main()
{
int t;
scanf("%d",&t);
while(t--){
memset(num,0,sizeof(num));
int q;
cnt=0;
scanf("%d",&q);
for(int i=1;i<=q;i++) fa[i]=i;
while(q--){
char s[3];
double x1,y1,x2,y2;
scanf("%s",s);
if(s[0]=='P'){
scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
l[++cnt]=line(point(x1,y1),point(x2,y2));
num[cnt]++;
for(int i=1;i<cnt;i++){
//i是父节点
if(cross_seg(l[i].p1,l[i].p2,l[cnt].p1,l[cnt].p2)){
int r1=find(i);
int r2=find(cnt);
if(r1!=r2){
fa[r2]=r1;
num[r1]+=num[r2];
}
}
}
}
else {
int k;
scanf("%d",&k);
int r=find(k);
printf("%d\n",num[r]);
}
}
if(t) printf("\n");
}
return 0;
}
Problem G. Interstellar Travel
Problem G. Interstellar Travel：
从(0, 0) —> (n, 0)，只要到了(n, 0)就行，所以目的是代价和最小，有时为负值。
解析：
- 第一点：排序规则，由题意以及代价最小得，x从小到大，y从大到小，id从小到大。
- 第二点：检验用check函数，要id较大的，叉积为负得。(利用凸包定理进行检验)
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 2e5 + 10;
typedef struct node{
LL x, y;
int id;
}Node;
Node q[N], ans[N];
bool cmp(node a, node b)
{
if(a.x == b.x)
{
if(a.y == b.y)
return a.id < b.id;
return a.y > b.y;//保证代价和最小
}
return a.x < b.x;
}
//check条件1：要id比我大的
//check条件2：要叉积为负的
bool check(node a, node b, node c)
{
if((b.x - a.x) * (c.y - a.y) == (b.y - a.y) * (c.x - a.x)) return c.id < b.id;
return ((b.x - a.x) * (c.y - a.y) > (b.y - a.y) * (c.x - a.x));
}
int main()
{
ios::sync_with_stdio(false);
int t;
cin >> t;
while(t -- )
{
int n;
cin >> n;
for(int i = 0; i < n; i ++ )
{
cin >> q[i].x >> q[i].y;
q[i].id = i + 1;
}
sort(q, q + n, cmp);
int k = 0;
for(int i = 0; i < n; i ++ )
{
//同类一个就行，并且要y大的。
if(i && q[i].x == q[i - 1].x && q[i].y == q[i - 1].y) continue;
while(k > 1 && check(ans[k - 2], ans[k - 1], q[i])) k -- ;
ans[k ++ ] = q[i];
}
for(int i = 0; i < k; i ++ )
if(i != k - 1) cout << ans[i].id << ' ';
else cout << ans[i].id << endl;
}
return 0;
}
Quoit Design：
最近点对问题。
第一步是划分，先以x从小到大排序，然后进行划分成s1和s2两个集合进行递归，递归剩余一个点或两个点的时候结束。
第二步，合并，如果最小点对在s1或者s2集合内部，那么直接合并即可，如果在s1和s2之间，也就是下面的图，我们用mid求他的中间值，这时用y来排序。(这里是剪枝)
//ECUST luoyongjun
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8;
const int MAXN = 100010;
const double INF = 1e20;
int sgn(double x){
if(fabs(x) < eps) return 0;
else return x<0?-1:1;
}
struct Point{
double x,y;
};
double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);}
bool cmpxy(Point A,Point B){ //排序：先对横坐标x排序，再对y排序
return sgn(A.x-B.x)<0 || (sgn(A.x-B.x)==0 && sgn(A.y-B.y)<0);
}
bool cmpy (Point A,Point B){return sgn(A.y-B.y)<0;}//只对y坐标排序
Point p[MAXN],tmp_p[MAXN];
double Closest_Pair(int left,int right){
double dis = INF;
if(left == right) return dis; //只剩1个点
if(left + 1 == right) //只剩2个点
return Distance(p[left], p[right]);
int mid = (left+right)/2; //分治
double d1 = Closest_Pair(left,mid); //求s1内的最近点对
double d2 = Closest_Pair(mid+1,right); //求s2内的最近点对
dis = min(d1,d2);
int k = 0;
for(int i=left;i<=right;i++) //在s1和s2中间附近找可能的最小点对
if(fabs(p[mid].x - p[i].x) <= dis) //按x坐标来找
tmp_p[k++] = p[i];
sort(tmp_p,tmp_p+k,cmpy);//按y坐标排序，用于剪枝。这里不能按x坐标排序
for(int i=0;i<k;i++)
for(int j=i+1;j<k;j++){
if(tmp_p[j].y - tmp_p[i].y >= dis) break; //剪枝
dis = min(dis,Distance(tmp_p[i],tmp_p[j]));
}
return dis; //返回最小距离
}
int main(){
int n;
while(~scanf("%d",&n) && n){
for(int i=0;i<n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
sort(p,p+n,cmpxy); //先排序
printf("%.2f\n",Closest_Pair(0,n-1)/2); //输出最短距离的一半
}
return 0;
}
Palace:
题意：也是最近点对问题，不过每次去掉一个点。
所以当去点最近点对的两个点时，特殊考虑一下即可。
#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
using namespace std;
typedef long long ll;
const int maxn=1e6;
const ll INF=1e18;
int X,Y;//最小两个点的下标
ll dis; //当前最小距离
struct node{
ll x,y;
int id;
}a[maxn],term[maxn];
bool cmpx(node a,node b)
{
return a.x<b.x;
}
bool cmpy(node a,node b)
{
return a.y<b.y;
}
ll dist(node a,node b)
{
return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
void minpairs(int left,int right,int opt)//opt表示删除点的序号
{
if(left==right)
return ;
if(left+1==right)
{
if(a[left].id!=opt&&a[right].id!=opt)
{
if(dist(a[left],a[right])<dis)
{
dis=dist(a[left],a[right]);
if(opt==-1)
{
X=a[left].id;
Y=a[right].id;
}
}
}
return ;
}
int mid=(left+right)/2;
minpairs(left,mid,opt);
minpairs(mid+1,right,opt);
int i,j,k=0;
for(i=left;i<=right;i++)
{
if(a[i].id!=opt&&abs(a[i].x-a[mid].x)<=dis)
term[k++]=a[i];
}
sort(term,term+k,cmpy);
for(i=0;i<k;i++)
{
for(j=i+1;j<k;j++)
{
if(term[j].y-term[i].y>dis)
break;
if(dist(term[j],term[i])<dis)
{
dis=dist(term[j],term[i]);
if(opt==-1)
{
X=term[i].id;
Y=term[j].id;
}
}
}
}
return ;
}
int main()
{
int t,n;
scanf("%d",&t);
while(t--)
{
ll ans=0;
scanf("%d",&n);
for(int i=0;i<n;i++)
{
scanf("%lld%lld",&a[i].x,&a[i].y);
a[i].id=i;
}
dis=INF;
sort(a,a+n,cmpx);
minpairs(0,n-1,-1);
ans+=dis*(n-2);
dis=INF;
minpairs(0,n-1,X);
ans+=dis;
dis=INF;
minpairs(0,n-1,Y);
ans+=dis;
printf("%lld\n",ans);
}
return 0;
}
最大三角形
//最大三角形的顶点一定是凸包的顶点，先求凸包再三重循环用向量算最大的三角形面积。
- 找到左下角的点
- 通过左下角的点计算角度，arctanα = y / x
- 排序以及找凸包顶点
#include<iostream>
#include<stdio.h>
#include<math.h>
#include<algorithm>
using namespace std;
#define N 50010
struct point
{
int x,y;
double angel;
} p[N],stack[N];
int top,n;
double dis(point a,point b)
{
return sqrt(double((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)));
}
bool mult(point p1,point p2,point p0)
{
return (p1.x-p0.x)*(p2.y-p0.y) >= (p2.x-p0.x)*(p1.y-p0.y);
}
bool cmp(point a,point b)
{
if(a.angel == b.angel)
{
if (a.x == b.x)
return a.y > b.y;
return a.x > b.x;
}
return a.angel < b.angel;
}
void graham()
{
int i,k=0;
for(i=0; i<n; i++)
if(p[i].y < p[k].y||((p[i].y == p[k].y)&&(p[i].x<p[k].x)))
k=i;
swap(p[0],p[k]);
for(i=1; i<n; i++)
p[i].angel=atan2(double(p[i].y-p[0].y),double(p[i].x-p[0].x));
sort(p+1,p+n,cmp);
stack[0]=p[0];
stack[1]=p[1];
stack[2]=p[2];
top=3;
for(i=3; i<n; i++)
{
while(top > 2 && mult(stack[top-2],stack[top-1],p[i])<=0)
top--;
stack[top++]=p[i];
}
}
int main()
{
int i,j,k;
double ans,t,s1,s2,s3,max;
while(scanf("%d",&n)!=EOF)
{
for(i=0; i<n; i++)
scanf("%d%d",&p[i].x,&p[i].y);
graham();
ans=0;
max=0;
for(i=0; i<top-2; i++)
for(j=i+1; j<top-1; j++)
for(k=j+1; k<top; k++)
{
if(mult(stack[i],stack[j],stack[k])==0)
continue;
s1=dis(stack[i],stack[j]);
s2=dis(stack[j],stack[k]);
s3=dis(stack[k],stack[i]);
t=(s1+s2+s3)/2;
ans=t*(t-s1)*(t-s2)*(t-s3);
if(ans>max)
max=ans;
}
printf("%.2lf\n",sqrt(max));
}
return 0;
}
The widest road
题意：
两凸包间的最短距离，这道题需要判断两个凸包是不是相离的，如果是，直接输出，否则相交输出零。
#include<cstdio>
#include<cmath>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<queue>
#include<map>
#include<stack>
#include<set>
#define e exp(1.0); //2.718281828
#define mod 1000000007
#define INF 0x7fffffff
#define inf 0x3f3f3f3f
typedef long long LL;
using namespace std;
#define zero(x) (((x)>0?(x):(-x))<eps)
const double eps=1e-8;
const double pi=acos(-1.0);
//判断数k的符号 -1负数 1正数 0零
int dcmp(double k) {
return k<-eps?-1:k>eps?1:0;
}
inline double sqr(double x) {
return x*x;
}
struct point {
double x,y;
point() {};
point(double a,double b):x(a),y(b) {};
void input() {
scanf("%lf %lf",&x,&y);
}
friend point operator + (const point &a,const point &b) {
return point(a.x+b.x,a.y+b.y);
}
friend point operator - (const point &a,const point &b) {
return point(a.x-b.x,a.y-b.y);
}
friend bool operator == (const point &a,const point &b) {
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
friend point operator * (const point &a,const double &b) {
return point(a.x*b,a.y*b);
}
friend point operator * (const double &a,const point &b) {
return point(a*b.x,a*b.y);
}
friend point operator / (const point &a,const double &b) {
return point(a.x/b,a.y/b);
}
friend bool operator < (const point &a, const point &b) {
return a.x < b.x || (a.x == b.x && a.y < b.y);
}
double norm() {
return sqrt(sqr(x)+sqr(y));
}
};
//计算两个向量的叉积
double cross(const point &a,const point &b) {
return a.x*b.y-a.y*b.x;
}
double cross3(point A,point B,point C) { //叉乘
return (B.x-A.x)*(C.y-A.y)-(B.y-A.y)*(C.x-A.x);
}
//计算两个点的点积
double dot(const point &a,const point &b) {
return a.x*b.x+a.y*b.y;
}
double dot3(point A,point B,point C) { //点乘
return (C.x-A.x)*(B.x-A.x)+(C.y-A.y)*(B.y-A.y);
}
//向量长度
double length(const point &a) {
return sqrt(dot(a,a));
}
//两个向量的角度
double angle(const point &a,const point &b) {
return acos(dot(a,b)/length(a)/length(b));
}
//计算两个点的距离
double dist(const point &a,const point &b) {
return (a-b).norm();
}
//op沿远点逆时针旋转角度A
point rotate_point(const point &p,double A) {
double tx=p.x,ty=p.y;
return point(tx*cos(A)-ty*sin(A),tx*sin(A)+ty*cos(A));
}
double TriArea(const point &a, const point &b, const point &c) {
return fabs( cross( b - a, c - a ) ) / 2;
}
point Normal(const point &a) {
double L = length(a);
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;
}
//求点p到直线ab的距离
double DistanceToLine(point p, point a, point b) {
point v1 = b - a, v2 = p - a;
return fabs(cross(v1,v2)) / length(v1);
}
//求点p到线段ab的距离
double DistanceToSegment(point p, point a, point b) {
if(a==b) return length(p - a);
point v1 = b - a, v2 = p - a, v3 = p - b;
if(dcmp(dot(v1,v2)) < 0) return length(v2);
else if(dcmp(dot(v1,v3)) > 0) return length(v3);
else return fabs(cross(v1,v2)) / length(v1);
}
//判断直线a1a2和直线b1b2是否规范相交
bool SegmentProperIntersection(point a1, point a2, point b1, point b2) {
double c1 = cross(a2-a1,b1-a1), c2 = cross(a2-a1, b2-a1);
double c3 = cross(b2-b1, a1-b1), c4 = cross(b2-b1, a2-b1);
return dcmp(c1) * dcmp(c2) <0 && dcmp(c3) * dcmp(c4) < 0;
}
//判断点p是否在直线a1a2上
bool OnSegment(point p, point a1, point a2) {
return dcmp(cross(a1-p,a2-p)) ==0 && dcmp(dot(a1-p,a2-p))<0;
}
//判断线段a1a2和线段b1b2是否相交，可以在端点处相交
bool SegmentIntersection(point a1, point a2, point b1, point b2) {
return SegmentProperIntersection(a1, a2, b1, b2) || OnSegment(a1, b1, b2) || OnSegment(a2, b1, b2);
}
double SegmentToSegment(point a1, point a2, point b1, point b2) {
//线段间的最短距离分为四种情况
double t1 = DistanceToSegment(b1, a1, a2);
double t2 = DistanceToSegment(b2, a1, a2);
double t3 = DistanceToSegment(a1, b1, b2);
double t4 = DistanceToSegment(a2, b1, b2);
return min(t1,min(t2,min(t3,t4)));
}
//使点集逆时针转
void antiClockSort(point *ch, int n) {
double res = cross(ch[1] - ch[0], ch[2] - ch[0]);
if(dcmp(res) >= 0) return;
reverse(ch, ch+n);
}
int ConvexHull(point* P, int cnt, point* res) {
sort(P, P + cnt);
cnt = unique(P, P + cnt) - P;
int m = 0;
for (int i = 0; i < cnt; i++) {
while (m > 1 && cross(res[m - 1] - res[m - 2], P[i] - res[m - 2]) <= 0)
m--;
res[m++] = P[i];
}
int k = m;
for (int i = cnt - 2; i >= 0; i--) {
while (m > k && cross(res[m - 1] - res[m - 2], P[i] - res[m - 2]) <= 0)
m--;
res[m++] = P[i];
}
if (cnt > 1) m--;
return m;
}
/**********************************************/
//判断点是否在多边形内
int isPointInPolygon(point p, point *a, int n) {
int cnt = 0;
for(int i=0; i<n; ++i) {
if(OnSegment(p, a[i], a[(i+1)%n])) return -1;
int k = cross(a[(i+1)%n]-a[i], p-a[i]);
int d1 = a[i].y - p.y;
int d2 = a[(i+1)].y - p.y;
if(k>0 &&d1<=0 &&d2>0)//点在线段的左侧
cnt++;
if(k<0 &&d2<=0 &&d1>0)//点在线段的右侧
cnt++;
//k==0，点和线段共线的情况不考虑
}
if(cnt&1)return 1;
return 0;
}
//判断凸包是否相离
bool two_getaway_ConvexHull(point *cha, int n1, point *chb, int m1) {
if(n1==1 && m1==1) {
if(cha[0]==chb[0])
return false;
} else if(n1==1 && m1==2) {
if(OnSegment(cha[0], chb[0], chb[1]))
return false;
} else if(n1==2 && m1==1) {
if(OnSegment(chb[0], cha[0], cha[1]))
return false;
} else if(n1==2 && m1==2) {
if(SegmentIntersection(cha[0], cha[1], chb[0], chb[1]))
return false;
} else if(n1==2) {
for(int i=0; i<n1; ++i)
if(isPointInPolygon(cha[i], chb, m1))
return false;
} else if(m1==2) {
for(int i=0; i<m1; ++i)
if(isPointInPolygon(chb[i], cha, n1))
return false;
} else {
for(int i=0; i<n1; ++i) {
for(int j=0; j<m1; ++j) {
if(SegmentIntersection(cha[i], cha[(i+1)%n1], chb[j], chb[(j+1)%m1]))
return false;
}
}
for(int i=0; i<n1; ++i)
if(isPointInPolygon(cha[i], chb, m1))
return false;
for(int i=0; i<m1; ++i)
if(isPointInPolygon(chb[i], cha, n1))
return false;
}
return true;
}
//旋转卡壳求两个凸包最近距离
double solve(point *P, point *Q, int n, int m) {
if(n==1 && m==1) {
return length(P[0] - Q[0]);
} else if(n==1 && m==2) {
return DistanceToSegment(P[0], Q[0], Q[1]);
} else if(n==2 && m==1) {
return DistanceToSegment(Q[0], P[0], P[1]);
} else if(n==2 && m==2) {
return SegmentToSegment(P[0], P[1], Q[0], Q[1]);
}
int yminP = 0, ymaxQ = 0;
for(int i=0; i<n; ++i) if(P[i].y < P[yminP].y) yminP = i;
for(int i=0; i<m; ++i) if(Q[i].y > Q[ymaxQ].y) ymaxQ = i;
P[n] = P[0];
Q[n] = Q[0];
double INF2 = 1e100;
double arg, ans = INF2;
for(int i=0; i<n; ++i) {
//当叉积负正转正时，说明点ymaxQ就是对踵点
while((arg=cross(P[yminP] - P[yminP+1],Q[ymaxQ+1] - Q[ymaxQ])) < -eps)
ymaxQ = (ymaxQ+1)%m;
double ret;
if(arg > eps) { //卡住第二个凸包上的点。
ret = DistanceToSegment(Q[ymaxQ], P[yminP], P[yminP+1]);
ans = min(ans,ret);
} else { //arg==0，卡住第二个凸包的边
ret = SegmentToSegment(P[yminP],P[yminP+1],Q[ymaxQ],Q[ymaxQ+1]);
ans = min(ans,ret);
}
yminP = (yminP+1)%n;
}
return ans;
}
double mindis_twotubao(point *P, point *Q, int n, int m){
//尼玛，hdu2823要判是否分离，poj3608不判
//return min(solve(P, Q, n, m),solve(Q,P,m,n));
//判断凸包是不是相离，如果不是，输出0
if(two_getaway_ConvexHull(P,n,Q,m)==true) return min(solve(P, Q, n, m),solve(Q,P,m,n));
else return 0.0;
}
const int N=10005;
point a[N],b[N];
point cha[N],chb[N];
int main() {
int n,m;
while(scanf("%d%d",&n,&m)!=EOF){
for(int i=0;i<n;++i) scanf("%lf%lf",&a[i].x,&a[i].y);
for(int i=0;i<m;++i) scanf("%lf%lf",&b[i].x,&b[i].y);
//先求凸包
int n1 = ConvexHull(a, n, cha);
int m1 = ConvexHull(b, m, chb);
printf("%.4f\n",mindis_twotubao(cha,chb,n1,m1));
}
return 0;
}
矩形面积：
给出n个矩形，求能覆盖所有矩形的最小的矩形的面积。
题解：对所有点求凸包，然后旋转卡壳，对没一条边求该边的最左最右和最上的三个点。
利用叉积面积求高，利用点积的性质求最左右点和长度，更新面积最小值即可。
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#define MAX 50010
using namespace std;
struct Point{
double x,y;
Point(double x=0,double y=0):x(x),y(y){}
};
Point P[MAX],ch[MAX];
typedef Point Vector;
typedef Point point;
Vector operator - (Point A,Point B)
{
return Vector(A.x-B.x,A.y-B.y);
}
bool operator <(const Point &a,const Point &b)
{
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
const double eps=1e-10;
int dcmp(double x)
{
if(fabs(x)<eps) return 0; else return x<0?-1:1;
}
bool operator ==(const Point &a,const Point &b)
{
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Cross(Vector A,Vector B)
{
return A.x*B.y-A.y*B.x;
}
double dot(Vector A,Vector B)
{
return A.x*B.x+A.y*B.y;
}
int ConvexHull(Point *p,int n)
{
sort(p,p+n);
n=unique(p,p+n)-p;
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;
}
double Length(Vector A)
{
return dot(A,A);
}
double rotating_calipers(Point *p,int n)
{
int l=1,r=1,w;
double ans=1e30;
p[n]=p[0];
for(int i=0;i<n;i++)
{
//注意这里等于0一定要算上
//这里debug了整整一个小时 - -|||||
//找到至高点
while(dcmp(Cross(p[i+1]-p[i],p[w+1]-p[i])-Cross(p[i+1]-p[i],p[w]-p[i]))>=0) //因为边平行的时候面积相等 虽然如此但还是要继续找下一个 横着爬不动的意思
w=(w+1)%n;
//找到最右的点 不可能向左的
while(dcmp(dot(p[i+1]-p[i],p[r+1]-p[i])-dot(p[i+1]-p[i],p[r]-p[i]))>0) //凸包不可能凹进去 所以不需要等号 加深对凸包求解过程的理解
r=(r+1)%n;
if(i==0) l=r;
while(dcmp(dot(p[i+1]-p[i],p[l+1]-p[i])-dot(p[i+1]-p[i],p[l]-p[i]))<=0) //必须加等号 否则凸包遇到直边的时候拐不过来 上不去 第二组样例可以完美说明问题
l=(l+1)%n;
double d=Length(p[i+1]-p[i]);
double area=fabs(Cross(p[i+1]-p[i],p[w]-p[i]))
*fabs(dot(p[i+1]-p[i],p[r]-p[i])-dot(p[i+1]-p[i],p[l]-p[i]))/d;
//cout<<fabs(Cross(p[i+1]-p[i],p[w]-p[i]))<<" ";
//cout<<"rrrrr "<<r<<" lll "<<l<<endl;
//cout<<dot(p[i+1]-p[i],p[r]-p[i])<<" "<<dot(p[i+1]-p[i],p[l]-p[i])<<endl;
ans=min(ans,area);
}
return ans;
}
int main()
{
int t,n,cas=1;
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
n*=4;
for(int i=0;i<n;i++)
scanf("%lf%lf",&P[i].x,&P[i].y);
int m=ConvexHull(P,n);
//for(int i=0;i<n;i++)
//cout<<ch[i].x<<" "<<ch[i].y<<endl;
double ans;
if(m<3) ans=0;
else ans=rotating_calipers(ch,m);
long long tmp = ans+0.5;
printf("Case #%d:\n%lld\n",cas++,tmp);
}
return 0;
}
Run ：
给定n个人的比赛状况，现在已经跑得距离，以及比赛速度，时间无限，问有多少人可以在某一段时间跑到第一名。
以时间为横轴，距离为纵轴，每一个人的跑步状况可以等效于已知起点，速度的一条射线我们可以加一个边界，对这些半平面求交，很明显可以选中的半平面代表的选手才可以跑到第一名。减去认为加的边界2，就可以得到答案。
#include <stdio.h>
#include <iostream>
#include <algorithm>
#include <sstream>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <string>
#include <time.h>
#include <math.h>
#include <queue>
#include <stack>
#include <set>
#include <map>
using namespace std;
#define INF 10000000000400.0
#define eps 1e-8
#define pi acos(-1.0)
typedef long long ll;
int dcmp(double x){
if(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 &a,const Point &b){
return Point(a.x+b.x,a.y+b.y);
}
Point operator - (const Point &a,const Point &b){
return Point(a.x-b.x,a.y-b.y);
}
Point operator * (const Point &a,const double &p){
return Point(a.x*p,a.y*p);
}
Point operator / (const Point &a,const double &p){
return Point(a.x/p,a.y/p);
}
bool operator < (const Point &a,const Point &b){
return a.x<b.x||(a.x==b.x&&a.y<b.y);
}
bool operator == (const Point &a,const Point &b){
return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;
}
double Dot(Point a,Point b){
return a.x*b.x+a.y*b.y;
}
double Length(Point a){
return sqrt(Dot(a,a));
}
double Angle(Point a,Point b){
return acos(Dot(a,b)/Length(a)/Length(b));
}
double angle(Point a){
return atan2(a.y,a.x);
}
double Cross(Point a,Point b){
return a.x*b.y-a.y*b.x;
}
Point vecunit(Point a){
return a/Length(a);
}
Point Normal(Point a){
return Point(-a.y,a.x)/Length(a);
}
Point Rotate(Point a,double rad){
return Point(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
double Area2(Point a,Point b,Point c){
return Length(Cross(b-a,c-a));
}
struct Line{
Point p,v;
double ang;
Line(){};
Line(Point p,Point v):p(p),v(v){
ang=atan2(v.y,v.x);
}
bool operator < (const Line &L) const {
return ang<L.ang;
}
};
bool OnLeft(const Line &L,const Point &p){
return dcmp(Cross(L.v,p-L.p))>0;
}
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;
}
Point GetLineIntersection(Line a,Line b){
return GetLineIntersection(a.p,a.v,b.p,b.v);
}
vector<Point> HPI(vector<Line> L){
int n=L.size();
sort(L.begin(),L.end());//将所有半平面按照极角排序。
/*for(int i=0;i<n;i++){
cout<<"han "<<i<<" ";
printf("%.2lf %.2lf %.2lf %.2lf %.2lf\n",L[i].p.x,L[i].p.y,L[i].v.x,L[i].v.y,L[i].ang);
}*/
int first,last;
vector<Point> p(n);
vector<Line> q(n);
vector<Point> ans;
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]=GetLineIntersection(q[last-1],q[last]);//计算队列顶部半平面交点。
}
while(first<last&&!OnLeft(q[first],p[last-1]))last--;//删除队列顶部的无用半平面。
if(last-first<=1)return ans;//半平面退化
p[last]=GetLineIntersection(q[last],q[first]);//计算队列顶部与首部的交点。
for(int i=first;i<=last;i++)ans.push_back(p[i]);//将队列中的点复制。
return ans;
}
double PolyArea(vector<Point> p){
int n=p.size();
double ans=0;
for(int i=1;i<n-1;i++)
ans+=Cross(p[i]-p[0],p[i+1]-p[0]);
return fabs(ans)/2;
}
int main()
{
//freopen("data.in","r",stdin);
//freopen("data.out","w",stdout);
int T,n;
cin>>T;
while(T--){
cin>>n;
vector<Line> L;
L.push_back(Line(Point(INF,INF),Point(-1,0)));
L.push_back(Line(Point(0,INF),Point(0,-1)));
while(n--){
double a,b;
scanf("%lf%lf",&a,&b);
L.push_back(Line(Point(0,a),Point(1,b)));
}
vector<Point> ans=HPI(L);
// for(int i=0;i<ans.size();i++)cout<<ans[i].x<<" "<<ans[i].y<<endl;
printf("%d\n",ans.size()-2);
}
return 0;
}
Mission Impossible ：
题目意思：有一块半径为r的圆形蛋糕，其中心在原点，有一个人在某一个点（x，y），现把蛋糕按两点所在直线切蛋糕，求切了n次后，这个人所在蛋糕的面积占总面积的百分比。
很容易想到先进行半平面交求出这个人所在位置的区域，再根据这个区域（多边形）求与圆的交，就是就个人得到的蛋糕。
//140MS
#include<iostream>
#include<cstdio>
#include<math.h>
#define eps 1e-8
#define maxn 10000
#define PI acos(-1.0)
struct point{double x,y;};
point p[maxn],save[maxn],temp[maxn];
point points[maxn][2];
point cen,people;
int n,ns,m;
double r;
double xmult(point p1,point p2,point p0)
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p1.y-p0.y)*(p2.x-p0.x);
}
double dmult(point p1,point p2,point p0)
{
return (p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y);
}
double distance(point p1,point p2)
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
void mycopy(point *s,int &ns,point *temp,int n)
{
int i;
ns=n;
for(i=0;i<n;i++)s[i]=temp[i];
}
void getline(point p1,point p2,double &a,double &b,double &c)
{
a=p2.y-p1.y;
b=p1.x-p2.x;
c=p2.x*p1.y-p1.x*p2.y;
}
point intersection(point p1,point p2,point p3,point p4)
{
point ret=p1;
double t=((p1.x-p3.x)*(p3.y-p4.y)-(p1.y-p3.y)*(p3.x-p4.x))
/((p1.x-p2.x)*(p3.y-p4.y)-(p1.y-p2.y)*(p3.x-p4.x));
ret.x+=(p2.x-p1.x)*t;
ret.y+=(p2.y-p1.y)*t;
return ret;
}
void cut(point side)
{
int i,j;
temp[0].x=-10000;temp[0].y=-10000;
temp[1].x=-10000;temp[1].y=10000;
temp[2].x=10000;temp[2].y=10000;
temp[3].x=10000;temp[3].y=-10000;
ns=m=4;
for(i=0;i<n;i++)
{
double a,b,c;
if(xmult(side,points[i][1],points[i][0])>=0)
getline(points[i][0],points[i][1],a,b,c);
else getline(points[i][1],points[i][0],a,b,c);
int cnt=0;
for(j=0;j<ns;j++)
{
if(a* temp[j].x+b* temp[j].y+c>=0)
{
save[cnt++]=temp[j];
}
else
{
point p1=temp[(j-1+ns)%ns],p2=temp[(j+1)%ns];
if(a*p1.x+b*p1.y+c>0)
save[cnt++]=intersection(points[i][0],points[i][1],p1,temp[j]);
if(a*p2.x+b*p2.y+c>0)
save[cnt++]=intersection(points[i][0],points[i][1],p2,temp[j]);
}
}
mycopy(temp,ns,save,cnt);
}
}
double cirtri(point pa,point pb,point po,double r)
{
double a,b,c,x,y;
double area=xmult(pa,pb,po)/2;
a=distance(po,pb);
b=distance(po,pa);
c=distance(pa,pb);
if(a<=r&&b<=r)//1
{
return area;
}
else if(a<r&&b>=r)//2
{
x=(dmult(pa,po,pb)+sqrt(c*c*r*r-xmult(pa,po,pb)*xmult(pa,po,pb)))/c;
return asin(area*(c-x)*2/c/b/r)*r*r/2+area*x/c;
}
else if(a>=r&&b<r)//3
{
y=(dmult(pb,po,pa)+sqrt(c*c*r*r-xmult(pb,po,pa)*xmult(pb,po,pa)))/c;
return asin(area*(c-y)*2/c/a/r)*r*r/2+area*y/c;
}
else if(fabs(2*area)>=r*c||dmult(pb,po,pa)<=0
||dmult(pa,po,pb)<=0)//4
{
if(dmult(pa,pb,po)<0)
{
if(xmult(pa,pb,po)<0)
{
return (-PI-asin(area*2/a/b))*r*r/2;
}
else return (PI-asin(area*2/a/b))*r*r/2;
}
else return asin(area*2/a/b)*r*r/2;
}
else //5
{
x=(dmult(pa,po,pb)+sqrt(c*c*r*r-xmult(pa,po,pb)*xmult(pa,po,pb)))/c;
y=(dmult(pb,po,pa)+sqrt(c*c*r*r-xmult(pb,po,pa)*xmult(pb,po,pa)))/c;
return (asin(area*(1-x/c)*2/r/b)
+asin(area*(1-y/c)*2/r/a))*r*r/2
+area*((y+x)/c-1);
}
}
int main()
{
int cas;
scanf("%d",&cas);
int i,j,k;
for(k=1;k<=cas;k++)
{
scanf("%lf%d",&r,&n);
for(i=0;i<n;i++)
scanf("%lf%lf%lf%lf",&points[i][0].x,&points[i][0].y,
&points[i][1].x,&points[i][1].y);
scanf("%lf%lf",&people.x,&people.y);
cut(people);
double res=0;
cen.x=0;cen.y=0;
for(i=0;i<ns;i++)
res+=cirtri(save[i],save[(i+1)%ns],cen,r);
double ans=fabs(100.0*res)/(PI*r*r);
printf("Case %d: %.5lf%%\n",k,ans);
}
return 0;
}
An Easy Physics Problem：
小球碰撞圆柱体是否会反弹到另一点需要验证。
#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-8; //本题如果设定eps=1e-10，会Wrong Answer
int sgn(double x){ //判断x是否等于0
if(fabs(x) < eps) return 0;
else return x<0?-1:1;
}
struct Point{ //定义点和基本运算
double x,y;
Point(){}
Point(double x,double y):x(x),y(y){}
Point operator + (Point B){return Point(x+B.x,y+B.y);}
Point operator - (Point B){return Point(x-B.x,y-B.y);}
Point operator * (double k){return Point(x*k,y*k);}
Point operator / (double k){return Point(x/k,y/k);}
};
typedef Point Vector; //定义向量
double Dot(Vector A,Vector B){return A.x*B.x + A.y*B.y;} //点积
double Len(Vector A){return sqrt(Dot(A,A));} //向量的长度
double Len2(Vector A){return Dot(A,A);} //向量长度的平方
double Cross(Vector A,Vector B){return A.x*B.y - A.y*B.x;} //叉积
double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);}
struct Line{
Point p1,p2;
Line(){}
Line(Point p1,Point p2):p1(p1),p2(p2){}
};
typedef Line Segment; //定义线段，两端点是p1,p2
int Point_line_relation(Point p,Line v){
int c = sgn(Cross(p-v.p1,v.p2-v.p1));
if(c < 0)return 1; //1：p在v的左边
if(c > 0)return 2; //2：p在v的右边
return 0; //0：p在v上
}
double Dis_point_line(Point p, Line v){ //点到直线的距离
return fabs(Cross(p-v.p1,v.p2-v.p1))/Distance(v.p1,v.p2);
}
//点到线段的距离
double Dis_point_seg(Point p, Segment v){
if(sgn(Dot(p- v.p1,v.p2-v.p1))<0 || sgn(Dot(p- v.p2,v.p1-v.p2))<0)
return min(Distance(p,v.p1),Distance(p,v.p2));
return Dis_point_line(p,v);
}
//点在直线上的投影
Point Point_line_proj(Point p, Line v){
double k=Dot(v.p2-v.p1,p-v.p1)/Len2(v.p2-v.p1);
return v.p1+(v.p2-v.p1)*k;
}
//点p对直线v的对称点
Point Point_line_symmetry(Point p, Line v){
Point q = Point_line_proj(p,v);
return Point(2*q.x-p.x,2*q.y-p.y);
}
struct Circle{
Point c; //圆心
double r; //半径
Circle(){}
Circle(Point c,double r):c(c),r(r){}
Circle(double x,double y,double _r){c=Point(x,y);r = _r;}
};
//线段和圆的关系：0 线段在圆内, 1 线段和圆相切, 2 线段在圆外
int Seg_circle_relation(Segment v,Circle C){
double dst = Dis_point_seg(C.c,v);
if(sgn(dst-C.r) < 0) return 0;
if(sgn(dst-C.r) ==0) return 1;
return 2;
}
//直线和圆的关系：0 直线在圆内, 1 直线和圆相切, 2 直线在圆外
int Line_circle_relation(Line v,Circle C){
double dst = Dis_point_line(C.c,v);
if(sgn(dst-C.r) < 0) return 0;
if(sgn(dst-C.r) ==0) return 1;
return 2;
}
//直线和圆的交点，pa, pb是交点。返回值是交点个数
int Line_cross_circle(Line v,Circle C,Point &pa,Point &pb){
if(Line_circle_relation(v, C)==2) return 0; //无交点
Point q = Point_line_proj(C.c,v); //圆心在直线上的投影点
double d = Dis_point_line(C.c,v); //圆心到直线的距离
double k = sqrt(C.r*C.r-d*d);
if(sgn(k) == 0){ //1个交点，直线和圆相切
pa = q; pb = q; return 1;
}
Point n=(v.p2-v.p1)/ Len(v.p2-v.p1); //单位向量
pa = q + n*k;
pb = q - n*k;
return 2; //2个交点
}
int main() {
int T; scanf("%d", &T);
for (int cas = 1; cas <= T; cas++) {
Circle O; Point A,B,V;
scanf("%lf%lf%lf", &O.c.x, &O.c.y, &O.r);
scanf("%lf%lf%lf%lf", &A.x, &A.y, &V.x, &V.y);
scanf("%lf%lf", &B.x, &B.y);
Line l(A, A+V); //射线
Line t(A, B);
//情况（1）直线和圆不相交，而且直线经过点
if(Point_line_relation(B,l)==0
&& Seg_circle_relation(t,O)>=1 && sgn(Cross(B-A,V))==0)
printf("Case #%d: Yes\n", cas);
else{
Point pa, pb; //直线和圆的交点
//情况（2）直线和圆相切，不经过点
if(Line_cross_circle(l,O,pa,pb) != 2)
printf("Case #%d: No\n",cas);
//情况（3）直线和圆相交
else{
Point cut; //直线和圆的碰撞点
if(Distance(pa,A) > Distance(pb,A)) cut = pb;
else cut = pa;
Line mid(cut, O.c); //圆心到碰撞点的的直线;
Point en = Point_line_symmetry(A,mid); //镜像点
Line light(cut, en); //反射线
if(Distance(light.p2,B) > Distance(light.p1,B))
swap(light.p1, light.p2); if(sgn(Cross(light.p2-light.p1,
Point(B.x-cut.x,B.y-cut.y)))== 0)
printf("Case #%d: Yes\n", cas);
else
printf("Case #%d: No\n", cas);
}
}
}
return 0;
}
Buried memory :
最小圆覆盖问题。
#include <bits/stdc++.h>
using namespace std;
#define eps 1e-8
const int maxn = 100000+5;
int sgn(double x)
{
if (fabs(x)<eps) return 0;
else return x<0? -1:1;
}
struct Point
{
double x,y;
};
double Distance(Point A, Point B){return hypot(A.x-B.x,A.y-B.y);}
//求三角形abc的外接圆圆心
Point circle_center(const Point a, const Point b, const Point c)
{
Point center;
double a1 = b.x-a.x,b1=b.y-a.y,c1=(a1*a1+b1*b1)/2;
double a2 = c.x-a.x,b2=c.y-a.y,c2=(a2*a2+b2*b2)/2;
double d=a1*b2-a2*b1;
center.x=a.x+(c1*b2-c2*b1)/d;
center.y=a.y+(a1*c2-a2*c1)/d;
return center;
}
//求最小圆覆盖，返回圆心c和半径r:
void min_cover_circle(Point *p, int n,Point &c, double &r)
{
random_shuffle(p,p+n); //打乱所有点
c = p[0];r = 0; //第一个点
for (int i = 1; i < n; ++i) //剩下所有点
{
if (sgn(Distance(p[i],c)-r)>0) //pi在圆外部
{
c=p[i];r=0; //将圆心设为pi半径为0
for (int j = 0; j < i; ++j) //重新检查前面的点
{
if (sgn(Distance(p[j],c)-r)>0)//两点定圆
{
c.x=(p[i].x+p[j].x)/2;
c.y=(p[i].y+p[j].y)/2;
r=Distance(p[j],c);
for (int k = 0; k < j; ++k)
{
if (sgn(Distance(p[k],c)-r)>0)
{
c=circle_center(p[i],p[j],p[k]);
r=Distance(p[i],c);
}
}
}
}
}
}
}
int main(int argc, char const *argv[])
{
int n;
Point p[maxn];
Point c;double r;
while(~scanf("%d",&n)&&n)
{
for (int i = 0; i < n; ++i)
{
scanf("%lf%lf",&p[i].x,&p[i].y);
}
min_cover_circle(p,n,c,r);
printf("%.2lf %.2lf %.2lf\n",c.x,c.y,r);
}
return 0;
}
Weapon：
两个异面直线的距离。
这道题是两个无限延伸的圆柱体，转化为两条直线间的距离，然后减去两个截面圆的半径，取min。
#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
using namespace std;
struct mq
{
double x; //x,y,z表示垂直于圆表面的向量
double y;
double z;
double a; //a,b,c圆心的坐标
double b;
double c;
double r; //r圆的半径
};
mq node[42];
double solve(mq p1,mq p2)
{
double a1,b1,c1,a2,b2,c2;
double s1,s2,s3; //s向量
double q1,q2,q3;
double ans1,ans2,ans;
a1=p1.x,b1=p1.y,c1=p1.z;
a2=p2.x,b2=p2.y,c2=p2.z;
s1=b1*c2-b2*c1,s2=c1*a2-c2*a1;
s3=a1*b2-a2*b1;
q1=p2.a-p1.a,q2=p2.b-p1.b;
q3=p2.c-p1.c;
ans1=fabs(q1*s1+q2*s2+q3*s3);
ans2=sqrt(s1*s1+s2*s2+s3*s3);
ans=ans1/ans2;
return ans;
}
int main()
{
int tes,n,i,j;
double x1,y1,z1,x2,y2,z2,x3,y3,z3;
double a1,b1,c1,a2,b2,c2;
scanf("%d",&tes);
while(tes--)
{
scanf("%d",&n);
for(i=0;i<n;i++)
{
scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",&x1,&y1,&z1,&x2,&y2,&z2,&x3,&y3,&z3);
node[i].a=x1,node[i].b=y1,node[i].c=z1,
a1=x2-x1,b1=y2-y1,c1=z2-z1;
a2=x3-x1,b2=y3-y1,c2=z3-z1;
node[i].r=sqrt(a1*a1+b1*b1+c1*c1); //半径
node[i].x=b1*c2-b2*c1,node[i].y=c1*a2-c2*a1;
node[i].z=a1*b2-a2*b1;
}
int flag=0;
double res=10000000;
double tmp;
for(i=0;i<n;i++)
{
for(j=i+1;j<n;j++)
{
tmp=solve(node[i],node[j]); //tmp返回的是两条中间的线的距离
tmp=tmp-node[i].r-node[j].r;
if(tmp<=0)
{
flag=1;
break;
}
if(tmp<res)
res=tmp;
}
if(flag)
break;
}
if(flag) puts("Lucky");
else printf("%.2f\n",res);
}
return 0;
}
tetrahedron
首先假设内切球球心为(x0,x1,x2)，可以用r=3V / (S1+S2+S3+S4)得出半径，
这样对于四个平面列出三个方程，解得
xn=∑3i=0Aixn⋅Si/(S1+S2+S3+S4)
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#define LD double
#define sqr(x) ((x)*(x))
#define eps 1e-13
using namespace std;
struct node
{
LD x,y,z;
void scan()
{
scanf("%lf%lf%lf",&x,&y,&z);
}
void print()
{
printf("%.4lf %.4lf %.4lf\n",x,y,z);
}
LD length()
{
return sqrt(sqr(x)+sqr(y)+sqr(z));
}
node operator+(const node &tmp)
{
return (node){x+tmp.x,y+tmp.y,z+tmp.z};
}
node operator-(const node &tmp)
{
return (node){x-tmp.x,y-tmp.y,z-tmp.z};
}
node operator/(LD tmp)
{
return (node){x/tmp,y/tmp,z/tmp};
}
node operator*(LD tmp)
{
return (node){x*tmp,y*tmp,z*tmp};
}
};
node cross(node a,node b)
{
node ans;
ans.x = a.y*b.z - b.y*a.z;
ans.y = b.x*a.z - a.x*b.z;
ans.z = a.x*b.y - b.x*a.y;
return ans;
}
LD dist(node a,node b)
{
return (b-a).length();
}
LD dot(node a,node b)
{
return a.x*b.x + a.y*b.y + a.z*b.z;
}
LD get_angle(node a,node b)
{
LD tmp = dot(a,b)/a.length()/b.length();
return acos(tmp);
}
node get_node(node A,node B,node C)
{
LD Lth = (B-A).length() + (C-A).length() + (C-B).length();
cout << sqr(Lth-4) << endl;
LD r = fabs(cross(B-A,C-A).length()) / Lth;
cout << r*r << endl;
node v1 = C-A;
node v2 = B-A;
node v = (v1+v2)/(v1+v2).length();
LD d = (C-A).length()/2;
LD L = sqrt(sqr(d)+sqr(r));
v = v*L;
return A+v;
}
int main()
{
node A,B,C,D;
while(~scanf("%lf%lf%lf",&A.x,&A.y,&A.z))
{
B.scan();
C.scan();
D.scan();
if(fabs(dot(cross(B-A,C-A),D-A)) < eps)
{
puts("O O O O");
continue;
}
LD S1 = fabs(cross(B-D,C-D).length())/2;
LD S2 = fabs(cross(D-A,C-A).length())/2;
LD S3 = fabs(cross(B-A,D-A).length())/2;
LD S4 = fabs(cross(B-A,C-A).length())/2;
LD Ve = fabs(dot(cross(B-A,C-A),D-A))/6;
//计算体积的叉积公式
LD R = 3*Ve / ((S1+S2+S3+S4));
node ans;
ans.x = (S1*A.x + S2*B.x + S3*C.x + S4*D.x)/(S1+S2+S3+S4);
ans.y = (S1*A.y + S2*B.y + S3*C.y + S4*D.y)/(S1+S2+S3+S4);
ans.z = (S1*A.z + S2*B.z + S3*C.z + S4*D.z)/(S1+S2+S3+S4);
printf("%.4lf %.4lf %.4lf %.4lf\n",ans.x,ans.y,ans.z,R);
}
return 0;
}
Super Star：
最小球覆盖
模拟退火的方法。
#include <bits/stdc++.h>
using namespace std;
struct Point{
double x,y,z;
}p[35];
double dist(Point a,Point b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)+(a.z-b.z)*(a.z-b.z));
}
double ac(int n)
{
double ans = 1e9;
Point tmp;
tmp.x=tmp.y=tmp.z=0;
int s=1; double step=1000;//超时的话可以换100
double eps=1e-8;
while(step>eps)//模拟退火
{
for(int i=1;i<=n;i++){if(dist(tmp,p[s])<dist(tmp,p[i])) s=i;}
double Dist=dist(tmp,p[s]); ans = min(ans, Dist);
tmp.x+=(p[s].x-tmp.x)/Dist*step;
tmp.y+=(p[s].y-tmp.y)/Dist*step;
tmp.z+=(p[s].z-tmp.z)/Dist*step;
step*=0.99;
}
return ans;
}
int main(){
int n;
while(~scanf("%d",&n)){
if(!n) break;
for(int i=1;i<=n;i++){
scanf("%lf %lf %lf",&p[i].x,&p[i].y,&p[i].z);
}
double ans=ac(n);
printf("%.5f\n",ans);
}
}
3D Convex Hull ：
多面体求面的个数。
#include<cstdio>
#include<vector>
#include<cstring>
#include<iostream>
#include<algorithm>
#define PR 1e-8
#define N 510
using namespace std;
struct TPoint
{
double x,y,z;
TPoint() {}
TPoint(double _x,double _y,double _z):x(_x),y(_y),z(_z) {}
TPoint operator - (const TPoint p)
{
return TPoint(x-p.x,y-p.y,z-p.z);
}
TPoint operator * (const TPoint p)
{
return TPoint(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
}
double operator ^ (const TPoint p)
{
return x*p.x+y*p.y+z*p.z;
}//点积
};
struct fac
{
int a,b,c;//凸包一个面上的三个点的编号
bool ok;//该面是否是最终凸包中的面
};
struct T3dull
{
int n;//初始点数
TPoint ply[N];//初始点
int trianglecnt;//凸包上三角形数
fac tri[N];//凸包三角形
int vis[N][N];//点i到点j是属于哪个面
double dist(TPoint a)
{
return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}//两点长度
double area(TPoint a,TPoint b,TPoint c)
{
return dist((b-a)*(c-a));
}//三角形面积*2
double volume(TPoint a,TPoint b,TPoint c,TPoint d)
{
return (b-a)*(c-a)^(d-a);
}//四面体有向体积*6
double ptoplane(TPoint &p,fac &f)//正：点在面同向
{
TPoint m=ply[f.b]-ply[f.a],n=ply[f.c]-ply[f.a],t=p-ply[f.a];
return (m*n)^t;
}
void deal(int p,int a,int b)
{
int f=vis[a][b];
fac add;
if(tri[f].ok)
{
if((ptoplane(ply[p],tri[f]))>PR) dfs(p,f);
else
{
add.a=b,add.b=a,add.c=p,add.ok=1;
vis[p][b]=vis[a][p]=vis[b][a]=trianglecnt;
tri[trianglecnt++]=add;
}
}
}
void dfs(int p,int cnt)//维护凸包，如果点p在凸包外更新凸包
{
tri[cnt].ok=0;
deal(p,tri[cnt].b,tri[cnt].a);
deal(p,tri[cnt].c,tri[cnt].b);
deal(p,tri[cnt].a,tri[cnt].c);
}
bool same(int s,int e)//判断两个面是否为同一面
{
TPoint a=ply[tri[s].a],b=ply[tri[s].b],c=ply[tri[s].c];
return fabs(volume(a,b,c,ply[tri[e].a]))<PR&&fabs(volume(a,b,c,ply[tri[e].b]))<PR&&fabs(volume(a,b,c,ply[tri[e].c]))<PR;
}
void construct()//构建凸包
{
int i,j;
trianglecnt=0;
if(n<4) return ;
bool tmp=1;
for(int i=1; i<n; i++) //前两点不共点
{
if((dist(ply[0]-ply[i]))>PR)
{
swap(ply[1],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
tmp=1;
for(i=2; i<n; i++) //前三点不共线
{
if((dist((ply[0]-ply[1])*(ply[1]-ply[i])))>PR)
{
swap(ply[2],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
tmp=1;
for(i=3; i<n; i++) //前四点不共面
{
if(fabs((ply[0]-ply[1])*(ply[1]-ply[2])^(ply[0]-ply[i]))>PR)
{
swap(ply[3],ply[i]);
tmp=0;
break;
}
}
if(tmp) return ;
fac add;
for(int i=0; i<4; i++) //构建初始四面体
{
add.a=(i+1)%4,add.b=(i+2)%4,add.c=(i+3)%4,add.ok=1;
if((ptoplane(ply[i],add))>0)
swap(add.b,add.c);
vis[add.a][add.b]=vis[add.b][add.c]=vis[add.c][add.a]=trianglecnt;
tri[trianglecnt++]=add;
}
for(int i=4; i<n; i++) //构建更新凸包
{
for(j=0; j<trianglecnt; j++)
{
if(tri[j].ok&&(ptoplane(ply[i],tri[j]))>PR)
{
dfs(i,j);
break;
}
}
}
int cnt=trianglecnt;
trianglecnt=0;
for(i=0; i<cnt; i++)
{
if(tri[i].ok)
tri[trianglecnt++]=tri[i];
}
}
double area()//表面积
{
double ret=0;
for(int i=0; i<trianglecnt; i++)
ret+=area(ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
return ret/2.0;
}
double volume()//体积
{
TPoint p(0,0,0);
double ret=0;
for(int i=0; i<trianglecnt; i++)
ret+=volume(p,ply[tri[i].a],ply[tri[i].b],ply[tri[i].c]);
return fabs(ret/6);
}
int facetri()//表面三角形数
{
return trianglecnt;
}
int facepolygon()//表面多边形数
{
int ans=0,i,j,k;
for(i=0; i<trianglecnt; i++)
{
for(j=0,k=1; j<i; j++)
{
if(same(i,j))
{
k=0;
break;
}
}
ans+=k;
}
return ans;
}
} hull;
int main()
{
while(~scanf("%d",&hull.n))
{
int i;
for(i=0; i<hull.n; i++)
scanf("%lf%lf%lf",&hull.ply[i].x,&hull.ply[i].y,&hull.ply[i].z);
hull.construct();
printf("%d\n",hull.facepolygon());
}
return 0;
}
Rescue ：
求空间几何的重心。
题目中这句话：为了实现最佳控制，整艘船的控制中心位于质量的中心。就代表了求重心到凸多面体的最短距离，先求出重心，也就是控制室，然后求出重心到每一个面的距离，然后取min。
#include <bits/stdc++.h>
#define FF(i, a, b) for(int i=a; i<b; i++)
#define FD(i, a, b) for(int i=a; i>=b; i--)
#define REP(i, n) for(int i=0; i<n; i++)
#define CLR(a, b) memset(a, b, sizeof(a))
#define debug puts("**debug**")
#define LL long long
#define PB push_back
#define MP make_pair
#define eps 1e-8
using namespace std;
int dcmp(double x) { if(fabs(x) < eps) return 0; else return x < 0 ? -1 : 1; }
struct Point3 {
double x, y, z;
Point3(double x=0, double y=0, double z=0):x(x),y(y),z(z) { }
};
typedef Point3 Vector3;
Vector3 operator + (const Vector3& A, const Vector3& B) { return Vector3(A.x+B.x, A.y+B.y, A.z+B.z); }
Vector3 operator - (const Point3& A, const Point3& B) { return Vector3(A.x-B.x, A.y-B.y, A.z-B.z); }
Vector3 operator * (const Vector3& A, double p) { return Vector3(A.x*p, A.y*p, A.z*p); }
Vector3 operator / (const Vector3& A, double p) { return Vector3(A.x/p, A.y/p, A.z/p); }
bool operator == (const Point3& a, const Point3& b) {
return dcmp(a.x-b.x) == 0 && dcmp(a.y-b.y) == 0 && dcmp(a.z-b.z) == 0;
}
Point3 read_point3() {
Point3 p;
scanf("%lf%lf%lf", &p.x, &p.y, &p.z);
return p;
}
double Dot(const Vector3& A, const Vector3& B) { return A.x*B.x + A.y*B.y + A.z*B.z; }
double Length(const Vector3& A) { return sqrt(Dot(A, A)); }
double Angle(const Vector3& A, const Vector3& B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
Vector3 Cross(const Vector3& A, const Vector3& B) { return Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x); }
double Area2(const Point3& A, const Point3& B, const Point3& C) { return Length(Cross(B-A, C-A)); }
double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D) { return Dot(D-A, Cross(B-A, C-A)); }
Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D) { return (A + B + C + D)/4.0; }
double rand01() { return rand() / (double)RAND_MAX; }
double randeps() { return (rand01() - 0.5) * eps; }
Point3 add_noise(const Point3& p) {
return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
}
struct Face {
int v[3];
Face(int a, int b, int c) { v[0] = a; v[1] = b; v[2] = c; }
Vector3 Normal(const vector<Point3>& P) const {
return Cross(P[v[1]]-P[v[0]], P[v[2]]-P[v[0]]);
}
// f是否能看见P[i]
int CanSee(const vector<Point3>& P, int i) const {
return Dot(P[i]-P[v[0]], Normal(P)) > 0;
}
};
// 增量法求三维凸包
// 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动
vector<Face> CH3D(const vector<Point3>& P) {
int n = P.size();
vector<vector<int> > vis(n);
for(int i = 0; i < n; i++) vis[i].resize(n);
vector<Face> cur;
cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动，前三个点不共线
cur.push_back(Face(2, 1, 0));
for(int i = 3; i < n; i++) {
vector<Face> next;
// 计算每条边的“左面”的可见性
for(int j = 0; j < cur.size(); j++) {
Face& f = cur[j];
int res = f.CanSee(P, i);
if(!res) next.push_back(f);
for(int k = 0; k < 3; k++) vis[f.v[k]][f.v[(k+1)%3]] = res;
}
for(int j = 0; j < cur.size(); j++)
for(int k = 0; k < 3; k++) {
int a = cur[j].v[k], b = cur[j].v[(k+1)%3];
if(vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见
next.push_back(Face(a, b, i));
}
cur = next;
}
return cur;
}
struct ConvexPolyhedron {
int n;
vector<Point3> P, P2;
vector<Face> faces;
bool read() {
if(scanf("%d", &n) != 1) return false;
P.resize(n);
P2.resize(n);
for(int i = 0; i < n; i++) { P[i] = read_point3(); P2[i] = add_noise(P[i]); }
faces = CH3D(P2);
return true;
}
Point3 centroid() {
Point3 C = P[0];
double totv = 0;
Point3 tot(0,0,0);
for(int i = 0; i < faces.size(); i++) {
Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
double v = -Volume6(p1, p2, p3, C);
totv += v;
tot = tot + Centroid(p1, p2, p3, C)*v;
}
return tot / totv;
}
double mindist(Point3 C) {
double ans = 1e30;
for(int i = 0; i < faces.size(); i++) {
Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
}
return ans;
}
}P1;
int main()
{
int n, m;
ConvexPolyhedron P1, P2;
while(P1.read()) {
Point3 C1 = P1.centroid();
double d1 = P1.mindist(C1);
printf("%.3f\n", d1);
}
return 0;
}
525
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
