ferinの競プロ帳

競プロについてのメモ

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

問題ページ

概要

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

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

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

n \leq 16a,b,c \leq 1000d \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) でできるので XO(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;  
}