计算几何练习

以前基本没有做过计算几何， 为了练习， 现在做了几道。

 3630: [JLOI2014]镜面通道

题目大意：给定一个二维的通道，通道内有一些正方形和圆形的零件(允许重叠)，求最少删掉多少零件可以使光线通过反射通过这个通道

首先我们大概可以感知到只要通道的两端是联通的就一定可以有一种可行方案 （想了想并不会证，，似乎有一个什么定理？）

然后不就是一个简单最小割了吗，，，通道的上面是S， 下面是T，两个元素相交的话就连一条边。就好了。

于是现在是一个判断 两个圆之间， 两个矩形之间， 一个圆和一个矩形之间 是否相交的问题。

两个圆很好判，有矩形的话因为这道题中所有的矩形都是平行于坐标轴的所以只有四个顶点是关键点， 圆也只有最上面最下面最左面最右面四个点是关键点，然后枚举一下随便判一判就好啦，，感觉很不优美。

  1 #include <iostream>
  2 #include <cstdio>
  3 #include <cstring>
  4 #include <algorithm>
  5 #include <cmath>
  6 #include <queue>
  7 #define MAXN 1005
  8 #define INF 1<<29
  9 #define T 2*n+1
 10 #define fi 2*i-1
 11 #define si 2*i
 12 #define fj 2*j-1
 13 #define sj 2*j
 14 using namespace std;
 15 int X, Y, n, head[MAXN], ee = -1, dis[MAXN];
 16 struct point{
 17     int kk, x1, y1, r1, x2, y2;
 18 }p[MAXN];
 19 struct Edge{
 20     int to, next, f;
 21 }edge[5000005];
 22 inline void addedge(int x, int y, int f){//printf("!! %d %d %d
", x, y, f);
 23 edge[++ ee].to = y; edge[ee].f = f; edge[ee].next = head[x]; head[x] = ee;}
 24 inline double dist(int x1, int y1, int x2, int y2){return sqrt((double)((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)));}
 25 inline void Add(int i, int j){//printf("***** %d %d
", i, j); 
 26 addedge(si, fj, INF); addedge(fj, si, 0); addedge(sj, fi, INF); addedge(fi, sj, 0);}
 27 bool inz(int x, int y, int i){
 28     if(x >= p[i].x1 && x <= p[i].x2 && y >= p[i].y1 && y <= p[i].y2) return 1; return 0;
 29 }
 30 int q[1000005];
 31 bool bfs(){
 32     for(int i = 0; i <= T; i ++) dis[i] = -1;
 33     int hh = 0, tt = 1;
 34     dis[0] = 1; q[0] = 0;
 35     while(hh != tt){
 36         int u = q[hh ++];
 37         for(int i = head[u]; i != -1; i = edge[i].next) if(dis[edge[i].to] == -1 && edge[i].f){
 38             dis[edge[i].to] = dis[u] + 1;
 39             q[tt ++] = edge[i].to;
 40             if(edge[i].to == T)return 1;
 41         }
 42     }return 0;
 43 }
 44 int dfs(int now, int maxf){
 45     int cnt = 0;
 46     if(now == T) return maxf;
 47     for(int i = head[now]; i != -1; i = edge[i].next) if(dis[edge[i].to] == dis[now] + 1 && edge[i].f){
 48         int f = dfs(edge[i].to, min(maxf - cnt, edge[i].f));
 49         edge[i].f -= f; edge[i^1].f += f; cnt += f;
 50         if(!f) dis[edge[i].to] = -1;
 51         if(cnt == maxf) return cnt;
 52     }return cnt;
 53 }
 54 int dinic(int s, int t){
 55     int ret = 0;
 56     while(bfs()){ ret += dfs(s, INF);}
 57     return ret;
 58 }
 59 int main(){
 60     scanf("%d%d%d", &X, &Y, &n);
 61     memset(head, -1 , sizeof(head));
 62     for(int i = 1; i <= n; i ++){
 63         scanf("%d", &p[i].kk);
 64         if(p[i].kk == 1) scanf("%d%d%d", &p[i].x1, &p[i].y1, &p[i].r1);
 65         else scanf("%d%d%d%d", &p[i].x1, &p[i].y1, &p[i].x2, &p[i].y2);
 66     }
 67     for(int i = 1; i <= n; i ++) addedge(fi, si, 1), addedge(si, fi, 0);
 68     for(int i = 1; i <= n; i ++) {
 69         if(p[i].kk == 1 && p[i].y1 - p[i].r1 <= 0) addedge(0, fi, INF), addedge(fi, 0, 0);
 70         if(p[i].kk == 2 && p[i].y1 <= 0) addedge(0, fi, INF), addedge(fi, 0, 0);
 71     }
 72     for(int i = 1; i <= n; i ++){
 73         if(p[i].kk == 1 && p[i].y1 + p[i].r1 >= Y) addedge(si, T, INF), addedge(T, si, 0);
 74         if(p[i].kk == 2 && p[i].y2 >= Y) addedge(si, T, INF), addedge(T, si, 0);
 75     }
 76     for(int i = 1; i <= n; i ++)
 77         for(int j = i+1; j <= n; j ++){
 78             if(p[i].kk == 1 && p[j].kk == 1)
 79                 if(dist(p[i].x1, p[i].y1, p[j].x1, p[j].y1) <= p[i].r1 + p[j].r1) Add(i, j);
 80             if(p[i].kk + p[j].kk == 3){
 81                 int ii = i, jj = j; if(p[i].kk == 1) swap(ii, jj);
 82                 if(dist(p[ii].x1, p[ii].y1, p[jj].x1, p[jj].y1) <= p[jj].r1) {Add(i,j); continue;}
 83                 if(dist(p[ii].x1, p[ii].y2, p[jj].x1, p[jj].y1) <= p[jj].r1) {Add(i,j); continue;}
 84                 if(dist(p[ii].x2, p[ii].y1, p[jj].x1, p[jj].y1) <= p[jj].r1) {Add(i,j); continue;}
 85                 if(dist(p[ii].x2, p[ii].y2, p[jj].x1, p[jj].y1) <= p[jj].r1) {Add(i,j); continue;}
 86                 int ys = p[jj].x1 + p[jj].r1, yx = p[jj].x1 - p[jj].r1, yz = p[jj].y1 - p[jj].r1, yo = p[jj].y1 + p[jj].r1;
 87                 if(inz(ys,p[jj].y1,ii)){Add(i, j); continue;}
 88                 if(inz(yx,p[jj].y1,ii)) {Add(i, j); continue;}
 89                 if(inz(p[jj].x1,yz,ii)) {Add(i, j); continue;}
 90                 if(inz(p[jj].x1,yo,ii)) {Add(i, j); continue;}
 91             }
 92             if(p[i].kk == 2 && p[j].kk == 2){
 93                 if(inz(p[i].x1, p[i].y1, j)) {Add(i, j); continue;}
 94                 if(inz(p[i].x2, p[i].y1, j)) {Add(i, j); continue;}
 95                 if(inz(p[i].x1, p[i].y2, j)) {Add(i, j); continue;}
 96                 if(inz(p[i].x2, p[i].y2, j)) {Add(i, j); continue;}
 97                 if(inz(p[j].x1, p[j].y1, i)) {Add(i, j); continue;}
 98                 if(inz(p[j].x2, p[j].y1, i)) {Add(i, j); continue;}
 99                 if(inz(p[j].x1, p[j].y2, i)) {Add(i, j); continue;}
100                 if(inz(p[j].x2, p[j].y2, i)) {Add(i, j); continue;}
101             }
102         }
103  
104     cout << dinic(0, T) << endl;
105 //  system("pause");
106     return 0;
107 }

1502: [NOI2005]月下柠檬树

题目大意：给定一棵由圆台和圆锥构成的柠檬树，月光以α的夹角平行射向地面，求阴影部分面积

 把这棵树投影到地上是一个很简单的问题，注意宽度没有变化，高度被压缩了就可以了。

下面的问题就变成了求一堆圆和一堆梯形的并集的大小。

因为有圆又有梯形比较不好求，于是强行上 自适应simpson。

/*

这里插播一句所有和圆有关的东西都显然是不能直接simpson积分的，，因为 (x0 + 4*xmid + x1) / 6.0 * l 这个式子只适用于 ax^3 + bx^2 + cx^1  + d 这样的式子啊！！！ 带根号的式子显然是不行的啊！！！！

窝前两天在赛场上竟然试图直接用simpson求弓形面积。。。。。真是傻爆了啊。。。

*/

/*

大家都知道simpson公式是Newton–Cotes formulas的一个特殊情况，去wiki看了一下，还有这样一个叫做Boole's rule的式子：

可以处理5次以内的函数哦！

*/

/*

以下几题窝连续写了好多次最后都因为浏览器崩了没有保存下来。。。。。可能是这几道题嫌我太弱了并不愿意让我来写题解。。所以就不仔细写了。

*/

 关于自适应simpson的正确性和复杂度都不知道怎么证明。。求大神科普

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#define eps 0.000001
using namespace std;
int n;
double alpha, H[505], R[505], ll = 1<<30, rr = -(1<<30);
struct Line{double x1, y1, x2, y2;}L[505];
double S(double x){
    double ret = 0;
    for(int i = 1; i <= n; i ++) if(H[i]-R[i] <= x && H[i]+R[i] >= x)
        ret = max(ret, sqrt(R[i]*R[i] - (x-H[i])*(x-H[i])));
    for(int i = 1; i < n; i ++) if(L[i].x1 <= x & L[i].x2 >= x)
        ret = max(ret, L[i].y1 - (L[i].y1-L[i].y2) * (x-L[i].x1) / (L[i].x2-L[i].x1));
    return ret;
}
double Simpson(double l, double r){
    double mid = (l + r) / 2.0;
    return (S(l) + 4*S(mid) + S(r))/6;
}
double Calc(double l, double r){
    double mid = (l + r) / 2.0, ret = Simpson(l, r)*(r-l);
    if(fabs(ret - Simpson(l, mid)*(mid-l) - Simpson(mid, r)*(r-mid)) <= eps) return ret;
    return Calc(l, mid) + Calc(mid, r);
}
int main(){
    cin >> n >> alpha;
    alpha = 1.0 / tan(alpha); n ++;
    for(int i = 1; i <= n; i ++) scanf("%lf", &H[i]), H[i] *= alpha, H[i] += H[i-1];
    for(int i = 1; i < n; i ++) scanf("%lf", &R[i]), ll = min(ll, H[i]-R[i]), rr = max(rr, H[i]+R[i]); 
    R[n] = 0.0; ll = min(ll, H[n]-R[n]); rr = max(rr, H[n] + R[n]);
    for(int i = 1; i < n; i ++){
        double sina = (R[i]-R[i+1])/(H[i+1]-H[i]);
        L[i].x1 = H[i] + R[i]*sina; L[i].x2 = H[i+1] + R[i+1]*sina;
        L[i].y1 = sqrt(R[i]*R[i] - (L[i].x1-H[i])*(L[i].x1-H[i]));
        L[i].y2 = sqrt(R[i+1]*R[i+1] - (L[i].x2-H[i+1])*(L[i].x2-H[i+1]));
    }
    printf("%.2lf", 2*Calc(ll, rr));
    return 0;
}