计算几何 1

点和向量

const double eps = 1e-10;
struct Point
{
	double x,y; 
	Point(){}
	Point(double _x,double _y):x(_x),y(_y){}
};
typedef Point Vector;
//直线
struct Line {
  Point s, e;
  Line() {}
  Line(Point a, Point b) {s = a, e = b;}
};

Vector operator+(const Vector& p, const Vector& q)//向量加法
{
    return Vector(p.x + q.x, p.y + q.y);
}
Vector operator-(const Vector& p, const Vector& q)
{
    return Vector(p.x-q.x, p.y-q.y);
}
Vector operator*(const Vector& p, double k)
{
    return Vector(p.x * k, p.y * k);
}
Vector operator/(const Vector& p, double k)
{
    return Vector(p.x / k, p.y / k);
}

double Dot(const Vector& p, const Vector& q) { return p.x * q.x + p.y * q.y; }//点乘
double Cross(const Vector& p, const Vector& q) { return p.x * q.y-p.y * q.x; }//叉乘
int dcmp(double x)//精度三态函数（>0,<0,=0）
{
    if (fabs(x) < eps)return 0;
    else if (x > 0)return 1;
    return -1;
}
bool operator == (const Point& a, const Point& b)//向量相等
{
    return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}
double Length(Vector a)//模
{
    return sqrt(Dot(a, a));
}
double Angle(Vector a, Vector b)//夹角，弧度制
{
    return acos(Dot(a, b) / Length(a) / Length(b));
}

//直线倾斜角
double getLineAngle(Vector a) {
	return atan2(a.y, a.x);
}
double getLineAngle(Line a) {
	return atan2(a.e.y - a.s.y, a.e.x - a.s.x);
}
double getLineAngle(Vector a) {
	return atan2(a.y, a.x);
}

Vector Rotate(Vector a, double rad)//逆时针旋转
{
    return Vector(a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
}
double Distance(Point a, Point b)//两点间距离
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double Area(Point a, Point b, Point c)//三角形面积
{
    return fabs(Cross(b - a, c - a) / 2);
}

//线段相交，不包括端点
bool Intersect(Point A, Point B, Point C, Point D)
{
    double t1 = Cross(C - A, D - A) * Cross(C - B, D - B);
    double t2 = Cross(A - C, B - C) * Cross(A - D, B - D);
    return dcmp(t1) < 0 && dcmp(t2) < 0;
}
//线段相交，包括端点
bool StrictIntersect(Point A, Point B, Point C, Point D) //线段相交，包括端点
{
    return  
        dcmp(max(A.x, B.x) - min(C.x, D.x)) >= 0
        && dcmp(max(C.x, D.x) - min(A.x, B.x)) >= 0
        && dcmp(max(A.y, B.y) - min(C.y, D.y)) >= 0
        && dcmp(max(C.y, D.y) - min(A.y, B.y)) >= 0
        && dcmp(Cross(C - A, D - A) * Cross(C - B, D - B)) <= 0
        && dcmp(Cross(A - C, B - C) * Cross(A - D, B - D)) <= 0;
}
//点到直线距离
double DistanceToLine(Point A, Point M, Point N) {//点到直线距离
    return fabs(Cross(A - M, A - N) / Distance(M, N));
}

//两直线交点
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w)
{
    Vector u = P - Q;
    double t = Cross(w, u) / Cross(v, w);
    return P + v * t;
}
Point GetLineIntersection(Line a, Line b)
{
    return GetLineIntersection(a.st, a.ed - a.st, b.st, b.ed - b.st);
}

//判断点是否在凸多边形内
Point P[1005];
int n;
bool InsidePolygon(Point A){
    double alpha=0;
    for(int i=0;i<n;i++)
        alpha+=fabs(Angle(P[i]-A,P[(i+1)%n]-A));
    return dcmp(alpha-2*pi)==0;
}
double PolygonArea()//求多边形面积(叉积和计算法)
{
    double sum = 0;
    Point Ox = Point(0, 0);
    for (int i = 0; i < n; i++)
        sum += Cross(P[i] - Ox, P[(i + 1) % n] - Ox);
    if (sum < 0)sum = -sum;
    return sum / 2;
}

random_shuffle(b+1,b+n+1);//随机排序函数 

判断点在直线左or右

//A, B:直线上一点,C:待判断关系的点
int relation(Point A, Point B, Point C)
{
    // 1 left -1 right 0 in
    int c = dcmp(Cross((B - A),(C - A)));
    if(c < 0) return 1;
    else if(c > 0) return -1;
    return 0;
}

半平面交

const double eps = 1e-10;
struct Point{
    double x, y;
    Point() {}
    Point(double _x, double _y) :x(_x), y(_y) {}
};
struct Line {
  Point st, ed;
  Line() {}
  Line(Point a, Point b) {st = a, ed = b;}
};

const int maxn = 1e3;
Point p[maxn],ans[maxn];;
Line L[maxn];
int cnt,q[maxn];

typedef Point Vector;
Vector operator+(const Vector& p, const Vector& q)//向量加法
	{return Vector(p.x + q.x, p.y + q.y);}
Vector operator-(const Vector& p, const Vector& q)
	{return Vector(p.x-q.x, p.y-q.y);}
Vector operator*(const Vector& p, double k)
	{return Vector(p.x * k, p.y * k);}
Vector operator/(const Vector& p, double k)
	{return Vector(p.x / k, p.y / k);}
double Cross(const Vector& p, const Vector& q) //叉乘
	{ return p.x * q.y-p.y * q.x; }
int dcmp(double x){//精度三态函数（>0,<0,=0）
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

double getLineAngle(Vector a) 
	{return atan2(a.y, a.x);}
double getLineAngle(Line a) 
	{return atan2(a.ed.y - a.st.y, a.ed.x - a.st.x);}

double Area(Point a, Point b, Point c)//三角形面积
   { return (Cross(b - a, c - a) / 2.0);}

Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){//两直线交点
    Vector u = P - Q;
    double t = Cross(w, u) / Cross(v, w);
    return P + v * t;
}
Point GetLineIntersection(Line a, Line b)
   { return GetLineIntersection(a.st, a.ed - a.st, b.st, b.ed - b.st);}

//排序：极角小的排前面，极角相同时，最左边的排在最后面，以便去重
bool cmp(Line a, Line b) {
 	double A =  getLineAngle(a), B = getLineAngle(b);
	if (!dcmp(A-B)) return Area(a.st, a.ed, b.ed) < 0;
	return A < B;
}
// bc的交点是否在a的右侧
bool on_right(Line& a, Line& b, Line& c){
    Point o = GetLineIntersection(b, c);
    return dcmp(Area(a.st, a.ed, o)) <= 0;
}

double half_plane_intersection(){
    sort(L, L + cnt, cmp);
    int hh = 0, tt = -1;
    for (int i = 0; i < cnt; i ++ ){
        if (i && !dcmp(getLineAngle(L[i])- getLineAngle(L[i - 1]))) continue;
        while (hh + 1 <= tt && on_right(L[i], L[q[tt - 1]], L[q[tt]])) tt -- ;
        while (hh + 1 <= tt && on_right(L[i], L[q[hh]], L[q[hh + 1]])) hh ++ ;
        q[ ++ tt] = i;
    }
    while (hh + 1 <= tt && on_right(L[q[hh]], L[q[tt - 1]], L[q[tt]])) tt -- ;
    while (hh + 1 <= tt && on_right(L[q[tt]], L[q[hh]], L[q[hh + 1]])) hh ++ ;
    q[ ++ tt] = q[hh];
    int k = 0;
    for (int i = hh; i < tt; i ++ )
        ans[k ++ ] = GetLineIntersection(L[q[i]], L[q[i + 1]]);
    double res = 0;
    for (int i = 1; i + 1 < k; i ++ )
        res += Area(ans[0], ans[i], ans[i + 1]);
    return res ;
}

凸包

namespace ConvexHull {
    bool cmp1(Point a, Point b) {
        if (fabs(a.x - b.x) < eps)return a.y < b.y;
        return a.x < b.x;
    }
    //从左下角开始逆时针排列，去除凸包边上的点
    vector<Point> Andrew_s_monotone_chain(vector<Point> P) {
        int n = P.size(), k = 0;
        vector<Point> H(2 * n);
        sort(P.begin(), P.end(), cmp1);
        for (int i = 0; i < n; i++) {
            while (k >= 2 && Cross(H[k - 1] - H[k - 2], P[i] - H[k - 2]) < eps)k--;
            H[k++] = P[i];
        }
        int t = k + 1;
        for (int i = n - 1; i > 0; i--) {
            while (k >= t && Cross(H[k - 1] - H[k - 2], P[i - 1] - H[k - 2]) < eps)k--;
            H[k++] = P[i - 1];
            H.resize(k - 1);
        }
        return H;
    }
}
using namespace ConvexHull;

另一种

int cmp(Point a,Point b){//排序，对于x，y坐标
	return (a.x!=b.x) ? (a.x<b.x):(a.y<b.y);
}

void andrew(){//先右到左遍历找凸包，在左到右找
    sort(p,p+n,cmp);
    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--;
}

int rotating_calipers(){//旋转卡壳同上
    int j=1,ans=0;
    for(int i=0;i<m;++i){
        while(Cross(ch[i+1]-ch[i],ch[j]-ch[i])<Cross(ch[i+1]-ch[i],ch[j+1]-ch[i])) j=(j+1)%m;
        ans=max(ans,max(Distance(ch[i],ch[j]),Distance(ch[i+1],ch[j])));
    }
    return ans;
}

旋转卡壳

double getmmx(){
	double ans=0;
	int siz=h.size();
	if(siz==2) return Distance(h[0],h[1]);  //只有两个点的时候不能构成三角形，所以要特判 
	int j=2;
	for(int i=0;i<siz;i++) {
		while(Cross(h[i+1]-h[i],h[j]-h[i])<Cross(h[i+1]-h[i],h[j+1]-h[i]) ) j=(j+1)%siz;  //实际上这也是在找面积，由于面积只有一个顶点，所以下一个j变小的时候就退出 
		ans=max(ans,max(Distance(h[i],h[j]),Distance(h[i+1],h[j]) ) ); //找到距离最远的一条线段，并且与当前最大值比较 
	}
	return ans;  //最后就是直径了 
}

最小圆覆盖

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define per(i,s,t) for(int i=s;i>=t;i--)
#define ms0(ar) memset(ar,0,sizeof ar)
#define ms1(ar) memset(ar,-1,sizeof ar)
#define Ri(x) scanf("%d",&x)
#define Rii(x,y) scanf("%d%d",&x,&y) 
const double eps = 1e-10;
const double pi=acos(-1.0);
struct Point{
	double x,y; 
	Point(){}
	Point(double _x,double _y):x(_x),y(_y){}
};
Point b[100010];
typedef Point Vector;
Vector operator+(const Vector& p, const Vector& q)//向量加法
{
    return Vector(p.x + q.x, p.y + q.y);
}
Vector operator-(const Vector& p, const Vector& q)
{
    return Vector(p.x-q.x, p.y-q.y);
}
Vector operator*(const Vector& p, double k)
{
    return Vector(p.x * k, p.y * k);
}
Vector operator/(const Vector& p, double k)
{
    return Vector(p.x / k, p.y / k);
}
double Dot(const Vector& p, const Vector& q) { return p.x * q.x + p.y * q.y; }//点乘
double Cross(const Vector& p, const Vector& q) { return p.x * q.y-p.y * q.x; }//叉乘
int dcmp(double x)//精度三态函数（>0,<0,=0）
{
    if (fabs(x) < eps)return 0;
    else if (x > 0)return 1;
    return -1;
}
bool operator == (const Point& a, const Point& b)//向量相等
{
    return dcmp(a.x - b.x) == 0 && dcmp(a.y - b.y) == 0;
}
double Length(Vector a)//模
{
    return sqrt(fabs(Dot(a, a)));
}
double Angle(Vector a, Vector b)//夹角，弧度制
{
    return acos(Dot(a, b) / Length(a) / Length(b));
}
Vector Rotate(Vector a, double rad)//逆时针旋转
{
    return Vector(a.x * cos(rad) - a.y * sin(rad), a.x * sin(rad) + a.y * cos(rad));
}
double Distance(Point a, Point b)//两点间距离
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double Area(Point a, Point b, Point c)//三角形面积
{
    return fabs(Cross(b - a, c - a) / 2);
}
bool Intersect(Point A, Point B, Point C, Point D)//线段相交，不包括端点
{
    double t1 = Cross(C - A, D - A) * Cross(C - B, D - B);
    double t2 = Cross(A - C, B - C) * Cross(A - D, B - D);
    return dcmp(t1) < 0 && dcmp(t2) < 0;
}

bool StrictIntersect(Point A, Point B, Point C, Point D) //线段相交，包括端点
{
    return  
        dcmp(max(A.x, B.x) - min(C.x, D.x)) >= 0
        && dcmp(max(C.x, D.x) - min(A.x, B.x)) >= 0
        && dcmp(max(A.y, B.y) - min(C.y, D.y)) >= 0
        && dcmp(max(C.y, D.y) - min(A.y, B.y)) >= 0
        && dcmp(Cross(C - A, D - A) * Cross(C - B, D - B)) <= 0
        && dcmp(Cross(A - C, B - C) * Cross(A - D, B - D)) <= 0;
}
double DistanceToLine(Point A, Point M, Point N) {//点到直线距离
    return fabs(Cross(A - M, A - N) / Distance(M, N));
}
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w)//两直线交点
{
    Vector u = P - Q;
    double t = Cross(w, u) / Cross(v, w);
    return P + v * t;
}
int n;
double R;
Point res;
double sqr(double x){return x*x;}
bool incircle(Point x){
	return (Distance(x,res)-R)<eps;
}
Point getres(Point a,Point b,Point c){
	Point p=(a+b)/2;    //ad中点 
	Point q=(a+c)/2;    //ac中点 
	Vector v=Rotate(b-a,pi/2.0),w=Rotate(c-a,pi/2.0);   //中垂线的方向向量 
	if (dcmp(Cross(v,w))==0){    //平行 
		if (dcmp(Distance(a,b)+Distance(b,c)-Distance(a,c))==0)
		   return (a+c)/2;
		if (dcmp(Distance(b,a)+Distance(a,c)-Distance(b,c))==0)
		   return (b+c)/2;
		if (dcmp(Distance(a,c)+Distance(c,b)-Distance(a,b))==0)
		   return (a+b)/2;
	}
	return GetLineIntersection(p,v,q,w);
}
int main(){
    scanf("%d",&n);
    int i,j,k;
    for(i=1;i<=n;i++)
        scanf("%lf%lf",&b[i].x,&b[i].y);
    random_shuffle(b+1,b+n+1);
    R=0;
    for(i=1;i<=n;i++)
        if(!incircle(b[i])){
            res=b[i];
            R=0.0;
            for(j=1;j<i;j++)
            if(!incircle(b[j])){
                res=(b[i]+b[j])/2;
                R=Distance(res,b[i]);
                for(k=1;k<j;k++)
                    if(!incircle(b[k])){
                        res=getres(b[i],b[j],b[k]);
                        R=Distance(b[i],res);
                    }
            }
    }
    printf("%.10lf\n%.10lf %.10lf",R,res.x,res.y);
}

自适应辛普森

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const double eps = 1e-12;

double f(double x)
{
    return sin(x) / x;//具体函数
}

double simpson(double l, double r)
{
    auto mid = (l + r) / 2;
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;
}

double asr(double l, double r, double s)
{
    auto mid = (l + r) / 2;
    auto left = simpson(l, mid), right = simpson(mid, r);
    if (fabs(left + right - s) < eps) return left + right;
    return asr(l, mid, left) + asr(mid, r, right);
}

int main()
{
    double l, r;
    scanf("%lf%lf", &l, &r);
    printf("%lf\n", asr(l, r, simpson(l, r)));
    return 0;
}

圆的面积并

#include<bits/stdc++.h> 
using namespace std;
typedef long long ll;
typedef double db;
#define rep(i,s,t) for(int i=s;i<=t;i++)
#define per(i,s,t) for(int i=s;i>=t;i--)
#define ms0(ar) memset(ar,0,sizeof ar)
#define ms1(ar) memset(ar,-1,sizeof ar)
#define Ri(x) scanf("%d",&x)
#define Rii(x,y) scanf("%d%d",&x,&y)

const int maxn=1005;

struct Point{
	double x,y; Point(){}
	Point(double x,double y):x(x),y(y){}
	bool operator < (const Point &p)const{return x<p.x;}
}q[maxn];

struct circle{
	double x,y,r;
	bool operator < (const circle &p)const{return r>p.r;}
}a[maxn];

const double eps = 1e-8, inf = 1e20, Pi = acos(-1);
int n;

#define sqr(x) (x)*(x)

double f(double x){
	int top=0;
	for(int i=1;i<=n;i++) if(fabs(x-a[i].x)<a[i].r){
		double len=sqrt(sqr(a[i].r)-sqr(x-a[i].x));
		q[++top]=Point(a[i].y-len,a[i].y+len);
	}
	double l=-inf,r=-inf,ret=0;
	sort(q+1,q+1+top);
	for(int i=1;i<=top;i++)
		if(q[i].x>r) ret+=r-l,l=q[i].x,r=q[i].y;
		else r=max(r,q[i].y);
	return ret+=r-l;
}

double Simpson(double L,double M,double R,double fL,double fM,double fR,int dep){
	double M1=(L+M)/2,M2=(M+R)/2,fM1=f(M1),fM2=f(M2);
	double g1=(M-L)*(fL+4*fM1+fM)/6, g2=(R-M)*(fM+4*fM2+fR)/6, g=(R-L)*(fL+4*fM+fR)/6;
	if(dep>11&&fabs(g-g1-g2)<1e-8) return g;
	return Simpson(L,M1,M,fL,fM1,fM,dep+1)+Simpson(M,M2,R,fM,fM2,fR,dep+1);
}

inline double dist(circle &a,circle &b){return sqrt(sqr(a.x-b.x)+sqr(a.y-b.y));}

void prework(){
	int cnt=0;
	sort(a+1,a+1+n);
	rep(i,1,n) {
		rep(j,1,cnt)  
			if(a[i].r+dist(a[i],a[j])<=a[j].r) 
				goto he;
		a[++cnt]=a[i]; he:;
	}
	n=cnt;
}

int main(){
	Ri(n);
	double l=inf,r=-inf,mid;
	rep(i,1,n) 
		scanf("%lf%lf%lf",&a[i].x,&a[i].y,&a[i].r),
		l=min(l,a[i].x-a[i].r),
		r=max(r,a[i].x+a[i].r);
	
	mid=(l+r)/2; prework();
	
	printf("%.3f",Simpson(l,mid,r,f(l),f(mid),f(r),0));
}