计算几何进阶

最新推荐文章于 2024-07-23 17:36:37 发布

胡耀文

最新推荐文章于 2024-07-23 17:36:37 发布

阅读量164

点赞数

文章标签：计算几何算法数据结构自适应希普森积分凸包

本文链接：https://blog.csdn.net/supreme567/article/details/120648197

版权

基础知识

2983. 玩具

在这里插入图片描述

#include<bits/stdc++.h>
using namespace std;
#define int long long
typedef pair<int, int> PII;
#define x first
#define y second
//#define d lld
int n;
PII a[5050],b[5050];//a在下面，b在上面
int ans[5050];
int cross(int x1,int y1,int x2,int y2){
    return x1*y2-x2*y1;
}
int area(PII a,PII b,PII c){
    return cross(b.x-a.x,b.y-a.y , c.x-a.x,c.y-a.y);
}
int find(int x,int y){
    int l=0,r=n;
    while(l<r){
        int mid=l+r>>1;
        if(area(a[mid],b[mid],{x,y})>0)r=mid;
        else l=mid+1;
    }
    return r;
}
signed main(){
    bool is_first=1;
    while(scanf("%lld",&n),n){
        int m,x1,y1,x2,y2;
        scanf("%lld%lld%lld%lld%lld",&m,&x1,&y1,&x2,&y2);//左上角，右下角
        for(int i=0;i<n;i++){
            int u,l;
            cin>>u>>l;
            b[i]={u,y1};//右上角
            a[i]={l,y2};//左下角
        }
        a[n]={x2,y2};
        b[n]={x2,y1};
        memset(ans,0,sizeof ans);
        while(m--){
            int x,y;
            cin>>x>>y;
            ans[find(x,y)]++;
        }
        if(is_first)is_first=0;
        else puts("");
        
        for(int i=0;i<=n;i++)printf("%lld: %lld\n",i,ans[i]);
        
    }
}

2984. 线段

在这里插入图片描述

#include<bits/stdc++.h>
using namespace std;
typedef pair<double, double> PDD;
#define x first
#define y second
int T,n,m;
const int N = 220;
const double eps = 1e-8;
PDD q[N],a[N],b[N];
bool cmp(double x,double y){
    if(fabs(x-y)<eps)return 0;
    if(x<y)return -1;
    else return 1;
}
int sign(double t){
    if(fabs(t)<eps)return 0;
    if(t>0)return 1;
    else return -1;
}
double cross(double x1,double y1,double x2,double y2){
    return x1*y2-x2*y1;
}
double area(PDD a,PDD b,PDD c){
    return cross(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
}
bool check(){
    for(int i=0;i<n*2;i++){
        for(int j=i+1;j<n*2;j++){
            if( !cmp(q[i].x,q[j].x) && !cmp(q[i].y,q[j].y) )continue;
            bool ok=1;
            for(int k=0;k<n;k++){
                int zuo=sign(area(q[j],q[i],a[k]));
                int you=sign(area(q[j],q[i],b[k]));
                if(zuo*you>0){ok=false;break;}
            }
            if(ok)return 1;
        }
    }
    return 0;
}

signed main(){
    cin>>T;
    while(T--){
        scanf("%d",&n);
        for(int i=0,k=0;i<n;i++){
            double x1,y1,x2,y2;
            scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
            q[k++]={x1,y1},q[k++]={x2,y2};
            a[i]  ={x1,y1},b[i]  ={x2,y2};
        }
        if(check())puts("Yes!");
        else puts("No!");
    }
}

在这里插入图片描述
思路就是判断遍历所有的端点，然后对于每一组端点，for循环遍历一下所有的的线段，看看是不是存在一组叉积ija ijb是同号，同号表示在同一边，向量ij没有穿过ab
如果是异号，表示向量ab分布在向量ij的两边，
只有存在两个端点能穿过所有的线段，才能check=true

凸包

1401. 围住奶牛

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 10010;

int n;
PDD q[N];
int stk[N];
bool used[N];

double get_dist(PDD a, PDD b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

PDD operator-(PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    return cross(b - a, c - a);
}

double andrew()
{
    sort(q, q + n);
    int top = 0;
    for (int i = 0; i < n; i ++ )
    {
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top -- ]] = false;
            else top -- ;
        }
        stk[ ++ top] = i;
        used[i] = true;
    }
    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }

    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    double res = andrew();
    printf("%.2lf\n", res);

    return 0;
}

上面是以下右为最大优先级的顺时针凸包
下面是从左下角为最大优先级的逆时针凸包

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 10010;

int n;
PDD q[N];
int stk[N];
bool used[N];

double get_dist(PDD a, PDD b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

PDD operator-(PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    return cross(b - a, c - a);
}

double andrew()
{
    sort(q, q + n);
    int top = 0;
    for (int i = 0; i < n; i ++ )
    {
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top -- ]] = false;
            else top -- ;
        }
        stk[ ++ top] = i;
        used[i] = true;
    }
    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
            top -- ;
        stk[ ++ top] = i;
    }

    double res = 0;
    
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    
    for(int i=2;i<=top;i++){
        auto t=q[stk[i]];
        cout<<t.x<<" "<<t.y<<endl;
    }
    
    return res;
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    double res = andrew();
    //printf("%.2lf\n", res);

    return 0;
}

2935. 信用卡凸包

在这里插入图片描述

多边形的外角和等于360°
也就是说其实算出凸包后，每个弧的周长等于2派R

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 40010;
const double pi = acos(-1);

int n, cnt;
PDD q[N];
int stk[N], top;
bool used[N];

PDD rotate(PDD a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

PDD operator- (PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    return cross(b - a, c - a);
}

double get_dist(PDD a, PDD b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

double andrew()
{
    sort(q, q + cnt);
    for (int i = 0; i < cnt; i ++ )
    {
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
        {
            // 凸包边界上的点即使被从栈中删掉，也不能删掉used上的标记
            if (area(q[stk[top - 1]], q[stk[top]], q[i]) < 0)
                used[stk[top -- ]] = false;
            else top -- ;
        }
        stk[ ++ top] = i;
        used[i] = true;
    }
    used[0] = 0;
    for (int i = cnt - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 1]], q[stk[top]], q[i]) <= 0)
            top -- ;
        stk[ ++ top] = i;
    }

    double res = 0;
    for (int i = 2; i <= top; i ++ )
        res += get_dist(q[stk[i - 1]], q[stk[i]]);
    return res;
}

int main()
{
    scanf("%d", &n);
    double a, b, r;
    scanf("%lf%lf%lf", &a, &b, &r);
    a = a / 2 - r, b = b / 2 - r;
    int dx[] = {1, 1, -1, -1}, dy[] = {1, -1, -1, 1};
    while (n -- )
    {
        double x, y, z;
        scanf("%lf%lf%lf", &x, &y, &z);
        for (int i = 0; i < 4; i ++ )
        {
            auto t = rotate({dx[i] * b, dy[i] * a}, -z);
            q[cnt ++ ] = {x + t.x, y + t.y};
        }
    }

    double res = andrew();
    printf("%.2lf\n", res + 2 * pi * r);

    return 0;
}

半平面交

2803. 凸多边形

在这里插入图片描述

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

#define x first
#define y second

using namespace std;

typedef pair<double,double> PDD;
const int N = 510;
const double eps = 1e-8;

int cnt;
struct Line
{
    PDD st, ed;
}line[N];
PDD pg[N], ans[N];
int q[N];

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

double get_angle(const Line& a)
{
    return atan2(a.ed.y - a.st.y, a.ed.x - a.st.x);
}

PDD operator-(PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

double cross(PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double area(PDD a, PDD b, PDD c)
{
    return cross(b - a, c - a);
}

bool cmp(const Line& a, const Line& b)
{
    double A = get_angle(a), B = get_angle(b);
    if (!dcmp(A, B)) return area(a.st, a.ed, b.ed) < 0;
    return A < B;
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    auto u = p - q;
    double t = cross(w, u) / cross(v, w);
    return {p.x + v.x * t, p.y + v.y * t};
}

PDD get_line_intersection(Line a, Line b)
{
    return get_line_intersection(a.st, a.ed - a.st, b.st, b.ed - b.st);
}

// bc的交点是否在a的右侧
bool on_right(Line& a, Line& b, Line& c)
{
    auto o = get_line_intersection(b, c);
    return sign(area(a.st, a.ed, o)) <= 0;
}

double half_plane_intersection()
{
    sort(line, line + cnt, cmp);
    int hh = 0, tt = -1;
    for (int i = 0; i < cnt; i ++ )
    {
        if (i && !dcmp(get_angle(line[i]), get_angle(line[i - 1]))) continue;
        while (hh + 1 <= tt && on_right(line[i], line[q[tt - 1]], line[q[tt]])) tt -- ;
        while (hh + 1 <= tt && on_right(line[i], line[q[hh]], line[q[hh + 1]])) hh ++ ;
        q[ ++ tt] = i;
    }
    while (hh + 1 <= tt && on_right(line[q[hh]], line[q[tt - 1]], line[q[tt]])) tt -- ;
    while (hh + 1 <= tt && on_right(line[q[tt]], line[q[hh]], line[q[hh + 1]])) hh ++ ;

    q[ ++ tt] = q[hh];//现在半平面交处理完毕，双端队列q【】存的是所有留下来的直线编号
    int k = 0;
    for (int i = hh; i < tt; i ++ )//现在要求所有线的交点，存进ans【】
        ans[k ++ ] = get_line_intersection(line[q[i]], line[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 / 2;
}

int main()
{
    int n, m;
    scanf("%d", &n);
    while (n -- )
    {
        scanf("%d", &m);
        for (int i = 0; i < m; i ++ ) scanf("%lf%lf", &pg[i].x, &pg[i].y);
        for (int i = 0; i < m; i ++ )
            line[cnt ++ ] = {pg[i], pg[(i + 1) % m]};
    }
    double res = half_plane_intersection();
    printf("%.3lf\n", res);

    return 0;
}

2957. 赛车

在这里插入图片描述

#include<bits/stdc++.h>
using namespace std;
#define x first
#define y second
typedef long double LD;
typedef pair<int, int> PII;
typedef pair< LD , LD> PDD;
const int N = 1e5+10;
const LD eps=1e-18;
int n,cnt;
struct Line{
    PDD st,ed;
    vector<int>ids;
}line[N];

int k[N];
int v[N];
int q[N];
int ans[N];
int sign(LD x)  // 符号函数
{
    if (fabs(x) < eps) return 0;  // x为0，则返回0
    if (x < 0) return -1;  // x为负数，则返回-1
    return 1;  // x为正数，则返回1
}

int dcmp(LD x, LD y)  // 比较两数大小
{
    if (fabs(x - y) < eps) return 0;  // x == y, 返回0
    if (x < y) return -1;  // x < y, 返回-1
    return 1;  // x > y, 返回1
}


PDD operator+ (PDD a, PDD b)  // 向量加法
{
    return {a.x + b.x, a.y + b.y};
}

PDD operator- (PDD a, PDD b)  //  向量减法
{
    return {a.x - b.x, a.y - b.y};
}

PDD operator* (PDD a, double t)  // 向量数乘
{
    return {a.x * t, a.y * t};
}

PDD operator/ (PDD a, double t)  // 向量除以常数
{
    return {a.x / t, a.y / t};
}

long double operator* (PDD a, PDD b)  // 外积、叉积
{
    return a.x * b.y - a.y * b.x;
}

long double operator& (PDD a, PDD b)  // 内积、点积
{
    return a.x * b.x + a.y * b.y;
}

long double area(PDD a, PDD b, PDD c)  // 以a, b, c为顶点的有向三角形面积
{
    return (b - a) * (c - a);
}

double get_len(PDD a)  // 求向量长度
{
    return sqrt(a & a);
}

double get_dist(PDD a, PDD b)  // 求两个点之间的距离
{
    return get_len(b - a);
}

double project(PDD a, PDD b, PDD c)  // 求向量ac在向量ab上的投影
{
    return ((c - a) & (b - a)) / get_len(b - a);
}

PDD rotate(PDD a, double b)  // 向量a逆时针旋转角度b
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

PDD norm(PDD a)  // 矩阵标准化（将长度变成1）
{
    return a / get_len(a);
}

bool on_segment(PDD p, PDD a, PDD b)  // 点p是否在线段ab上（包含端点a、b）
{
    return !sign((p - a) * (p - b)) && sign((p - a) & (p - b)) <= 0;
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)  // 求两直线交点：p + vt, q + wt
{
    auto u = p - q;
    auto t = (w * u) / (v * w);
    return {p.x+v.x*t , p.y+v.y*t};
}
PDD get_line_intersection(Line&a , Line&b){
    return get_line_intersection(a.st,a.ed-a.st,b.st,b.ed-b.st);
}

LD get_angle(const Line&a){
    return atan2(a.ed.y-a.st.y , a.ed.x-a.st.x);
}

bool on_right(Line&a,Line&b,Line&c){//判断bc的交点是否在a的右侧
    auto o=get_line_intersection(b,c);
    return area(a.st,a.ed,o)<0;
}

bool cmp(const Line&a,const Line&b){
    LD A=get_angle(a);
    LD B=get_angle(b);
    if(!dcmp(A,B))return area(a.st,a.ed,b.ed)<0;
    return A<B;
}

void half_plane_intersection(){
    sort(line,line+cnt,cmp);
    int hh=0,tt=-1;
    for(int i=0;i<cnt;i++){
        if(i&&!dcmp(get_angle(line[i-1]),get_angle(line[i])) )continue;
        while(hh+1<=tt&&on_right(line[i],line[q[tt-1]],line[q[tt]]))tt--;
        while(hh+1<=tt&&on_right(line[i],line[q[hh]],line[q[hh+1]]))hh++;
        q[++tt]=i;
    }
    while(hh+1<=tt&&on_right(line[q[hh]],line[q[tt-1]],line[q[tt]]))tt--;
    while(hh+1<=tt&&on_right(line[q[tt]],line[q[hh]],line[q[hh+1]]))hh++;
    
    int k=0;
    for(int i=hh;i<=tt;i++){
        for(auto id:line[q[i]].ids)ans[k++]=id;
    }
    sort(ans,ans+k);
    printf("%d\n",k);
    for(int i=0;i<k;i++)printf("%d ",ans[i]);
}

signed main(){
    //考虑到某些赛车会重复的问题（比如2和8号赛车），但是输出的时候2和8都要输出
    map<PII,vector<int>>ids;//为什么要用PII而不用PDD呢，
    //因为PDD在存储和判断的时候会有一些问题，而整数int完全不会有这种问题
    scanf("%d", &n);
    for(int i=1;i<=n;i++)scanf("%d",&k[i]);
    for(int i=1;i<=n;i++)scanf("%d",&v[i]);
    for(int i=1;i<=n;i++)ids[{k[i] , v[i]}].push_back(i);
    line[cnt++]={{0,10000},{0,0}};
    line[cnt++]={{0,0},{10000,0}};
    for(auto&[k,v]:ids)line[cnt++]={{0,k.x},{1,k.x+k.y} ,v};
    
    half_plane_intersection();
}

最小圆覆盖

3028. 最小圆覆盖

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 100010;
const double eps = 1e-12;
const double PI = acos(-1);

int n;
PDD q[N];
struct Circle
{
    PDD p;
    double r;
};

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

PDD operator- (PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

PDD operator+ (PDD a, PDD b)
{
    return {a.x + b.x, a.y + b.y};
}

PDD operator* (PDD a, double t)
{
    return {a.x * t, a.y * t};
}

PDD operator/ (PDD a, double t)
{
    return {a.x / t, a.y / t};
}

double operator* (PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

PDD rotate(PDD a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

double get_dist(PDD a, PDD b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    auto u = p - q;
    double t = w * u / (v * w);
    return p + v * t;
}

pair<PDD, PDD> get_line(PDD a, PDD b)
{
    return {(a + b) / 2, rotate(b - a, PI / 2)};
}

Circle get_circle(PDD a, PDD b, PDD c)
{
    auto u = get_line(a, b), v = get_line(a, c);
    auto p = get_line_intersection(u.x, u.y, v.x, v.y);
    return {p, get_dist(p, a)};
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    random_shuffle(q, q + n);

    Circle c({q[0], 0});
    for (int i = 1; i < n; i ++ )
        if (dcmp(c.r, get_dist(c.p, q[i])) < 0)
        {
            c = {q[i], 0};
            for (int j = 0; j < i; j ++ )
                if (dcmp(c.r, get_dist(c.p, q[j])) < 0)
                {
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for (int k = 0; k < j; k ++ )
                        if (dcmp(c.r, get_dist(c.p, q[k])) < 0)
                            c = get_circle(q[i], q[j], q[k]);
                }
        }

    printf("%.10lf\n", c.r);
    printf("%.10lf %.10lf\n", c.p.x, c.p.y);
    return 0;
}

2785. 信号增幅仪

在这里插入图片描述

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 50010;
const double eps = 1e-12;
const double PI = acos(-1);

int n;
PDD q[N];
struct Circle
{
    PDD p;
    double r;
};

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

PDD operator+ (PDD a, PDD b)
{
    return {a.x + b.x, a.y + b.y};
}

PDD operator- (PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

PDD operator* (PDD a, double t)
{
    return {a.x * t, a.y * t};
}

PDD operator/ (PDD a, double t)
{
    return {a.x / t, a.y / t};
}

double operator* (PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

PDD rotate(PDD a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

double get_dist(PDD a, PDD b)
{
    double dx = a.x - b.x;
    double dy = a.y - b.y;
    return sqrt(dx * dx + dy * dy);
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    auto u = p - q;
    double t = w * u / (v * w);
    return p + v * t;
}

pair<PDD, PDD> get_line(PDD a, PDD b)
{
    return {(a + b) / 2, rotate(b - a, PI / 2)};
}

Circle get_circle(PDD a, PDD b, PDD c)
{
    auto u = get_line(a, b), v = get_line(a, c);
    auto p = get_line_intersection(u.x, u.y, v.x, v.y);
    return {p, get_dist(p, a)};
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    double a, p;
    scanf("%lf%lf", &a, &p);
    for (int i = 0; i < n; i ++ )
    {
        q[i] = rotate(q[i], a / 180 * PI);
        q[i].x /= p;
    }

    random_shuffle(q, q + n);
    Circle c({q[0], 0});
    for (int i = 1; i < n; i ++ )
        if (dcmp(c.r, get_dist(c.p, q[i])) < 0)
        {
            c = {q[i], 0};
            for (int j = 0; j < i; j ++ )
                if (dcmp(c.r, get_dist(c.p, q[j])) < 0)
                {
                    c = {(q[i] + q[j]) / 2, get_dist(q[i], q[j]) / 2};
                    for (int k = 0; k < j; k ++ )
                        if (dcmp(c.r, get_dist(c.p, q[k])) < 0)
                            c = get_circle(q[i], q[j], q[k]);
                }
        }

    printf("%.3lf\n", c.r);

    return 0;
}

三维计算几何基础与三维凸包

2119. 最佳包裹

在这里插入图片描述

1. 三维向量表示(x, y, z)
2. 向量加减法、数乘运算，与二维相同
3. 模长 |A| = sqrt(x * x + y * y + z * z)
4. 点积
    (1) 几何意义：A·B = |A| * |B| * cos(C)
    (2) 代数求解：(x1, y1, z1) · (x2, y2, z2) = (x1x2, y1y2, z1z2);
5. 叉积
    (1) 几何意义：AxB = |A| * |B| * sin(C)，方向：右手定则
    (2) 代数求解：AxB = (y1z2 - z1y2, z1x2 - x1z2, x1y2 - x2y1)
6. 如何求平面法向量
    任取平面上两个不共线的向量A、B：AxB
7. 判断点D是否在平面里
    任取平面上两个不共线的向量A、B：先求法向量C = AxB，然后求平面上任意一点到D的向量E与C的点积，判断点积是否为0。
8. 求点D到平面的距离
    任取平面上两个不共线的向量A、B：先求法向量C = AxB。然后求平面上任意一点到D的向量E在C上的投影长度即可。即：E·C / |C|
9. 多面体欧拉定理
    顶点数 - 棱长数 + 表面数 = 2
10. 三维凸包

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

using namespace std;

const int N = 210;
const double eps = 1e-10;

int n, m;
bool g[N][N];

double rand_eps()
{
    return ((double)rand() / RAND_MAX - 0.5) * eps;
}

struct Point
{
    double x, y, z;
    void shake()
    {
        x += rand_eps(), y += rand_eps(), z += rand_eps();
    }
    Point operator- (Point t)
    {
        return {x - t.x, y - t.y, z - t.z};
    }
    double operator& (Point t)
    {
        return x * t.x + y * t.y + z * t.z;
    }
    Point operator* (Point t)
    {
        return {y * t.z - t.y * z, z * t.x - x * t.z, x * t.y - y * t.x};
    }
    double len()
    {
        return sqrt(x * x + y * y + z * z);
    }
}q[N];
struct Plane
{
    int v[3];
    Point norm()  // 法向量
    {
        return (q[v[1]] - q[v[0]]) * (q[v[2]] - q[v[0]]);
    }
    double area()  // 求面积
    {
        return norm().len() / 2;
    }
    bool above(Point a)
    {
        return ((a - q[v[0]]) & norm()) >= 0;
    }
}plane[N], np[N];

void get_convex_3d()
{
    plane[m ++ ] = {0, 1, 2};
    plane[m ++ ] = {2, 1, 0};
    for (int i = 3; i < n; i ++ )
    {
        int cnt = 0;
        for (int j = 0; j < m; j ++ )
        {
            bool t = plane[j].above(q[i]);
            if (!t) np[cnt ++ ] = plane[j];
            for (int k = 0; k < 3; k ++ )
                g[plane[j].v[k]][plane[j].v[(k + 1) % 3]] = t;
        }
        for (int j = 0; j < m; j ++ )
            for (int k = 0; k < 3; k ++ )
            {
                int a = plane[j].v[k], b = plane[j].v[(k + 1) % 3];
                if (g[a][b] && !g[b][a])
                    np[cnt ++ ] = {a, b, i};
            }
        m = cnt;
        for (int j = 0; j < m; j ++ ) plane[j] = np[j];
    }
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ )
    {
        scanf("%lf%lf%lf", &q[i].x, &q[i].y, &q[i].z);
        q[i].shake();
    }
    get_convex_3d();

    double res = 0;
    for (int i = 0; i < m; i ++ )
        res += plane[i].area();
    printf("%lf\n", res);
    return 0;
}

y总的代码过不了了现在，要1e-12改成1e-10才可以过
哎，计算几何的精度问题真的好难调，出题人本来不想在这里卡我们的，但是c++的浮点数精度运算的时候却会有一些偏差

旋转卡壳

2938. 周游世界

在这里插入图片描述

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

#define x first
#define y second

using namespace std;

typedef pair<int, int> PII;
const int N = 50010;


int n;
PII q[N];
int stk[N], top;
bool used[N];

PII operator- (PII a, PII b)
{
    return {a.x - b.x, a.y - b.y};
}

int operator* (PII a, PII b)
{
    return a.x * b.y - a.y * b.x;
}

int area(PII a, PII b, PII c)
{
    return (b - a) * (c - a);
}

int get_dist(PII a, PII b)
{
    int dx = a.x - b.x;
    int dy = a.y - b.y;
    return dx * dx + dy * dy;
}

void get_convex()
{
    sort(q, q + n);
    for (int i = 0; i < n; i ++ )
    {
        while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
        {
            if (area(q[stk[top - 2]], q[stk[top - 1]], q[i]) < 0)
                used[stk[ -- top]] = false;
            else top -- ;
        }
        stk[top ++ ] = i;
        used[i] = true;
    }

    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && area(q[stk[top - 2]], q[stk[top - 1]], q[i]) <= 0)
            top -- ;
        stk[top ++ ] = i;
    }
    top -- ;
}

int rotating_calipers()
{
    if (top <= 2) return get_dist(q[0], q[n - 1]);

    int res = 0;
    for (int i = 0, j = 2; i < top; i ++ )
    {
        auto d = q[stk[i]], e = q[stk[i + 1]];
        while (area(d, e, q[stk[j]]) < area(d, e, q[stk[j + 1]])) j = (j + 1) % top;
        res = max(res, max(get_dist(d, q[stk[j]]), get_dist(e, q[stk[j]])));
    }
    return res;
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%d%d", &q[i].x, &q[i].y);
    get_convex();
    printf("%d\n", rotating_calipers());

    return 0;
}

注意，当所有点共线时，求得的凸包只包含起点，但并不影响求最远点对，因为代码第70行特判了所有点共线的情况。另外如果打算求所有点共线的凸包，可以将代码第60行的小于等于号变成小于号。

2142. 最小矩形覆盖

在这里插入图片描述

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 50010;
const double eps = 1e-12, INF = 1e20;
const double PI = acos(-1);

int n;
PDD q[N];
PDD ans[N];
double min_area = INF;
int stk[N], top;
bool used[N];

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

PDD operator+ (PDD a, PDD b)
{
    return {a.x + b.x, a.y + b.y};
}

PDD operator- (PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

PDD operator* (PDD a, double t)
{
    return {a.x * t, a.y * t};
}

PDD operator/ (PDD a, double t)
{
    return {a.x / t, a.y / t};
}

double operator* (PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double operator& (PDD a, PDD b)
{
    return a.x * b.x + a.y * b.y;
}

double area(PDD a, PDD b, PDD c)
{
    return (b - a) * (c - a);
}

double get_len(PDD a)
{
    return sqrt(a & a);
}

double project(PDD a, PDD b, PDD c)
{
    return ((b - a) & (c - a)) / get_len(b - a);
}

PDD norm(PDD a)
{
    return a / get_len(a);
}

PDD rotate(PDD a, double b)
{
    return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)};
}

void get_convex()
{
    sort(q, q + n);
    for (int i = 0; i < n; i ++ )
    {
        while (top >= 2 && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) >= 0)
            used[stk[ -- top]] = false;
        stk[top ++ ] = i;
        used[i] = true;
    }
    used[0] = false;
    for (int i = n - 1; i >= 0; i -- )
    {
        if (used[i]) continue;
        while (top >= 2 && sign(area(q[stk[top - 2]], q[stk[top - 1]], q[i])) >= 0)
            top -- ;
        stk[top ++ ] = i;
    }
    reverse(stk, stk + top);
    top -- ;
}

void rotating_calipers()
{
    for (int i = 0, a = 2, b = 1, c = 2; i < top; i ++ )
    {
        auto d = q[stk[i]], e = q[stk[i + 1]];
        while (dcmp(area(d, e, q[stk[a]]), area(d, e, q[stk[a + 1]])) < 0) a = (a + 1) % top;
        while (dcmp(project(d, e, q[stk[b]]), project(d, e, q[stk[b + 1]])) < 0) b = (b + 1) % top;
        if (!i) c = a;
        while (dcmp(project(d, e, q[stk[c]]), project(d, e, q[stk[c + 1]])) > 0) c = (c + 1) % top;
        auto x = q[stk[a]], y = q[stk[b]], z = q[stk[c]];
        auto h = area(d, e, x) / get_len(e - d);
        auto w = ((y - z) & (e - d)) / get_len(e - d);
        if (h * w < min_area)
        {
            min_area = h * w;
            ans[0] = d + norm(e - d) * project(d, e, y);
            ans[3] = d + norm(e - d) * project(d, e, z);
            auto u = norm(rotate(e - d, -PI / 2));
            ans[1] = ans[0] + u * h;
            ans[2] = ans[3] + u * h;
        }
    }
}

int main()
{
    scanf("%d", &n);
    for (int i = 0; i < n; i ++ ) scanf("%lf%lf", &q[i].x, &q[i].y);
    get_convex();
    rotating_calipers();

    int k = 0;
    for (int i = 1; i < 4; i ++ )
        if (dcmp(ans[i].y, ans[k].y) < 0 || !dcmp(ans[i].y, ans[k].y) && dcmp(ans[i].x, ans[k].x))
            k = i;

    printf("%.5lf\n", min_area);
    for (int i = 0; i < 4; i ++, k ++ )
    {
        auto x = ans[k % 4].x, y = ans[k % 4].y;
        if (!sign(x)) x = 0;
        if (!sign(y)) y = 0;
        printf("%.5lf %.5lf\n", x, y);
    }

    return 0;
}

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define x first
#define y second
typedef pair<double,double>PDD;
const int N = 50010;
const double eps = 1e-12,INF=1e20;
const double PI = acos(-1);
int n;
PDD q[N];
PDD ans[N];
double min_area = INF;
int stk[N],top;
bool used[N];
int sign(double x){
    if(fabs(x)<eps)return 0;
    if(x<0)return -1;
    else return 1;
}

int dcmp(double a,double b){
    if(fabs(a-b)<eps)return 0;
    if(a<b)return -1;
    return 1;
}

PDD operator+(PDD a,PDD b){
    return {a.x+b.x,a.y+b.y};
}
PDD operator-(PDD a,PDD b){
    return {a.x-b.x,a.y-b.y};
}
PDD operator*(PDD a,double t){
    return {a.x*t,a.y*t};
}
PDD operator*(double t,PDD a){
    return {a.x*t,a.y*t};
}
PDD operator/(PDD a,double t){
    return {a.x/t,a.y/t};
}
double operator*(PDD a,PDD b){
    return a.x*b.y-a.y*b.x;
}
double operator&(PDD a,PDD b){
    return a.x*b.x+a.y*b.y;
}
double area(PDD a,PDD b,PDD c){
    return (b-a)*(c-a);
}
double get_len(PDD a){
    return sqrt(a & a);
}
double project(PDD a,PDD b,PDD c){//求向量ac在向量ab上的投影
    return ((b-a)&(c-a))/get_len(b-a);
}
PDD norm(PDD a){
    return a/get_len(a);
}
PDD rotate(PDD a,double b){//将向量a顺时针旋转b的角度
    return {a.x*cos(b)+a.y*sin(b) , -a.x*sin(b)+a.y*cos(b)};
}
void get_convex(){
    sort(q,q+n);
    for(int i=0;i<n;i++){
        while(top>=2 && sign(area(q[stk[top-2]] , q[stk[top-1]] , q[i])) >=0){
            used[stk[--top]]=0;
        }
        stk[top++]=i;
        used[i]=1;
    }
    used[0]=0;
    for(int i=n-1;i>=0;i--){
        if(used[i])continue;
        while(top>=2 && sign(area(q[stk[top-2]] , q[stk[top-1]] , q[i])) >=0){
            top--;
        }
        stk[top++]=i;
    }
    reverse(stk,stk+top);
    top--;
    //for(auto t:stk)cout<<t<<endl;
}

void rotating_calipers(){
    for(int i=0,a=2,b=1,c=2;i<top;i++){
        auto d=q[stk[i]] , e=q[stk[i+1]];
        while(dcmp(area(d,e,q[stk[a]]) , area(d,e,q[stk[a+1]]) )<0)a=(a+1)%top;
        while(dcmp(project(d,e,q[stk[b]]),project(d,e,q[stk[b+1]]) )<0)b=(b+1)%top;
        if(!i)c=a;
        while(dcmp(project(d,e,q[stk[c]]),project(d,e,q[stk[c+1]]) )>0)c=(c+1)%top;
        auto x = q[stk[a]];
        auto y = q[stk[b]];
        auto z = q[stk[c]];
        auto h = area(d,e,x)/get_len(e-d);
        //auto w = ( (y-z)&(e-d) ) / get_len(e-d);//其实就是求bc在ed上的投影
        auto w = -project(d,e,z) + project(d,e,y);
        if(h * w <min_area){
            min_area = h*w;
            ans[0]=d+norm(e-d)*project(d,e,y);
            ans[3]=d+norm(e-d)*project(d,e,z);
            auto u = norm(rotate(e - d , -PI/2) );
            ans[1]=ans[0] + u*h;
            ans[2]=ans[3] + u*h;
        }
    }
}

signed main(){
    scanf("%d",&n);
    for(int i=0;i<n;i++)scanf("%lf %lf",&q[i].x,&q[i].y);
    get_convex();
    rotating_calipers();
    int k=0;
    for(int i=1;i<4;i++){
        if(dcmp(ans[i].y,ans[k].y)<0 || !dcmp(ans[i].y,ans[k].y)&&dcmp(ans[i].x,ans[k].x)<0)
        k=i;
    }
    printf("%.5lf\n",min_area);
    for(int i=0;i<4;i++,k++){
        auto x=ans[k%4].x;
        auto y=ans[k%4].y;
        if(!sign(x))x=0;
        if(!sign(y))y=0;
        printf("%.5lf %.5lf\n",x,y);
    }
}

三角剖分

3034. 望远镜

在这里插入图片描述

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cmath>
#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 55;
const double eps = 1e-8;
const double PI = acos(-1);

double R;
int n;
PDD q[N], r;

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    return x < 0 ? -1 : 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    return x < y ? -1 : 1;
}

PDD operator+ (PDD a, PDD b) { return {a.x + b.x, a.y + b.y}; }
PDD operator- (PDD a, PDD b) { return {a.x - b.x, a.y - b.y}; }
PDD operator* (PDD a, double b) { return {a.x * b, a.y * b}; }
PDD operator/ (PDD a, double b) { return {a.x / b, a.y / b}; }
double operator* (PDD a, PDD b) { return {a.x * b.y - a.y * b.x}; }
double operator& (PDD a, PDD b) { return {a.x * b.x + a.y * b.y}; }
double area(PDD a, PDD b, PDD c) { return (b - a) * (c - a); }
double get_len(PDD a) { return sqrt(a & a); }
double get_dist(PDD a, PDD b) { return get_len(b - a); }
double project(PDD a, PDD b, PDD c) { return ((c - a) & (b - a)) / get_len(b - a); }
PDD rotate(PDD a, double b) { return {a.x * cos(b) + a.y * sin(b), -a.x * sin(b) + a.y * cos(b)}; }
PDD norm(PDD a) { return a / get_len(a); }
bool on_segment(PDD p, PDD a, PDD b) { return !sign((p - a) * (p - b)) && sign((p - a) & (p - b)) <= 0; }

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w) // 求两直线交点
{
    auto u = p - q;
    double t = (w * u) / (v * w);
    return p + v * t;
}

double get_circle_line_intersection(PDD a, PDD b, PDD& pa, PDD& pb) // 求圆与线段ab的交点，并返回点到线段距离
{
    auto e = get_line_intersection(a, b - a, r, rotate(b - a, PI / 2));
    auto mind = get_dist(r, e);
    if (!on_segment(e, a, b)) mind = min(get_dist(r, a), get_dist(r, b));
    if (dcmp(R, mind) <= 0) return mind;
    auto len = sqrt(R * R - get_dist(r, e) * get_dist(r, e));
    pa = e + norm(a - b) * len;
    pb = e + norm(b - a) * len;
    return mind;
}

double get_sector(PDD a, PDD b) //求扇形面积
{
    auto angle = acos((a & b) / get_len(a) / get_len(b));
     if (sign(a * b) < 0) angle *= -1;
     return R * R * angle / 2;
}

double get_circle_triangle_area(PDD a, PDD b) // 求圆与三角形相交部分面积
{
    auto da = get_dist(r, a), db = get_dist(r, b);
    if (dcmp(R, da) >= 0 && dcmp(R, db) >= 0) return a * b / 2;
    if (!sign(a * b)) return 0;
    PDD pa, pb;  //  直线与圆的交点
    auto mind = get_circle_line_intersection(a, b, pa, pb);
    if (dcmp(R, mind) <= 0) return get_sector(a, b);
    if (dcmp(R, da) >= 0) return a * pb / 2 + get_sector(pb, b);
    if (dcmp(R, db) >= 0) return pa * b / 2 + get_sector(a, pa);
    return get_sector(pb, b) + pa * pb / 2 + get_sector(a, pa);
}

double work()
{
    double res = 0;
    for (int i = 0; i < n; i++) // 求出多边形所有相邻两点与原点构成的三角形与圆的交的部分的面积的和
        res += get_circle_triangle_area(q[i], q[(i + 1) % n]);
    return fabs(res); // 因为顶点给定顺序是不确定的，因此需要输出绝对值
}

int main()
{
    while (cin >> R >> n)
    {
        for (int i = 0; i < n; i++)  cin >> q[i].x >> q[i].y;
        printf("%.2lf\n", work());
    }
}

扫描线

3068. 扫描线

在这里插入图片描述

#include<bits/stdc++.h>
#define x first
#define y second
#define int long long
using namespace std;
const int N=1e3+10;
typedef pair<int,int>PII;
PII l[N],r[N];
vector<int>xs;
int ans;
int n;
int range_area(int a,int b){
    vector<PII>v;
    for(int i=0;i<n;i++){
        if(l[i].x<=a&&r[i].x>=b)
        v.push_back({l[i].y,r[i].y});
    }
    if(v.empty())return 0;
    sort(v.begin(),v.end());
    int res=0;
    int st=v[0].x;
    int ed=v[0].y;
    for(int i=1;i<v.size();i++){
        if(ed>=v[i].x)ed=max(ed,v[i].y);
        else{
            res+=ed-st;
            st=v[i].x;
            ed=v[i].y;
        }
    }
    res+=ed-st;
    return res*(b-a);
}

signed main(){
    cin>>n;
    for(int i=0;i<n;i++){
        cin>>l[i].x>>l[i].y>>r[i].x>>r[i].y;
        xs.push_back(l[i].x);
        xs.push_back(r[i].x);
    }
    sort(xs.begin(),xs.end());
    xs.erase((xs.begin(),xs.end()),xs.end());
    for(int i=0;i+1<xs.size();i++){
        ans+=abs(range_area(xs[i],xs[i+1]));
    }
    cout<<ans;
}

2801. 三角形面积并

在这里插入图片描述

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

#define x first
#define y second

using namespace std;

typedef pair<double, double> PDD;
const int N = 110;
const double eps = 1e-8, INF = 1e6;

int n;
PDD tr[N][3];
PDD q[N];

int sign(double x)
{
    if (fabs(x) < eps) return 0;
    if (x < 0) return -1;
    return 1;
}

int dcmp(double x, double y)
{
    if (fabs(x - y) < eps) return 0;
    if (x < y) return -1;
    return 1;
}

PDD operator+ (PDD a, PDD b)
{
    return {a.x + b.x, a.y + b.y};
}

PDD operator- (PDD a, PDD b)
{
    return {a.x - b.x, a.y - b.y};
}

PDD operator* (PDD a, double t)
{
    return {a.x * t, a.y * t};
}

double operator* (PDD a, PDD b)
{
    return a.x * b.y - a.y * b.x;
}

double operator& (PDD a, PDD b)
{
    return a.x * b.x + a.y * b.y;
}

bool on_segment(PDD p, PDD a, PDD b)
{
    return sign((p - a) & (p - b)) <= 0;
}

PDD get_line_intersection(PDD p, PDD v, PDD q, PDD w)
{
    if (!sign(v * w)) return {INF, INF};
    auto u = p - q;
    auto t = w * u / (v * w);
    auto o = p + v * t;
    if (!on_segment(o, p, p + v) || !on_segment(o, q, q + w))
        return {INF, INF};
    return o;
}

double line_area(double a, int side)
{
    int cnt = 0;
    for (int i = 0; i < n; i ++ )
    {
        auto t = tr[i];
        if (dcmp(t[0].x, a) > 0 || dcmp(t[2].x, a) < 0) continue;
        if (!dcmp(t[0].x, a) && !dcmp(t[1].x, a))
        {
            if (side) q[cnt ++ ] = {t[0].y, t[1].y};
        }
        else if (!dcmp(t[2].x, a) && !dcmp(t[1].x, a))
        {
            if (!side) q[cnt ++ ] = {t[2].y, t[1].y};
        }
        else
        {
            double d[3];
            int u = 0;
            for (int j = 0; j < 3; j ++ )
            {
                auto o = get_line_intersection(t[j], t[(j + 1) % 3] - t[j], {a, -INF}, {0, INF * 2});
                if (dcmp(o.x, INF))
                    d[u ++ ] = o.y;
            }
            if (u)
            {
                sort(d, d + u);
                q[cnt ++ ] = {d[0], d[u - 1]};
            }
        }
    }
    if (!cnt) return 0;
    for (int i = 0; i < cnt; i ++ )
        if (q[i].x > q[i].y)
            swap(q[i].x, q[i].y);
    sort(q, q + cnt);
    double res = 0, st = q[0].x, ed = q[0].y;
    for (int i = 1; i < cnt; i ++ )
        if (q[i].x <= ed) ed = max(ed, q[i].y);
        else
        {
            res += ed - st;
            st = q[i].x, ed = q[i].y;
        }
    res += ed - st;
    return res;
}

double range_area(double a, double b)
{
    return (line_area(a, 1) + line_area(b, 0)) * (b - a) / 2;
}

int main()
{
    scanf("%d", &n);
    vector<double> xs;
    for (int i = 0; i < n; i ++ )
    {
        for (int j = 0; j < 3; j ++ )
        {
            scanf("%lf%lf", &tr[i][j].x, &tr[i][j].y);
            xs.push_back(tr[i][j].x);
        }
        sort(tr[i], tr[i] + 3);
    }
    for (int i = 0; i < n; i ++ )
        for (int j = i + 1; j < n; j ++ )
            for (int x = 0; x < 3; x ++ )
                for (int y = 0; y < 3; y ++ )
                {
                    auto o = get_line_intersection(tr[i][x], tr[i][(x + 1) % 3] - tr[i][x],
                                                    tr[j][y], tr[j][(y + 1) % 3] - tr[j][y]);
                    if (dcmp(o.x, INF))
                        xs.push_back(o.x);
                }
    sort(xs.begin(), xs.end());
    double res = 0;
    for (int i = 0; i + 1 < xs.size(); i ++ )
        if (dcmp(xs[i], xs[i + 1]))
            res += range_area(xs[i], xs[i + 1]);
    printf("%.2lf\n", res);
    return 0;
}

自适应辛普森积分

3074. 自适应辛普森积分

在这里插入图片描述

取小细条的左右端点和中间点：[a,f(a)],[(a+b)/2,f((a+b)/2)],[b,f(b)][a,f(a)],[(a+b)/2,f((a+b)/2)],[b,f(b)]，以这三个点做一条抛物线，这条抛物线与xx轴围成的面积SS就是小细条的面积.但是这样做的误差可能会较大，于是有了自适应辛普森积分法，做法如下：
然后再分别求出以[a,f(a)],[(a+b)/2,f((a+b)/2)][a,f(a)],[(a+b)/2,f((a+b)/2)]和[(a+b)/2,f((a+b)/2)],[b,f(b)][(a+b)/2,f((a+b)/2)],[b,f(b)]为两端点的小细条的面积S1S1和S2S2.
判断S−(S1+S2)S−(S1+S2)是否足够小，如果足够小就取SS，否则递归求两个更小细条的面积。

#include<bits/stdc++.h>
using namespace std;
const double eps=1e-4;
//eps开大了不精确可能会wa
//eps开小了递归次数增加可能会超时

double f(double x){
    return sin(x)/x;
}
double simpson(double l,double r){//辛普森积分公式，用于求l，r，mid为顶点的抛物线面积
    double mid=(l+r)/2;
    return (r-l)*(f(l)+4*f(mid)+f(r))/6;
}
double asr(double l,double r,double s){
    double mid=(l+r)/2;
    double left=simpson(l,mid);
    double right=simpson(mid,r);
    //if(fabs(left+right-s)<eps)return s;
    if(fabs(l-mid)<eps)return s;//原来开成两端足够小也是对的，这时eps取小一点就好了
    else return asr(l,mid,left)+asr(mid,r,right);
}
signed main(){
    double a,b;
    cin>>a>>b;
    printf("%.6f",asr(a,b,simpson(a,b)) );
}

圆的面积并

在这里插入图片描述

胡耀文

关注

0
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
计算几何进阶

基础知识2983. 玩具#include<bits/stdc++.h>using namespace std;#define int long longtypedef pair<int, int> PII;#define x first#define y second//#define d lldint n;PII a[5050],b[5050];//a在下面，b在上面int ans[5050];int cross(int x1,int y1,int x2,
复制链接

扫一扫