QOJ.ac

QOJ

IDProblemSubmitterResultTimeMemoryLanguageFile sizeSubmit timeJudge time
#552362#9227. Henry the Plumberucup-team4435WA 0ms3696kbC++2036.0kb2024-09-07 22:19:212024-09-07 22:19:21

Judging History

你现在查看的是最新测评结果

  • [2024-09-07 22:19:21]
  • 评测
  • 测评结果:WA
  • 用时:0ms
  • 内存:3696kb
  • [2024-09-07 22:19:21]
  • 提交

answer

#pragma GCC optimize("Ofast")

#include "bits/stdc++.h"

#define rep(i, n) for (int i = 0; i < (n); ++i)
#define rep1(i, n) for (int i = 1; i < (n); ++i)
#define rep1n(i, n) for (int i = 1; i <= (n); ++i)
#define repr(i, n) for (int i = (n) - 1; i >= 0; --i)
#define pb push_back
#define eb emplace_back
#define all(a) (a).begin(), (a).end()
#define rall(a) (a).rbegin(), (a).rend()
#define each(x, a) for (auto &x : a)
#define ar array
#define vec vector
#define range(i, n) rep(i, n)

using namespace std;


using ll = long long;
using ull = unsigned long long;
using ld = long double;
using str = string;
using pi = pair<int, int>;
using pl = pair<ll, ll>;

using vi = vector<int>;
using vl = vector<ll>;
using vpi = vector<pair<int, int>>;
using vvi = vector<vi>;

int Bit(int mask, int b) { return (mask >> b) & 1; }

template<class T>
bool ckmin(T &a, const T &b) {
    if (b < a) {
        a = b;
        return true;
    }
    return false;
}

template<class T>
bool ckmax(T &a, const T &b) {
    if (b > a) {
        a = b;
        return true;
    }
    return false;
}

// [l, r)
template<typename T, typename F>
T FindFirstTrue(T l, T r, const F &predicat) {
    --l;
    while (r - l > 1) {
        T mid = l + (r - l) / 2;
        if (predicat(mid)) {
            r = mid;
        } else {
            l = mid;
        }
    }
    return r;
}


template<typename T, typename F>
T FindLastFalse(T l, T r, const F &predicat) {
    return FindFirstTrue(l, r, predicat) - 1;
}

const int INFi = 2e9;
const ll INF = 2e18;
const int LG = 31;
const int N = 1e4;

const double eps = 1e-2;

namespace geom {
#define double ld
    using namespace std;

// https://victorlecomte.com/cp-geo.pdf
    const double inf = 1e100;
    const double PI = acos((double) -1.0);

    int sign(double x) { return (x > eps) - (x < -eps); }

    struct PT {
        double x, y;

        PT() { x = 0, y = 0; }

        PT(double x, double y) : x(x), y(y) {}

        PT(const PT &p) : x(p.x), y(p.y) {}

        void scan() { cin >> x >> y; }

        PT operator+(const PT &a) const { return PT(x + a.x, y + a.y); }

        PT operator-(const PT &a) const { return PT(x - a.x, y - a.y); }

        PT operator*(const double a) const { return PT(x * a, y * a); }

        friend PT operator*(const double &a, const PT &b) { return PT(a * b.x, a * b.y); }

        PT operator/(const double a) const { return PT(x / a, y / a); }

        bool operator==(PT a) const { return sign(a.x - x) == 0 && sign(a.y - y) == 0; }

        bool operator!=(PT a) const { return !(*this == a); }

        bool operator<(PT a) const { return sign(a.x - x) == 0 ? y < a.y : x < a.x; }

        bool operator>(PT a) const { return sign(a.x - x) == 0 ? y > a.y : x > a.x; }

        double norm() { return sqrt(x * x + y * y); }

        double norm2() { return x * x + y * y; }

        PT perp() { return PT(-y, x); }

        double arg() { return atan2(y, x); }

        PT truncate(double r) { // returns a vector with norm r and having same direction
            double k = norm();
            if (!sign(k)) return *this;
            r /= k;
            return PT(x * r, y * r);
        }
    };

    inline double dot(PT a, PT b) { return a.x * b.x + a.y * b.y; }

    inline double dist2(PT a, PT b) { return dot(a - b, a - b); }

    inline double dist(PT a, PT b) { return sqrt(dot(a - b, a - b)); }

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

    inline int orientation(PT a, PT b, PT c) { return sign(cross(b - a, c - a)); }

    PT perp(PT a) { return PT(-a.y, a.x); }

    PT rotateccw90(PT a) { return PT(-a.y, a.x); }

    PT rotatecw90(PT a) { return PT(a.y, -a.x); }

    PT rotateccw(PT a, double t) { return PT(a.x * cos(t) - a.y * sin(t), a.x * sin(t) + a.y * cos(t)); }

    PT rotatecw(PT a, double t) { return PT(a.x * cos(t) + a.y * sin(t), -a.x * sin(t) + a.y * cos(t)); }

    double SQ(double x) { return x * x; }

    double rad_to_deg(double r) { return (r * 180.0 / PI); }

    double deg_to_rad(double d) { return (d * PI / 180.0); }

    double get_angle(PT a, PT b) {
        double costheta = dot(a, b) / a.norm() / b.norm();
        return acos(max((double) -1.0, min((double) 1.0, costheta)));
    }

    struct p3 {
        double x, y, z;

        p3() {
            x = 0, y = 0;
            z = 0;
        }

        p3(double x, double y, double z) : x(x), y(y), z(z) {}

        p3(const p3 &p) : x(p.x), y(p.y), z(p.z) {}

        void scan() { cin >> x >> y >> z; }

        p3 operator+(const p3 &a) const { return p3(x + a.x, y + a.y, z + a.z); }

        p3 operator-(const p3 &a) const { return p3(x - a.x, y - a.y, z - a.z); }

        p3 operator*(const double a) const { return p3(x * a, y * a, z * a); }

        friend p3 operator*(const double &a, const p3 &b) { return p3(a * b.x, a * b.y, a * b.z); }

        p3 operator/(const double a) const { return p3(x / a, y / a, z / a); }

        bool operator==(p3 a) const { return sign(a.x - x) == 0 && sign(a.y - y) == 0 && sign(a.z - z) == 0; }

        bool operator!=(p3 a) const { return !(*this == a); }

        double abs() { return sqrt(x * x + y * y + z * z); }

        double sq() { return x * x + y * y + z * z; }

        p3 unit() { return *this / abs(); }
    } zero(0, 0, 0);

    double operator|(p3 v, p3 w) { //dot product
        return v.x * w.x + v.y * w.y + v.z * w.z;
    }

    p3 operator*(p3 v, p3 w) { //cross product
        return {v.y * w.z - v.z * w.y, v.z * w.x - v.x * w.z, v.x * w.y - v.y * w.x};
    }

    double sq(p3 v) { return v | v; }

    double abs(p3 v) { return sqrt(sq(v)); }

    p3 unit(p3 v) { return v / abs(v); }

    inline double dot(p3 a, p3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }

    inline double dist2(p3 a, p3 b) { return dot(a - b, a - b); }

    inline double dist(p3 a, p3 b) { return sqrt(dot(a - b, a - b)); }

    inline p3 cross(p3 a, p3 b) { return p3(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x); }

// if s is on the same side of the plane pqr as the vector pq * pr then it will be positive
// otherwise negative or 0 if on the plane
    double orient(p3 p, p3 q, p3 r, p3 s) { return (q - p) * (r - p) | (s - p); }

// returns orientation of p to q to r on the plane perpendicular to n
// assuming p, q, r are on the plane
    double orient_by_normal(p3 p, p3 q, p3 r, p3 n) { return (q - p) * (r - p) | n; }

    double get_angle(p3 a, p3 b) {
        double costheta = dot(a, b) / a.abs() / b.abs();
        return acos(max((double) -1.0, min((double) 1.0, costheta)));
    }

    double small_angle(p3 v, p3 w) {
        return acos(min(fabs(v | w) / abs(v) / abs(w), (double) 1.0));
    }

    struct plane {
        // n is the perpendicular normal vector to the plane
        p3 n;
        double d; // (n | p) = d
        // From normal n and offset d
        plane(p3 n, double d) : n(n), d(d) {}

        // From normal n and point P
        plane(p3 n, p3 p) : n(n), d(n | p) {}

        // From three non-collinear points P,Q,R
        plane(p3 p, p3 q, p3 r) : plane((q - p) * (r - p), p) {}

        // positive if on the same side as the normal, negative if on the opposite side, 0 if on the plane
        double side(p3 p) { return (n | p) - d; }

        // distance from point p to plane
        double dist(p3 p) { return fabs(side(p)) / abs(n); }

        // translate the plane by vector t
        plane translate(p3 t) { return {n, d + (n | t)}; }

        // shift the plane perpendicular to n by distance dist
        plane shiftUp(double dist) { return {n, d + dist * abs(n)}; }

        // orthogonal projection of point p onto plane
        p3 proj(p3 p) { return p - n * side(p) / sq(n); }

        // orthogonal reflection of point p onto plane
        p3 refl(p3 p) { return p - n * 2 * side(p) / sq(n); }

        pair<p3, p3> get_two_points_on_plane() {
            assert(sign(n.x) != 0 || sign(n.y) != 0 || sign(n.z) != 0);
            if (sign(n.x) == 0 && sign(n.y) == 0) return {p3(1, 0, d / n.z), p3(0, 1, d / n.z)};
            if (sign(n.y) == 0 && sign(n.z) == 0) return {p3(d / n.x, 1, 0), p3(d / n.x, 0, 1)};
            if (sign(n.z) == 0 && sign(n.x) == 0) return {p3(1, d / n.y, 0), p3(0, d / n.y, 1)};
            if (sign(n.x) == 0) return {p3(1, d / n.y, 0), p3(0, 0, d / n.z)};
            if (sign(n.y) == 0) return {p3(0, 1, d / n.z), p3(d / n.x, 0, 0)};
            if (sign(n.z) == 0) return {p3(d / n.x, 0, 1), p3(0, d / n.y, 0)};
            if (sign(d) != 0) return {p3(d / n.x, 0, 0), p3(0, d / n.y, 0)};
            return {p3(n.y, -n.x, 0), p3(-n.y, n.x, 0)};
        }
    };

    struct coords {
        // coordinate system for coplanar points
        // o is the origin, dx, dy, dz are unit vectors similar to normal 3D system
        // but dx and dy are on the plane
        p3 o, dx, dy, dz;

        // From three points P, Q, R on the plane
        coords(p3 p, p3 q, p3 r) : o(p) {
            dx = unit(q - p);
            dz = unit(dx * (r - p));
            dy = dz * dx;
        }

        // From four points P,Q,R,S: take directions PQ, PR, PS as is
        // it allows us to keep using integer coordinates but has some pitfalls
        // e.g. distances and angles are not preserved but relative positions are (convex hull works)
        coords(p3 p, p3 q, p3 r, p3 s) :
                o(p), dx(q - p), dy(r - p), dz(s - p) {}

        // 2D position vector of point p in this coordinate system centered at o
        // p must be on the plane
        PT pos2d(p3 p) {
            return {(p - o) | dx, (p - o) | dy};
        }

        // returns the 3D position vector of point p in this new coordinate system
        // p can be outside the plane
        p3 pos3d(p3 p) {
            return {(p - o) | dx, (p - o) | dy, (p - o) | dz};
        }

        // given 2D position vector p centered at o, return the original 3D position vector
        p3 pos3d(PT p) {
            return o + dx * p.x + dy * p.y;
        }
    };

    struct line3d {
        // d is the direction vector of the line
        p3 d, o; // p = o + k * d (k is a real parameter)
        line3d() {}

        // From two points P, Q
        line3d(p3 p, p3 q) : d(q - p), o(p) {}

        // From two planes p1, p2
        // assuming they are not parallel
        line3d(plane p1, plane p2) {
            d = p1.n * p2.n;
            o = (p2.n * p1.d - p1.n * p2.d) * d / sq(d);
            // o is actually the closest point on the line to the origin
        }

        double dist2(p3 p) { return sq(d * (p - o)) / sq(d); }

        double dist(p3 p) { return sqrt(dist2(p)); }

        // compare points by their projection on the line
        // so you can sort points on the line using this
        bool cmp_proj(p3 p, p3 q) { return (d | p) < (d | q); }

        // orthogonal projection of point p onto line
        p3 proj(p3 p) { return o + d * (d | (p - o)) / sq(d); }

        // orthogonal reflection of point p onto line
        p3 refl(p3 p) { return proj(p) * 2 - p; }

        // returns the intersection point of the line with plane p
        // assuming plane and line are not parallel
        p3 inter(plane p) {
            // assert((d | p.n) != 0); // no intersection if parallel
            return o - d * p.side(o) / (d | p.n);
        }
    };

// smallest distance between two lines
    double dist(line3d l1, line3d l2) {
        p3 n = l1.d * l2.d;
        if (n == zero) return l1.dist(l2.o); // parallel
        return fabs((l2.o - l1.o) | n) / abs(n);
    }

// closest point from line l2 to line l1
    p3 closest_on_l1(line3d l1, line3d l2) {
        p3 n2 = l2.d * (l1.d * l2.d);
        return l1.o + l1.d * ((l2.o - l1.o) | n2) / (l1.d | n2);
    }

// small angle between direction vectors of two lines
    double get_angle(line3d l1, line3d l2) {
        return small_angle(l1.d, l2.d);
    }

    bool is_parallel(line3d l1, line3d l2) {
        return l1.d * l2.d == zero;
    }

    bool is_perpendicular(line3d l1, line3d l2) {
        return sign((l1.d | l2.d)) == 0;
    }

// small angle between normal vectors of two planes
    double get_angle(plane p1, plane p2) {
        return small_angle(p1.n, p2.n);
    }

    bool is_parallel(plane p1, plane p2) {
        return p1.n * p2.n == zero;
    }

    bool is_perpendicular(plane p1, plane p2) {
        return sign((p1.n | p2.n)) == 0;
    }

    double get_angle(plane p, line3d l) {
        return PI / 2 - small_angle(p.n, l.d);
    }

    bool is_parallel(plane p, line3d l) {
        return sign((p.n | l.d)) == 0;
    }

    bool is_perpendicular(plane p, line3d l) {
        return p.n * l.d == zero;
    }

// returns the line perpendicular to plane p and passing through point o
    line3d perp_through(plane p, p3 o) { return line3d(o, o + p.n); }

// returns the plane perpendicular to line l and passing through point o
    plane perp_through(line3d l, p3 o) { return plane(l.d, o); }

// returns two points on intesection line of two planes formed by points
// a1, b1, c1 and a2, b2, c2 respectively
    pair<p3, p3> plane_plane_intersection(p3 a1, p3 b1, p3 c1, p3 a2, p3 b2, p3 c2) {
        p3 n1 = (b1 - a1) * (c1 - a1);
        p3 n2 = (b2 - a2) * (c2 - a2);
        double d1 = n1 | a1, d2 = n2 | a2;
        p3 d = n1 * n2;
        if (d == zero) return make_pair(zero, zero);
        p3 o = (n2 * d1 - n1 * d2) * d / (d | d);
        return make_pair(o, o + d);
    }

// returns center of circle passing through three
// non-colinear and co-planer points a, b and c
    p3 circle_center(p3 a, p3 b, p3 c) {
        p3 v1 = b - a, v2 = c - a;
        double v1v1 = v1 | v1, v2v2 = v2 | v2, v1v2 = v1 | v2;
        double base = 0.5 / (v1v1 * v2v2 - v1v2 * v1v2);
        double k1 = base * v2v2 * (v1v1 - v1v2);
        double k2 = base * v1v1 * (v2v2 - v1v2);
        return a + v1 * k1 + v2 * k2;
    }

// segment ab to point c
    double distance_from_segment_to_point(p3 a, p3 b, p3 c) {
        if (sign(dot(b - a, c - a)) < 0) return dist(a, c);
        if (sign(dot(a - b, c - b)) < 0) return dist(b, c);
        return fabs(cross((b - a).unit(), c - a).abs());
    }

    double distance_from_triangle_to_point(p3 a, p3 b, p3 c, p3 d) {
        plane P(a, b, c);
        p3 proj = P.proj(d);
        double dis = min(distance_from_segment_to_point(a, b, d),
                         min(distance_from_segment_to_point(b, c, d), distance_from_segment_to_point(c, a, d)));
        int o = sign(orient_by_normal(a, b, proj, P.n));
        int inside = o == sign(orient_by_normal(b, c, proj, P.n));
        inside &= o == sign(orient_by_normal(c, a, proj, P.n));
        if (inside) return (d - proj).abs();
        return dis;
    }

    double distance_from_triangle_to_segment(p3 a, p3 b, p3 c, p3 d, p3 e) {
        double l = 0.0, r = 1.0;
        int cnt = 100;
        double ret = inf;
        while (cnt--) {
            double mid1 = l + (r - l) / 3.0, mid2 = r - (r - l) / 3.0;
            double x = distance_from_triangle_to_point(a, b, c, d + (e - d) * mid1);
            double y = distance_from_triangle_to_point(a, b, c, d + (e - d) * mid2);
            if (x < y) {
                r = mid2;
                ret = x;
            } else {
                ret = y;
                l = mid1;
            }
        }
        return ret;
    }

// triangles are solid
    double distance_from_triangle_to_triangle(p3 a, p3 b, p3 c, p3 d, p3 e, p3 f) {
        double ret = inf;
        ret = min(ret, distance_from_triangle_to_segment(a, b, c, d, e));
        ret = min(ret, distance_from_triangle_to_segment(a, b, c, e, f));
        ret = min(ret, distance_from_triangle_to_segment(a, b, c, f, d));
        ret = min(ret, distance_from_triangle_to_segment(d, e, f, a, b));
        ret = min(ret, distance_from_triangle_to_segment(d, e, f, b, c));
        ret = min(ret, distance_from_triangle_to_segment(d, e, f, c, a));
        return ret;
    }


    bool operator<(p3 p, p3 q) {
        return tie(p.x, p.y, p.z) < tie(q.x, q.y, q.z);
    }

    struct edge {
        int v;
        bool same; // is the common edge between two faces in the same order?
    };

// Given a series of faces (lists of points) of a polyhedron, reverse some of them
// so that their orientations are consistent (all area vectors of the faces either pointing outwards or inwards)
// just compute the area vector of one face to see if it’s pointing outwards or inwards
    vector<vector<p3>> reorient(vector<vector<p3>> fs) {
        int n = fs.size();
        // Find the common edges and create the resulting graph
        vector<vector<edge>> g(n);
        map<pair<p3, p3>, int> es;
        for (int u = 0; u < n; u++) {
            for (int i = 0, m = fs[u].size(); i < m; i++) {
                p3 a = fs[u][i], b = fs[u][(i + 1) % m];
                // Let’s look at edge a-b
                if (es.count({a, b})) { // seen in same order
                    int v = es[{a, b}];
                    g[u].push_back({v, true});
                    g[v].push_back({u, true});
                } else if (es.count({b, a})) { // seen in different order
                    int v = es[{b, a}];
                    g[u].push_back({v, false});
                    g[v].push_back({u, false});
                } else { // not seen yet
                    es[{a, b}] = u;
                }
            }
        }
        // Perform BFS to find which faces should be flipped
        vector<bool> vis(n, false), flip(n);
        flip[0] = false;
        queue<int> q;
        q.push(0);
        while (!q.empty()) {
            int u = q.front();
            q.pop();
            for (edge e: g[u]) {
                if (!vis[e.v]) {
                    vis[e.v] = true;
                    // If the edge was in the same order,
                    // exactly one of the two should be flipped
                    flip[e.v] = (flip[u] ^ e.same);
                    q.push(e.v);
                }
            }
        }
        // Actually perform the flips
        for (int u = 0; u < n; u++)
            if (flip[u]) {
                reverse(fs[u].begin(), fs[u].end());
            }
        return fs;
    }

// O(n^2), O(n) faces in the hull
    struct CH3D {
        struct face {
            int a, b, c; // the number of three points on one face of the convex hull
            bool ok;  // whether the face belongs to the face on the final convex hull
        };
        int n; // initial vertex number
        vector<p3> P;
        int num; // convex hull surface triangle number
        vector<face> F; // convex surface triangles
        vector<vector<int>> g;

        void init(vector<p3> p) {
            P = p;
            n = p.size();
            F.resize(8 * n + 1);
            g.resize(n + 1, vector<int>(n + 1));
        }

        double len(p3 a) {
            return sqrt(a | a);
        }

        p3 cross(const p3 &a, const p3 &b, const p3 &c) {
            return (b - a) * (c - a);
        }

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

        double volume(p3 a, p3 b, p3 c, p3 d) {
            return (b - a) * (c - a) | (d - a);
        }

        // positive: p3 in the same direction
        double dblcmp(p3 &p, face &f) {
            p3 m = P[f.b] - P[f.a];
            p3 n = P[f.c] - P[f.a];
            p3 t = p - P[f.a];
            return (m * n) | t;
        }

        void deal(int p, int a, int b) {
            int f = g[a][b]; // search for another plane adjacent to the edge
            face add;
            if (F[f].ok) {
                if (dblcmp(P[p], F[f]) > eps) dfs(p, f);
                else {
                    add.a = b;
                    add.b = a;
                    add.c = p; // pay attention to the order here, to be right-handed
                    add.ok = true;
                    g[p][b] = g[a][p] = g[b][a] = num;
                    F[num++] = add;
                }
            }
        }

        // recursively search all faces that should be removed from the convex hull
        void dfs(int p, int now) {
            F[now].ok = 0;
            deal(p, F[now].b, F[now].a);
            deal(p, F[now].c, F[now].b);
            deal(p, F[now].a, F[now].c);
        }

        bool same(int s, int t) {
            p3 &a = P[F[s].a];
            p3 &b = P[F[s].b];
            p3 &c = P[F[s].c];
            return fabs(volume(a, b, c, P[F[t].a])) < eps &&
                   fabs(volume(a, b, c, P[F[t].b])) < eps &&
                   fabs(volume(a, b, c, P[F[t].c])) < eps;
        }

        // building a 3D convex hull
        void create_hull() {
            int i, j, tmp;
            face add;
            num = 0;
            if (n < 4)return;

            // ensure that the first four points are not coplanar
            bool flag = true;
            for (i = 1; i < n; i++) {
                if (len(P[0] - P[i]) > eps) {
                    swap(P[1], P[i]);
                    flag = false;
                    break;
                }
            }
            if (flag) return;
            flag = true;
            // make the first three points not collinear
            for (i = 2; i < n; i++) {
                if (len((P[0] - P[1]) * (P[1] - P[i])) > eps) {
                    swap(P[2], P[i]);
                    flag = false;
                    break;
                }
            }
            if (flag) return;
            flag = true;
            // make the first four points not coplanar
            for (int i = 3; i < n; i++) {
                if (fabs((P[0] - P[1]) * (P[1] - P[2]) | (P[0] - P[i])) > eps) {
                    swap(P[3], P[i]);
                    flag = false;
                    break;
                }
            }
            if (flag) return;

            for (i = 0; i < 4; i++) {
                add.a = (i + 1) % 4;
                add.b = (i + 2) % 4;
                add.c = (i + 3) % 4;
                add.ok = true;
                if (dblcmp(P[i], add) > 0)swap(add.b, add.c);
                g[add.a][add.b] = g[add.b][add.c] = g[add.c][add.a] = num;
                F[num++] = add;
            }
            for (i = 4; i < n; i++) {
                for (j = 0; j < num; j++) {
                    if (F[j].ok && dblcmp(P[i], F[j]) > eps) {
                        dfs(i, j);
                        break;
                    }
                }
            }
            tmp = num;
            for (i = num = 0; i < tmp; i++)
                if (F[i].ok) F[num++] = F[i];

        }

        double surface_area() {
            double res = 0;
            if (n == 3) {
                p3 p = cross(P[0], P[1], P[2]);
                res = len(p) / 2.0;
                return res;
            }
            for (int i = 0; i < num; i++) {
                res += area(P[F[i].a], P[F[i].b], P[F[i].c]);
            }
            return res / 2.0;
        }

        double volume() {
            double res = 0;
            p3 tmp(0, 0, 0);
            for (int i = 0; i < num; i++) {
                res += volume(tmp, P[F[i].a], P[F[i].b], P[F[i].c]);
            }
            return fabs(res / 6.0);
        }

        int number_of_triangles() { // number of surface triangles
            return num;
        }

        int number_of_polygons() {  // number of surface polygons
            int i, j, res, flag;
            for (i = res = 0; i < num; i++) {
                flag = 1;
                for (j = 0; j < i; j++) {
                    if (same(i, j)) {
                        flag = 0;
                        break;
                    }
                }
                res += flag;
            }
            return res;
        }

        p3 centroid() { // center of gravity
            p3 ans(0, 0, 0), o(0, 0, 0);
            double all = 0;
            for (int i = 0; i < num; i++) {
                double vol = volume(o, P[F[i].a], P[F[i].b], P[F[i].c]);
                ans = ans + (o + P[F[i].a] + P[F[i].b] + P[F[i].c]) / 4.0 * vol;
                all += vol;
            }
            ans = ans / all;
            return ans;
        }

        double point_to_face_distance(p3 p, int i) {
            return fabs(volume(P[F[i].a], P[F[i].b], P[F[i].c], p) /
                        len((P[F[i].b] - P[F[i].a]) * (P[F[i].c] - P[F[i].a])));
        }
    };

// https://desktop.arcgis.com/en/arcmap/10.7/map/projections/geographic-coordinate-system.htm
// given the radius of the sphere, latitude and longitude of a point in degrees
// return the 3D coordinates of the point on the sphere assuming the sphere is centered at the origin
    p3 get_sphere(double r, double lat, double lon) {
        lat *= PI / 180, lon *= PI / 180;
        return {r * cos(lat) * cos(lon), r * cos(lat) * sin(lon), r * sin(lat)};
    }

    int sphere_line_intersection(p3 o, double r, line3d l, pair<p3, p3> &out) {
        double h2 = r * r - l.dist2(o);
        if (h2 < 1e-8) return 0; // the line doesn’t touch the sphere
        p3 p = l.proj(o);
        p3 h = l.d * sqrt(h2) / abs(l.d); // vector parallel to l, of length h
        out = {p - h, p + h};
        return 1 + (h2 > 0);
    }

// The shortest distance between two points A and B on a sphere (O, r) is
// given by travelling along plane OAB and on the surface of the sphere. It is called the great-circle distance
// if a and b are outside the sphere, then it will give the distance between their projections on the sphere
    double great_circle_dist(p3 o, double r, p3 a, p3 b) {
        // s = r * theta
        return r * get_angle(a - o, b - o);
    }

// Assume that the sphere is centered at the origin
// We will call a segment [AB] valid if A and B are not
// opposite each other on the sphere
    bool validSegment(p3 a, p3 b) {
        return a * b != zero || (a | b) > 0;
    }

    bool proper_intersection(p3 a, p3 b, p3 c, p3 d, p3 &out) {
        p3 ab = a * b, cd = c * d; // normals of planes OAB and OCD
        int oa = sign(cd | a),
                ob = sign(cd | b),
                oc = sign(ab | c),
                od = sign(ab | d);
        out = ab * cd * od; // four multiplications => careful with overflow!
        return (oa != ob && oc != od && oa != oc);
    }

// Assume that the sphere is centered at the origin
    bool point_on_sphere_segment(p3 a, p3 b, p3 p) {
        p3 n = a * b;
        if (n == zero)
            return a * p == zero && (a | p) > 0;
        return (n | p) == 0 && (n | a * p) >= 0 && (n | b * p) <= 0;
    }

    struct DirectionSet : vector<p3> {
        using vector::vector; // import constructors
        void insert(p3 p) {
            for (p3 q: *this) if (p * q == zero) return;
            push_back(p);
        }
    };

// Assume that the sphere is centered at the origin
// it returns the direction vectors of the intersection points
// to get the actual points, scale the direction vectors to the radius of the sphere
    DirectionSet segment_segment_intersection_on_sphere(p3 a, p3 b, p3 c, p3 d) {
        assert(validSegment(a, b) && validSegment(c, d));
        p3 out;
        if (proper_intersection(a, b, c, d, out)) return {out};
        DirectionSet s;
        if (point_on_sphere_segment(c, d, a)) s.insert(a);
        if (point_on_sphere_segment(c, d, b)) s.insert(b);
        if (point_on_sphere_segment(a, b, c)) s.insert(c);
        if (point_on_sphere_segment(a, b, d)) s.insert(d);
        return s;
    }

// small angle between spherical segments ab and ac
// assume that the sphere is centered at the origin
// all points a, b, c are on the sphere
    double angle_on_sphere(p3 a, p3 b, p3 c) {
        return get_angle(a * b, a * c);
    }

// oriented angle between spherical segments ab and ac
// that is how much we rotate counterclockwise to get from ab to ac
// assume that the sphere is centered at the origin
// all points a, b, c are on the sphere
    double oriented_angle_on_sphere(p3 a, p3 b, p3 c) {
        if ((a * b | c) >= 0) return angle_on_sphere(a, b, c);
        else return 2 * PI - angle_on_sphere(a, b, c);
    }

// Assume that the sphere is centered at the origin
// the polygon is simple and given in counterclockwise order
// for each consecutive pair of points, the counterclockwise left
// part of the segment is considered to be inside the surface area that the polygon encloses
// if the polygon is outside the sphere, the projection of the polygon on the sphere will be considered
    double area_on_the_surface_of_the_sphere(double r, vector<p3> p) {
        int n = p.size();
        double sum = -(n - 2) * PI;
        for (int i = 0; i < n; i++) {
            sum += oriented_angle_on_sphere(p[(i + 1) % n], p[(i + 2) % n], p[i]);
        }
        return r * r * sum;
    }

// Assume that O is the origin
// it returns 0 if O is outside the polyhedron
// 1 if O is inside the polyhedron, and the vector areas of the faces are oriented towards the outside;
// −1 if O is inside the polyhedron, and the vector areas of the faces are oriented towards the inside.
    int winding_number_3D(vector<vector<p3>> fs) {
        double sum = 0;
        for (vector<p3> f: fs) {
            sum += remainder(area_on_the_surface_of_the_sphere(1, f), 4 * PI);
        }
        return round(sum / (4 * PI));
    }

    struct sphere {
        p3 c;
        double r;

        sphere() {}

        sphere(p3 c, double r) : c(c), r(r) {}
    };

// spherical cap is a portion of a sphere cut off by a plane
// https://en.wikipedia.org/wiki/Spherical_cap
    struct spherical_cap {
        p3 c;
        double r;

        spherical_cap() {}

        spherical_cap(p3 c, double r) : c(c), r(r) {}

        // angle th is the polar angle between the rays from the center of the sphere to one edge of the cap
        // and orthogonal line from the center of the sphere to the plane of the cap

        // height of the cap (just like real world cap)
        double height(double th) {
            return r * (1 - cos(th));
        }

        // radius of the base of the cap
        double base_radius(double th) {
            return r * sin(th);
        }

        // volume of the cap
        double volume(double th) {
            double h = height(th);
            return PI * h * h * (3 * r - h) / 3.0;
        }

        // surface area of the cap
        double surface_area(double th) {
            double h = height(th);
            return 2 * PI * r * h;
        }
    };

// returns the sphere passing through four points
    sphere circumscribed_sphere(p3 a, p3 b, p3 c, p3 d) {
        assert(sign(plane(a, b, c).side(d)) != 0);

        plane u = plane(a - b, (a + b) / 2);
        plane v = plane(b - c, (b + c) / 2);
        plane w = plane(c - d, (c + d) / 2);

        assert(!is_parallel(u, v));
        assert(!is_parallel(v, w));
        line3d l1(u, v), l2(v, w);
        assert(sign(dist(l1, l2)) == 0);

        p3 C = closest_on_l1(l1, l2);
        return sphere(C, abs(C - a));
    }

// https://mathworld.wolfram.com/Sphere-SphereIntersection.html
// it won't work if one sphere is totally inside the other sphere
// handle that case separately
// returns the surface area and volume of the intersection
    pair<double, double> sphere_sphere_intersection(sphere s1, sphere s2) {
        double d = abs(s1.c - s2.c);
        if (sign(d - s1.r - s2.r) >= 0) return {0, 0}; // not intersecting
        // only the distance matters, so we will now consider the centers
        // of the big sphere to be (0, 0, 0) and (d, 0, 0) for the small sphere
        // we can transform the results back to w.r.t the real centers if we want

        double R = max(s1.r, s2.r);
        double r = min(s1.r, s2.r);
        double y = R + r - d;
        double x = (R * R - r * r + d * d) / (2 * d);
        // the intersecting plane is parallel to the yz plane
        // with the above x value as its x coordinate
        double w = d * d - r * r + R * R;
        double a = sqrt(4 * d * d * R * R - w * w) / (2.0 * d);
        // a is the radius of the intersecting circle on the intersecting plane
        // with center (x, 0)
        double h1 = R - x;
        double h2 = y - h1;
        // h1 is the height of the intersecting spherical cap of the big sphere
        // h2 is for the small sphere

        // total volume of the whole intersection = sum of the volumes of the spherical caps
        double volume = PI * h1 * h1 * (3 * R - h1) / 3.0 + PI * h2 * h2 * (3 * r - h2) / 3.0;
        // total surface area of the intersecting spherical caps
        double surface_area = 2 * PI * R * h1 + 2 * PI * r * h2;
        return make_pair(surface_area, volume);
    }

    sphere smallest_enclosing_sphere(vector<p3> p) {
        int n = p.size();
        p3 c(0, 0, 0);
        for (int i = 0; i < n; i++) c = c + p[i];
        c = c / n;

        double ratio = 0.1;
        int pos = 0;
        int it = 100000;
        while (it--) {
            pos = 0;
            for (int i = 1; i < n; i++) {
                if (sq(c - p[i]) > sq(c - p[pos])) pos = i;
            }
            c = c + (p[pos] - c) * ratio;
            ratio *= 0.998;
        }
        return sphere(c, abs(c - p[pos]));
    }

// it returns the angle of the spherical cap that is formed by the intersection of all tangents
    double tangent_from_point_to_sphere(p3 p, sphere s) {
        double d = abs(p - s.c);
        if (sign(d - s.r) < 0) return -1; // inside the sphere, so no tangent
        if (sign(d - s.r) == 0) return -2; // on the sphere, handle separately
        double tangent_length = sqrt(d * d - s.r * s.r);
        double th = acos(s.r / d);
        return th;
    }


    struct pyramid {
        int n;     // number of sides of the pyramid
        double l;  // length of each side
        double ang;

        pyramid(int _n, double _l) {
            n = _n;
            l = _l;
            ang = PI / n;
        }

        double base_area() {
            return l * l * n / (4 * tan(ang));
        }

        double height() {
            return l * sqrt(1 - 1 / (4 * sin(ang) * sin(ang)));
        }

        double volume() {
            return base_area() * height() / 3;
        }
    };

    struct cylinder {
        double r, h; // radius and height
        cylinder(double _r, double _h) {
            r = _r;
            h = _h;
        }

        double volume() {
            return PI * r * r * h;
        }

        double surface_area() {
            return 2 * PI * r * h + 2 * PI * r * r;
        }
    };

    struct cone {
        double r, h; // radius and height
        cone(double _r, double _h) {
            r = _r;
            h = _h;
        }

        double volume() {
            return PI * r * r * h / 3.0;
        }

        double surface_area() {
            return PI * r * (r + sqrt(h * h + r * r));
        }
    };
};

using geom::p3;
using geom::plane;

const ld EPS = eps;

void solve() {
    p3 p1, d1, p2, d2;
    cin >> p1.x >> p1.y >> p1.z;
    cin >> d1.x >> d1.y;
    d1.z = 0;
    cin >> p2.x >> p2.y >> p2.z;
    cin >> d2.x >> d2.y;
    d2.z = 0;

    plane a1(d1, p1);
    plane a2(d2, p2);
    if (a1.dist(p2) <= EPS && a2.dist(p1) <= EPS) {
        cout << "2\n";
        return;
    }

    if (0 && geom::is_parallel(a1, a2)) {
        cout << "4\n";
        return;
    }

    auto line = geom::line3d(a1, a2);


    p3 O = (p1 + p2) / 2;
    ld r = geom::dist(p1, O);

    pair<p3, p3> out;
    int res = geom::sphere_line_intersection(O, r, line, out);

    if (res != 0) {
        cout << "3\n";
		return;
    }

    cout << "4\n";

}

signed main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0);
    cout << setprecision(12) << fixed;
    int t = 1;
    cin >> t;
    rep(i, t) {
        solve();
    }
    return 0;
}

Details

Tip: Click on the bar to expand more detailed information

Test #1:

score: 0
Wrong Answer
time: 0ms
memory: 3696kb

input:

2
-1 -1 3
1 1
2 2 3
2 2
5 5 1
3 0
7 6 -2
1 -2

output:

3
3

result:

wrong answer 1st numbers differ - expected: '4', found: '3'