整理一些计算几何的东西供自己学习。
文章目录
旋转卡壳
这篇是本魔芋不负责口胡的，纯属鬼扯。
传送门
半平面交
看代码吧，实在不知道该说啥。
bzoj2618: [Cqoi2006]凸多边形
求凸多边形面积交，半平面交模板题。
关于atan2
atan2(y,x):返回复数x+y*i的辐角,取值(-pi,pi].
用于极角排序,比atan稳定.
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<queue>
#include<set>
#define For(i,a,b) for(int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(int i=(a);i>=(b);i--)
#define Formylove return 0
const int N=20007;
typedef long long LL;
typedef double db;
using namespace std;
int n,k,cnt;
db ans;
template<typename T>void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator <(const pt &A,const pt &B) {
return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0);
}
}p[N];
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
struct Line {
pt a,b; db slop;
Line(){}
Line(pt a,pt b,db slop):a(a),b(b),slop(slop){}
}L[N],que[N];
bool cmp(const Line&A,const Line&B) {
return dcmp(A.slop-B.slop)<0||(dcmp(A.slop-B.slop)==0&&dcmp(cross(A.b-A.a,B.b-A.a))<0);
}
pt inter(Line A,Line B) {
db s1=cross(B.b-A.a,A.b-A.a);
db s2=cross(A.b-A.a,B.a-A.a);
db t=s1/(s1+s2);
return B.b+(B.a-B.b)*t;
}
int ck(Line A,Line B,Line C) {
pt P=inter(A,B);
return dcmp(cross(P-C.a,C.b-C.a))>=0;
}
db area(int n) {
db rs=0; p[n+1]=p[1];
For(i,1,n) rs+=cross(p[i],p[i+1]);
return fabs(rs)/2.0;
}
db solve(int cnt) {
int ql=1,qr=0;
sort(L+1,L+cnt+1,cmp);
que[++qr]=L[1]; que[++qr]=L[2];
For(i,3,cnt) {
while(ql<qr&&ck(que[qr-1],que[qr],L[i])) qr--;
while(ql<qr&&ck(que[ql+1],que[ql],L[i])) ql++;
que[++qr]=L[i];
}
while(ql+1<qr&&ck(que[qr-1],que[qr],que[ql])) qr--;
while(ql+1<qr&&ck(que[ql+1],que[ql],que[qr])) ql++;
que[qr+1]=que[ql];
int tot=0;
For(i,ql,qr) p[++tot]=inter(que[i],que[i+1]);
return area(tot);
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
read(n);
For(i,1,n) {
read(k);
For(j,1,k) { read(p[j].x); read(p[j].y); }
p[k+1]=p[1];
For(j,1,k) L[++cnt]=Line(p[j],p[j+1],atan2(p[j+1].y-p[j].y,p[j+1].x-p[j].x));
}
printf("%.3lf\n",solve(cnt));
Formylove;
}
最小圆覆盖
一个写得非常好的博客
感谢之前汪神教我的这个三点求圆心。
bzoj2823: [AHOI2012]信号塔
模板题，求圆心和半径。
upd:
如果j到i,k到j而不是j到i-1,k到j-1会WA，精度炸了ORZ。。。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<queue>
#include<set>
#define For(i,a,b) for(int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(int i=(a);i>=(b);i--)
#define Formylove return 0
const int N=500007;
typedef long long LL;
typedef double db;
using namespace std;
int n;
db R;
template<typename T>void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-15
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator <(const pt &A,const pt &B) {
return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0);
}
}p[N],O;
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return dot(A,A); }
pt get_circle_center(pt A,pt B,pt C) {
pt c=B-A,b=C-A,t;
db l1=lenth(c)/2,l2=lenth(b)/2;
t.x=(l1*b.y-l2*c.y)/(c.x*b.y-b.x*c.y);
t.y=(l1*b.x-l2*c.x)/(c.y*b.x-b.y*c.x);
return A+t;
}
void min_circle_cover(int n,pt &c,db &r) {
random_shuffle(p+1,p+n+1);
c=p[1],r=0;
For(i,2,n) if(dcmp(lenth(p[i]-c)-r)>0) {
c=p[i],r=0;
For(j,1,i-1) if(dcmp(lenth(p[j]-c)-r)>0) {
c=(p[i]+p[j])/2,r=lenth(p[i]-c);
For(k,1,j-1) if(dcmp(lenth(p[k]-c)-r)>0) {
c=get_circle_center(p[i],p[j],p[k]);
r=lenth(p[i]-c);
}
}
}
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
srand(332748118);
read(n);
For(i,1,n) scanf("%lf%lf",&p[i].x,&p[i].y);
min_circle_cover(n,O,R);
printf("%.2lf %.2lf %.2lf\n",O.x,O.y,sqrt(R));
Formylove;
}
平面最近点对
HDU-1007 Quoit Design
模板题，求平面最近点对，分治。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<queue>
#include<cmath>
const int N=1e5+7;
typedef long long LL;
typedef double db;
using namespace std;
int T,n;
template<typename T>void read(T &x) {
char ch=getchar(); x=0; T f=1;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
struct pt {
db x,y;
}p[N],tp[N];
bool cmp(const pt&A,const pt&B) { return A.x<B.x||(A.x==B.x&&A.y<B.y); }
bool cmpy(const pt&A,const pt&B) { return A.y<B.y; }
db sqr(db x) {return x*x;}
db dis(pt A,pt B) { return sqrt(sqr(A.x-B.x)+sqr(A.y-B.y)); }
db solve(int l,int r) {
db d=1e18;
if(l>=r) return d;
if(l+1==r) return dis(p[l],p[r]);
int mid=((l+r)>>1);
d=min(d,solve(l,mid));
d=min(d,solve(mid+1,r));
int k=0;
for(int i=l;i<=r;i++)
if(fabs(p[i].x-p[mid].x)<=d) tp[++k]=p[i];
sort(tp+1,tp+k+1,cmpy);
for(int i=1;i<=k;i++)
for(int j=i+1;j<=k&&tp[j].y-tp[i].y<d;j++)
d=min(d,dis(tp[i],tp[j]));
return d;
}
int main() {
#ifdef DEBUG
freopen(".in","r",stdin);
freopen(".out","w",stdout);
#endif
for(;;) {
read(n);
if(!n) break;
for(int i=1;i<=n;i++) scanf("%lf%lf",&p[i].x,&p[i].y);
sort(p+1,p+n+1,cmp);
printf("%.2lf\n",solve(1,n)/2);
}
return 0;
}
三维凸包
//Achen
#include<bits/stdc++.h>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
const int N=2019;
typedef long long LL;
typedef double db;
using namespace std;
int n,tot;
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x>0?1:-1); }
db Rand() { return (rand()/(db)RAND_MAX-0.5)*eps; }
struct vc {
db x,y,z;
vc(){}
vc(db x,db y,db z):x(x),y(y),z(z){}
void shake() { x+=Rand(); y+=Rand(); z+=Rand(); }
}v[N];
vc operator +(const vc&A,const vc&B) { return vc(A.x+B.x,A.y+B.y,A.z+B.z); }
vc operator -(const vc&A,const vc&B) { return vc(A.x-B.x,A.y-B.y,A.z-B.z); }
vc cross(vc A,vc B) { return vc(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x); }
db dot(vc A,vc B) { return A.x*B.x+A.y*B.y+A.z*B.z; }
db lenth(vc A) { return sqrt(dot(A,A)); }
struct face {
int a[3];
vc normal() { return cross(v[a[1]]-v[a[0]],v[a[2]]-v[a[0]]); }
db area() { return lenth(normal())/2.0; }
}f[N],tpf[N];
bool cansee(face fc,vc b) { return dot(b-v[fc.a[0]],fc.normal())>0; }
bool vis[N][N],ok;
void get_3D_ham(int n) {
f[++tot]=(face){1,2,3};
f[++tot]=(face){3,2,1};
For(i,4,n) {
int cc=0;
For(j,1,tot) {
if(!(ok=cansee(f[j],v[i]))) tpf[++cc]=f[j];
For(k,0,2) vis[f[j].a[k]][f[j].a[(k+1)%3]]=ok;
}
For(j,1,tot) For(k,0,2) {
int a=f[j].a[k],b=f[j].a[(k+1)%3];
if(vis[a][b]&&!vis[b][a]) tpf[++cc]=(face){a,b,i};
}
For(j,1,cc) f[j]=tpf[j]; tot=cc;
}
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
read(n);
For(i,1,n) {
scanf("%lf %lf %lf",&v[i].x,&v[i].y,&v[i].z);
v[i].shake();
}
get_3D_ham(n);
db ans=0;
For(i,1,tot)
ans+=f[i].area();
printf("%.3lf\n",ans);
//cerr<<clock()<<endl;
Formylove;
}
Pick定理、欧拉公式和圆的反演
我是可爱的传送门
poj2954Triangle(pick定理)
关于圆的小模板
两圆面积交
poj2546Circular Area
vjudge传送门
求两圆的面积交。余弦定理求角度，正弦定理求面积。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
#define pi acos(-1.0)
const int N=2019;
typedef long long LL;
typedef double db;
using namespace std;
int n,k;
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator <(const pt &A,const pt &B) {
return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0);
}
}p[N];
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }
struct circle {
pt o; db r;
friend bool operator <(const circle&A,const circle&B) {
return A.r<B.r;
}
void read() { scanf("%lf%lf%lf",&o.x,&o.y,&r); }
}C1,C2;
db get_inter_area(circle A,circle B) {
db d=lenth(A.o-B.o);
if(dcmp(d)==0||dcmp(d-(A.r-B.r))<=0) return B.r*B.r*pi;
if(dcmp(d-(A.r+B.r))>=0) return 0;
db theta1=acos((A.r*A.r+d*d-B.r*B.r)/(2.0*A.r*d));
db theta2=acos((B.r*B.r+d*d-A.r*A.r)/(2.0*B.r*d));
return A.r*A.r*theta1+B.r*B.r*theta2-A.r*A.r*sin(theta1*2.0)/2.0-B.r*B.r*sin(theta2*2.0)/2.0;
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
C1.read(); C2.read();
if(C1<C2) swap(C1,C2);
printf("%.3lf\n",get_inter_area(C1,C2));
Formylove;
}
圆的交点
UVA - 10969 Sweet Dream
读入n个圆r,x,y,后面的圆会覆盖前面的，求最终可以看到的边长。
求出每个圆和它后面的圆的交点，排序，考虑每两个点之间的弧能否看到，若这段弧的中点被它之后的圆覆盖，则看不到，否则能看到。注意特判圆内部，还有弧度统一到(0,2pi].
求交点用的余弦公式和反三角函数，误差较大但好写。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
#define pi acos(-1.0)
const int N=2019;
typedef long long LL;
typedef double db;
using namespace std;
int T,n,no[N];
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator <(const pt &A,const pt &B) {
return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0);
}
}p[N];
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }
struct circle {
pt o; db r;
friend bool operator <(const circle&A,const circle&B) {
return A.r<B.r;
}
void read() { scanf("%lf%lf%lf",&r,&o.x,&o.y); }
}C[N];
vector<db>vc[N],yvc[N];
db mo(db x) { return x<0?x+pi*2.0:(x>=pi*2.0?x-pi*2:x); }
void get_circle_inter(int id,circle A,circle B) {
db d=lenth(A.o-B.o);
if(dcmp(d)==0||dcmp(d-fabs(A.r-B.r))<=0) {
if(dcmp(B.r-A.r)>=0) no[id]=1;
return;
}
if(dcmp(d-(A.r+B.r))>=0) return ;
db theta=acos((A.r*A.r+d*d-B.r*B.r)/(2.0*A.r*d));
db bs=atan2(B.o.y-A.o.y,B.o.x-A.o.x);
vc[id].push_back(mo(bs-theta));
vc[id].push_back(mo(bs+theta));
//return A.r*A.r*theta1+B.r*B.r*theta2-A.r*A.r*sin(theta1*2.0)/2.0-B.r*B.r*sin(theta2*2.0)/2.0;
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
read(T);
while(T--) {
db ans=0;
read(n);
For(i,1,n) C[i].read(),no[i]=0,vc[i].clear(),yvc[i].clear();
Rep(i,n,1) {
For(j,i+1,n) get_circle_inter(i,C[i],C[j]);
int up=vc[i].size();
if(no[i]) continue;
if(up==0) ans+=pi*C[i].r*2.0;
else {
For(j,0,up-1)
yvc[i].push_back(vc[i][j]);
sort(vc[i].begin(),vc[i].end());
For(j,0,up-1) {
db l=vc[i][j],r=vc[i][(j+1)%up];
if(j==up-1) r+=2.0*pi;
db md=mo((l+r)/2.0);
int fl=1;
for(int k=0;k<up;k+=2) {
db ll=yvc[i][k],rr=yvc[i][k+1];
if(rr<ll) rr+=2.0*pi;
if((dcmp(md-ll)>=0&&dcmp(md-rr)<=0)||(dcmp(md+2.0*pi-ll)>=0&&dcmp(md+2.0*pi-rr)<=0)) {
fl=0; break;
}
}
if(fl) ans+=2.0*C[i].r*pi*(r-l)/(2.0*pi);
}
}
}
printf("%.3lf\n",ans);
}
Formylove;
}
两圆的公切线
UVA - 10674 Tangents
没搞清楚输出格式wa了好久。。重合时输出-1,没有时输出0。。。
注释在代码里，画图易看出。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
#define pi acos(-1.0)
const int N=2019;
typedef long long LL;
typedef double db;
using namespace std;
int T,cnt;
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator ==(const pt &A,const pt &B) { return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0; }
friend bool operator <(const pt &A,const pt &B) { return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0); }
}p[N];
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }
struct Line {
pt a,b;
Line(){}
Line(pt a,pt b):a(a),b(b){}
friend bool operator <(const Line&A,const Line&B) { return A.a<B.a||(A.a==B.a&&A.b<B.b); }
void out() {
printf("%.5lf %.5lf %.5lf %.5lf %.5lf\n",a.x,a.y,b.x,b.y,lenth(a-b));
}
}L[N];
struct circle {
pt o; db r;
friend bool operator <(const circle&A,const circle&B) { return A.r<B.r; }
void read() { scanf("%lf%lf%lf",&o.x,&o.y,&r); }
pt get(db t) { return pt(o.x+r*cos(t),o.y+r*sin(t)); }
}C1,C2;
void get_circle_inter(circle A,circle B) {
int fl=1;
if(A<B) swap(A,B),fl=0;
db d=lenth(A.o-B.o);
if(dcmp(d)==0&&dcmp(A.r-B.r)==0) { cnt=-1; return ; }
if(dcmp(d-(A.r-B.r))<0) return; //内含，无切线
db bs=atan2(B.o.y-A.o.y,B.o.x-A.o.x);
if(dcmp(d-(A.r-B.r))==0) { //内切，一条切线
pt x=A.get(bs),y=B.get(bs);
L[++cnt]=fl?Line(x,y):Line(y,x); return;
}
db th=acos((A.r-B.r)/d); //两条外切线
pt x=A.get(bs+th),y=B.get(bs+th);
L[++cnt]=fl?Line(x,y):Line(y,x);
x=A.get(bs-th),y=B.get(bs-th);
L[++cnt]=fl?Line(x,y):Line(y,x);
if(dcmp(d-(A.r+B.r))==0) { //外切，里面一条切线
x=A.get(bs),y=B.get(bs+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
}
else if(dcmp(d-(A.r+B.r))>0) { //外离,里面两条切线
th=acos((A.r+B.r)/d);
x=A.get(bs+th),y=B.get(bs+th+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
x=A.get(bs-th),y=B.get(bs-th+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
}
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
for(;;) {
C1.read(); C2.read();
if(C1.o.x==0&&C1.o.y==0&&C1.r==0&&C2.o.x==0&&C2.o.y==0&&C2.r==0) break;
cnt=0;
get_circle_inter(C1,C2);
printf("%d\n",cnt);
if(!cnt||cnt==-1) continue;
sort(L+1,L+cnt+1);
For(i,1,cnt) L[i].out();
}
Formylove;
}
圆的反演
具体见上方“Pick定理、欧拉公式和圆的反演”
HDU - 4773 Problem of Apollonius
给定两相离圆和圆外一点，求过该点的与两外两圆相切的圆。
反演不改变相切关系，不过点的圆反演后为圆，不过点的直线反演后为过点的圆，那么把两圆按该点以任意半径反演后，求两圆的切线，再把切线反演回圆即得答案。
圆关于圆的反演，直线关于圆的反演和公切线的模板。
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
#define pi acos(-1.0)
const int N=2019;
typedef long long LL;
typedef double db;
using namespace std;
int T,cnt,tot;
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-10
int dcmp(db x) { return fabs(x)<eps?0:(x<0?-1:1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator ==(const pt &A,const pt &B) { return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0; }
friend bool operator <(const pt &A,const pt &B) { return dcmp(A.x-B.x)<0||(dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)<0); }
void read() { scanf("%lf%lf",&x,&y); }
}sp;
pt operator +(const pt &A,const pt &B) { return pt(A.x+B.x,A.y+B.y); }
pt operator -(const pt &A,const pt &B) { return pt(A.x-B.x,A.y-B.y); }
pt operator *(const pt &A,const db &B) { return pt(A.x*B,A.y*B); }
pt operator /(const pt &A,const db &B) { return pt(A.x/B,A.y/B); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
db dot(pt A,pt B) { return A.x*B.x+A.y*B.y; }
db lenth(pt A) { return sqrt(dot(A,A)); }
struct Line {
pt a,b;
Line(){}
Line(pt a,pt b):a(a),b(b){}
friend bool operator <(const Line&A,const Line&B) { return A.a<B.a||(A.a==B.a&&A.b<B.b); }
void out() {
printf("%.5lf %.5lf %.5lf %.5lf %.5lf\n",a.x,a.y,b.x,b.y,lenth(a-b));
}
}L[N];
db node_to_line(pt A,Line B) {
return fabs(cross(A-B.a,B.b-B.a))/lenth(B.b-B.a);
}
pt foot(pt A,Line B) {
if(B.a.x==B.b.x) return pt(B.a.x,A.y);
if(B.a.y==B.b.y) return pt(A.x,B.a.y);
db k1=(B.a.y-B.b.y)/(B.a.x-B.b.x),k2=-1.0/k1,b=B.a.y-k1*B.a.x;
db x=(k2*A.x-A.y+b)/(k2-k1),y=k1*x+b;
return pt(x,y);
}
struct circle {
pt o; db r;
friend bool operator <(const circle&A,const circle&B) { return A.r<B.r; }
pt get(db t) { return pt(o.x+r*cos(t),o.y+r*sin(t)); }
void read() { scanf("%lf%lf%lf",&o.x,&o.y,&r); }
void out() { printf("%lf %lf %lf\n",o.x,o.y,r); }
}C1,C2,C1f,C2f,ansC[N];
void get_circle_inter(circle A,circle B) {
int fl=1;
if(A<B) swap(A,B),fl=0;
db d=lenth(A.o-B.o);
if(dcmp(d)==0&&dcmp(A.r-B.r)==0) { cnt=-1; return ; }
if(dcmp(d-(A.r-B.r))<0) return; //内含，无切线
db bs=atan2(B.o.y-A.o.y,B.o.x-A.o.x);
if(dcmp(d-(A.r-B.r))==0) { //内切，一条切线
pt x=A.get(bs),y=B.get(bs);
L[++cnt]=fl?Line(x,y):Line(y,x); return;
}
db th=acos((A.r-B.r)/d); //两条外切线
pt x=A.get(bs+th),y=B.get(bs+th);
L[++cnt]=fl?Line(x,y):Line(y,x);
x=A.get(bs-th),y=B.get(bs-th);
L[++cnt]=fl?Line(x,y):Line(y,x);
if(dcmp(d-(A.r+B.r))==0) { //外切，里面一条切线
x=A.get(bs),y=B.get(bs+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
}
else if(dcmp(d-(A.r+B.r))>0) { //外离,里面两条切线
th=acos((A.r+B.r)/d);
x=A.get(bs+th),y=B.get(bs+th+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
x=A.get(bs-th),y=B.get(bs-th+pi);
L[++cnt]=fl?Line(x,y):Line(y,x);
}
}
circle get_circle_inv(circle A,circle P) {
circle rs;
db d=lenth(A.o-P.o),dmi=d+A.r,dmx=d-A.r;
dmi=P.r*P.r/dmi; dmx=P.r*P.r/dmx;
rs.o=P.o+(A.o-P.o)*((dmx+dmi)/2.0/d);
rs.r=(dmx-dmi)/2.0;
return rs;
}
circle get_line_inv(Line A,circle P) {
circle rs;
pt ft=foot(P.o,A);
if(ft==P.o) { rs.r=0; return rs; }
db d=lenth(P.o-ft),dmx=P.r*P.r/d;
rs.o=P.o+(ft-P.o)*(dmx/d);
rs.o=(rs.o+P.o)/2;
rs.r=lenth(P.o-rs.o);
return rs;
}
void ck(circle C) {
if(C.r==0) return ;
if(dcmp(lenth(sp-C.o)-C.r)!=0) return ;
if(dcmp(lenth(C.o-C1.o)-C.r-C1.r)!=0) return ;
if(dcmp(lenth(C.o-C2.o)-C.r-C2.r)!=0) return ;
ansC[++tot]=C;
}
int main() {
//freopen("1.in","r",stdin);
//freopen("1.out","w",stdout);
read(T);
while(T--) {
C1.read(); C2.read(); sp.read();
circle P; P.o=sp; P.r=5.0;
C1f=get_circle_inv(C1,P);
C2f=get_circle_inv(C2,P);
cnt=tot=0;
get_circle_inter(C1f,C2f);
For(i,1,cnt) ck(get_line_inv(L[i],P));
printf("%d\n",tot);
For(i,1,tot) ansC[i].out();
}
Formylove;
}
自适应辛普森
floyd传递闭包
bzoj1027: [JSOI2007]合金
因为x+y+z=1，直接把每个三元组看成前两元为坐标的点
两块合金能合成的点在他们的线段上
n块合金能合成的点在他们构成的凸包内。
题目转换为给定m个点，求一个点数最少的凸包，使n个指定点都在凸包内。
若所有的n个点都在向量i,j的左侧，则
f
[
i
]
[
j
]
=
1
f[i][j]=1
f[i][j]=1，floyd求得最小的
f
[
i
]
[
i
]
f[i][i]
f[i][i]即答案。
特判一个点和两个点的凸包。
当只需要判断连通性的时候可以用bitset优化floyd传递闭包。
For(k,1,n) For(i,1,n) For(j,1,n)
f[i][j]|=(f[i][k]&f[k][j]);
bitset即：
For(k,1,n) For(i,1,n) if(f[i][k])
f[i]|=f[k];
合金代码：
//Achen
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdlib>
#include<cstdio>
#include<vector>
#include<cmath>
#define For(i,a,b) for(register int i=(a);i<=(b);i++)
#define Rep(i,a,b) for(register int i=(a);i>=(b);i--)
#define Formylove return 0
const int N=507;
typedef long long LL;
typedef double db;
using namespace std;
int n,m;
template<typename T> void read(T &x) {
char ch=getchar(); T f=1; x=0;
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=-1,ch=getchar();
for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0'; x*=f;
}
#define eps 1e-7
int dcmp(db x) { return fabs(x)<eps?0:(x>0?1:-1); }
struct pt {
db x,y;
pt(){}
pt(db x,db y):x(x),y(y){}
friend bool operator <(const pt&A,const pt&B) {
return A.x<B.x||(A.x==B.x&&A.y<B.y);
}
void read() { scanf("%lf%lf",&x,&y); }
}p[N],q[N];
pt operator -(const pt&A,const pt&B) { return pt(A.x-B.x,A.y-B.y); }
db cross(pt A,pt B) { return A.x*B.y-A.y*B.x; }
int f[N][N];
int main() {
freopen("1027.in","r",stdin);
freopen("1027.out","w",stdout);
read(m); read(n);
db nouse;
For(i,1,m) { p[i].read(); cin>>nouse; }
For(i,1,n) { q[i].read(); cin>>nouse; }
if(m==1) {
For(i,1,n) if(q[i].x!=p[1].x||q[i].y!=p[1].y) {
puts("-1");
Formylove;
}
puts("1"); Formylove;
}
For(i,1,m) For(j,1,m) {
if(i==j) { f[i][j]=m+1; continue; }
int ok=1;
For(k,1,n) if(dcmp(cross(p[j]-p[i],q[k]-p[i]))<0) {
ok=0; break;
}
f[i][j]=ok?1:m+1;
}
For(k,1,m) For(i,1,m) For(j,1,m)
f[i][j]=min(f[i][j],f[i][k]+f[k][j]);
int ans=m+1;
For(i,1,m) ans=min(ans,f[i][i]);
if(ans==m+1) ans=-1;
if(ans==2) {
db l=1e9,r=-1e9;
For(i,1,m) For(j,1,m) if(f[i][j]==1&&f[j][i]==1) {
l=min(min(l,p[i].x),p[j].x);
r=max(max(r,p[i].x),p[j].x);
}
For(i,1,n) if(dcmp(q[i].x-l)<0||dcmp(q[i].x-r)>0) {
puts("-1");
Formylove;
}
}
printf("%d\n",ans);
Formylove;
}
521
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
