ICPC 2024 - H. Maxwell’s Demon

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

namespace {

struct TimeFrac {
    int scaled = 0;
    int denom = 1;
};

bool operator<(const TimeFrac& a, const TimeFrac& b) {
    return 1LL * a.scaled * b.denom < 1LL * b.scaled * a.denom;
}

bool same_time(const TimeFrac& a, const TimeFrac& b) {
    return 1LL * a.scaled * b.denom == 1LL * b.scaled * a.denom;
}

struct XorVector {
    array<unsigned long long, 4> words{};

    bool any() const {
        return words[0] || words[1] || words[2] || words[3];
    }
};

void xor_with(XorVector& a, const XorVector& b) {
    for (int i = 0; i < 4; ++i) {
        a.words[i] ^= b.words[i];
    }
}

int highest_bit(const XorVector& v) {
    for (int block = 3; block >= 0; --block) {
        if (v.words[block] == 0) {
            continue;
        }
        return block * 64 + 63 - __builtin_clzll(v.words[block]);
    }
    return -1;
}

long long extended_gcd(long long a, long long b, long long& x, long long& y) {
    if (b == 0) {
        x = (a >= 0 ? 1 : -1);
        y = 0;
        return llabs(a);
    }
    long long x1 = 0;
    long long y1 = 0;
    const long long g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

long long gcd_ll(long long a, long long b) {
    a = llabs(a);
    b = llabs(b);
    while (b != 0) {
        const long long r = a % b;
        a = b;
        b = r;
    }
    return a;
}

long long mod_inverse(long long a, long long mod) {
    long long x = 0;
    long long y = 0;
    const long long g = extended_gcd(a, mod, x, y);
    if (g != 1) {
        return -1;
    }
    x %= mod;
    if (x < 0) {
        x += mod;
    }
    return x;
}

long long positive_mod(long long value, long long mod) {
    value %= mod;
    if (value < 0) {
        value += mod;
    }
    return value;
}

vector<int> generate_times(const int w,
                           const int h,
                           const int d,
                           const int px,
                           const int py,
                           const int vx,
                           const int vy) {
    const int q0 = abs(px);
    const int qv = (px > 0 ? vx : -vx);
    const int ax = abs(qv);
    if (ax == 0) {
        return {};
    }

    const int s0 = (qv < 0 ? q0 : 2 * w - q0);
    const long long limit_scaled = 1LL * 2 * w * h * ax;
    vector<int> hits;

    vector<int> targets = {d};
    if (2 * h - d != d) {
        targets.push_back(2 * h - d);
    }

    for (const int target_y : targets) {
        const long long A = 1LL * 2 * w * vy;
        const long long M = 1LL * 2 * h * ax;
        const long long B =
            1LL * ax * (target_y - py) - 1LL * vy * s0;

        long long first_k = -1;
        long long step_k = -1;

        if (A == 0) {
            if (positive_mod(B, M) == 0) {
                first_k = 0;
                step_k = 1;
            }
        } else {
            const long long g = gcd_ll(A, M);
            if (B % g == 0) {
                const long long A2 = A / g;
                const long long B2 = B / g;
                const long long M2 = M / g;
                const long long inv = mod_inverse(positive_mod(A2, M2), M2);
                first_k = positive_mod(inv * positive_mod(B2, M2), M2);
                step_k = M2;
            }
        }

        if (first_k == -1) {
            continue;
        }

        for (long long k = first_k;; k += step_k) {
            const long long scaled = 1LL * s0 + 1LL * 2 * w * k;
            if (scaled >= limit_scaled) {
                break;
            }
            hits.push_back(static_cast<int>(scaled));
            if (step_k == 0) {
                break;
            }
        }
    }

    sort(hits.begin(), hits.end());
    hits.erase(unique(hits.begin(), hits.end()), hits.end());
    return hits;
}

bool insert_vector(XorVector value, array<XorVector, 256>& basis) {
    while (value.any()) {
        const int bit = highest_bit(value);
        if (!basis[bit].any()) {
            basis[bit] = value;
            return true;
        }
        xor_with(value, basis[bit]);
    }
    return false;
}

bool in_span(XorVector value, const array<XorVector, 256>& basis) {
    while (value.any()) {
        const int bit = highest_bit(value);
        if (!basis[bit].any()) {
            return false;
        }
        xor_with(value, basis[bit]);
    }
    return true;
}

void solve() {
    int w;
    int h;
    int d;
    int red_count;
    int blue_count;
    cin >> w >> h >> d >> red_count >> blue_count;

    const int particle_count = red_count + blue_count;
    vector<vector<int>> hit_times(particle_count);
    vector<int> denominators(particle_count, 1);
    XorVector target;

    for (int i = 0; i < particle_count; ++i) {
        int px;
        int py;
        int vx;
        int vy;
        cin >> px >> py >> vx >> vy;

        const int initial_side = (px > 0 ? 1 : 0);
        const int desired_side = (i < red_count ? 0 : 1);
        if (initial_side != desired_side) {
            target.words[i / 64] |= 1ULL << (i % 64);
        }

        const int qv = (px > 0 ? vx : -vx);
        denominators[i] = abs(qv);
        hit_times[i] = generate_times(w, h, d, px, py, vx, vy);
    }

    if (!target.any()) {
        cout << "0.0\n";
        return;
    }

    struct HeapNode {
        int particle = 0;
        int index = 0;
    };
    auto cmp = [&](const HeapNode& lhs, const HeapNode& rhs) {
        const TimeFrac a{hit_times[lhs.particle][lhs.index],
                         denominators[lhs.particle]};
        const TimeFrac b{hit_times[rhs.particle][rhs.index],
                         denominators[rhs.particle]};
        return b < a;
    };

    priority_queue<HeapNode, vector<HeapNode>, decltype(cmp)> pq(cmp);
    for (int i = 0; i < particle_count; ++i) {
        if (!hit_times[i].empty()) {
            pq.push({i, 0});
        }
    }

    array<XorVector, 256> basis;
    while (!pq.empty()) {
        const HeapNode first = pq.top();
        pq.pop();
        const TimeFrac current{hit_times[first.particle][first.index],
                               denominators[first.particle]};

        XorVector event;
        event.words[first.particle / 64] |= 1ULL << (first.particle % 64);

        if (first.index + 1 < static_cast<int>(hit_times[first.particle].size())) {
            pq.push({first.particle, first.index + 1});
        }

        while (!pq.empty()) {
            const HeapNode next = pq.top();
            const TimeFrac next_time{hit_times[next.particle][next.index],
                                     denominators[next.particle]};
            if (!same_time(current, next_time)) {
                break;
            }
            pq.pop();
            event.words[next.particle / 64] |= 1ULL << (next.particle % 64);
            if (next.index + 1 <
                static_cast<int>(hit_times[next.particle].size())) {
                pq.push({next.particle, next.index + 1});
            }
        }

        if (insert_vector(event, basis) && in_span(target, basis)) {
            cout << fixed << setprecision(10)
                 << static_cast<long double>(current.scaled) /
                        static_cast<long double>(current.denom)
                 << '\n';
            return;
        }
    }

    cout << "impossible\n";
}

}  // namespace

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

    solve();
    return 0;
}

\documentclass[11pt]{article}
\usepackage[margin=1in]{geometry}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{enumitem}

\title{ICPC World Finals 2024\\H. Maxwell’s Demon}
\author{}
\date{}

\begin{document}
\maketitle

\section*{Problem Summary}

Particles move inside two mirrored chambers. Whenever one or more particles hit the demon's position on the center wall, the demon may either let all of them pass or reflect all of them. Passing flips the side of every particle in that simultaneous hit set. We must find the earliest time when the demon can make all red particles end up on the left and all blue particles on the right.

\section*{Key Observations}

\begin{itemize}[leftmargin=*]
    \item Whether a particle reflects at the center wall or passes through it does not change \emph{when} it returns to the demon. It only flips which chamber the particle is in.
    \item Therefore each time instant when some particles are at the demon gives one XOR operation on the side vector: for every particle present at that time, its left/right bit is toggled.
    \item The whole problem becomes linear algebra over $\mathbb{F}_2$:
    \begin{itemize}[leftmargin=*]
        \item the initial state is the current side of every particle;
        \item the target state is ``red = left, blue = right'';
        \item each usable demon time contributes one bit vector whose $i$th bit is $1$ iff particle $i$ is at the demon at that time.
    \end{itemize}
    \item Because all positions and velocities are integers, every particle state repeats after time $2wh$: both unfolded coordinates advance by multiples of the chamber periods. So it is enough to examine hit times in $[0,2wh)$.
\end{itemize}

\section*{Algorithm}

\begin{enumerate}[leftmargin=*]
    \item For each particle, replace the chamber side by its distance to the center wall. This gives a 1D reflected motion on $[0,w]$ that is independent of the demon's previous choices.

    \item Let $q_0$ be the initial distance to the center wall and let $u$ be the signed velocity in this 1D model. The particle hits the center wall at times
    \[
    t = t_0 + k \cdot \frac{2w}{|u|},
    \]
    where $t_0 = \frac{q_0}{|u|}$ if it is moving toward the center and
    \[
    t_0 = \frac{2w-q_0}{|u|}
    \]
    otherwise.

    \item For the vertical coordinate, use the standard unfolded reflection model on period $2h$. A particle is at height $d$ exactly when the unfolded coordinate is congruent to $d$ or $2h-d$ modulo $2h$.

    \item Combining the horizontal hit progression with the vertical congruence gives a linear congruence in the hit index $k$. Solve it with extended Euclid. This yields at most two arithmetic progressions of valid demon-hit times for each particle. Generate all such times in $[0,2wh)$.

    \item Merge the per-particle time lists in increasing order. For each distinct time, build one bit vector describing the set of particles present at the demon.

    \item Maintain a Gaussian-elimination basis over $\mathbb{F}_2$ for all event vectors seen so far. After each new independent event vector is inserted, test whether the target side-change vector lies in the current span. The first time this happens is the answer.
\end{enumerate}

\section*{Correctness Proof}

We prove that the algorithm returns the correct answer.

\paragraph{Lemma 1.}
For every particle, the set of times at which it is located at the demon is independent of all previous pass/reflect choices made by the demon.

\paragraph{Proof.}
Passing through the center wall and reflecting at the center wall only change which chamber the particle occupies; they do not change its absolute distance from the center wall or its vertical motion. Hence the future times at which it returns to the point $(0,d)$ are fixed in advance. \qed

\paragraph{Lemma 2.}
At any fixed time $t$, allowing particles through the demon toggles exactly the side bits of the particles present at time $t$, and leaves all other bits unchanged.

\paragraph{Proof.}
By the problem statement, at time $t$ either all particles at the demon pass through or all reflect. Passing switches a particle from one chamber to the other, i.e. toggles its side bit. Particles not at the demon are unaffected. \qed

\paragraph{Lemma 3.}
After processing all demon times up to time $T$, the set of reachable side vectors is exactly the affine space obtained from the initial side vector by adding the span of all event vectors up to time $T$.

\paragraph{Proof.}
By Lemma 2, each demon action adds one event vector over $\mathbb{F}_2$, and actions at different times commute because XOR is commutative. Therefore the reachable side changes are exactly all XOR sums of available event vectors, i.e. their span. \qed

\paragraph{Theorem.}
The algorithm outputs the correct answer for every valid input.

\paragraph{Proof.}
By Lemma 1, the algorithm enumerates the correct event vectors. By Lemma 3, the target arrangement is reachable by time $T$ if and only if the target difference vector is in the span of event vectors up to $T$. The algorithm processes times in increasing order and stops at the first time this condition becomes true, so the reported time is exactly the earliest feasible one. If it never becomes true during one full period $2wh$, then no new event vectors will ever appear later, so the arrangement is impossible. \qed

\section*{Complexity Analysis}

Let $n=r+b \le 200$. Each particle has $O(wh)$ hits in one period in the worst case, and we only examine one full period. The total number of generated events is therefore manageable. Gaussian elimination uses dimension at most $200$, while the basis rank is also at most $200$. The memory usage is dominated by the per-particle hit lists.

\section*{Implementation Notes}

\begin{itemize}[leftmargin=*]
    \item The code stores times as exact fractions using a particle-specific denominator, and compares them by cross multiplication to avoid floating-point ordering errors.
    \item The answer is printed as a floating-point number only after the exact earliest rational time has been found.
    \item If the initial side vector already equals the target vector, the answer is immediately $0$.
\end{itemize}

\end{document}