点和向量
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));
}