被卡精度了,大佬们帮调调。
查看原帖
被卡精度了,大佬们帮调调。
460020
mmr123楼主2024/9/27 16:23
#include <bits/stdc++.h>
using namespace std;

// #define double long double

const double eps = numeric_limits<double>::epsilon();

double _fabs(double x) {
    if (x > -eps && x < eps) {
        return 0;
    } else {
        return (x < 0) ? (-x) : x;
    }
}

bool epse(double a, double b) {
    return (_fabs(a - b) < eps) ? true : false;
}

struct Segment {
    double x, y;
};

struct Coordinate {
    double x, y;

    double operator-(Coordinate const& a) const {
        return sqrt(pow(this->x - a.x, 2) + pow(this->y - a.y, 2));
    }
};

struct Circle {
    Coordinate pos;
    double     r;

    int relation(Circle x) {
        double d = this->pos - x.pos;

        if (epse(d, this->r + x.r)) {
            return 0;
        } else if (d > this->r + x.r) {
            return 1;
        } else {
            return -1;
        }
    }

    bool in(Coordinate p) {
        double d = this->pos - p;

        if (d < this->r) {
            return true;
        } else {
            return false;
        }
    }
};

int    n;
int    cnt;
double ans;

vector<Circle> c;
vector<Circle> c_tmp;

bool cover[1005];

int             st, ed;
vector<Segment> seg;

double f(double x) {
    seg.clear();

    for (int i = st; i <= ed; i++) {
        if (x > c[i].pos.x + c[i].r || x < c[i].pos.x - c[i].r) {
            continue;
        }
        double l = x - c[i].pos.x;
        l        = sqrt(c[i].r * c[i].r - l * l);
        seg.push_back({c[i].pos.y - l, c[i].pos.y + l});
    }

    if (seg.empty()) {
        return 0;
    }

    sort(seg.begin(), seg.end(), [](Segment a, Segment b) {
        if (a.x == b.x) return a.y < b.y;
        return a.x < b.y;
    });

    double ans = 0, last = seg[0].x;

    for (int i = 0; i < seg.size(); ++i) {
        if (seg[i].y > last) {
            ans += seg[i].y - max(last, seg[i].x);
            last = seg[i].y;
        }
    }

    return ans;
}

double calc(double l, double r, double sl, double sr) {
    double mid = (l + r) / 2.0;
    return (r - l) * (sl + 4.0 * f(mid) + sr) / 6.0;
}

double simpson(double l, double r, double S, double sl, double sr) {
    double mid = (l + r) / 2.0, s = f(mid);

    double t1 = calc(l, mid, sl, s), t2 = calc(mid, r, s, sr);

    if (_fabs(t1 + t2 - S) < eps) { return S; }
    ++cnt;
    return simpson(l, mid, t1, sl, s) + simpson(mid, r, t2, s, sr);
}

int main() {

    // freopen("SP8073.in", "r", stdin);
    // freopen("SP8073.out", "w", stdout);

    cin >> n;

    for (int i = 1; i <= n; i++) {
        double x, y, z;
        cin >> x >> y >> z;

        c_tmp.push_back({{x, y}, z});
    }

    sort(c_tmp.begin(), c_tmp.end(), [](Circle a, Circle b) {
        return a.r > b.r;
    });

    for (auto i = 0; i < c_tmp.size(); i++) {
        if (!cover[i]) {
            c.push_back(c_tmp[i]);
        }
        for (int j = i + 1; j < c_tmp.size(); j++) {
            if (c_tmp[i].relation(c_tmp[j]) == -1) {
                cover[j] = true;
            }
        }
    }

    sort(c.begin(), c.end(), [](Circle a, Circle b) {
        return a.pos.x - a.r < b.pos.x - b.r;
    });

    double l, r;
    int    i = 0, j;

    while (i < c.size()) {
        l = c[i].pos.x - c[i].r;
        r = c[i].pos.x + c[i].r;
        for (j = i + 1; j < c.size(); j++) {
            if (c[j].pos.x - c[j].r > r) {
                break;
            }
            r = max(r, c[j].pos.x + c[j].r);
        }
        st = i, ed = j - 1, i = j;
        double mid = (l + r) / 2.0;
        ans += simpson(l, r, f(mid), f(l), f(r));
    }

    printf("%lf\n", ans);

    // fclose(stdin);
    // fclose(stdout);

    // system("pause");

    return 0;
}

自适应辛普森法,SPOJ 上在 14 - 16 测试点之间 WA。

2024/9/27 16:23
加载中...