【NOI2009】 描边
题目描述
描边 【问题描述】 小Z自幼就酷爱数学。聪明的他特别喜欢研究一些数学小问题。 有一天，小Z在一张纸上选择了n个点，并用铅笔将它们两两连接起来，构成n(n-1)/2条线段。由于铅笔很细，可以认为这些线段的宽度为0。 望着这些线段，小Z陷入了冥想中。他认为这些线段中的一部分比较重要，需要进行强调。因此小Z拿出了毛笔，将它们重新进行了描边。毛笔画在纸上，会形成一个半径为r的圆。在对一条线段进行描边时，毛笔的中心(即圆心)将从线段的一个端点开始，沿着该线段描向另一个端点。下图即为在一张4个点的图中，对其中一条线段进行描边强调后的情况。 现在，小Z非常想知道在描边之后纸面上共有多大面积的区域被强调，你能帮助他解答这个问题么？ 【输入文件】 这是一道提交答案型试题，所有的输入文件 path1.in~path10.in 已在相应目录下。 输入文件 path*.in 第一行包含一个正整数n，表示选择的点的数目。 第2至第n+1行，第i+1行有两个实数xi, yi，表示点i的坐标为(xi, yi)。 第n+2行有一个正整数m，表示小Z认为比较重要的线段的条数。 第n+3至第n+m+2行，每行有两个正整数a, b表示一条线段。a, b两个数分别表示该线段的两个端点的编号。 第n+m+3行，有一个实数r，表示毛笔在纸上形成的圆的半径。 第n+m+4行，有四个实数p1, p2, p3, p4，为评分使用的参数。 【输出文件】 输出文件path*.out仅包含一行，即为描边后被强调区域的总面积。 【输入样例】 2 1 1 1 2 1 1 2 1 0.00001 0.001 0.1 1 【输出样例】 5.1415927 【样例说明】 如下图所示。 【评分标准】 每个测试点单独评分。 本题设有4个评分参数p1,p2,p3,p4 (p1< p2 < p3 < p4)，已在输入文件中给出。你的得分将按照如下规则给出： 若你的答案与标准答案相差不超过p1，则该测试点你将得到满分； 否则，若你的答案与标准答案相差不超过p2，则你将得到该测试点70%的分数； 否则，若你的答案与标准答案相差不超过p3，则你将得到该测试点40%的分数； 否则，若你的答案与标准答案相差不超过p4，则你将得到该测试点10%的分数； 否则该测试点你的得分为0。
辛普森(pascal)
1 program path; 2 uses math; 3 const zero=1e-10; 4 //pi=3.141592653589793238; 5 maxn=5000; 6 type 7 circle=record x,y:extended; end; 8 rect=record 9 x1,y1,x2,y2,x3,y3,x4,y4:extended; 10 end; 11 interval=record l,r:extended; end; 12 con=record x,y:longint; end; 13 ty=array[0..maxn] of interval; 14 var 15 n,m,tot:longint; 16 r:extended; 17 d:array[0..maxn] of con; 18 b:array[0..maxn] of rect; 19 a:array[0..maxn] of circle; 20 t,c:ty; 21 get:array[0..maxn] of boolean; 22 //================= 23 procedure sort(l,r:longint; var a:ty); 24 var 25 i,j:longint; 26 m:extended; 27 t:interval; 28 begin 29 i:=l; j:=r; m:=a[(i+j) >> 1].l; 30 repeat 31 while a[i].l<m do inc(i); 32 while a[j].l>m do dec(j); 33 if i<=j then 34 begin 35 t:=a[i]; a[i]:=a[j]; a[j]:=t; 36 inc(i); dec(j); 37 end; 38 until i>j; 39 if i<r then sort(i,r,a); 40 if j>l then sort(l,j,a); 41 end; 42 //================= 43 procedure up(var l:extended; var r:extended; x,x1,y1,x2,y2:extended); 44 var 45 t:extended; 46 begin 47 if (x-x1>=-zero)and(x-x2<=zero)or(x-x2>=-zero)and(x-x1<=zero) then 48 begin 49 t:=(y2-y1)/(x2-x1)*(x-x1)+y1; 50 l:=min(l,t); r:=max(r,t); 51 end; 52 end; 53 //================= 54 function f(x:extended):extended; 55 var 56 i,t:longint; 57 len,sl,sr:extended; 58 begin 59 t:=0; 60 for i:=1 to n do 61 if abs(a[i].x-x)<r then 62 begin 63 len:=sqrt(sqr(r)-sqr(a[i].x-x)); 64 inc(t); 65 c[t].l:=a[i].y-len; c[t].r:=a[i].y+len; 66 end; 67 for i:=1 to m do 68 begin 69 sl:=1e10; sr:=-1e10; 70 up(sl,sr,x,b[i].x1,b[i].y1,b[i].x2,b[i].y2); 71 up(sl,sr,x,b[i].x2,b[i].y2,b[i].x3,b[i].y3); 72 up(sl,sr,x,b[i].x3,b[i].y3,b[i].x4,b[i].y4); 73 up(sl,sr,x,b[i].x4,b[i].y4,b[i].x1,b[i].y1); 74 if (sl<1e9)and(sr>-1e9) then 75 begin inc(t); c[t].l:=sl; c[t].r:=sr; end; 76 end; 77 if t=0 then exit(0); 78 sort(1,t,c); 79 sl:=c[1].l; sr:=c[1].r; 80 len:=0; 81 for i:=2 to t do 82 if sr<c[i].l then 83 begin 84 len:=len+sr-sl; 85 sl:=c[i].l; sr:=c[i].r; 86 end else if sr<c[i].r then sr:=c[i].r; 87 len:=len+sr-sl; 88 exit(len); 89 end; 90 //================= 91 function simpson(l,r,fl,fm,fr:extended):extended; inline; 92 begin exit((fl+4*fm+fr)*(r-l)/6); end; 93 //================= 94 function area(l,m,r,fl,fm,fr,pre:extended):extended; inline; 95 var 96 lm,rm,left,right,flm,frm:extended; 97 begin 98 lm:=(l+m)/2; rm:=(r+m)/2; 99 flm:=f(lm); frm:=f(rm); 100 left:=simpson(l,m,fl,flm,fm); 101 right:=simpson(m,r,fm,frm,fr); 102 if abs(left+right-pre)<zero then exit(pre) else 103 exit(area(l,lm,m,fl,flm,fm,left)+area(m,rm,r,fm,frm,fr,right)); 104 end; 105 //================= 106 function find(l,r:extended):extended; inline; 107 var 108 m,fl,fm,fr,pre:extended; 109 begin 110 m:=(l+r)/2; 111 fl:=f(l); fm:=f(m); fr:=f(r); 112 pre:=simpson(l,r,fl,fm,fr); 113 exit(area(l,m,r,fl,fm,fr,pre)); 114 end; 115 //================= 116 procedure init; 117 var 118 i:longint; 119 dx,dy,alp,x,y:extended; 120 begin 121 alp:=random(maxlongint)/maxlongint*2*pi; 122 dx:=random(maxlongint)/maxlongint*1000; 123 dy:=random(maxlongint)/maxlongint*1000; 124 read(n); 125 for i:=1 to n do 126 begin 127 read(x,y); 128 a[i].x:=cos(alp)*(x+dx)-sin(alp)*(y+dy); 129 a[i].y:=sin(alp)*(x+dx)+cos(alp)*(y+dy); 130 end; 131 read(m); 132 for i:=1 to m do 133 begin 134 read(d[i].x,d[i].y); 135 get[d[i].x]:=true; get[d[i].y]:=true; 136 end; 137 read(r); 138 end; 139 //================= 140 procedure prepare; 141 var 142 i,x,y:longint; 143 dis,dx,dy:extended; 144 begin 145 for i:=1 to m do 146 begin 147 x:=d[i].x; y:=d[i].y; 148 dis:=sqrt(sqr(a[x].x-a[y].x)+sqr(a[x].y-a[y].y)); 149 dx:=-r*(a[y].y-a[x].y)/dis; 150 dy:= r*(a[y].x-a[x].x)/dis; 151 b[i].x1:=a[x].x+dx; b[i].y1:=a[x].y+dy; 152 b[i].x2:=a[x].x-dx; b[i].y2:=a[x].y-dy; 153 b[i].x3:=a[y].x-dx; b[i].y3:=a[y].y-dy; 154 b[i].x4:=a[y].x+dx; b[i].y4:=a[y].y+dy; 155 end; 156 x:=0; 157 for i:=1 to n do 158 if get[i] then 159 begin inc(x); a[x]:=a[i]; end; 160 n:=x; x:=0; 161 for i:=1 to n do 162 begin 163 inc(x); 164 t[x].l:=a[i].x-r; t[x].r:=a[i].x+r; 165 end; 166 for i:=1 to m do 167 begin 168 inc(x); 169 t[x].l:=min(min(b[i].x1,b[i].x2),min(b[i].x3,b[i].x4)); 170 t[x].r:=max(max(b[i].x1,b[i].x2),max(b[i].x3,b[i].x4)); 171 end; 172 tot:=x; 173 end; 174 //================= 175 procedure main; 176 var 177 i:longint; 178 ans,l,r:extended; 179 begin 180 sort(1,tot,t); 181 ans:=0; 182 l:=t[1].l; r:=t[1].r; 183 for i:=2 to tot do 184 if r<t[i].l then 185 begin 186 ans:=ans+find(l,r); 187 l:=t[i].l; r:=t[i].r; 188 end else 189 if r<t[i].r then r:=t[i].r; 190 ans:=ans+find(l,r); 191 writeln(ans:0:10); 192 end; 193 //================= 194 begin 195 assign(input,'path1.in'); reset(input); 196 assign(output,'path1.out'); rewrite(output); 197 randomize; 198 init; 199 prepare; 200 main; 201 close(input); close(output); 202 end.
辛普森（c++）
1 #include <cstdio> 2 #include <cstdlib> 3 #include <cmath> 4 #include <algorithm> 5 #include <string> 6 #include <cstring> 7 #include <ctime> 8 #define rep(i, n) for (int i=1; i<n+1; ++i) 9 #define sqr(x) ((x) * (x)) 10 typedef long double LD; 11 using namespace std; 12 const int maxn=5000; 13 const LD pi=3.141592654, zero=1e-8; 14 int n,m,tot; LD R; 15 bool get[maxn]; 16 17 struct point { int x,y;} d[maxn]; 18 struct circle { LD x,y;} a[maxn]; 19 struct rect { LD x1,y1,x2,y2,x3,y3,x4,y4;} rec[maxn]; 20 struct inter { LD l,r;} t[maxn], c[maxn]; 21 22 inline bool cmp(const inter &a,const inter &b) {return a.l < b.l;} 23 24 inline void up(LD &l, LD &r, LD x, LD x1, LD y1, LD x2, LD y2) 25 { 26 if ((x1<x+zero && x-zero<x2) || (x2<x+zero && x-zero<x1)) 27 { 28 LD y = (y2 - y1)/(x2 - x1)*(x-x1)+y1; 29 l = (l > y? y:l) ; 30 r = (r < y? y:r) ; 31 } 32 } 33 34 LD f(LD x) 35 { 36 int j = 0; 37 rep(i,n) if (fabs(a[i].x-x)<R) 38 { 39 LD len = sqrt(sqr(R)-sqr(a[i].x-x)); 40 c[++j].r = a[i].y + len; 41 c[j].l = a[i].y - len; 42 } 43 rep(i,m) 44 { 45 LD l = 1e10, r = -1e10; 46 up(l, r, x, rec[i].x1, rec[i].y1, rec[i].x2, rec[i].y2); 47 up(l, r, x, rec[i].x2, rec[i].y2, rec[i].x3, rec[i].y3); 48 up(l, r, x, rec[i].x3, rec[i].y3, rec[i].x4, rec[i].y4); 49 up(l, r, x, rec[i].x4, rec[i].y4, rec[i].x1, rec[i].y1); 50 if (l<1e9 && r>-1e9) { c[++j].l = l, c[j].r = r;} 51 } 52 if (j == 0) {return 0;} 53 sort(c+1, c+j+1, cmp); 54 LD l = c[1].l, r = c[1].r, s = 0; 55 for (int i=2; i<j+1; ++i) if (c[i].l>r+zero) s += r - l, l = c[i].l, r = c[i].r; 56 else r = (r<c[i].r? c[i].r:r); 57 return s+r-l; 58 } 59 60 inline LD simpson(LD l, LD r, LD fl, LD fm, LD fr) 61 {return (r - l) / 6 * (fl + fr + 4 * fm);} 62 63 inline LD area(LD l, LD fl, LD m, LD fm, LD r, LD fr, LD pre) 64 { 65 LD ls = (l + m)/2, rs = (m + r)/2; 66 LD fls = f(ls), frs = f(rs); 67 LD la = simpson(l,m,fl,fls,fm), ra = simpson(m,r,fm,frs,fr); 68 return fabs(la+ra-pre) < zero ? pre:area(l,fl,ls,fls,m,fm,la)+area(m,fm,rs,frs,r,fr,ra); 69 } 70 71 inline LD calc(LD l, LD r) 72 { 73 LD fl = f(l), fr = f(r), m = (l+r)/2, fm = f(m); 74 LD pre = simpson(l,r,fl,fm,fr); 75 return area(l,fl,m,fm,r,fr,pre); 76 } 77 78 void init() 79 { 80 srand(time(0)); 81 LD alp = (LD)rand() / RAND_MAX * 2 * pi, 82 dx = (LD)rand() / RAND_MAX * 1000, 83 dy = (LD)rand() / RAND_MAX * 1000; 84 scanf("%d",&n); 85 rep(i,n) 86 { 87 LD x,y; 88 scanf("%Lf%Lf", &x, &y); 89 a[i].x = cos(alp)*(dx+x) - sin(alp)*(dy+y); 90 a[i].y = sin(alp)*(dx+x) + cos(alp)*(dy+y); 91 } 92 scanf("%d",&m); 93 rep(i,m) scanf("%d%d", &d[i].x, &d[i].y) , get[d[i].x]=get[d[i].y]=1; 94 scanf("%Lf",&R); 95 } 96 97 void prepare() 98 { 99 int j=0; 100 rep(i,n) if (get[i]) a[++j]=a[i]; 101 n=j; 102 rep(i,m) 103 { 104 int x = d[i].x , y = d[i].y; 105 LD dis = sqrt(sqr(a[x].x - a[y].x)+sqr(a[x].y - a[y].y)); 106 LD dx =-R*(a[y].y - a[x].y) / dis , dy =R*(a[y].x - a[x].x) / dis; 107 rec[i].x1 = a[x].x + dx , rec[i].y1 = a[x].y + dy, 108 rec[i].x2 = a[x].x - dx , rec[i].y2 = a[x].y - dy, 109 rec[i].x3 = a[y].x - dx , rec[i].y3 = a[y].y - dy, 110 rec[i].x4 = a[y].x + dx , rec[i].y4 = a[y].y + dy; 111 } 112 j = 0; 113 rep(i,n) t[++j].l = a[i].x - R , t[j].r = a[i].x + R; 114 rep(i,m) 115 { 116 t[++j].l = min(min(rec[i].x1, rec[i].x2), min(rec[i].x3, rec[i].x4)); 117 t[j].r = max(max(rec[i].x1, rec[i].x2), max(rec[i].x3, rec[i].x4)); 118 } 119 tot = j; 120 sort(t+1,t+tot+1,cmp); 121 } 122 123 void work() 124 { 125 LD ans = 0, l = t[1].l, r = t[1].r; 126 for (int i = 2; i<tot+1; ++i) 127 if (t[i].l>r+zero) ans += calc(l,r), l = t[i].l, r = t[i].r; 128 else r = (r<t[i].r? t[i].r:r); 129 ans += calc(l,r); 130 printf("%.8Lf", ans); 131 } 132 133 int main() 134 { 135 freopen("path10.in","r",stdin); 136 freopen("path10.out","w",stdout); 137 init(); 138 prepare(); 139 work(); 140 return 0; 141 }
