POJ 3675 裸的题了。直接上模板就行了。注意的是,求出的是有向面积,有可能是负数。
#include <iostream> #include <cstdio> #include <cmath> #include <vector> #include <cstring> #include <algorithm> #include <string> #include <set> #include <functional> #include <numeric> #include <sstream> #include <stack> #define CL(arr, val) memset(arr, val, sizeof(arr)) #define REP(i, n)for((i) = 0; (i) < (n); ++(i)) #define FOR(i, l, h)for((i) = (l); (i) <= (h); ++(i)) #define FORD(i, h, l)for((i) = (h); (i) >= (l); --(i)) #define L(x)(x) << 1 #define R(x)(x) << 1 | 1 #define MID(l, r)(l + r) >> 1 #define Min(x, y)x < y ? x : y #define Max(x, y)x < y ? y : x #define E(x)(1 << (x)) #define iabs(x) (x) < 0 ? -(x) : (x) #define OUT(x)printf("%I64d ", x) #define lowbit(x)(x)&(-x) #define Read()freopen("data.in", "r", stdin) #define Write()freopen("data.out", "w", stdout); const double eps = 1e-8; typedef long long LL; const int inf = ~0u>>2; using namespace std; inline double max (double a, double b) { if (a > b) return a; else return b; } inline double min (double a, double b) { if (a < b) return a; else return b; } inline int fi (double a){ if (a > eps) return 1; else if (a >= -eps) return 0; else return -1; } class Vector{ public: double x, y; Vector (void) {} Vector (double x0, double y0) : x(x0), y(y0) {} double operator * (const Vector& a) const { return x * a.y - y * a.x; } double operator % (const Vector& a) const { return x * a.x + y * a.y; } Vector verti (void) const { return Vector(-y, x); } double length (void) const { return sqrt(x * x + y * y); } Vector adjust (double len) { double ol = len / length(); return Vector(x * ol, y * ol); } Vector oppose (void) { return Vector(-x, -y); } }; class point{ public: double x, y; point (void) {} point (double x0, double y0) : x(x0), y(y0) {} Vector operator - (const point& a) const { return Vector(x - a.x, y - a.y); } point operator + (const Vector& a) const { return point(x + a.x, y + a.y); } }; class segment{ public: point a, b; segment (void) {} segment (point a0, point b0) : a(a0), b(b0) {} point intersect (const segment& s) const { Vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b; double s1 = v1 * v2, s2 = v3 * v4; double se = s1 + s2; s1 /= se, s2 /= se; return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1); } point pverti (const point& p) const { Vector t = (b - a).verti(); segment uv(p, p + t); return intersect(uv); } bool on_segment (const point& p) const { if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 && fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true; else return false; } }; double radius; point polygon[110]; double kuras_area (point a, point b, double cx, double cy){ point ori(cx, cy); int sgn = fi((b - ori) * (a - ori)); double da = (a - ori).length(), db = (b - ori).length(); int ra = fi(da - radius), rb = fi(db - radius); double angle = acos(((b - ori) % (a - ori)) / (da * db)); segment t(a, b); point h, u; Vector seg; double ans, dlt, mov, tangle; if (fi(da) == 0 || fi(db) == 0) return 0; else if (sgn == 0) return 0; else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn; else if (ra >= 0 && rb >= 0){ h = t.pverti(ori); dlt = (h - ori).length(); if (!t.on_segment(h) || fi(dlt - radius) >= 0) return radius * radius * (angle / 2) * sgn; else{ ans = radius * radius * (angle / 2); tangle = acos(dlt / radius); ans -= radius * radius * tangle; ans += radius * sin(tangle) * dlt; return ans * sgn; } } else{ h = t.pverti(ori); dlt = (h - ori).length(); seg = b - a; mov = sqrt(radius * radius - dlt * dlt); seg = seg.adjust(mov); if (t.on_segment(h + seg)) u = h + seg; else u = h + seg.oppose(); if (ra == 1) swap(a, b); ans = fabs((a - ori) * (u - ori)) / 2; tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length())); ans += radius * radius * (tangle / 2); return ans * sgn; } } const double pi = acos(-1.0); int main (){ int n; while(scanf("%lf",&radius)!=EOF){ scanf("%d",&n); for(int i=0;i<n;i++){ scanf("%lf%lf",&polygon[i].x,&polygon[i].y); } double ans=0; for(int i=0;i<n;i++){ ans+=kuras_area(polygon[i],polygon[(i+1)%n],0,0); } printf("%.2lf ",fabs(ans)); } }
HDU 3982。很明显,求出切割的半平面交之后就可以得出一个多边形。然后求多边形与圆面积相交部分即可。
#include <iostream> #include <cstdio> #include <cstring> #include <algorithm> #include <cmath> using namespace std; const double eps = 1e-6; typedef long long LL; const int inf = ~0u>>2; const int MAXN=4222; inline double max (double a, double b) { if (a > b) return a; else return b; } inline double min (double a, double b) { if (a < b) return a; else return b; } inline int fi (double a){ if (a > eps) return 1; else if (a >= -eps) return 0; else return -1; } struct Point { double x,y; Point(double xx,double yy){x=xx; y=yy;} Point(){} }convex[MAXN],cherry; class Vector{ public: double x, y; Vector (void) {} Vector (double x0, double y0) : x(x0), y(y0) {} double operator * (const Vector& a) const { return x * a.y - y * a.x; } //向量叉乘 double operator % (const Vector& a) const { return x * a.x + y * a.y; } //向量点乘 Vector verti (void) const { return Vector(-y, x); } double length (void) const { return sqrt(x * x + y * y); } Vector adjust (double len) { double ol = len / length(); return Vector(x * ol, y * ol); } Vector oppose (void) { return Vector(-x, -y); } }; class point{ public: double x, y; point (void) {} point (double x0, double y0) : x(x0), y(y0) {} Vector operator - (const point& a) const { return Vector(x - a.x, y - a.y); } point operator + (const Vector& a) const { return point(x + a.x, y + a.y); } void _copy(Point t){ x=t.x; y=t.y; } }; class segment{ public: point a, b; segment (void) {} segment (point a0, point b0) : a(a0), b(b0) {} point intersect (const segment& s) const { Vector v1 = s.a - a, v2 = s.b - a, v3 = s.b - b, v4 = s.a - b; double s1 = v1 * v2, s2 = v3 * v4; double se = s1 + s2; s1 /= se, s2 /= se; return point(a.x * s2 + b.x * s1, a.y * s2 + b.y * s1); } point pverti (const point& p) const { Vector t = (b - a).verti(); segment uv(p, p + t); return intersect(uv); } bool on_segment (const point& p) const { if (fi(min(a.x, b.x) - p.x) <= 0 && fi(p.x - max(a.x, b.x)) <= 0 && fi(min(a.y, b.y) - p.y) <= 0 && fi(p.y - max(a.y, b.y)) <= 0) return true; else return false; } }; double radius; point polygon[MAXN];//多边形 double kuras_area (point a, point b, double cx, double cy){ point ori(cx, cy); int sgn = fi((b - ori) * (a - ori)); double da = (a - ori).length(), db = (b - ori).length(); int ra = fi(da - radius), rb = fi(db - radius); double angle = acos(((b - ori) % (a - ori)) / (da * db)); segment t(a, b); point h, u; Vector seg; double ans, dlt, mov, tangle; if (fi(da) == 0 || fi(db) == 0) return 0; else if (sgn == 0) return 0; else if (ra <= 0 && rb <= 0) return fabs((b - ori) * (a - ori)) / 2 * sgn; else if (ra >= 0 && rb >= 0){ h = t.pverti(ori); dlt = (h - ori).length(); if (!t.on_segment(h) || fi(dlt - radius) >= 0) return radius * radius * (angle / 2) * sgn; else{ ans = radius * radius * (angle / 2); tangle = acos(dlt / radius); ans -= radius * radius * tangle; ans += radius * sin(tangle) * dlt; return ans * sgn; } } else{ h = t.pverti(ori); dlt = (h - ori).length(); seg = b - a; mov = sqrt(radius * radius - dlt * dlt); seg = seg.adjust(mov); if (t.on_segment(h + seg)) u = h + seg; else u = h + seg.oppose(); if (ra == 1) swap(a, b); ans = fabs((a - ori) * (u - ori)) / 2; tangle = acos(((u - ori) % (b - ori)) / ((u - ori).length() * (b - ori).length())); ans += radius * radius * (tangle / 2); return ans * sgn; } } const double pi = acos(-1.0); struct Line { Point a,b; double ang; }line[MAXN/2],st[MAXN/2]; int n,cnt; double operator *(const Point x,const Point y){ return x.x*y.y-x.y*y.x; } Point operator -(Point x,const Point y){ x.x-=y.x; x.y-=y.y; return x; } Point operator * (const Line &x,const Line &y){ //交点 double a1=(y.b-x.a)*(y.a-x.a),a2=(y.a-x.b)*(y.b-x.b); Point r; r.x=(x.a.x*a2+x.b.x*a1)/(a2+a1); r.y=(x.a.y*a2+x.b.y*a1)/(a2+a1); return r; } bool operator ==(const Point &a,const Point &b){ return fabs(a.x-b.x)<eps&&fabs(a.y-b.y)<eps; } bool operator <(const Line &x,const Line &y){ if(fabs(x.ang-y.ang)<eps) return (y.b-x.a)*(x.b-y.a)>eps; return x.ang <y.ang; } bool JudgeOut(const Line &x,const Point &p){ return (p-x.a)*(x.b-x.a)>eps; } bool Parellel(const Line &x,const Line &y){ return fabs((x.b-x.a)*(y.b-y.a))<eps; } void HplaneI(){ sort(line,line+n); int top=1,bot=0,tmp=1; for(int i=1;i<n;i++){ if(line[i].ang-line[i-1].ang>eps) line[tmp++]=line[i]; } n=tmp; st[0]=line[0],st[1]=line[1]; for(int i=2;i<n;i++){ if(Parellel(st[top],st[top-1])||Parellel(st[bot],st[bot+1]))return ; while(bot<top&&JudgeOut(line[i],st[top]*st[top-1])) top--; while(bot<top&&JudgeOut(line[i],st[bot]*st[bot+1])) bot++; st[++top]=line[i]; } while(bot<top&&JudgeOut(st[bot],st[top]*st[top-1])) top--; while(bot<top&&JudgeOut(st[top],st[bot]*st[bot+1])) bot++; if(top<=bot+1) return ; st[++top]=st[bot]; cnt=0; for(int i=bot;i<top;i++){ polygon[cnt++]._copy(st[i]*st[i+1]); } } int main(){ int T,icase=0; scanf("%d",&T); while(T--){ scanf("%lf%d",&radius,&n); for(int i=0;i<n;i++){ scanf("%lf%lf%lf%lf",&line[i].a.x,&line[i].a.y,&line[i].b.x,&line[i].b.y); } scanf("%lf%lf",&cherry.x,&cherry.y); for(int i=0;i<n;i++){ if((line[i].a-cherry)*(line[i].b-cherry)<0){ swap(line[i].a,line[i].b); } line[i].ang=atan2(line[i].b.y-line[i].a.y,line[i].b.x-line[i].a.x); } line[n].a=Point(-radius,-radius); line[n].b=Point(radius,-radius); line[n++].ang=0; line[n].a=Point(radius,-radius); line[n].b=Point(radius,radius); line[n++].ang=pi/2; line[n].a=Point(radius,radius); line[n].b=Point(-radius,radius); line[n++].ang=pi; line[n].a=Point(-radius,radius); line[n].b=Point(-radius,-radius); line[n++].ang=-pi/2; cnt=0; HplaneI(); double ans=0; for(int i=0;i<cnt;i++){ // printf("%lf %lf ",polygon[i].x,polygon[i].y); ans+=kuras_area(polygon[i],polygon[(i+1)%cnt],0,0); } printf("Case %d: ",++icase); printf("%.5lf",fabs(fabs(ans)-pi*radius*radius)<eps?100:fabs(ans)/(pi*radius*radius)*100); puts("%"); } return 0; }