Mst oriented in O(NM)

.gitignore

@@ -0,0 +1,5 @@
+polyup
+url
+url_internal
+*.pdf
+*.out

algo/mst_oriented/01.in

@@ -0,0 +1,6 @@
+5 5 0
+0 1 5
+0 2 3
+2 3 2
+3 4 2
+4 2 2

algo/mst_oriented/02.in

@@ -0,0 +1,12 @@
+9 11 0
+0 1 4
+2 6 4
+1 2 3
+2 3 3
+3 4 3
+4 5 3
+5 1 3
+4 7 1
+7 5 1
+5 4 1
+7 8 2

algo/mst_oriented/mst_oriented.cpp

@@ -0,0 +1,202 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+
+using namespace std;
+using lli = long long int;
+
+using pii   = tuple<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;
+
+/** WARNING **
+ *************
+
+ This algorithm considers the graph to have no self-pointing edge, ie. no edge
+ (u, u)
+*/
+
+const int INFTY = 1000*1000*1000; // FIXME large enough for this problem?
+
+struct Edge {
+    int from, to;
+    int weight;
+    Edge() {}
+    Edge(int from, int to, int weight): from(from), to(to), weight(weight) {}
+};
+
+struct Vert {
+    int selfId;
+    int ufFather;
+    int ufHeight;
+    vector<int> out;
+    vector<int> in;
+    int parentEdge;
+    bool suppressed;
+
+    Vert(int id):
+        selfId(id), ufFather(id), ufHeight(1),
+        parentEdge(-1), suppressed(false)
+    {}
+};
+
+int root;
+vector<Edge> edges;
+vector<Vert> vertices;
+int cost = 0;
+
+int ufFind(int vert) {
+    Vert& v = vertices[vert];
+    if(v.ufFather == vert)
+        return vert;
+    v.ufFather = ufFind(v.ufFather);
+    return v.ufFather;
+}
+
+void ufUnion(int iv1, int iv2) {
+    iv1 = ufFind(iv1);
+    iv2 = ufFind(iv2);
+    Vert& v1 = vertices[iv1];
+    Vert& v2 = vertices[iv2];
+    if(v1.ufHeight < v2.ufHeight)
+        v1.ufFather = iv2;
+    else if(v2.ufHeight < v1.ufHeight)
+        v2.ufFather = iv1;
+    else {
+        v1.ufFather = iv2;
+        v2.ufHeight++;
+    }
+}
+
+void buildMst(int vertId) {
+    printf("> %d\n", vertId);
+    Vert& vert = vertices[vertId];
+    int minIncWeight = INFTY;
+    int minId = -1;
+
+    for(int incEdge: vert.in) {
+        Edge& edge = edges[incEdge];
+        if(minIncWeight > edge.weight) {
+            minIncWeight = edge.weight;
+            minId = incEdge;
+        }
+    }
+
+    if(minId == -1) {
+        assert(vertId == root);
+        return;
+    }
+
+    for(int incEdge: vert.in) {
+        Edge& edge = edges[incEdge];
+        edge.weight -= minIncWeight;
+    }
+    cost += minIncWeight;
+
+    vert.parentEdge = minId;
+    int parent = edges[minId].from;
+    int parentCmp = ufFind(parent), selfCmp = ufFind(vertId);
+    if(parentCmp == selfCmp) {
+        // Build current cycle
+        vector<int> cycle;
+        cycle.push_back(parent);
+        while(cycle.back() != vertId)
+            cycle.push_back(edges[vertices[cycle.back()].parentEdge].from);
+
+        vertices.push_back(Vert(vertices.size()));
+        Vert& condensed = vertices.back();
+        printf("Compressing into %d: ", condensed.selfId);
+        for(int cyc: cycle)
+            printf("%d ", cyc);
+        puts("");
+        vector<int> nInEdges, nOutEdges;
+        for(int cycleVert: cycle) {
+            Vert& v = vertices[cycleVert];
+            v.suppressed = true;
+            for(int inEdge: v.in) {
+                edges[inEdge].to = condensed.selfId;
+                nInEdges.push_back(inEdge);
+            }
+            for(int outEdge: v.out) {
+                edges[outEdge].from = condensed.selfId;
+                if(edges[outEdge].to != condensed.selfId)
+                    condensed.out.push_back(outEdge);
+            }
+        }
+
+        for(int inEdge: nInEdges) {
+            if(edges[inEdge].from != condensed.selfId)
+                condensed.in.push_back(inEdge);
+        }
+
+        ufUnion(vertId, condensed.selfId);
+        buildMst(condensed.selfId);
+    }
+    else {
+        ufUnion(vertId, parent);
+        if(ufFind(vertId) != ufFind(root))
+            buildMst(parent);
+    }
+}
+
+void dfs(int pos, vector<bool>& seen) {
+    if(seen[pos])
+        return;
+    seen[pos] = true;
+    for(int outEdge: vertices[pos].out)
+        dfs(edges[outEdge].to, seen);
+}
+
+bool allReachable() {
+    vector<bool> seen(vertices.size(), false);
+    dfs(root, seen);
+    for(size_t pos=0; pos < seen.size(); ++pos) {
+        if(!seen[pos]) {
+            fprintf(stderr, "NOT REACHABLE: %lu\n", pos);
+            return false;
+        }
+    }
+    return true;
+}
+
+int main() {
+    int nbVert, nbEdge;
+    scanf("%d%d%d", &nbVert, &nbEdge, &root);
+    for(int vert=0; vert < nbVert; ++vert)
+        vertices.push_back(Vert(vert));
+    for(int edgeId=0; edgeId < nbEdge; ++edgeId) {
+        Edge edge;
+        scanf("%d%d%d", &edge.from, &edge.to, &edge.weight);
+        edges.push_back(edge);
+        vertices[edge.from].out.push_back(edges.size() - 1);
+        vertices[edge.to].in.push_back(edges.size() - 1);
+    }
+
+    if(!allReachable()) {
+        printf("-1\n");
+        puts("BLEH");
+        return 0;
+    }
+
+    for(int vert=0; vert < nbVert; ++vert) {
+        if(vertices[vert].suppressed || ufFind(vert) == root)
+            continue;
+        buildMst(vert);
+    }
+
+    printf("%d\n", cost);
+}

day2/j/01.in

@@ -0,0 +1,5 @@
+3 2 11
+1
+3
+1
+

day2/j/j.cpp

@@ -0,0 +1,175 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+#define FORU(i, a, b) for(lli i = (lli)(a); i < (lli)(b); ++i)
+#define FORD(i, a, b) for(lli i = (lli)(b)-1; i >= (lli)(a); --i)
+#define ALL(x) (x).begin(), (x).end()
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+#define pb push_back
+
+
+
+using namespace std;
+using lli = long long int;
+
+using pii   = tuple<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;
+
+const lli INFTY = ((1ll << 60) - 1) * 2;
+
+// =====
+
+struct Evt {
+    enum {
+        JACKPOT, ENTRY
+    } type;
+
+    lli date;
+    lli nbPers;
+
+    bool operator<(const Evt& e) const {
+        if(date == e.date)
+            return type == Evt::ENTRY;
+        return date < e.date;
+    }
+};
+
+lli gcd(lli a, lli b) {
+    if(b == 0)
+        return a;
+    return gcd(b, a % b);
+}
+
+struct Frac {
+    Frac() : num(0), den(1) {}
+    Frac(lli n, lli d): num(n), den(d) { simpl(); }
+    lli num, den;
+
+    void simpl() {
+        lli g = gcd(num, den);
+        num /= g;
+        den /= g;
+    }
+
+    /*
+    Frac operator+(const Frac& b) {
+        return Frac(b.den * num + den * b.num, den * b.den);
+    }
+    */
+
+    double val() const {
+        return ((double)num / (double)den);
+    }
+
+    Frac& operator*=(const Frac& b) {
+        num *= b.num;
+        den *= b.den;
+        simpl();
+        return *this;
+    }
+
+    double operator*(double b) const {
+        return b * val();
+    }
+};
+
+vector<lli> joinTimes;
+vector<Evt> jackpots;
+lli ticketPrice, jackpotVal;
+map<lli, double> expec;
+
+void getJackpots() {
+    lli slope = 0, cMoney = 0;
+    vector<lli> orderedJoins(joinTimes);
+    sort(orderedJoins.begin(), orderedJoins.end());
+    orderedJoins.push_back(INFTY);
+
+    lli pEvent = 0;
+
+    for(size_t j = 0; j < orderedJoins.size() - 1; ++j) {
+        lli join = orderedJoins[j];
+        lli nEvent = orderedJoins[j+1];
+        cMoney += (join - pEvent) * slope;
+        slope += ticketPrice;
+        pEvent = join;
+        while(cMoney + (nEvent - pEvent - 1) * slope >= jackpotVal) {
+            jackpots.push_back(Evt({Evt::JACKPOT,
+                    (jackpotVal - cMoney + slope - 1) / slope + pEvent,
+                    0}));
+            pEvent = jackpots.back().date;
+            slope -= ticketPrice;
+            cMoney = 0;
+        }
+    }
+}
+
+void fillExpec() {
+    vector<Evt> events(jackpots);
+    for(lli join: joinTimes)
+        events.push_back(Evt({Evt::ENTRY, join, 0}));
+    sort(events.begin(), events.end());
+    int cPers = 0;
+    for(Evt& evt: events) {
+        cPers += (evt.type == Evt::ENTRY) ? 1 : -1;
+        evt.nbPers = cPers + ((evt.type == Evt::JACKPOT)?1:0);
+    }
+
+    double cExpec = 0;
+    Frac cProd;
+    for(auto evtIt = events.begin(); evtIt != events.end(); ++evtIt) {
+        const Evt& evt = *evtIt;
+        if(evt.nbPers == 1 && evt.type == Evt::ENTRY) {
+            cProd = Frac(1, 1);
+            cExpec = 0;
+            for(auto cJack = evtIt; cJack != events.end(); ++cJack) {
+                if(cJack->type != Evt::JACKPOT)
+                    continue;
+
+                cExpec +=
+                    cProd * ((double)cJack->date) / ((double)cJack->nbPers);
+                cProd *= Frac(cJack->nbPers - 1, cJack->nbPers);
+
+                if(cJack->nbPers == 1)
+                    break;
+            }
+        }
+
+        if(evt.type == Evt::JACKPOT) {
+            cExpec -= ((double)evt.date) / ((double)evt.nbPers);
+            cExpec *= ((double)evt.nbPers) / ((double)evt.nbPers - 1.);
+        }
+        else {
+            expec.insert({evt.date, ticketPrice * (cExpec - evt.date)});
+        }
+    }
+}
+
+int main() {
+    int nbPlayers;
+    scanf("%d%lld%lld", &nbPlayers, &ticketPrice, &jackpotVal);
+    ticketPrice /= 2;
+    joinTimes.resize(nbPlayers);
+    FOR(i, nbPlayers)
+        scanf("%lld", &(joinTimes[i]));
+    getJackpots();
+    fillExpec();
+
+    for(const lli& joinTime: joinTimes) {
+        printf("%.15lf\n", expec[joinTime]);
+    }
+}
+

day3/g/01.in

@@ -0,0 +1,3 @@
+2 3
+WWB
+BBB

day3/g/03.in

@@ -0,0 +1,4 @@
+3 3
+BBB
+WBW
+WBW

day3/g/2.in

@@ -0,0 +1,4 @@
+3 3
+WWW
+WBW
+WWW

day3/g/g.cpp

@@ -0,0 +1,135 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+#define FORU(i, a, b) for(lli i = (lli)(a); i < (lli)(b); ++i)
+#define FORD(i, a, b) for(lli i = (lli)(b)-1; i >= (lli)(a); --i)
+#define ALL(x) (x).begin(), (x).end()
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+#define pb push_back
+
+
+
+using namespace std;
+using lli = long long int;
+
+using pii   = tuple<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;
+// =====
+
+const int DIM = 50 * 50 + 1;
+int nbRows, nbCols, nbDim;
+
+typedef bitset<DIM> Matr;
+
+struct PermVect {
+    vector<Matr> v;
+    vector<int> lineId;
+
+    PermVect(){}
+    PermVect(size_t n) {
+        v.resize(n);
+        lineId.resize(n);
+        FOR(i, lineId.size()) lineId[i] = i;
+    }
+
+    Matr& operator[](size_t pos) { return v[lineId[pos]]; }
+    const Matr& operator[](size_t pos) const { return v[lineId[pos]]; }
+    void swapRow(size_t i, size_t j) {
+        swap(lineId[i], lineId[j]);
+    }
+};
+
+PermVect crosses;
+
+inline size_t at(int row, int col) {
+    return row * nbCols + col;
+}
+
+void dump() {
+    FOR(rowId, nbDim) {
+        const Matr& row = crosses[rowId];
+        FOR(pos, nbDim)
+            printf("%c", row[pos]?'1':'0');
+        printf("  %c\n", row[nbDim]?'1':'0');
+    }
+    puts("===");
+}
+
+bool solve() {
+    for(int trigAt=0; trigAt < nbDim; ++trigAt) {
+        bool found = false;
+        for(int dim=trigAt; dim < nbDim; ++dim) {
+            if(!crosses[dim][trigAt])
+                continue;
+            found = true;
+            crosses.swapRow(trigAt, dim);
+            for(int nDim=trigAt+1; nDim < nbDim; ++nDim) {
+                if(crosses[nDim][trigAt])
+                    crosses[nDim] ^= crosses[trigAt];
+            }
+            break;
+        }
+        if(!found)
+            return false;
+    }
+    for(int dim=nbDim-1; dim >= 0; --dim) {
+        for(int nDim=dim-1; nDim >= 0; --nDim) {
+            if(crosses[nDim][dim])
+                crosses[nDim] ^= crosses[dim];
+        }
+    }
+    return true;
+}
+
+int main(void) {
+    scanf("%d%d", &nbRows, &nbCols);
+    nbDim = nbRows * nbCols - 1;
+    crosses = PermVect(nbDim);
+    FOR(row, nbRows) {
+        FOR(col, nbCols) {
+            FOR(i, nbCols)
+                crosses[at(row, col)].set(at(row, i));
+            FOR(i, nbRows)
+                crosses[at(row, col)].set(at(i, col));
+        }
+    }
+
+    FOR(row, nbRows) {
+        getchar();
+        FOR(col, nbCols)
+            crosses[at(row, col)][nbDim] = getchar() == 'B';
+    }
+
+
+    if(!solve()) {
+        dump();
+        puts("No");
+    }
+    else {
+        puts("Yes");
+        set<int> ids;
+        FOR(row, nbDim)
+            if(crosses[row][nbDim])
+                ids.insert(row);
+        printf("%lu\n", ids.size());
+        FOR(i, nbDim) {
+            if(ids.find(crosses.lineId[i]) != ids.end())
+                printf("%lld %lld\n", (i/nbCols) + 1, (i%nbCols) + 1);
+        }
+    }
+    return 0;
+}

day3/g/g.orig.cpp

@@ -0,0 +1,212 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+#define mp make_pair
+#define pb push_back
+
+#define fst first
+#define snd second
+
+
+using namespace std;
+using lli = long long int;
+using ll = long long int;
+
+using pii   = tuple<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;
+
+
+template<int N> // N = dimension
+struct gauss {
+  bitset<N> A[N];     // upper diagonal matrix (free family)
+  int V[N]; int nv=0; // identify the columns
+  bitset<N> B[N];     // how elements of A are generated ?
+
+  gauss() { reset(); }
+  void reset() {
+    // memset(this,0,sizeof(gauss<N>)); should work as well
+    FOR(i,N) A[i].reset();
+    nv=0;
+    FOR(i,N) B[i].reset();
+  }
+
+  // add a matrix column (a vector in the free family)
+  void add(int id, bitset<N> x) {
+    if(nv == N) return;
+    bitset<N> y; y[nv]=1;
+    FOR(i,N) if(x[i]) {
+      if(A[i][i]) { x ^= A[i]; y ^= B[i]; }
+      else { A[i] = x; V[nv++] = id; B[i] = y; return; }
+    }
+  }
+
+  // reach a vector x
+  // returns true when possible, false otherwise
+  // y contains the vectors of V used
+  bool solve(bitset<N> x, bitset<N> &y){
+    y.reset();
+    FOR(i,N) if(x[i]&&A[i][i]) {
+      x ^= A[i];
+      y ^= B[i];
+    }
+    if(x.any()) return false;
+    return true;
+  }
+};
+
+// example usage : solve a linear system in Z/2Z
+// M is a matrix (list of columns)
+// X is a vector
+// solves M Y = X
+// returns in Y the indices of columns of M with xor X
+const int DIM=2500;
+bool solve(vector<bitset<DIM> > M, bitset<DIM> X, vi &Y) {
+  gauss<DIM> G;
+  FOR(i,M.size()) {
+    G.add(i, M[i]);
+  }
+  bitset<DIM> y;
+  if (!G.solve(X, y))
+      return false;
+  for (int i = 0; i < DIM; ++i)
+      if (y[i])
+        Y.pb(G.V[i]);
+  return true;
+}
+
+
+// Inspire de https://github.com/stjepang/snippets/blob/master/gauss.cpp
+// Elimination de Gauss
+// Resout un systeme d'equations lineaires
+// Complexite : O(nb_lins * nb_cols^2)
+// Si le systeme a au moins une solution, value contiendra une solution possible
+const int MAX_NB_COLS = 2500;
+const double eps = 1e-8;
+
+struct Gauss {
+  // posnz[i] = -1 si la i-eme composante est libre
+  int posnz[MAX_NB_COLS];
+  // value[i] = la valeur de X(i) verifiant l'equation ci-dessous 
+  int value[MAX_NB_COLS];
+  // vrai ssi. le systeme a >= 1 solution
+  bool has_solution;
+
+  // mat[0 .. nb_lins-1][0 .. nb_cols-1] * X = mat[0 .. nb_lins-1][nb_cols]
+  Gauss(int mat[][MAX_NB_COLS + 1], int nb_lins, int nb_cols) {
+    fill(posnz, posnz + nb_cols, -1);
+
+
+    for (int i = 0; i < nb_lins; ++i, puts(""))
+        for (int j = 0; j < nb_cols; ++j)
+        printf("%d ", mat[i][j]);
+    int posnz_cur = 0; 
+    for (int col = 0; col < nb_cols; ++col) {
+      int max_lin = posnz_cur;
+
+      for (int lin = max_lin + 1; lin < nb_lins; ++lin)
+        if (fabs(mat[lin][col]) > fabs(mat[max_lin][col]))
+          max_lin = lin;
+
+      // La colonne est nulle
+      // Condition de la forme == 0 si on est dans les entiers
+      if (mat[max_lin][col] == 0) continue;
+
+      for (int i = 0; i <= nb_cols; ++i)
+        swap(mat[max_lin][i], mat[posnz_cur][i]);
+
+      for (int lin = 0; lin < nb_lins; ++lin) {
+        if (lin == posnz_cur) continue;
+        // Pour Gauss modulaire : remplacer par l'inverse de mat[posnz_cur][col]
+        int factor = mat[lin][col] * (1 ^ mat[posnz_cur][col]);
+        for (int i = 0; i <= nb_cols; ++i) {
+          mat[lin][i] -= factor * mat[posnz_cur][i];
+          mat[lin][i] %= 2;
+          mat[lin][i] += 2;
+          mat[lin][i] %= 2;
+        }
+          // Gauss mod : rajouter le modulo
+      }
+
+      posnz[col] = posnz_cur++;
+    }
+
+    // Genere une solution valide
+    for (int col = 0; col < nb_cols; ++col) {
+      if (posnz[col] != -1)
+        value[col] = mat[posnz[col]][nb_cols] * (1 ^ mat[posnz[col]][col]);
+      // Gauss mod
+      else 
+        value[col] = 0;
+    }
+
+    for (int i = 0; i < nb_lins; ++i, puts(""))
+        for (int j = 0; j < nb_cols; ++j)
+        printf("%d ", mat[i][j]);
+
+    // Verifie que la solution generee est valide
+    has_solution = true;
+    for (int lin = 0; lin < nb_lins; ++lin) {
+      int sum = 0;
+      for (int col = 0; col < nb_cols; ++col) {
+        sum += mat[lin][col] * value[col]; // Gauss mod
+        sum %= 2;
+      }
+      if (sum - mat[lin][nb_cols] != 0) // Gauss mod
+        has_solution = false;
+    }
+  }
+};
+
+
+const int MAXN = 50, MAXM = 50;
+int grid[MAXN][MAXM];
+int mat[MAXN * MAXM][MAXN * MAXM + 1];
+int N, M;
+
+int linof(int i) { return i / M; }
+int colof(int i) { return i % M; }
+
+int main() {
+    scanf("%d%d", &N, &M); 
+    for (int i = 0; i < N; ++i)
+        for (int j = 0; j < M; ++j) {
+           char c;
+           scanf(" %c", &c);
+           grid[i][j] = c == 'W';
+        }
+    vector<bitset<DIM>> mat(DIM);
+    bitset<DIM> X;
+    for (int i = 0; i < N * M; ++i)
+        for (int j = 0; j < N * M; ++j)
+            mat[i][j] = linof(i) == linof(j) || colof(i) == colof(j);
+    for (int i = 0; i < N * M; ++i)
+        X[i] = grid[linof(i)][colof(i)] ^ 1;
+    /*
+    for (int i = 0; i < N * M; ++i)
+        mat[i][N * M] = grid[linof(i)][colof(i)] ^ 1;
+    */
+    vi Y;
+    if (solve(mat, X, Y)) {
+        puts("Yes");
+        printf("%d\n", (int)Y.size());
+        for (int x: Y)
+            printf("%d %d\n", linof(x) + 1, colof(x) + 1);
+    } else
+        puts("No");
+    return 0;
+}

day3/g/gen.py

@@ -0,0 +1,7 @@
+from random import randint
+
+print("50 50")
+for row in range(50):
+    for col in range(50):
+        print('B' if randint(0, 1) else 'W', end='')
+    print('')

day3/g/large.in

@@ -0,0 +1,51 @@
+50 50
+BWBBWWBWWBBBBWWWWBWWBBWBBWWBWBBWBWBWWBWWBBWBBWBBBW
+WBWWWWWBBWBWWBBBBBWWBWWWBWBWBBBBWWBBWWWBWBBBBBBBWW
+WWWBWBBBWWBBBBWBBBWBBWBWWWWBBBWBWWWBWWBBWWBWWWWBWB
+BBBWWWWWWBBWWBWWBBWWBWBBWWWBWWWBWWBBBBBWBWBBBBBBBB
+BBBWBWWWBBBBWBWWBBBBWWBWBBBBWBWBWBBWWWWBBBBWBWBBBB
+WWBBWWWWBBWBBBBWBBBWBBWWWBWBWBBWWWWBBWBWWWBBBWBBWW
+BBBWBBBWWWWBWWBWBWBBBBWBBWWBWBBBWBBBWWWWBBBBWBBWWW
+WBBWBWWBWWBWBWBWWBBWBWWWWBWBBBBBBBWBWWWWWWBBBWBWBB
+BWBBWWBWBWWBBBBBBWWWBWWWBWWBBBBBBWBWWWWBWWBWWBBWWW
+WWWBWWWWWWWWBWWBWBWWBWBWWWWWWBWBBWBWWBWWBWBBBBWWBB
+BBWWBWBWBWWWWBBWWWBBBBBBBBWWBBBBWWBWWBBWWWWWWWBWWW
+BWBWWWWBBBWWWBWWWWBWWBBBBWWBWWBBBWWWBWBBBBBWWBBBWB
+BWWWWWWWBWWBBBBBBWWWWBWBWWWWWBWBWWBBBWBBWBWWBBWWWB
+WBBBWWBWWBWWBWBWBBBBBBBWWWBWBWBWWWWBWBBWBBBWWBWBWW
+WBBWWBWWWBBBBWBBWBBBBBWBWWBBBBBBWWWBWBBWBBBBWBBBWW
+WWBWWBBWWWWWWWBWBBWBBWBWWWWBWWBBWWWBWBBWWBWWWWWBWW
+BBBWBWWWWBWBBWWWBWWBWWBWWBWBWBBBBBWBWBWBWBBWBWBBWW
+WBWBWWWBBWBWBBBBWBBBWBWBWBBBBWBBWBWBWBBBWWWWBBBBWB
+WWWBWBWBWBWWWBBBWBBBBBWWBBBWWWWWWWWBWWBWWBBWBBWBWB
+BBWWWWWBBBWBBWWBBWWBBBWBWBBWBBWWWBBWWBWWWBWBWWWWBB
+WWBWBBBBBWBBBWBWWWWWWBWWBWBWWWWWWBBBWWWBWBBBWBBBWB
+BBBWWBWBWBBWWBBWBWBWBWBWWBWBWWWWBBWWBWBBWWBBWWWBWW
+BBWWBBWWWBBWBWBWBWWBWWBWBBWWBBWWWBBWWBWWWWBBWWBWWB
+WWBBBBBBWWBBBWBBWBBWWBBWWBWBWWBWBBBBWWWWWWBBBBBBBW
+WWBWWWWWBWBBWBWWBBBBWWWBWWWBWBBWBBBWWBBWWBBBWWBWBB
+WBBWBBBBBWBBWBBWWBBBBBBWBWBBBWWBWBWBWBWBWWWBBBBWBB
+BWWBBWBWBBWWWWWBBBWWBWWBBBBBBWWWBWWWBBWBBWWWWWBBWW
+WWBBBWBWBWWWWBBWWWBWWBWWWBWWBBBWBBWWBWWBBWBBWWWBBW
+WBWWBWBWBWWBBBBWBWWWBWWWBWBBWWWBWBBBWBBBBWBWWBWWWB
+WBBWBWBBBBBWWWBWBWBWBBWBBBBBWWWBBWBBWWWBWWBBWWWBBB
+BBWBWWBWBBWBWWWBBWWWBWBWWWWWWBBBWWBWBBBWWBWBWWBWWB
+WWBBBBWBWWWBBWWBWWWBWBBWWWWBBBWBBWWBWBBBWBWWWBBWWW
+WWWBWWWBWBWWBWBBBWWWWBBWBWBBWBBBWWBBWBWWWBBBBWWWBB
+WBWWWBWWWBBBBWBWBBBWWBBWBWWWBWWWBBWBWWWBWWBBBBWWWB
+BBWBBWBWWBWWWWBWBWWWWWBBWBBBBWWWWBWBBBBWWWBBWWBBBB
+BBWBWWBBWBBWBBBWBWBWBBWWWWBBBWWWBWWBBBBBWWBBWBWBBB
+WBWBWWWBBWBBBBBWWBWWWBWBBWWBWBWBBWWBBBWBBWWBBBWWBB
+BBWWBWWBBBWBBBWBWBBWBWBBBBBWWBWBBWBWWBBBWWWWBBBBBW
+BBBWWWBBWWBWWBWWWBBWBWWBWWBWBBWBBBWWWBWBWWBWWWBWWB
+BWWWBBBWWWBBWWWWBBBWBWWBWWBBWWBWBWWWWWBBBWWBBWWBBW
+WWBBWBBBWBWBWWWBWWWWBBBWBBBBBBWWBBBBWWBBBWBWBBWBWB
+WBBBBWWWWBBWWWWWBWBBBWBWWBWBWWBBBWWWWWBWWBBWWWWWBW
+BBWWBWBWBBBBWWBWWWBWWBWWWBWBWBWWWWWBBWBWBBBWWBBBWW
+BWBWWBBWWBWBWBWWBWBWBWBBWBWBWWBBWWWBWBWWWBWBWBBBBB
+WWBWWBBBBBWBWBWBWBBWBWWWWWWBBWBWBWWWBWWBWBWWWWWBBB
+BBWBWBWBWWBBBBBBBWBWBBBBWBBWWBBWWWBBBBBBBWWWWWBBBW
+BWBWBWBWBWBWBBWBBWWBWBBWBWWWWBWWBWBWWBWBBWBBWWBBBW
+WWBWWBBWBWBWBBBBBBWBWBBWWBWBBBWBWBBBBBBWWBBWBWBWBW
+WWBBBBBWBWBWWWBWBBWWBWBWWBWBWWWBWWWBBWBBBWBWBWBBBB
+BWBBBBBWWBWWWBWWWBBBWWBWBWBBBBBBBWBWBBWWWBWWBBWWBB

day4/8.cpp

@@ -0,0 +1,136 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+#define ALL(x) (x).begin(), (x).end()
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+#define mp make_pair
+
+using namespace std;
+using lli = long long int;
+
+using pii   = pair<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;
+
+const int BOUND = 1000*1000;
+//const int BOUND = 5;
+int nbChairs;
+
+int ask(int x, int y) {
+    printf("? %d %d\n", x, y);
+    fflush(stdout);
+    int out;
+    scanf("%d", &out);
+    return out;
+}
+
+struct Query {
+    unordered_map<int, int> memo;
+
+    int query(int v) {
+        if(memo.find(v) != memo.end())
+            return memo[v];
+        memo[v] = int_query(v);
+        return memo[v];
+    }
+    virtual int int_query(int v) = 0;
+};
+
+struct RowQuery: public Query {
+    int int_query(int v) {
+        return ask(BOUND + 42, v);
+    }
+};
+
+struct ColQuery: public Query {
+    int int_query(int v) {
+        return ask(v, BOUND + 42);
+    }
+};
+
+struct FindQuery: public Query {
+    vector<int> rows, cols;
+    vii prev;
+    vector<int> numInRow;
+    int cRow, foundThisRow;
+    FindQuery(): cRow(0), foundThisRow(0) {}
+
+    int int_query(int v) {
+        assert(v < (int)cols.size());
+        int ans = ask(cols[v], rows[cRow]);
+        for(const pii& p: prev)
+            if(p.second <= v && p.first < cRow)
+                ans--;
+        return ans;
+    }
+    void found(int col) {
+        prev.push_back(mp(cRow, col));
+        memo.clear();
+        foundThisRow++;
+        if(foundThisRow >= numInRow[cRow]) {
+            cRow++;
+            foundThisRow = 0;
+        }
+    }
+    void addRow(int v) {
+        if(rows.empty() || rows.back() != v) {
+            rows.push_back(v);
+            numInRow.push_back(1);
+        }
+        else
+            numInRow.back()++;
+        //printf("* Row %d (%d th)\n", v, numInRow.back());
+    }
+    void addCol(int v) {
+        if(cols.empty() || cols.back() != v)
+            cols.push_back(v);
+        //printf("* Col %d\n", v);
+    }
+};
+
+int dicho(Query* qu, int beg, int end, int seek) {
+    if(end - beg == 1)
+        return beg;
+    int m = (beg + end) / 2 - 1;
+    int mVal = qu->query(m);
+    if(mVal >= seek)
+        return dicho(qu, beg, m+1, seek);
+    return dicho(qu, m+1, end, seek);
+}
+
+int main(void) {
+    scanf("%d", &nbChairs);
+
+    FindQuery qFind;
+    RowQuery qRow;
+    ColQuery qCol;
+
+    FOR(i, nbChairs) {
+        qFind.addRow(dicho(&qRow, -BOUND, BOUND+1, i+1));
+        qFind.addCol(dicho(&qCol, -BOUND, BOUND+1, i+1));
+    }
+
+    FOR(i, nbChairs) {
+        int col = dicho(&qFind, 0, qFind.cols.size(), qFind.foundThisRow + 1);
+        qFind.found(col);
+    }
+    printf("!");
+    for(const pii& chair: qFind.prev)
+        printf(" %d %d", qFind.cols[chair.second], qFind.rows[chair.first]);
+    printf("\n");
+    fflush(stdout);
+    return 0;
+}

day4/8_int.py

@@ -0,0 +1,67 @@
+from random import randint
+from sys import stdin, stderr, stdout
+BOUND = 10**6
+#BOUND = 5
+
+pts = []
+
+while len(pts) < 100:
+    x = randint(-BOUND, BOUND)
+    y = randint(-BOUND, BOUND)
+    if (x, y) in pts:
+        continue
+    pts.append((x, y))
+
+print(len(pts))
+stdout.flush()
+
+nbQuestions = 0
+
+
+def question(l):
+    global nbQuestions
+    x, y = [int(x) for x in l.strip().split()]
+    nbQuestions += 1
+    out = 0
+    for (xi, yi) in pts:
+        if xi <= x and yi <= y:
+            out += 1
+    return out
+
+
+def check(l):
+    global nbQuestions
+    gPts = set()
+    print("{} questions, {} allowed".format(
+        nbQuestions,
+        50 * len(pts)), file=stderr)
+    for pt in l:
+        gPts.add(pt)
+    if len(gPts) != len(pts):
+        print('Fail', file=stderr)
+        return
+    for pt in gPts:
+        if pt not in pts:
+            print('Fail (not in)', file=stderr)
+            return
+    print('Success', file=stderr)
+
+
+while True:
+    line = input()
+    if line[0] == '!':
+        spl = line[2:].strip().split()
+        l = []
+        for i in range(0, len(spl)//2):
+            l.append((int(spl[2*i]), int(spl[2*i+1])))
+        check(l)
+        print(l, file=stderr)
+        print(pts, file=stderr)
+        break
+    elif line[0] == '?':
+        print(question(line[2:]))
+        stdout.flush()
+    elif line[0] == '*':
+        print(line[1:], file=stderr)
+    else:
+        print("OMG", file=stderr)

day4/run.sh

@@ -0,0 +1,3 @@
+mkfifo plug.fifo
+./a.out < plug.fifo | python 8_int.py > plug.fifo
+rm plug.fifo

misc/voronoi_twal.py

@@ -0,0 +1,380 @@
+import heapq
+import sys
+from math import sqrt
+from functools import reduce
+from random import randint
+
+# J'ai utilisé dans cet exo l'algorithme de Fortune pour calculer le diagramme de Voronoï des points d'accès wifi
+# J'ai utilisé une skiplist à la place d'un arbre binaire balancé car j'ai un peu la flemme de coder un rb-tree ^^
+# Avant d'utiliser une skiplist, j'avais une liste doublement chaînée utilisant la même interface, car c'est plus facile à
+# coder et du coup moins de chance d'avoir des bugs. Avec j'avais une complexitée de O(n*(n+h)) où h est le nombre de points
+# sur l'enveloppe convexe. Avec une skiplist, j'ai une complexité de O(n*(log(n)+h)), mais j'ai été déçu car la constante
+# de temps fait que la skiplist devient avantageuse autour de N=20000, et l'énoncé donne N<=2000...
+# On m'a néanmoins conseillé de donner l'algorithme avec la meilleure complexité, donc la voici.
+
+# Une priority queue car elles sont nécessairement synchronisées dans python... grr.
+class PriorityQueue:
+    def __init__(self, elems):
+        self.queue = [[prio, i, elem] for (i, (elem, prio)) in enumerate(elems)]
+        self.entryFinder = {}
+        for entry in self.queue:
+            self.entryFinder[entry[2]] = entry
+        heapq.heapify(self.queue)
+        self.counter = len(self.queue)
+
+    def push(self, elem, priority):
+        entry = [priority, self.counter, elem]
+        self.entryFinder[elem] = entry
+        heapq.heappush(self.queue, entry)
+        self.counter += 1
+
+    def top(self):
+        while self.queue:
+            candidate = self.queue[0][2]
+            if candidate != None:
+                return candidate
+            heapq.heappop(self.queue)
+        return None
+
+    def pop(self):
+        while self.queue:
+            candidate = heapq.heappop(self.queue)[2]
+            if candidate != None:
+                del self.entryFinder[candidate]
+                return candidate
+        return None
+
+    def contains(self, elem):
+        return elem in self.entryFinder
+
+    def remove(self, elem):
+        entry = self.entryFinder.pop(elem)
+        entry[2] = None
+
+    def empty(self):
+        return self.top() == None
+
+# Une classe représentant l'enveloppe convexe d'un ensemble de points
+class ConvexHull:
+    def __init__(self, sites):
+        # On la calcule avec un parcours de graham
+        sortedSites = sorted(sites, key=lambda s:s.x)
+        up = reduce(ConvexHull.keepLeft, sortedSites, [])
+        down = reduce(ConvexHull.keepLeft, reversed(sortedSites), [])
+        self.hull = up + down[1:-1]
+
+    def det(a, b, n):
+        return (b.x-a.x)*(n.y-a.y) - (b.y-a.y)*(n.x-a.x)
+
+    def keepLeft(hull, n):
+        while len(hull) > 1 and ConvexHull.det(hull[-2], hull[-1], n) <= 0:
+            hull.pop()
+        hull.append(n)
+        return hull
+
+    def isLeftTo(a, b, s):
+        return (s.x-a.x)*(b.y-a.y)-(s.y-a.y)*(b.x-a.x) < 0
+
+    def isInsideHull(self, s):
+        for i in range(len(self.hull)):
+            if not ConvexHull.isLeftTo(self.hull[i-1], self.hull[i], s):
+                return False
+        return True
+
+# Classe représentant un point dans le plan. C'est nommé "Site" car c'est comme ça que sont nommés les points
+# dont il faut calculer le diagramme de voronoi dans les papiers de Fortune
+class Site:
+    def __init__(self, x, y):
+        self.x = x
+        self.y = y
+
+    def dist(self, other):
+        return sqrt((self.x-other.x)**2 + (self.y-other.y)**2)
+
+# Une parabole de la beachline
+class Parabola:
+    def __init__(self, ctx, site):
+        self.ctx = ctx
+        self.site = site
+        self.ledge = None
+        self.redge = None
+
+    def getEquation(self):
+        x = self.site.x
+        y = self.site.y
+        l = self.ctx.sweepline
+        div = 2*(y-l)
+        if div == 0: #FIXME
+            return None
+        # Retourne (a, b, c) tels que y = a*x^2 + b*x + c
+        # Obtenus en résolvant en y sqrt((x-xa)^2 + (y-ya)^2) = y-l avec xa = self.site.x, ya = self.site.y
+        return (1/div, -2*x/div, (x*x + y*y - l*l)/div)
+
+    def intersect(self, other):
+        if self.getEquation() == None:
+            return self.site.x
+        if other.getEquation() == None:
+            return other.site.x
+        (a, b, c) = tuple(map(lambda x, y: x-y, self.getEquation(), other.getEquation()))
+        if a == 0: #FIXME
+            # b*x + c = 0
+            return -c/b
+        else:
+            # a*x^2 + b*x + c = 0
+            discr = sqrt(b*b-4*a*c)
+            x1 = (-b-discr)/(2*a)
+            x2 = (-b+discr)/(2*a)
+            if self.site.y < other.site.y:
+                return min(x1, x2)
+            else:
+                return max(x1, x2)
+
+    # Donne le point où la parabole va être "écrasée" par ses deux voisines
+    def circleCenter(self):
+        return self.ledge.intersect(self.redge)
+
+    def hasCircleCenter(self):
+        return self.ledge != None and self.redge != None and self.circleCenter() != None
+
+    # Renvoie la position de la sweepline quand la parabole va disparaitre
+    def circleEventY(self):
+        voronoiVertex = self.circleCenter()
+        return voronoiVertex.y + voronoiVertex.dist(self.site)
+
+    # Renvoie une copie, sans ledge et redge car elles sont recalculées après une copie
+    def copy(self):
+        return Parabola(self.ctx, self.site)
+
+# Représente une arête du diagramme de voronoï
+class Edge:
+    def __init__(self, leftSite, rightSite):
+        self.ls = leftSite
+        self.rs = rightSite
+        # a*x + b*y = c
+        # Obtenus en résolvant en x et y sqrt((x-ls.x)^2 + (y-ls.y)^2) = sqrt((x-lr.x)^2 + (y-lr.y)^2)
+        self.a = 2*(rightSite.x - leftSite.x)
+        self.b = 2*(rightSite.y - leftSite.y)
+        self.c = rightSite.x**2 + rightSite.y**2 - (leftSite.x**2 + leftSite.y**2)
+
+    def intersect(self, other):
+        # Résolution via la méthode de Cramer
+        det = self.a*other.b - self.b*other.a
+        if det == 0: #FIXME
+            return None
+        return Site((self.c*other.b - self.b*other.c)/det, (self.a*other.c - self.c*other.a)/det)
+
+# Représente un noeud d'une skip-list doublement chaînée
+class SkipListNode:
+    def __init__(self, height, parabola):
+        self.prv = [None]*height
+        self.nxt = [None]*height
+        self.parabola = parabola
+
+# Représente une skiplist
+class SkipList:
+    def __init__(self, elems):
+        self.maxHeight = SkipList.randomHeight()
+        self.head = SkipListNode(self.maxHeight, elems[0])
+        updateList = [self.head]*self.maxHeight
+        for i in range(1, len(elems)):
+            self.insertAfter(elems[i], updateList)
+
+    # Renvoie une hauteur aléatoire
+    def randomHeight():
+        height = 1
+        while randint(1, 2) != 1:
+            height += 1
+        return height
+
+    # Insertion d'un élément en fournissant le résultat de getElemByX() dans updateList
+    # Ça met aussi à jour updateList pour insérer à nouveau un élément après celui qu'on vient d'insérer
+    def insertAfter(self, elem, updateList):
+        height = SkipList.randomHeight()
+        self.maxHeight = max(height, self.maxHeight)
+        newNode = SkipListNode(height, elem)
+        # On ajuste la hauteur de updateList et self.head.nxt
+        while len(self.head.nxt) < self.maxHeight:
+            self.head.nxt.append(None)
+        while len(updateList) < self.maxHeight:
+            updateList.append(self.head)
+
+        for i in range(height):
+            # Là c'est comme dans une liste doublement chaînée sauf qu'il y a des [i] partout
+            newNode.prv[i] = updateList[i]
+            newNode.nxt[i] = updateList[i].nxt[i]
+            if updateList[i].nxt[i] != None:
+                updateList[i].nxt[i].prv[i] = newNode
+            updateList[i].nxt[i] = newNode
+            updateList[i] = newNode
+        return newNode
+
+    # Enlève un élément
+    def remove(self, node):
+        for i in range(len(node.nxt)-1, -1, -1):
+            # Là c'est aussi comme pour une liste doublement chaînée mais avec des [i] partout
+            node.prv[i].nxt[i] = node.nxt[i]
+            if node.nxt[i] != None:
+                node.nxt[i].prv[i] = node.prv[i]
+            if self.head.nxt[i] == None:
+                self.maxHeight -= 1
+
+    # On donne un X, et ça dit quelle est la parabole présente à cet endroit
+    # Ça renvoie les éléments à modifier si on insère une autre parabole après celle cherchée
+    def getElemByX(self, x):
+        updateList = [None]*self.maxHeight
+        node = self.head
+        for i in range(self.maxHeight-1, -1, -1):
+            while node.nxt[i] != None and self.getX(node) < x:
+                node = node.nxt[i]
+            # Souvent on "overshoot" la parabole et il faut revenir un peu en arrière
+            if node.prv[0] != None and self.getX(node.prv[0]) > x:
+                node = node.prv[i]
+            updateList[i] = node
+        return updateList
+
+    # Renvoie le x où node.parabola se termine
+    def getX(self, node):
+        return node.parabola.intersect(node.nxt[0].parabola)
+
+# Représente un évènement, circle event ou site event
+class Event:
+    def __init__(self, isSiteEvent, elem):
+        self.isSiteEvent = isSiteEvent
+        self.elem = elem
+
+    # Pour pouvoir enlever des éléments de la hashtable de la priority queue
+    def __hash__(self):
+        return hash(self.isSiteEvent) ^ hash(self.elem)
+
+    def __eq__(self, other):
+        return other != None and self.isSiteEvent == other.isSiteEvent and self.elem == other.elem
+
+class Context:
+    def __init__(self, sites):
+        self.sweepline = 0
+        self.sites = [Site(x, y) for (x, y) in sites]
+        minY = min([site.y for site in self.sites])
+        # Les sites qui seront dans la beachline
+        sitesBL = [site for site in self.sites if site.y == minY] #FIXME
+        # Les sites qui seront dans la priority queue
+        sitesQ = [site for site in self.sites if site.y != minY] #FIXME
+        self.queue = PriorityQueue([(Event(True, site), site.y) for site in sitesQ])
+        self.beachline = SkipList(list(map(lambda site: Parabola(self, site), sorted(sitesBL, key=lambda site: site.x))))
+        # On calcule les arêtes associées aux paraboles de la beachline
+        node = self.beachline.head
+        while node.nxt[0] != None:
+            self.updateREdge(node)
+            node = node.nxt[0]
+        self.sweepline = minY
+        self.maxDist = 0
+        self.hull = ConvexHull(self.sites)
+
+    # Enlève un circle event de la queue
+    def removeCircleEvent(self, node):
+        if self.queue.contains(Event(False, node)):
+            self.queue.remove(Event(False, node))
+
+    # Ajoute un circle event dans la queue après quelques vérifications
+    def addCircleEvent(self, node):
+        def isRight(a, b, s):
+            return (b.x-a.x)*(s.y-a.y)-(b.y-a.y)*(s.x-a.x) >= 0
+
+        parabola = node.parabola
+        if parabola.hasCircleCenter():
+            pos = parabola.circleEventY()
+            center = parabola.circleCenter()
+            # On ajoute l'event si il est derrière la sweepline et si la parabole va en effet disparaitre
+            # (elle disparait si les points de la parabole de gauche, du milieu, de droite tourne dans le sens horaire, ce que calcule isRight)
+            if pos > self.sweepline and isRight(node.prv[0].parabola.site, node.parabola.site, node.nxt[0].parabola.site):
+                self.queue.push(Event(False, node), pos)
+
+    # Calcule l'arête à droite d'une parabole
+    def updateREdge(self, node):
+        if node.nxt[0] != None:
+            edge = Edge(node.parabola.site, node.nxt[0].parabola.site)
+            node.parabola.redge = edge
+            node.nxt[0].parabola.ledge = edge
+
+    # On a trouvé un point du diagramme de voronoi ! C'est là qu'il risque d'y avoir des jaloux
+    # Si le point du diagramme est en dehors de l'enveloppe convexe de tous les points alors
+    # Il n'y aura pas de jaloux là, il faut donc ne pas mettre à jour self.maxDist
+    def handleResult(self, v, p1, p2, p3):
+        d = v.dist(p1)
+        if d > self.maxDist and self.hull.isInsideHull(v):
+            self.maxDist = d
+
+    # La partie intéressante de l'algorithme, le traitement d'un event
+    def processEvent(self):
+        event = self.queue.pop()
+        # Si c'est un site event
+        if event.isSiteEvent:
+            # On récupère quelques infos
+            site = event.elem
+            self.sweepline = site.y
+            updateList = self.beachline.getElemByX(site.x)
+            leftNode = updateList[0]
+            xpos = self.beachline.getX(leftNode) if leftNode.nxt[0] != None else float("nan")
+            # On ajoute notre parabole
+            middleNode = self.beachline.insertAfter(Parabola(self, site), updateList)
+            # On enlève un éventuel circle event de la parabole qu'on "coupe"
+            self.removeCircleEvent(leftNode)
+            # Si on coupe une parabole
+            if site.x != xpos: #FIXME
+                # On ajoute la parabole qu'on coupe car pour l'instant il n'y a que la partie de gauche
+                rightNode = self.beachline.insertAfter(leftNode.parabola.copy(), updateList)
+                self.updateREdge(rightNode)
+            # Si on s'intercale entre deux paraboles
+            else:
+                # On a trouvé un point du diagramme de voronoï, youhou !
+                self.handleResult(Edge(site, leftNode.parabola.site).intersect(Edge(site, middleNode.nxt[0].parabola.site)), site, leftNode.parabola.site, middleNode.nxt[0].parabola.site)
+
+            # On re-calcule les arêtes
+            self.updateREdge(leftNode)
+            self.updateREdge(middleNode)
+
+            # On ajoute les éventuels circle events
+            self.addCircleEvent(leftNode)
+            if site.x != xpos: #FIXME
+                self.addCircleEvent(rightNode)
+
+        # Si c'est un circle event
+        else:
+            # On récupère quelques infos
+            middleNode = event.elem
+            leftNode = middleNode.prv[0]
+            rightNode = middleNode.nxt[0]
+            self.sweepline = middleNode.parabola.circleEventY()
+
+            # On a trouvé un point du diagramme de voronoï, youhou !
+            self.handleResult(middleNode.parabola.circleCenter(), leftNode.parabola.site, middleNode.parabola.site, rightNode.parabola.site)
+
+            # On enlève les éventuels circle events
+            self.removeCircleEvent(leftNode)
+            self.removeCircleEvent(rightNode)
+            # On enlève la parabole qui disparait
+            self.beachline.remove(middleNode)
+            # On ajoute les nouveaux circle events
+            self.updateREdge(leftNode)
+            self.addCircleEvent(leftNode)
+            self.addCircleEvent(rightNode)
+
+    def process(self):
+        # Tant qu'il y a du boulot, on le traite !
+        while not self.queue.empty():
+            self.processEvent()
+
+def wifi(N, coords):
+    # La fonction "wifi" est maintenant très compliquée à coder
+    ctx = Context(coords)
+    ctx.process()
+    return ctx.maxDist
+
+if __name__ == '__main__':
+    N = int(input())
+    coords = [tuple(map(int, input().split())) for _ in range(N)]
+    f = wifi(N, coords)
+    if f == 0: #FIXME
+        print("0")
+    else:
+        print("%.3f" % f)
+

template.cpp

@@ -0,0 +1,25 @@
+#include <bits/stdc++.h>
+
+#define FOR(i, n) for(lli i = 0; i < (lli)(n); ++i)
+#define ALL(x) (x).begin(), (x).end()
+
+#define X(A) get<0>(A)
+#define Y(A) get<1>(A)
+#define Z(A) get<2>(A)
+#define W(A) get<3>(A)
+
+#define mt make_tuple
+
+using namespace std;
+using lli = long long int;
+
+using pii   = tuple<lli, lli>;
+using piii  = tuple<lli, lli, lli>;
+using vi    = vector<lli>;
+using vii   = vector<pii>;
+using viii  = vector<piii>;
+using vvi   = vector<vi>;
+using vvii  = vector<vii>;
+using vviii = vector<viii>;
+using vb    = vector<bool>;
+using vvb   = vector<vb>;