Editorial for Expected Bounding Box (Stretch)

Remember to use this editorial only when stuck, and not to copy-paste code from it. Please be respectful to the problem author and editorialist.
Submitting an official solution before solving the problem yourself is a bannable offence.

Approach

A location $(x, y)$ ~(x, y)~ lies in the selected bounding rectangle if and only if there is at least one selected point to its left, one to its right, one above it, and one below it. Integrating that probability over the plane gives the expected area.

The probability is constant on every open rectangle between consecutive distinct $x$ ~x~-coordinates and consecutive distinct $y$ ~y~-coordinates, so it is enough to work with compressed strips. For a strip, write $A, B, C, D$ ~A, B, C, D~ for the products of $(1 - p_i)$ ~(1 - p_i)~ over the four strict quadrants (left-down, left-up, right-down, right-up). Inclusion-exclusion on the four bad events (no selected point on a given side) yields

$P(\text{inside}) = 1 - AB - AC - BD - CD + ABC + ABD + ACD + BCD - ABCD$ .

Sweep the distinct $x$ ~x~-coordinates from left to right. Between two consecutive $x$ ~x~-values the left/right partition is fixed, and the expected vertical length is a weighted sum of the displayed formula over all $y$ ~y~-gaps. Crossing a point moves it from the right set to the left set and multiplies a $y$ ~y~-prefix or $y$ ~y~-suffix of those products. A lazy segment tree stores the ten monomials of the inclusion-exclusion, supports range multiplies on $A, B, C, D$ ~A, B, C, D~, and returns the current expected vertical length in $O(\log n)$ ~O(\log n)~.

Process every point with the same $x$ ~x~-coordinate before advancing to the next strip. The whole sweep is $O(n \log n)$ ~O(n \log n)~.

The four side-events are strongly dependent because a single point belongs to two sides at once. Multiplying the four marginal success probabilities is therefore wrong; inclusion-exclusion is required.

Solution (C++)

#include <bits/stdc++.h>
using namespace std;

const int MOD = 1000000007;

int addmod(int a, int b) {
    int s = a + b;
    if (s >= MOD) {
        s -= MOD;
    }
    return s;
}

int submod(int a, int b) {
    int s = a - b;
    if (s < 0) {
        s += MOD;
    }
    return s;
}

int mulmod(long long a, long long b) {
    return (int)(a * b % MOD);
}

int modpow(long long a, long long e) {
    long long r = 1;
    a %= MOD;
    while (e) {
        if (e & 1) {
            r = r * a % MOD;
        }
        a = a * a % MOD;
        e >>= 1;
    }
    return (int)r;
}

struct Node {
    int s1;
    int sAB, sAC, sBD, sCD;
    int sABC, sABD, sACD, sBCD;
    int sABCD;
    int mulA, mulB, mulC, mulD;
};

vector<Node> st;
int n_gaps;

void apply_one(Node &n, int mulA, int mulB, int mulC, int mulD) {
    if (mulA != 1) {
        n.sAB = mulmod(n.sAB, mulA);
        n.sAC = mulmod(n.sAC, mulA);
        n.sABC = mulmod(n.sABC, mulA);
        n.sABD = mulmod(n.sABD, mulA);
        n.sACD = mulmod(n.sACD, mulA);
        n.sABCD = mulmod(n.sABCD, mulA);
        n.mulA = mulmod(n.mulA, mulA);
    }
    if (mulB != 1) {
        n.sAB = mulmod(n.sAB, mulB);
        n.sBD = mulmod(n.sBD, mulB);
        n.sABC = mulmod(n.sABC, mulB);
        n.sABD = mulmod(n.sABD, mulB);
        n.sBCD = mulmod(n.sBCD, mulB);
        n.sABCD = mulmod(n.sABCD, mulB);
        n.mulB = mulmod(n.mulB, mulB);
    }
    if (mulC != 1) {
        n.sAC = mulmod(n.sAC, mulC);
        n.sCD = mulmod(n.sCD, mulC);
        n.sABC = mulmod(n.sABC, mulC);
        n.sACD = mulmod(n.sACD, mulC);
        n.sBCD = mulmod(n.sBCD, mulC);
        n.sABCD = mulmod(n.sABCD, mulC);
        n.mulC = mulmod(n.mulC, mulC);
    }
    if (mulD != 1) {
        n.sBD = mulmod(n.sBD, mulD);
        n.sCD = mulmod(n.sCD, mulD);
        n.sABD = mulmod(n.sABD, mulD);
        n.sACD = mulmod(n.sACD, mulD);
        n.sBCD = mulmod(n.sBCD, mulD);
        n.sABCD = mulmod(n.sABCD, mulD);
        n.mulD = mulmod(n.mulD, mulD);
    }
}

void push(int p) {
    Node &n = st[p];
    if (n.mulA == 1 && n.mulB == 1 && n.mulC == 1 && n.mulD == 1) {
        return;
    }
    apply_one(st[p << 1], n.mulA, n.mulB, n.mulC, n.mulD);
    apply_one(st[p << 1 | 1], n.mulA, n.mulB, n.mulC, n.mulD);
    n.mulA = n.mulB = n.mulC = n.mulD = 1;
}

void pull(int p) {
    Node &n = st[p];
    const Node &L = st[p << 1];
    const Node &R = st[p << 1 | 1];
    n.s1 = addmod(L.s1, R.s1);
    n.sAB = addmod(L.sAB, R.sAB);
    n.sAC = addmod(L.sAC, R.sAC);
    n.sBD = addmod(L.sBD, R.sBD);
    n.sCD = addmod(L.sCD, R.sCD);
    n.sABC = addmod(L.sABC, R.sABC);
    n.sABD = addmod(L.sABD, R.sABD);
    n.sACD = addmod(L.sACD, R.sACD);
    n.sBCD = addmod(L.sBCD, R.sBCD);
    n.sABCD = addmod(L.sABCD, R.sABCD);
}

void build(int p, int l, int r, const vector<int> &dy) {
    st[p].mulA = st[p].mulB = st[p].mulC = st[p].mulD = 1;
    if (l == r) {
        int w = dy[l];
        st[p].s1 = w;
        st[p].sAB = w;
        st[p].sAC = w;
        st[p].sBD = w;
        st[p].sCD = w;
        st[p].sABC = w;
        st[p].sABD = w;
        st[p].sACD = w;
        st[p].sBCD = w;
        st[p].sABCD = w;
        return;
    }
    int mid = (l + r) >> 1;
    build(p << 1, l, mid, dy);
    build(p << 1 | 1, mid + 1, r, dy);
    pull(p);
}

void update(int p, int l, int r, int ql, int qr, int mA, int mB, int mC, int mD) {
    if (ql > qr || qr < l || r < ql) {
        return;
    }
    if (ql <= l && r <= qr) {
        apply_one(st[p], mA, mB, mC, mD);
        return;
    }
    push(p);
    int mid = (l + r) >> 1;
    update(p << 1, l, mid, ql, qr, mA, mB, mC, mD);
    update(p << 1 | 1, mid + 1, r, ql, qr, mA, mB, mC, mD);
    pull(p);
}

void range_mul(int ql, int qr, int mA, int mB, int mC, int mD) {
    if (n_gaps <= 0 || ql > qr) {
        return;
    }
    update(1, 0, n_gaps - 1, ql, qr, mA, mB, mC, mD);
}

int vertical_expectation() {
    if (n_gaps <= 0) {
        return 0;
    }
    const Node &n = st[1];
    int ans = n.s1;
    ans = submod(ans, n.sAB);
    ans = submod(ans, n.sAC);
    ans = submod(ans, n.sBD);
    ans = submod(ans, n.sCD);
    ans = addmod(ans, n.sABC);
    ans = addmod(ans, n.sABD);
    ans = addmod(ans, n.sACD);
    ans = addmod(ans, n.sBCD);
    ans = submod(ans, n.sABCD);
    return ans;
}

struct Point {
    int x, y, q, yid;
};

bool cmp_x(const Point &a, const Point &b) {
    return a.x < b.x;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    int n;
    cin >> n;
    vector<Point> pts(n);
    vector<int> ys;
    ys.reserve(n);
    for (int i = 0; i < n; i++) {
        int x, y, a, b;
        cin >> x >> y >> a >> b;
        int invb = modpow(b, MOD - 2);
        int q = mulmod((b - a) % MOD, invb);
        pts[i].x = x;
        pts[i].y = y;
        pts[i].q = q;
        ys.push_back(y);
    }

    sort(ys.begin(), ys.end());
    ys.erase(unique(ys.begin(), ys.end()), ys.end());
    n_gaps = (int)ys.size() - 1;
    if (n_gaps <= 0) {
        cout << 0 << "\n";
        return 0;
    }

    vector<int> dy(n_gaps);
    for (int i = 0; i < n_gaps; i++) {
        long long d = (long long)ys[i + 1] - ys[i];
        d %= MOD;
        if (d < 0) {
            d += MOD;
        }
        dy[i] = (int)d;
    }

    for (int i = 0; i < n; i++) {
        pts[i].yid = (int)(lower_bound(ys.begin(), ys.end(), pts[i].y) - ys.begin());
    }

    st.assign(4 * n_gaps + 4, Node());
    build(1, 0, n_gaps - 1, dy);

    for (int i = 0; i < n; i++) {
        int yid = pts[i].yid;
        int q = pts[i].q;
        range_mul(yid, n_gaps - 1, 1, 1, q, 1);
        range_mul(0, yid - 1, 1, 1, 1, q);
    }

    sort(pts.begin(), pts.end(), cmp_x);

    int ans = 0;
    bool has_prev = false;
    int prev_x = 0;
    int i = 0;
    while (i < n) {
        int j = i;
        while (j < n && pts[j].x == pts[i].x) {
            j++;
        }
        if (has_prev) {
            long long dx = (long long)pts[i].x - prev_x;
            dx %= MOD;
            if (dx < 0) {
                dx += MOD;
            }
            ans = addmod(ans, mulmod((int)dx, vertical_expectation()));
        }
        has_prev = true;
        prev_x = pts[i].x;
        for (int k = i; k < j; k++) {
            int yid = pts[k].yid;
            int q = pts[k].q;
            int iq = modpow(q, MOD - 2);
            range_mul(yid, n_gaps - 1, q, 1, iq, 1);
            range_mul(0, yid - 1, 1, q, 1, iq);
        }
        i = j;
    }

    cout << ans << "\n";
}

Comments

There are no comments at the moment.