TrafficSimulation

a.cpp

#if 1
#include <array>
#include <bitset>
#include <deque>
#include <map>
#include <queue>
#include <set>
#include <unordered_set>
#include <unordered_map>
#include <vector>
#include <algorithm>
#include <functional>
#include <string>
#include <fstream>
#include <iostream>
#include <sstream>
#include <random>
#include <cmath>
#include <cstring>
using namespace std;
using i8 = int8_t; using i16 = int16_t; using i32 = int32_t; using i64 = int64_t; using u8 = uint8_t; using u16 = uint16_t; using u32 = uint32_t; using u64 = uint64_t; using f32 = float; using f64 = double;
template <class T> inline T sq(T a) { return a * a; }
template <class T> inline bool chmin(T& a, const T& b) { if (b < a) { a = b; return true; } return false; }
template <class T> inline bool chmax(T& a, const T& b) { if (a < b) { a = b; return true; } return false; }
void _u() { cerr << endl; } template <class H, class... T> void _u(H&& h, T&&... t) { cerr << h << ", "; _u(move(t)...); }
#define _polyfill(v) v
#define _overload3(a, b, c, d, ...) d
#define _rep3(i, a, b) for (i32 i = i32(a); i < i32(b); ++i)
#define _rep2(i, b) _rep3(i, 0, b)
#define _rep1(i) for (i32 i = 0; ; ++i)
#define rep(...) _polyfill(_overload3(__VA_ARGS__, _rep3, _rep2, _rep1)(__VA_ARGS__))
#define A(a) begin(a), end(a)
#define S(a) i32((a).size())
#define U(...) { cerr << #__VA_ARGS__ << ": "; _u(__VA_ARGS__); }
#define lin U(__LINE__)
#define exi exit(0);
#endif

class timer {
    vector<timer> timers;
    int n = 0;
public:
    double limit = 9.8;
    double t = 0;
    timer() {}
    timer(int size) : timers(size) {}
    bool elapses() const {
        return time() - t > limit;
    }
    void measure() {
        t = time() - t;
        ++n;
    }
    void measure(char id) {
        timers[id].measure();
    }
    void print() {
        if (n % 2)
            measure();
        for (int i = 0; i < 128; ++i) {
            if (timers[i].n)
                cerr << (char)i << ' ' << timers[i].t << 's' << endl;
        }
        cerr << "  " << t << 's' << endl;
    }
    static double time() {
#ifdef LOCAL
        return __rdtsc() / 3.3e9;
#else
        unsigned long long a, d;
        __asm__ volatile ("rdtsc" : "=a" (a), "=d" (d));
        return (d << 32 | a) / 2.8e9;
#endif
    }
} timer(128);

constexpr bool Deterministic = 0;

class rando {
    unsigned y;
public:
    rando(unsigned y) : y(y) {}
    unsigned next() {
        return y ^= (y ^= (y ^= y << 13) >> 17) << 5;
    }
    int next(int b) {
        return next() % b;
    }
    int next(int a, int b) {
        return next(b - a) + a;
    }
    double nextDouble(double b = 1) {
        return b * next() / 4294967296.0;
    }
    double nextDouble(double a, double b) {
        return nextDouble(b - a) + a;
    }
    int operator() (int b) {
        return next(b);
    }
} rando(Deterministic ? 2463534242 : random_device()());

namespace Tercel
{
  struct Vector {
    double x;
    double y;

    bool operator==(const Vector& v) const
    {
      return (x == v.x && y == v.y);
    }

    bool operator<(const Vector& v) const
    {
      return x != v.x ? x < v.x : y < v.y;
    }
  };

  struct Circle
  {
    Vector center;
    double radius;
  };

  class Triangle {
  public:
    const Vector* p1, * p2, * p3;

  private:
    inline const Vector* getMinVertex() const
    {
      return *p1 < *p2 && *p1 < *p3 ? p1 : *p2 < *p3 ? p2 : p3;
    }

    inline const Vector* getMidVertex() const
    {
      return *p1 < *p2 && *p2 < *p3 || *p3 < *p2 && *p2 < *p1 ? p2 :
             *p2 < *p3 && *p3 < *p1 || *p1 < *p3 && *p3 < *p2 ? p3 : p1;
    }

    inline const Vector* getMaxVertex() const
    {
      return *p2 < *p1 && *p3 < *p1 ? p1 : *p3 < *p2 ? p2 : p3;
    }

  public:
    bool operator==(const Triangle& t) const
    {
      return *p1 == *t.p1 && *p2 == *t.p2 && *p3 == *t.p3 ||
             *p1 == *t.p2 && *p2 == *t.p3 && *p3 == *t.p1 ||
             *p1 == *t.p3 && *p2 == *t.p1 && *p3 == *t.p2 ||
             *p1 == *t.p3 && *p2 == *t.p2 && *p3 == *t.p1 ||
             *p1 == *t.p2 && *p2 == *t.p1 && *p3 == *t.p3 ||
             *p1 == *t.p1 && *p2 == *t.p3 && *p3 == *t.p2;
    }

    bool operator<(const Triangle& t) const
    {
      return !(*getMinVertex() == *t.getMinVertex()) ?
                 *getMinVertex() < *t.getMinVertex() :
             !(*getMidVertex() == *t.getMidVertex()) ?
                 *getMidVertex() < *t.getMidVertex() :
                 *getMaxVertex() < *t.getMaxVertex();
    }

    bool hasCommonPoints(const Triangle& t) const
    {
      return *p1 == *t.p1 || *p1 == *t.p2 || *p1 == *t.p3 ||
             *p2 == *t.p1 || *p2 == *t.p2 || *p2 == *t.p3 ||
             *p3 == *t.p1 || *p3 == *t.p2 || *p3 == *t.p3;
    }
  };

  class Delaunay2d
  {
    typedef const std::set<Vector>         ConstVertexSet;
    typedef ConstVertexSet::const_iterator ConstVertexIter;
    typedef std::set<Triangle>             TriangleSet;
    typedef std::set<Triangle>::iterator   TriangleIter;
    typedef std::map<Triangle, bool>       TriangleMap;

  private:
    static inline void addElementToRedundanciesMap(TriangleMap* pRddcMap,
                                                   Triangle& t)
    {
      TriangleMap::iterator it = pRddcMap->find(t);
      if(it != pRddcMap->end() && it->second)
      {
        pRddcMap->erase(it);
        pRddcMap->insert(TriangleMap::value_type(t, false));
      }
      else
      {
        pRddcMap->insert(TriangleMap::value_type(t, true));
      }
    }

  public:
    static void getDelaunayTriangles(ConstVertexSet& pVertexSet,
                                     TriangleSet* triangleSet)
    {
      Triangle hugeTriangle;
      {
        double maxX, maxY; maxX = maxY = -DBL_MAX;
        double minX, minY; minX = minY = DBL_MAX;
        for(ConstVertexIter it = pVertexSet.begin();
            it != pVertexSet.end(); ++it)
        {
          double x = it->x;
          double y = it->y;
          if(maxX < x) maxX = x; if(minX > x) minX = x;
          if(maxY < y) maxY = y; if(minY > y) minY = y;
        }
        double centerX  = (maxX - minX) * 0.5;
        double centerY  = (maxY - minY) * 0.5;
        double dX   = maxX - centerX;
        double dY   = maxY - centerY;
        double radius = sqrt(dX * dX + dY * dY) + 1.0;
        Vector* p1 = new Vector;
        p1->x    = centerX - sqrt(3.0) * radius;
        p1->y    = centerY - radius;
        Vector* p2 = new Vector;
        p2->x    = centerX + sqrt(3.0) * radius;
        p2->y    = centerY - radius;
        Vector* p3 = new Vector;
        p3->x    = centerX;
        p3->y    = centerY + 2.0 * radius;
        hugeTriangle.p1 = p1;
        hugeTriangle.p2 = p2;
        hugeTriangle.p3 = p3;
      }
      triangleSet->insert(hugeTriangle);
      for(ConstVertexIter vIter = pVertexSet.begin();
          vIter != pVertexSet.end(); ++vIter)
      {
        const Vector* p = &*vIter;
        TriangleMap rddcMap;
        for(TriangleIter tIter = triangleSet->begin();
            tIter != triangleSet->end();)
        {
          Triangle t = *tIter;
          Circle   c;
          {
            double x1 = t.p1->x;  double y1 = t.p1->y;
            double x2 = t.p2->x;  double y2 = t.p2->y;
            double x3 = t.p3->x;  double y3 = t.p3->y;
            double m = 2.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
            double x = ((y3-y1)*(x2*x2-x1*x1+y2*y2-y1*y1)
                       +(y1-y2)*(x3*x3-x1*x1+y3*y3-y1*y1))/m;
            double y = ((x1-x3)*(x2*x2-x1*x1+y2*y2-y1*y1)
                       +(x2-x1)*(x3*x3-x1*x1+y3*y3-y1*y1))/m;
            c.center.x = x;
            c.center.y = y;
            double dx   = t.p1->x - x;
            double dy   = t.p1->y - y;
            double radius = sqrt(dx * dx + dy * dy);
            c.radius = radius;
          }
          double dx = c.center.x - p->x;
          double dy = c.center.y - p->y;
          double dist = sqrt(dx * dx + dy * dy);
          if(dist < c.radius)
          {
            Triangle t1;
            t1.p1 = p; t1.p2 = t.p1; t1.p3 = t.p2;
            addElementToRedundanciesMap(&rddcMap, t1);
            Triangle t2;
            t2.p1 = p; t2.p2 = t.p2; t2.p3 = t.p3;
            addElementToRedundanciesMap(&rddcMap, t2);
            Triangle t3;
            t3.p1 = p; t3.p2 = t.p3; t3.p3 = t.p1;
            addElementToRedundanciesMap(&rddcMap, t3);
            triangleSet->erase(tIter++);
          }
          else ++tIter;
        }
        for(TriangleMap::iterator iter = rddcMap.begin();
          iter != rddcMap.end(); ++iter)
        {
          if(iter->second) triangleSet->insert(iter->first);
        }
      }
      for(TriangleIter tIter = triangleSet->begin();
          tIter != triangleSet->end(); )
      {
        if(hugeTriangle.hasCommonPoints(*tIter))
          triangleSet->erase(tIter++);
        else ++tIter;
      }
      delete hugeTriangle.p1;
      delete hugeTriangle.p2;
      delete hugeTriangle.p3;
    }
  };
}

struct P {
    i32 x, y;
};

struct HouseToHouse {
    i32 from, to, size;
    f64 distance;
    vector<i32> path;
    i32 rails;
    f64 time;
};

struct Node {
    i32 v;
    f64 cost, dist;
    f64 score;
    bool operator < (const Node& o) const {
        return score > o.score;
    }
};

f64 RoadCost, RailCost;
vector<P> ps, ps2;
vector<i32> to;
i32 i2r[1024];
vector<array<i32, 2>> es;
vector<i32> esPerP[1024];
f64 distances[1024][1024];
f64 cost, dist;
f64 curr, best;
vector<HouseToHouse> h2hs, bh2hs;
vector<i32> kk[1024][1024];
vector<pair<i32, i32>> roads2;
vector<i8> isRails;
vector<f64> memo;
vector<i32> pre;
i32 ct[1024][1024];
vector<pair<f64, i32>> heap;
vector<Node> q;
vector<array<i32, 3>> roads;

inline pair<f64, f64> flip(i32 i) {
    f64 costdiff = 0;
    f64 distancediff = 0;
    if (isRails[i]) {
        costdiff = (RoadCost - RailCost) * distances[roads2[i].first][roads2[i].second];
        for (i32 j : kk[roads2[i].first][roads2[i].second]) {
            auto& h2h = h2hs[j];
            distancediff -= sq(h2h.time / (h2h.rails ? 1 : 2)) * h2h.size;
            h2h.time -= distances[roads2[i].first][roads2[i].second] / 2;
            --h2h.rails;
            h2h.time += distances[roads2[i].first][roads2[i].second] * sq(1 + 0.02 * ct[roads2[i].first][roads2[i].second]);
            distancediff += sq(h2h.time / (h2h.rails ? 1 : 2)) * h2h.size;
        }
    } else {
        costdiff = (RailCost - RoadCost) * distances[roads2[i].first][roads2[i].second];
        for (i32 j : kk[roads2[i].first][roads2[i].second]) {
            auto& h2h = h2hs[j];
            distancediff -= sq(h2h.time / (h2h.rails ? 1 : 2)) * h2h.size;
            h2h.time -= distances[roads2[i].first][roads2[i].second] * sq(1 + 0.02 * ct[roads2[i].first][roads2[i].second]);
            ++h2h.rails;
            h2h.time += distances[roads2[i].first][roads2[i].second] / 2;
            distancediff += sq(h2h.time / (h2h.rails ? 1 : 2)) * h2h.size;
        }
    }
    isRails[i] = !isRails[i];
    return { costdiff, distancediff };
}

inline void hc() {
    roads2.clear();
    rep (i, S(h2hs)) {
        auto& h2h = h2hs[i];
        rep(j, 1, S(h2h.path)) {
            roads2.emplace_back(min(h2h.path[j - 1], h2h.path[j]), max(h2h.path[j - 1], h2h.path[j]));
            kk[roads2.back().first][roads2.back().second].emplace_back(i);
        }
    }
    sort(A(roads2));
    roads2.erase(unique(A(roads2)), roads2.end());
    isRails.assign(S(roads2), 1);
    f64 cost = RailCost * ::cost;
    f64 dist = ::dist;
    rep (i, S(h2hs)) {
        auto& h2h = h2hs[i];
        h2h.time = h2h.distance / 2;
        h2h.rails = S(h2h.path) - 1;
    }
    curr = cost * sqrt(dist);
    while (true) {
        bool updated = false;
        rep (i, S(roads2)) {
            auto p = flip(i);
            f64 next = (cost + p.first) * sqrt(dist + p.second);
            if (curr > next) {
                cost += p.first;
                dist += p.second;
                curr = next;
                updated = true;
            } else flip(i);
        }
        if (!updated) break;
    }
    if (best > curr) {
        best = curr;
        bh2hs = h2hs;
        roads.resize(S(roads2));
        rep (i, S(roads)) {
            roads[i][0] = i2r[roads2[i].first];
            roads[i][1] = i2r[roads2[i].second];
            roads[i][2] = isRails[i];
        }
    }
    for (auto& p : roads2) {
        kk[p.first][p.second].clear();
    }
}

i32 main() {
    timer.measure();
    cin >> RoadCost >> RailCost;
    {
        vector<vector<i32>> a(100, vector<i32>(100, -1));
        vector<vector<i32>> b(1000, vector<i32>(1000, 0));
        vector<i32> i2i(1000);
        to.assign(1000, 0);
        rep (i, 1000) {
            f64 fx, fy;
            cin >> to[i] >> fx >> fy;
            i32 x = i32(fx);
            i32 y = i32(fy);
            P p;
            p.x = x;
            p.y = y;
            if (a[y][x] == -1) {
                a[y][x] = S(ps);
                i2r[a[y][x]] = S(ps2);
                ps.emplace_back(p);
            }
            ps2.emplace_back(p);
            i2i[i] = a[y][x];
        }
        rep (i, 1000) {
            i32 i1 = min(i2i[i], i2i[to[i]]);
            i32 i2 = max(i2i[i], i2i[to[i]]);
            if (i1 != i2) {
                if (b[i1][i2] == 0) {
                    HouseToHouse h2h;
                    h2h.from = i1;
                    h2h.to = i2;
                    h2hs.emplace_back(h2h);
                }
                ++b[i1][i2];
            }
        }
        for (auto& h2h : h2hs) {
            h2h.size = b[h2h.from][h2h.to];
        }
    }
    {
        vector<i32> p2i(10000);
        set<Tercel::Vector> vs;
        set<Tercel::Triangle> ts;
        rep (i, S(ps)) {
            auto& p = ps[i];
            p2i[p.y * 100 + p.x] = i;
            Tercel::Vector v;
            v.x = p.x;
            v.y = p.y;
            vs.emplace(v);
        }
        Tercel::Delaunay2d::getDelaunayTriangles(vs, &ts);
        for (auto& t : ts) {
            i32 i1 = p2i[i32(t.p1->y) * 100 + i32(t.p1->x)];
            i32 i2 = p2i[i32(t.p2->y) * 100 + i32(t.p2->x)];
            i32 i3 = p2i[i32(t.p3->y) * 100 + i32(t.p3->x)];
            es.emplace_back(array<i32, 2>{ min(i1, i2), max(i1, i2) });
            es.emplace_back(array<i32, 2>{ min(i1, i3), max(i1, i3) });
            es.emplace_back(array<i32, 2>{ min(i2, i3), max(i2, i3) });
        }
    }
    sort(A(es));
    es.erase(unique(A(es)), es.end());
    for (auto& e : es) {
        esPerP[e[0]].emplace_back(e[1]);
        esPerP[e[1]].emplace_back(e[0]);
    }
    rep (i, 1, S(ps)) {
        rep (j, i) {
            distances[i][j] = distances[j][i] = sqrt(sq(ps[i].x - ps[j].x) + sq(ps[i].y - ps[j].y));
        }
    }
    best = 1e9;
    while (true) {
        random_shuffle(A(h2hs), rando);
        cost = dist = 0;
        memset(ct, 0, sizeof(ct));
        for (auto& h2h : h2hs) {
            h2h.distance = 0;
            h2h.path.clear();
        }
        i32 la = -1;
        rep (iter) {
            f64 pt[32];
            rep (i, 32) {
                pt[i] = max(0.08, pow(iter == 0 ? 0.92 : iter == 1 ? 0.79 : 0.73, i));
            }
            rep (h, S(h2hs)) {
                if (la == h) {
                    la = -1;
                    break;
                }
                auto& h2h = h2hs[h];
                curr = cost * sqrt(dist);
                dist -= sq(h2h.distance / 2) * h2h.size;
                rep (i, 1, S(h2h.path)) {
                    ct[h2h.path[i - 1]][h2h.path[i]] -= h2h.size;
                    ct[h2h.path[i]][h2h.path[i - 1]] -= h2h.size;
                    if (ct[h2h.path[i - 1]][h2h.path[i]] == 0) {
                        cost -= distances[h2h.path[i - 1]][h2h.path[i]];
                    }
                }
                memo.assign(1000, 1e9);
                pre.assign(1000, -1);
                heap.clear();
                heap.emplace_back(distances[h2h.from][h2h.to] * 0.08, h2h.from);
                memo[h2h.from] = distances[h2h.from][h2h.to] * 0.08;
                while (heap.size()) {
                    i32 c = heap[0].second;
                    f64 d = heap[0].first;
                    pop_heap(A(heap), greater<pair<f64, i32>>());
                    heap.pop_back();
                    if (c == h2h.to) break;
                    if (memo[c] != d) continue;
                    d -= distances[h2h.to][c] * 0.08;
                    for (i32 n : esPerP[c]) {
                        f64 nd = d + distances[c][n] * (ct[c][n] < 32 ? pt[ct[c][n]] : 0.08) + distances[h2h.to][n] * 0.08;
                        if (memo[n] > nd) {
                            memo[n] = nd;
                            pre[n] = c;
                            heap.emplace_back(nd, n);
                            push_heap(A(heap), greater<pair<f64, i32>>());
                        }
                    }
                }
                h2h.path.clear();
                h2h.path.emplace_back(h2h.to);
                h2h.distance = 0;
                for (i32 i = pre[h2h.to]; i != -1; i = pre[i]) {
                    if (ct[h2h.path.back()][i] == 0) {
                        cost += distances[h2h.path.back()][i];
                    }
                    h2h.distance += distances[h2h.path.back()][i];
                    ct[h2h.path.back()][i] += h2h.size;
                    ct[i][h2h.path.back()] += h2h.size;
                    h2h.path.emplace_back(i);
                }
                reverse(A(h2h.path));
                dist += sq(h2h.distance / 2) * h2h.size;
                if (abs(curr - cost * sqrt(dist)) > 1e-9) {
                    la = h;
                }
            }
            hc();
            if (la == -1 || timer.elapses()) break;
        }
        if (timer.elapses()) break;
        for (auto& h2h : h2hs) {
            dist -= sq(h2h.distance / 2) * h2h.size;
            rep (i, 1, S(h2h.path)) {
                ct[h2h.path[i - 1]][h2h.path[i]] -= h2h.size;
                ct[h2h.path[i]][h2h.path[i - 1]] -= h2h.size;
                if (ct[h2h.path[i - 1]][h2h.path[i]] == 0) {
                    cost -= distances[h2h.path[i - 1]][h2h.path[i]];
                }
            }
            memo.assign(1000, 1e9);
            pre.assign(1000, -1);
            q.clear();
            q.emplace_back();
            q.back().v = h2h.from;
            q.back().cost = q.back().dist = 0;
            memo[h2h.from] = q.back().score = cost * sqrt(dist + sq(distances[h2h.from][h2h.to] / 2) * h2h.size);
            while (q.size()) {
                auto n = q[0];
                pop_heap(A(q));
                q.pop_back();
                if (n.v == h2h.to) {
                    cost += n.cost;
                    dist += sq(n.dist / 2) * h2h.size;
                    break;
                }
                if (memo[n.v] != n.score) continue;
                for (i32 ne : esPerP[n.v]) {
                    f64 co = n.cost + (ct[n.v][ne] ? 0 : distances[n.v][ne]);
                    f64 di = n.dist + distances[n.v][ne];
                    f64 score = (cost + co) * sqrt(dist + sq((di + distances[h2h.to][ne]) / 2) * h2h.size);
                    if (memo[ne] > score) {
                        memo[ne] = score;
                        pre[ne] = n.v;
                        q.emplace_back();
                        q.back().cost = co;
                        q.back().dist = di;
                        q.back().score = score;
                        q.back().v = ne;
                        push_heap(A(q));
                    }
                }
            }
            h2h.path.clear();
            h2h.path.emplace_back(h2h.to);
            h2h.distance = 0;
            for (i32 i = pre[h2h.to]; i != -1; i = pre[i]) {
                ct[h2h.path.back()][i] += h2h.size;
                ct[i][h2h.path.back()] += h2h.size;
                h2h.distance += distances[h2h.path.back()][i];
                h2h.path.emplace_back(i);
            }
            reverse(A(h2h.path));
        }
        hc();
    }
    h2hs = bh2hs;
    rep (i, 1, 1000) {
        rep (j, i) {
            auto& p1 = ps2[i];
            auto& p2 = ps2[j];
            if (p1.x == p2.x && p1.y == p2.y) {
                array<i32, 3> q;
                q[0] = j;
                q[1] = i;
                q[2] = 0;
                roads.emplace_back(q);
            }
        }
    }
    cout << 0 << endl;
    cout << roads.size() << endl;
    for (auto& r : roads) {
        cout << r[0] << ' ' << r[1] << ' ' << r[2] << endl;
    }
    vector<i32> a;
    rep (i, 1000) {
        i32 be = -1, en = -1;
        rep (j, S(ps)) {
            if (be == -1 && ps[j].x == ps2[i].x && ps[j].y == ps2[i].y) {
                be = j;
            }
            if (en == -1 && ps[j].x == ps2[to[i]].x && ps[j].y == ps2[to[i]].y) {
                en = j;
            }
        }
        a.clear();
        rep (j, S(h2hs)) {
            auto& h2h = h2hs[j];
            if (h2h.from == be && h2h.to == en) {
                for (auto& p : h2h.path) {
                    a.emplace_back(i2r[p]);
                }
            } else if (h2h.from == en && h2h.to == be) {
                for (auto& p : h2h.path) {
                    a.emplace_back(i2r[p]);
                }
                reverse(A(a));
            }
        }
        if (a[0] != i) a.insert(a.begin(), i);
        if (a.back() != to[i]) a.emplace_back(to[i]);
        a.erase(a.begin());
        cout << a.size();
        for (i32 j : a) cout << ' ' << j;
        cout << endl;
    }
    return 0;
}