2020-05-03

AtCoder Regular Contest 037 D - Chaotic Polygons

AtCoder

レベルが $L$ のときの答え、レベルが $L$ のときに2点を結ぶ方法、レベルが $L$ のときに同じ点を2回通るパターンを各レベルについて計算していく。それぞれ $r_1(L), r_2(L), r_3(L)$ と表す。

$r1(L) =$ 各三角形内で完結するパターン + 三角形をまたぐパターン $= r1(L-1) \times 3 + r2(L-1)^3$
$r2(L) =$ 三角形を2つ使うパターン + 三角形を3つ使うパターン $= r2(L-1)^2 + (r2(L-1)^3 - r3(L-1)^2 \times r2(L-1))$
$r3(L) = r2(L-1)^2 \times r3(L-1) - r3(L-1)^3$

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for(ll i=(ll)a; i<(ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
constexpr ll INF = 1LL<<60;  
  
template<ll MOD>  
struct modint {  
    ll x;  
    modint(): x(0) {}  
    modint(ll y) : x(y>=0 ? y%MOD : y%MOD+MOD) {}  
    static constexpr ll mod() { return MOD; }  
    // e乗  
    modint pow(ll e) {  
        ll a = 1, p = x;  
        while(e > 0) {  
            if(e%2 == 0) {p = (p*p) % MOD; e /= 2;}  
            else {a = (a*p) % MOD; e--;}  
        }  
        return modint(a);  
    }  
    modint inv() const {  
        ll a=x, b=MOD, u=1, y=1, v=0, z=0;  
        while(a) {  
            ll q = b/a;  
            swap(z -= q*u, u);  
            swap(y -= q*v, v);  
            swap(b -= q*a, a);  
        }  
        return z;  
    }  
    // Comparators  
    bool operator <(modint b) { return x < b.x; }  
    bool operator >(modint b) { return x > b.x; }  
    bool operator<=(modint b) { return x <= b.x; }  
    bool operator>=(modint b) { return x >= b.x; }  
    bool operator!=(modint b) { return x != b.x; }  
    bool operator==(modint b) { return x == b.x; }  
    // Basic Operations  
    modint operator+(modint r) const { return modint(*this) += r; }  
    modint operator-(modint r) const { return modint(*this) -= r; }  
    modint operator*(modint r) const { return modint(*this) *= r; }  
    modint operator/(modint r) const { return modint(*this) /= r; }  
    modint &operator+=(modint r) {  
        if((x += r.x) >= MOD) x -= MOD;  
        return *this;  
    }  
    modint &operator-=(modint r) {  
        if((x -= r.x) < 0) x += MOD;  
        return *this;  
    }  
    modint &operator*=(modint r) {  
    #if !defined(_WIN32) || defined(_WIN64)  
        x = (ll)x * r.x % MOD; return *this;  
    #endif  
        unsigned long long y = (unsigned long long)x * r.x;  
        unsigned xh = (unsigned) (y >> 32), xl = (unsigned) y, d, m;  
        asm(  
            "divl %4; \n\t"  
            : "=a" (d), "=d" (m)  
            : "d" (xh), "a" (xl), "r" (MOD)  
        );  
        x = m;  
        return *this;  
    }  
    modint &operator/=(modint r) { return *this *= r.inv(); }  
    // increment, decrement  
    modint operator++() { x++; return *this; }  
    modint operator++(signed) { modint t = *this; x++; return t; }  
    modint operator--() { x--; return *this; }  
    modint operator--(signed) { modint t = *this; x--; return t; }  
    // 平方剰余のうち一つを返す なければ-1  
    friend modint sqrt(modint a) {  
        if(a == 0) return 0;  
        ll q = MOD-1, s = 0;  
        while((q&1)==0) q>>=1, s++;  
        modint z=2;  
        while(1) {  
            if(z.pow((MOD-1)/2) == MOD-1) break;  
            z++;  
        }  
        modint c = z.pow(q), r = a.pow((q+1)/2), t = a.pow(q);  
        ll m = s;  
        while(t.x>1) {  
            modint tp=t;  
            ll k=-1;  
            FOR(i, 1, m) {  
                tp *= tp;  
                if(tp == 1) { k=i; break; }  
            }  
            if(k==-1) return -1;  
            modint cp=c;  
            REP(i, m-k-1) cp *= cp;  
            c = cp*cp, t = c*t, r = cp*r, m = k;  
        }  
        return r.x;  
    }  
  
    template<class T>  
    friend modint operator*(T l, modint r) { return modint(l) *= r; }  
    template<class T>  
    friend modint operator+(T l, modint r) { return modint(l) += r; }  
    template<class T>  
    friend modint operator-(T l, modint r) { return modint(l) -= r; }  
    template<class T>  
    friend modint operator/(T l, modint r) { return modint(l) /= r; }  
    template<class T>  
    friend bool operator==(T l, modint r) { return modint(l) == r; }  
    template<class T>  
    friend bool operator!=(T l, modint r) { return modint(l) != r; }  
    // Input/Output  
    friend ostream &operator<<(ostream& os, modint a) { return os << a.x; }  
    friend istream &operator>>(istream& is, modint &a) {   
        is >> a.x;  
        a.x = ((a.x%MOD)+MOD)%MOD;  
        return is;  
    }  
    friend string to_frac(modint v) {  
        static map<ll, PII> mp;  
        if(mp.empty()) {  
            mp[0] = mp[MOD] = {0, 1};  
            FOR(i, 2, 1001) FOR(j, 1, i) if(__gcd(i, j) == 1) {  
                mp[(modint(i) / j).x] = {i, j};  
            }  
        }  
        auto itr = mp.lower_bound(v.x);  
        if(itr != mp.begin() && v.x - prev(itr)->first < itr->first - v.x) --itr;  
        string ret = to_string(itr->second.first + itr->second.second * ((int)v.x - itr->first));  
        if(itr->second.second > 1) {  
            ret += '/';  
            ret += to_string(itr->second.second);  
        }  
        return ret;  
    }  
};  
using mint = modint<1000000007>;  
  
int main() {  
    ll l;  
    cin >> l;  
  
    struct state { mint r1, r2, r3; };  
    auto dfs = [&](auto &&self, ll v) -> state {  
        if(v == 0) {  
            state ret({1, 2, 1});  
            // dump(v, ret);  
            return ret;  
        }  
        state st = self(self, v-1);  
        state ret({  
            3 * st.r1 + st.r2*st.r2*st.r2,  
            st.r2*st.r2 + st.r2*st.r2*st.r2 - st.r3*st.r2*st.r3,  
            st.r2*st.r3*st.r2 - st.r3*st.r3*st.r3  
        });  
        // dump(v, ret);  
        return ret;  
    };   
  
    cout << dfs(dfs, l).r1 << endl;  
  
    return 0;  
}

2020-05-02

チーム練習 (04/29)

codeforces

2020 Petrozavodsk Winter Camp, Jagiellonian U Contest をやった。本番51位相当、参加者のレベル帯めっちゃ高くないか

(開始) 手分けして読む。L,Bが解かれているので考える。
(0:22) Lは貪欲でいいことがわかるのでolpheが書きはじめる。
(0:25) Lが通る。
(0:39) Bは高速ゼータ変換でいいことにdivくんが気づいて書く。
(0:59) Iが貪欲にできる限り大きいものから取っていいことにolpheが気づいて書き始める。
(1:37) Gについて色々考えた結果、マッチング全探索して2点間を直接線分で結ぶと一つは解が存在するに反例が見つからなかったので書くことになった。
(1:55) Gを出すとn<=6なのにMLE(？？？) 未定義を踏んでないか確認したけどなさそう。なんかやばそうなバグらせかたをしていそうなので実装をolpheにかわる。二人がHの考察を詰めてて、実装できる段階になった。
(2:29) Gが通る。divくんがHの実装をはじめる。Eを考えるけどつらそうなのでFを考える。隣接swapとreverseだけで揃える場合は転倒数で計算できそうなので、シャッフルをする場合の期待値を求めたい。頑張ればできそうな気もするけど遷移がよくわからない。
(3:09) Hが通る。Jを考える。
(3:37) Jは状態をうまくもってUFでつないで管理するとできるっぽいのでolpheが実装を始める。EJあたりを考えたけど通らなくて終了。

解説を読むとEが誤読だったことが判明マジでごめん
Jは連結にするべきパターンが少し抜けてたらしい
cin高速化を複数回実行していたのを1回になるようにしたらGのMLEが取れた。標準入出力のあとにcin高速化のいつもの関数を呼ぶのは処理系定義らしく、気をつけた方がよさそう。
Gはもうちょっと早く実装に入ってもよかったかなあ

A 2020 Petrozavodsk Winter Camp, Jagiellonian U Contest A問題 Bags of Candies - ferinの競プロ帳
E 2020 Petrozavodsk Winter Camp, Jagiellonian U Contest E問題 Contamination - ferinの競プロ帳
F 2020 Petrozavodsk Winter Camp, Jagiellonian U Contest F問題 The Halfwitters - ferinの競プロ帳
J 2020 Petrozavodsk Winter Camp, Jagiellonian U Contest J問題 Space Gophers - ferinの競プロ帳

C 解説の $x_i \geq 0$ の条件を消せると書いてある部分がよくわからない

2020-05-02

2020 Petrozavodsk Winter Camp, Jagiellonian U Contest A問題 Bags of Candies

codeforces

問題ページ

概要

$1 \sim N$ の $N$ 種類の味のキャンディがある
バッグに1か2種類を入れる
ただし2種類の場合、味の $\text{gcd}$ が2以上でなければならない
バッグの数を最小化しろ
$N \leq 10^{11}$

解法

2種類のペアをつくる数を最大化すればバッグの数を最小化できる。味1のキャンディと $N/2$ より大きいの素数の味のキャンディをペアの片方にすることは不可能である。これらの味の集合を $D$ とすると、ペアの個数は $\lfloor \frac{N-|D|}{2} \rfloor$ 以下である。

この上限を必ず達成できることを示す。 $f(x)=(x$ の素因数のうち最大のもの $)$ とする。 $D$ に属さない味のうち、 $f(x)$ が同じ味のものをまとめた集合を考える。要素数が偶数の集合についてはこれらの集合内でペアをつくればよい。要素数が奇数の集合については $2f(x)$ 以外の要素でペアをつくり、余った $2f(x)$ については別の集合同士でペアをつくればよい。したがって、 $D$ に属さない味のキャンディでペアを構成することは常に可能である。

あとは $N/2$ より大きい素数の個数を高速に求められれば良い。Counting Primes の高速な提出を写経するとこれはできる。

想定解に $10^7$ ごとに素数の数を数えておいて埋め込めばよいと書かれていたので計算したらこどふぉのファイル長の制限に引っかかったの何？

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for(ll i=(ll)a; i<(ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
constexpr ll INF = 1LL<<60;  
  
// N以下の素数の数を数える  
// https://judge.yosupo.jp/submission/7976 の写経でN<=10^11 が25ms  
__attribute__((target("avx"), optimize("O3", "unroll-loops")))  
ll prime_count(const ll N) {  
    if(N <= 1) return 0;  
    if(N == 2) return 1;  
    const int v = sqrtl(N);  
    int s = (v+1)/2;  
    vector<int> smalls(s), roughs(s);   
    for(int i=1; i<s; ++i) smalls[i] = i;  
    for(int i=0; i<s; ++i) roughs[i] = 2*i+1;  
    vector<ll> larges(s);   
    for(int i=0; i<s; ++i) larges[i] = (N/(2*i+1)-1)/2;  
    vector<bool> skip(v+1);  
    const auto divide = [](ll n, ll d) -> int { return double(n) / d; };  
    const auto half = [](int n) -> int { return (n - 1) >> 1; };  
    int pc = 0;  
    for(int p=3; p<=v; p+=2) if(!skip[p]) {  
        int q = p*p;  
        if (ll(q)*q > N) break;  
        skip[p] = true;  
        for(int i=q; i<=v; i+=2*p) skip[i] = true;  
        int ns=0;  
        for(int k=0; k<s; ++k) {  
            int i = roughs[k];  
            if(skip[i]) continue;  
            ll d = ll(i) * p;  
            larges[ns] = larges[k] - (d<=v ? larges[smalls[d>>1]-pc] : smalls[half(divide(N, d))]) + pc;  
            roughs[ns++] = i;  
        }  
        s = ns;  
        for(int i=half(v), j=((v/p)-1)|1; j>=p; j-=2) {  
            int c = smalls[j>>1]-pc;  
            for(int e=(j*p)>>1; i>=e; --i) smalls[i] -= c;  
        }  
        ++pc;  
    }  
    larges[0] += ll(s+2*(pc-1))*(s-1)/2;  
    for(int k=1; k<s; ++k) larges[0] -= larges[k];  
    for(int l=1; l<s; ++l) {  
        int q = roughs[l];  
        ll M = N/q;  
        int e = smalls[half(M / q)]-pc;  
        if(e<l+1) break;  
        ll t = 0;  
        for(int k = l + 1; k <= e; ++k) t += smalls[half(divide(M, roughs[k]))];  
        larges[0] += t-ll(e-l)*(pc+l-1);  
    }  
    return larges[0]+1;  
}  
  
int main() {  
    ll t;  
    cin >> t;  
    REP(i, t) {  
        ll n;  
        cin >> n;  
        ll pair = (n - prime_count(n) + prime_count(n/2) - 1) / 2;  
        cout << n-pair << "\n";  
    }  
  
    return 0;  
}

2020-05-02

2020 Petrozavodsk Winter Camp, Jagiellonian U Contest J問題 Space Gophers

codeforces

問題ページ

トンネルが連結なパターンを2通りに分ける

並行な場合
隣接している4パターンのみ考えられる。mapかsetで存在するトンネルを保持しておけばこの判定はできる。
垂直な場合
x軸に垂直なトンネルとy軸に垂直なトンネルが連結となるのはz座標の差が1以下であるときのみである。x軸に垂直でz座標がzのトンネルとx軸に垂直でz座標がz-1,z,zのトンネルを連結にする。連結なトンネルについてすべての集合を持つ必要はないので、トンネルの集合を圧縮する。これをすべてのzについて行う。垂直な方向が違う場合についても対称に同様に計算できる。

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for(ll i=(ll)a; i<(ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
constexpr ll INF = 1LL<<60;  
  
struct UnionFind {  
    vector<ll> par, s;  
    UnionFind(ll n) : par(n), s(n, 1) { iota(ALL(par), 0); }  
    ll find(ll x) {  
        if(par[x] == x) return x;  
        return par[x] = find(par[x]);  
    }  
    void unite(ll x, ll y) {  
        x = find(x), y = find(y);  
        if(x == y) return;  
        if(s[x] < s[y]) par[x] = y, s[y] += s[x];  
        else par[y] = x, s[x] += s[y];  
    }  
    bool same(int x, int y) { return find(x) == find(y); }  
    ll size(int x) { return s[find(x)]; }  
};  
  
void solve() {  
    ll n;  
    cin >> n;  
    vector<ll> x(n), y(n), z(n);  
    map<tuple<ll,ll,ll>, ll> mp;  
    REP(i, n) {  
        cin >> x[i] >> y[i] >> z[i];  
        mp[{x[i], y[i], z[i]}] = i;  
    }  
  
    // 並行で連結な分をつなぐ  
    UnionFind uf(n);  
    REP(i, n) {  
        if(x[i]==-1) {  
            if(mp.find({-1, y[i]-1, z[i]}) != mp.end()) uf.unite(i, mp[{-1, y[i]-1, z[i]}]);  
            if(mp.find({-1, y[i]+1, z[i]}) != mp.end()) uf.unite(i, mp[{-1, y[i]+1, z[i]}]);  
            if(mp.find({-1, y[i], z[i]-1}) != mp.end()) uf.unite(i, mp[{-1, y[i], z[i]-1}]);  
            if(mp.find({-1, y[i], z[i]+1}) != mp.end()) uf.unite(i, mp[{-1, y[i], z[i]+1}]);  
        } else if(y[i]==-1) {  
            if(mp.find({x[i]-1, -1, z[i]}) != mp.end()) uf.unite(i, mp[{x[i]-1, -1, z[i]}]);  
            if(mp.find({x[i]+1, -1, z[i]}) != mp.end()) uf.unite(i, mp[{x[i]+1, -1, z[i]}]);  
            if(mp.find({x[i], -1, z[i]-1}) != mp.end()) uf.unite(i, mp[{x[i], -1, z[i]-1}]);  
            if(mp.find({x[i], -1, z[i]+1}) != mp.end()) uf.unite(i, mp[{x[i], -1, z[i]+1}]);  
        } else {  
            if(mp.find({x[i]-1, y[i], -1}) != mp.end()) uf.unite(i, mp[{x[i]-1, y[i], -1}]);  
            if(mp.find({x[i]+1, y[i], -1}) != mp.end()) uf.unite(i, mp[{x[i]+1, y[i], -1}]);  
            if(mp.find({x[i], y[i]-1, -1}) != mp.end()) uf.unite(i, mp[{x[i], y[i]-1, -1}]);  
            if(mp.find({x[i], y[i]+1, -1}) != mp.end()) uf.unite(i, mp[{x[i], y[i]+1, -1}]);  
        }  
    }  
  
    // 垂直で連結な分をつなぐ  
    {  
        vector<ll> vx, vy, vz;  
        REP(i, n) {  
            if(x[i]==-1) vx.push_back(i);  
            else if(y[i]==-1) vy.push_back(i);  
            else vz.push_back(i);  
        }  
  
        map<ll,vector<ll>> xz, yz, xy, zy, yx, zx;  
        for(auto i: vx) xz[z[i]].push_back(i);  
        for(auto i: vy) yz[z[i]].push_back(i);  
        for(auto i: vx) xy[y[i]].push_back(i);  
        for(auto i: vz) zy[y[i]].push_back(i);  
        for(auto i: vy) yx[x[i]].push_back(i);   
        for(auto i: vz) zx[x[i]].push_back(i);  
  
        auto unite = [&](vector<ll> &a, vector<ll> &b) {  
            if(a.size() && b.size()) {  
                for(auto i: b) uf.unite(a[0], i);  
                for(auto i: a) uf.unite(i, b[0]);  
                a = {uf.find(a[0])};  
                b = {uf.find(b[0])};  
            }  
        };  
  
        // (-1, ya, za) (xb, -1, zb)  
        for(auto i: vx) {  
            if(yz.find(z[i]-1) != yz.end()) unite(xz[z[i]], yz[z[i]-1]);  
            if(yz.find(z[i]) != yz.end()) unite(xz[z[i]], yz[z[i]]);  
            if(yz.find(z[i]+1) != yz.end()) unite(xz[z[i]], yz[z[i]+1]);  
        }  
        // (-1, ya, za) (xb, yb, -1)  
        for(auto i: vx) {  
            if(zy.find(y[i]-1) != zy.end()) unite(xy[y[i]], zy[y[i]-1]);  
            if(zy.find(y[i]) != zy.end()) unite(xy[y[i]], zy[y[i]]);  
            if(zy.find(y[i]+1) != zy.end()) unite(xy[y[i]], zy[y[i]+1]);  
        }  
        // (xa, -1, za) (xb, yb, -1)  
        for(auto i: vy) {  
            if(zx.find(x[i]-1) != zx.end()) unite(yx[x[i]], zx[x[i]-1]);  
            if(zx.find(x[i]) != zx.end()) unite(yx[x[i]], zx[x[i]]);  
            if(zx.find(x[i]+1) != zx.end()) unite(yx[x[i]], zx[x[i]+1]);  
        }  
    }  
  
    map<PII,ll> xy, xz, yz;  
    REP(i, n) {  
        if(x[i]==-1) yz[{y[i], z[i]}] = i;  
        else if(y[i]==-1) xz[{x[i], z[i]}] = i;  
        else xy[{x[i], y[i]}] = i;  
    }  
  
    ll q;  
    cin >> q;  
    REP(i, q) {  
        ll sx, sy, sz, gx, gy, gz;  
        cin >> sx >> sy >> sz >> gx >> gy >> gz;  
  
        ll vs = -1;  
        if(yz.find({sy, sz}) != yz.end()) vs = yz[{sy, sz}];  
        if(xy.find({sx, sy}) != xy.end()) vs = xy[{sx, sy}];  
        if(xz.find({sx, sz}) != xz.end()) vs = xz[{sx, sz}];  
  
        ll vg = -1;  
        if(yz.find({gy, gz}) != yz.end()) vg = yz[{gy, gz}];  
        if(xy.find({gx, gy}) != xy.end()) vg = xy[{gx, gy}];  
        if(xz.find({gx, gz}) != xz.end()) vg = xz[{gx, gz}];  
  
        if(vs!=-1 && vg!=-1 && uf.same(vs, vg)) cout << "YES\n";  
        else cout << "NO\n";  
    }  
}  
  
int main() {  
    ll test;  
    cin >> test;  
    REP(i, test) solve();  
  
    return 0;  
}

2020-05-02

2020 Petrozavodsk Winter Camp, Jagiellonian U Contest F問題 The Halfwitters

codeforces

問題ページ

概要

$n$ 要素の順列が与えられる
以下の操作を繰り返して順列を昇順に並べる

隣接要素のswap： $a$ 分かかる
順列のreverse： $b$ 分かかる
ランダムシャッフル： $c$ 分かかる

初期状態が $d$ 個与えられるのでそれぞれについて、最適に操作したときにかかる時間の期待値を求めろ
出力は分数で行うこと

$n \leq 16$ 、 $a,b,c \leq 1000$ 、 $d \leq 10000$
複数テストケースで $\sum{d} \leq 10^5$

解法

隣接swapとreverseのみを使って順列を昇順に並べるときにかかる時間は転倒数を用いて計算することができる。 $\text{dp}_{\text{inv}} = ($ 転倒数が $\text{inv}$ のときにかかる時間 $)$ とすると、 $\text{dp}_{\text{inv}} = \min(a \times \text{inv}, (n(n-1)/2-a) \times \text{inv} + b)$ となる。

ランダムな順列を整列するのにかかる時間の期待値を $X$ とすると、転倒数が $\text{inv}$ のときの答えは $\text{ans}_{\text{inv}} = \min(\text{dp}_{\text{inv}}, X+c)$ となる。あとは $X$ を求められればよい。 $\text{cnt}_{i} = ($ 転倒数が $i$ となる順列の個数 $)$ とすると、 $X = \sum_{\text{inv}} \text{cnt}_{\text{inv}} * \text{ans}_{\text{inv}}$ となる。 $\text{ans}_{\text{inv}}$ のうち、 $X+c$ になるものと $\text{dp}_{\text{inv}}$ になるものがそれぞれ何かがわかれば $X$ が求められる。

$\text{dp}$ をソートしておくと、ある $j$ について $\text{ans}_{i} = \text{dp}_{i}\ (i \leq j)$ 、 $\text{ans}_{i} = X+c\ (i \geq j+1)$ となる。 $\displaystyle S = \sum_{l=0}^j \text{cnt}_l \cdot \text{dp}_l$ で $\displaystyle M = \sum_{l=j+1}^k \text{cnt}_l$ とする。ただし、 $\text{cnt}$ は $\text{dp}$ のソートに合わせて並べ替えておくものとする。このとき $X \cdot n! = S + M \cdot (X+c)$ であるから、 $X$ を計算することができる。

$j$ を一旦決め打って $X$ を計算し、 $j$ に関する制約を満たすのであれば、そのときの $X$ が正しい値である。

$\text{cnt}$ に関してはDPで計算することができる。 $A \lbrack i \rbrack \lbrack j \rbrack =($ 長さ $i$ で転倒数 $j$ の順列の個数 $)$ とすると、遷移は $A \lbrack i \rbrack \lbrack j \rbrack += A \lbrack i-1 \rbrack \lbrack j-k \rbrack \ (0 \leq k \lt i)$ となる。

計算量について

まず最初に $\text{cnt}$ をDPで計算しておく。これは $O(n^4)$ でできる。

各テストケースについて $\text{dp}$ を計算する。転倒数は高々 $n(n-1)/2$ なので $O(n^2)$ でできる。固定した各 $j$ について $X$ の計算は $O(1)$ でできるので $X$ は $O(n^2)$ でできる。よってテストケース数を $T$ とすると、この部分の計算量は $O(Tn^2)$ である。

与えられた順列について $O(n^2)$ で転倒数を求めることで答えを計算できる。この部分の計算量は $O((\sum d)\cdot n^2)$ である。

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for(ll i=(ll)a; i<(ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
constexpr ll INF = 1LL<<60;  
  
// 分数ライブラリ  
// 常に約分されているとする  
// 負のときは常にaを負にする  
struct fraction {  
    ll a, b;  
    fraction(ll x=0, ll y=1) : a(x), b(y) {  
        ll g = __gcd(a, b);  
        a /= g, b /= g;  
        if(b < 0) a *= -1, b *= -1;  
    }  
    // comparator  
    bool operator<(fraction r) const { return a*r.b < b*r.a; }  
    bool operator>(fraction r) const { return a*r.b > b*r.a; }  
    bool operator<=(fraction r) const { return a*r.b <= b*r.a; }  
    bool operator>=(fraction r) const { return a*r.b >= b*r.a; }  
    bool operator==(fraction r) const { return a*r.b == b*r.a; }  
    bool operator!=(fraction r) const { return a*r.b != b*r.a; }  
    // Basic Operations  
    fraction operator+(fraction r) const { return fraction(*this) += r; }  
    fraction operator-(fraction r) const { return fraction(*this) -= r; }  
    fraction operator*(fraction r) const { return fraction(*this) *= r; }  
    fraction operator/(fraction r) const { return fraction(*this) /= r; }  
    fraction &operator+=(fraction r) {  
        ll g = __gcd(abs(a*r.b+b*r.a), b*r.b);  
        ll na = (a*r.b+b*r.a)/g, nb = (b*r.b)/g;  
        a = na, b = nb;  
        return *this;  
    }  
    fraction &operator-=(fraction r) {  
        ll g = __gcd(abs(a*r.b-b*r.a), b*r.b);  
        ll na = (a*r.b-b*r.a)/g, nb = (b*r.b)/g;  
        a = na, b = nb;  
        return *this;  
    }  
    fraction &operator*=(fraction r) {  
        ll g1 = __gcd(a, r.a), g2 = __gcd(b, r.b);  
        a = a/g1*r.a, b = b/g2*r.b;  
        return *this;  
    }  
    fraction &operator/=(fraction r) {  
        ll g1 = __gcd(a, r.b), g2 = __gcd(b, r.a);  
        a = a/g1*r.b, b = b/g2*r.a;  
        if(b < 0) a *= -1, b *= -1;  
        return *this;  
    }  
    friend fraction abs(fraction a) {  
        a.a = abs(a.a);  
        return a;  
    }  
    // output  
    friend ostream &operator<<(ostream& os, fraction a) {  
        return os << a.a << "/" << a.b;  
    }  
};  
  
int main() {  
    // cnt[i][j] = (長さiで転倒数がjの順列の個数)  
    vector<vector<ll>> cnt(17);  
    FOR(i, 1, 17) {  
        cnt[i].resize(i*(i-1)/2+1);  
        if(i==1) {  
            cnt[i][0] = 1;  
            continue;  
        }  
  
        REP(j, i*(i-1)/2+1) REP(k, i) {  
            if(0<=j-k && j-k<=(i-1)*(i-2)/2) {  
                cnt[i][j] += cnt[i-1][j-k];  
            }  
        }  
    }  
  
    vector<ll> fact(17);  
    fact[0] = 1;  
    FOR(i, 1, 17) fact[i] = fact[i-1] * i;  
  
    auto solve = [&] {  
        ll n, a, b, c, d;  
        cin >> n >> a >> b >> c >> d;  
  
        // dp[i] = (転倒数がiの順列を隣接swap,reverseで並べる時間)  
        vector<ll> dp(n*(n-1)/2+1);  
        REP(i, n*(n-1)/2+1) dp[i] = min(a*i, a*(n*(n-1)/2-i) + b);  
  
        // x = (randomな順列から整列するのにかかる時間の期待値)  
        auto calc = [&] {  
            vector<PII> sdp(n*(n-1)/2+1);  
            REP(i, n*(n-1)/2+1) sdp[i] = PII(dp[i], i);  
            sort(ALL(sdp));  
            vector<ll> scnt(n*(n-1)/2+1);  
            REP(i, n*(n-1)/2+1) scnt[i] = cnt[n][sdp[i].second];  
  
            ll s = 0, m = 0;  
            REP(i, n*(n-1)/2+1) m += scnt[i];  
            fraction tx(s+m*c, fact[n]-m);  
            if(sdp[0].first*tx.b >= tx.a+c*tx.b) return tx;  
            REP(i, n*(n-1)/2+1) {  
                s += scnt[i] * sdp[i].first;  
                m -= scnt[i];  
  
                fraction tx(s+m*c, fact[n]-m);  
                bool cond1 = sdp[i].first*tx.b <= tx.a+c*tx.b;  
                bool cond2 = i+1>n*(n-1)/2 || sdp[i+1].first*tx.b >= tx.a+c*tx.b;   
                if(cond1 && cond2) return tx;  
            }  
            assert(false);  
        };  
        fraction x = calc();  
  
        // ans[i] = (転倒数がiのときの答え)  
        vector<fraction> ans(n*(n-1)/2+1);  
        REP(i, n*(n-1)/2+1) {  
            if(dp[i]*x.b < x.a + c*x.b) ans[i] = fraction(dp[i], 1);  
            else ans[i] = x + fraction(c, 1);  
        }  
  
        REP(i, d) {  
            vector<ll> p(n);  
            REP(j, n) cin >> p[j], p[j]--;  
  
            ll inv = 0;  
            REP(j, n) REP(k, j) if(p[k] > p[j]) inv++;  
            cout << ans[inv] << "\n";  
        }  
    };  
  
    ll test;  
    cin >> test;  
    REP(i, test) solve();  
  
    return 0;  
}

2020-05-02

2020 Petrozavodsk Winter Camp, Jagiellonian U Contest E問題 Contamination

codeforces

問題ページ

概要

$n$ 個の円(中心が $(c_x,c_y)$ で半径 $r$ ) が存在
$q$ 個のクエリが与えられる
2点 $(p_x,p_y),(q_x,q_y)$ を以下の条件を満たしつつ行き来できるか？

円の内部を通らない
$y_{\min}\leq y \leq y_{\max}$ となる範囲のみ

$n,q \leq 10^6$
円は重ならない

解法

$p_x \leq q_x$ と仮定しても一般性を失わない。クエリが'NO'となる条件は「 $p_x \leq c_x \leq q_x$ かつ $y_{\min} \leq c_y-r \leq c_y+r \leq y_{\max}$ となる円が存在する」ことである。

この判定は平面走査でできる。 $y=\text{const}$ のときに $\text{seg} \lbrack i \rbrack = (x=i$ を含む円の上端 $)$ となるようにセグ木を保持して管理する。 $\text{const}$ を小さい方から動かしていく。このときセグ木を更新するタイミングは以下の3パターンである。

円の下端
$\text{seg} \lbrack c_x \rbrack = c_y+r$ と更新
円の上端
$\text{seg} \lbrack c_x \rbrack = -\infty$ と更新
$y_{\min}$ に達したタイミング
区間 $\lbrack p_x, q_x \rbrack$ の $\max \geq y_{\max}$ であれば'NO'

したがって $O( (q+n)\log(q+n) )$ で解けた。

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for(ll i=(ll)a; i<(ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
constexpr ll INF = 1LL<<60;  
  
// 根が1  
template<typename Monoid>  
struct segtree {  
    using T = typename Monoid::T;  
    int n;  
    vector<T> dat;  
  
    segtree(int n_) {  
        n = 1;  
        while(n < n_) n <<= 1;  
        dat.assign(n*2, Monoid::id());  
    }  
    void build(vector<T> v) {  
        REP(i, v.size()) dat[i+n] = v[i];  
        for(int i=n-1; i>0; --i) dat[i] = Monoid::op(dat[i*2], dat[i*2+1]);  
    }  
  
    T query(int a, int b) {  
        if(a >= b) return Monoid::id();  
        T l = Monoid::id(), r = Monoid::id();  
        for(a+=n, b+=n; a<b; a>>=1, b>>=1) {  
            if(a&1) l = Monoid::op(l, dat[a++]);  
            if(b&1) r = Monoid::op(dat[--b], r);  
        }  
        return Monoid::op(l, r);  
    }  
    void update(int i, T v) {  
        i += n;  
        dat[i] = v;  
        while(i >>= 1) {  
            dat[i] = Monoid::op(dat[i*2], dat[i*2+1]);  
        }  
    }  
  
    friend ostream &operator <<(ostream& out,const segtree<Monoid>& seg){  
        out << "---------------------" << endl;  
        int cnt = 1;  
        for(int i=1; i<=seg.n; i*=2) {  
            REP(j, i) {  
                out << seg.dat[cnt] << " ";  
                cnt++;  
            }  
            out << endl;  
        }  
        out << "---------------------" << endl;  
        return out;  
    }  
};  
  
template<typename Tp>  
struct max_monoid {  
    using T = Tp;  
    static constexpr Tp id() { return -INF; }  
    static Tp op(const Tp &a, const Tp &b) { return max(a, b); }  
};  
  
int main() {  
    ll n, q;  
    cin >> n >> q;  
    vector<PII> ord(q+2*n);  
    vector<ll> cx(n), cy(n), r(n);  
    REP(i, n) {  
        cin >> cx[i] >> cy[i] >> r[i];  
        ord[2*i] = {i, 0};  
        ord[2*i+1] = {i, 2};  
    }  
    vector<ll> px(q), py(q), qx(q), qy(q), ymin(q), ymax(q);  
    REP(i, q) {  
        cin >> px[i] >> py[i] >> qx[i] >> qy[i] >> ymin[i] >> ymax[i];  
        if(px[i] > qx[i]) swap(px[i], qx[i]), swap(py[i], qy[i]);  
        ord[2*n + i] = {i, 1};  
    }  
    sort(ALL(ord), [&](PII a, PII b){  
        ll vl, vr;  
        if(a.second == 1) vl = ymin[a.first];  
        else if(a.second == 0) vl = cy[a.first] - r[a.first];  
        else if(a.second == 2) vl = cy[a.first] + r[a.first];  
        if(b.second == 1) vr = ymin[b.first];  
        else if(b.second == 0) vr = cy[b.first] - r[b.first];  
        else if(b.second == 2) vr = cy[b.first] + r[b.first];  
        return PII(vl, a.second) < PII(vr, b.second);  
    });  
  
    vector<ll> vs(2*q+n);  
    REP(i, q) vs[2*i] = px[i], vs[2*i+1] = qx[i];  
    REP(i, n) vs[2*q+i] = cx[i];  
    sort(ALL(vs));  
    vs.erase(unique(ALL(vs)), vs.end());  
    REP(i, q) {  
        px[i] = lower_bound(ALL(vs), px[i]) - vs.begin();  
        qx[i] = lower_bound(ALL(vs), qx[i]) - vs.begin();   
    }  
    REP(i, n) cx[i] = lower_bound(ALL(vs), cx[i]) - vs.begin();  
  
    vector<bool> ans(q, true);  
    segtree<max_monoid<ll>> seg(vs.size());  
    seg.build(vector<ll>(vs.size(), -INF));  
    for(auto [i, type]: ord) {  
        if(type == 1) {  
            if(ymax[i] <= seg.query(px[i], qx[i]+1)) ans[i] = false;  
        } else if(type == 0) {  
            seg.update(cx[i], cy[i]+r[i]);  
        } else if(type == 2) {  
            seg.update(cx[i], -INF);  
        }  
    }  
  
    REP(i, q) cout << (ans[i] ? "YES" : "NO") << "\n";  
  
    return 0;  
}

2020-04-27

天下一プログラマーコンテスト2014予選B C - 天下一王国の歴史

AtCoder

問題ページ

まず答えの一番上の行を適当に定める。このとき入力の1~n-1行目までが正しく生成されるように、答えの2~n行目までを定めると一意に定まる。入力のn行目についても正しく生成されるような、一番上の行を探したい。 $O(2^n)$ 通りを全列挙することは不可能なので、一番上の行を変化させたときに、生成されるn行目がどのように変化するか考える。

$n=10$ のときに $(0,3)$ を変化させた場合に、答えのどこが変化するのかを以下に示す。規則的な変化になっている。

diff:  
    0123456789  
i=0 ...x......  
i=1 ..x.x.....  
i=2 .x.x.x....  
i=3 x.x.x.x...  
i=4 .x.x.x.x..  
i=5 ..x.x.x.x.  
i=6 ...x.x.x.x  
i=7 ....x.x.x.  
i=8 .....x.x..  
i=9 ......x...

ここで一番下の行に注目すると、隣接する変化しているマスの個数は必ず偶数個である。一番上の行をどのように変化させたとしても同様のことが示せ、 $n$ 行目は変化しないことがわかる。入力で与えられる文字列は必ず生成できる制約なので、一番上の行をどのように固定したとしても条件を満たす答えを得ることができる。

#include <bits/stdc++.h>  
using namespace std;  
using ll = long long;  
using PII = pair<ll, ll>;  
#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)  
#define REP(i, n) FOR(i, 0, n)  
#define ALL(x) x.begin(), x.end()  
template<typename T> void chmin(T &a, const T &b) { a = min(a, b); }  
template<typename T> void chmax(T &a, const T &b) { a = max(a, b); }  
struct FastIO {FastIO() { cin.tie(0); ios::sync_with_stdio(0); }}fastiofastio;  
#ifdef DEBUG  
#include "../program_contest_library/memo/dump.hpp"  
#else  
#define dump(...)  
#endif  
const ll INF = 1LL<<60;  
  
int main(void) {  
    ll n;  
    cin >> n;  
    vector<string> v(n);  
    REP(i, n) cin >> v[i];  
  
    ll dx[] = {-1, 0, 1}, dy[] = {0, -1, 0};  
    auto f = [&](string s) {  
        vector<string> ret(n, string(n, '.'));  
        ret[0] = s;  
        FOR(i, 1, n) REP(j, n) {  
            ll cnt = 0;  
            REP(k, 3) {  
                ll ni = i-1+dy[k], nj = j+dx[k];  
                if(0<=ni && ni<n && 0<=nj && nj<n) {  
                    cnt += ret[ni][nj]=='#';  
                }  
            }  
            if(v[i-1][j]=='#') {  
                if(cnt%2==0) ret[i][j] = '#';  
            } else {  
                if(cnt%2) ret[i][j] = '#';  
            }  
        }  
        return ret;  
    };  
  
    auto ret = f(string(n, '.'));  
    REP(i, n) cout << ret[i] << "\n";  
  
    return 0;  
}

結構すき