AOJ 2606 Perm Query

競プロ AOJ

問題文

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2606

解法

まず，制約を読むと，ある項についてたかだか40回の遷移で元に戻ることがわかるので，愚直にこれを求めます．
このとき，それぞれの項が c[i] 回 で元に戻るとします．
すると，[l, r] が与えられた時，この区間がちょうど最初の状態に戻るのは，L = lcm(c[l], c[l+1], ..., c[r]) 回後であることは明らかです．
そして，各項は L / c[i] 回サイクルを回ることになります．1サイクルでの和s[i] を求めておけば，これは s[i] * L / c[i] を l から r まで足した和が答えになります．

区間の LCM は segment tree で求められるのでそうします．

このままだと O(QNlogN) なのですが，sum[k][i] := （第 i 項までの，サイクルが k であるようなものの s[i] の和）という累積和を持っておけば，k が L の約数であるようなものについてのみ (L / k) * (sum[k][r] - sum[k][l]) を考えて総和を取れば，求めたい答えになるのでそうします．k はたかだか 40 なので間に合います．

ソースコード

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

using ll = long long;

constexpr ll M = 1e9 + 7;
constexpr int MAX_K = 41;

// セグツリは略

int main() {
    int N, Q;
    cin >> N >> Q;
    vector<ll> p(N);
    for(int i = 0; i < N; ++i) {
        cin >> p[i];
    }

    vector<ll> cycle(N);
    vector<vector<ll>> sum(MAX_K, vector<ll>(N + 1));
    for(int i = 0; i < N; ++i) {
        ll s = p[i];
        ll now = p[i];
        int k = 1;
        for(k = 1; k < MAX_K; ++k) {
            if(now != i + 1) {
                now = p[now - 1];
                s = (s + now) % M;
            } else {
                break;
            }
        }
        cycle[i] = k;
        sum[k][i + 1] = s;
        for(int k = 0; k < MAX_K; ++k) {
            sum[k][i + 1] = (sum[k][i] + sum[k][i + 1]) % M;
        }
    }
    segment_tree<LCM> lcm_seg(cycle);
    while(Q--) {
        int l, r;
        cin >> l >> r;
        l--;
        ll L = lcm_seg.query(l, r);
        ll res = 0;
        for(int k = 1; k < MAX_K; ++k) {
            if(L % k == 0) {
                (res += ((L / k) % M) * (sum[k][r] - sum[k][l] + M)) %= M;
            }
        }
        cout << res << endl;
    }
}

AOJ 2617 Air Pollution

競プロ AOJ

問題文

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2617

解法

この問題を解く上でポイントとなるのは，大気を循環させたときに，(p[i-1], p[i], p[i+1]) の累積和が不変であるということです．

今大気の状態が (a, b, c, d) であるとしましょう．累積和を取ると，(a, a + b, a + b + c, a + b + c + d) です．
b を選んで大気を循環させてみます．すると，(a + b, -b, b + c, d) となります．累積和を取ると，(a + b, a, a + b + c, a + b + c + d) になっています．
よく見ると，累積和で最初の2つがswapされているだけだということに気づきます．つまり大気の循環は，累積和を取った上での隣接2項間のswapに対応します．

さて，最終的に l[i] が達成すべき大気の綺麗さですが，これは正の値なので，累積和を取ると単調増加になっているはずです．
つまり，もし達成できるなら，これは最初の累積和の転倒数を求める問題と同じです．
達成できない場合というのは，累積和をソートした初項が負であるか，隣接二項の差が l[i] より小さいかのどちらかですから，これを確認すれば良いです．

ソースコード

#include <bits/stdc++.h>
using namespace std;
 
using pii = pair<int, int>;
using ll = long long;
 
template <typename T>
class fenwick_tree {
public:
    fenwick_tree(int n) : n(n), dat(n, 0) {}
 
    void add(int i, T value) {
        for(; i<n; i |= i + 1) {
            dat[i] += value;
        }
    }
 
    T sum(int i) const {
        T res = 0;
        for(; i>=0; i = (i & (i+1)) - 1) {
            res += dat[i];
        }
        return res;
    }
    T sum(int l, int r) const {
        return sum(r-1) - sum(l-1);
    }
 
private:
    const int n;
    std::vector<T> dat;
};
 
 
int main() {
    int n;
    scanf("%d", &n);
    vector<pair<ll, int>> v(n);
    vector<int> l(n);
    for(int i = 0; i < n; ++i) {
        int p;
        scanf("%d", &p);
        if(i == 0) {
            v[i] = make_pair(p, i);
        } else {
            v[i] = make_pair(v[i - 1].first + p, i);
        }
    }
    for(int i = 0; i < n; ++i) {
        cin >> l[i];
    }
    sort(v.begin(), v.end() - 1);
    bool check = v[0].first >= l[0];
    for(int i = 1; i < n; ++i) {
        check &= v[i].first - v[i - 1].first >= l[i];
    }
    if(!check) {
        cout << -1 << endl;
        return 0;
    }
 
    fenwick_tree<int> bit(n);
    ll res = 0;
    for(int i = 0; i < n - 1; ++i) {
        res += i - bit.sum(v[i].second);
        bit.add(v[i].second, 1);
    }
    printf("%lld\n", res);
}

感想

自力で気づけず，教えてもらって解きました．
結構面白い問題だとは思いましたが，AOJ-ICPCで900もあるかと言われると微妙…．

ARC058 E - 和風いろはちゃん

競プロ AtCoder

問題文

http://arc058.contest.atcoder.jp/tasks/arc058_c

解法

X, Y, Z が小さいので bitDP でやります．
dp[i][S] := i 番目まで見たときに，直前の和が X + Y + Z 以下になるところまでの状態が S であるような場合の数
とします．
たとえば，X = 5, Y = 7, Z = 5 のときに，"575" という状態を "10000 1000000 10000" という (X + Y + Z) ビットからなる列として表します．
XYZ であるか，という状態は，よく考えると下から Z ビット目，(Z + Y) ビット目，(X + Y + Z) ビット目が 1 になっているか，という問と同じなので，そういうマスクを作っておけば判定は簡単です．
過去に XYZ を含んだ数列は次からは任意の数を取って良いので，その数列を特別な場所に保存しておくと単に 10 倍するだけで済みます．

ソースコード

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

using ll = long long;

constexpr int mod = 1e9 + 7;

int main() {
    int N, X, Y, Z;
    cin >> N >> X >> Y >> Z;
    int const M = 1 << (X + Y + Z);
    vector<vector<ll>> dp(N + 1, vector<ll>(M + 1));
    int const mask = (1 << (Z - 1)) | (1 << (Z + Y - 1)) | (1 << (Z + Y + X - 1));
    dp[0][0] = 1;

    for(int i = 0; i < N; ++i) {
        for(int S = 0; S < M; ++S) {
            if(dp[i][S] == 0) {
                continue;
            }
            for(int j = 1; j <= 10; ++j) {
                int nS = ((S << j) | (1 << (j - 1))) & (M - 1);
                if((nS & mask) == mask) {
                    (dp[i + 1][M] += dp[i][S]) %= mod;
                } else {
                    (dp[i + 1][nS] += dp[i][S]) %= mod;
                }
            }
        }
        (dp[i + 1][M] += dp[i][M] * 10) %= mod;
    }
    cout << dp[N].back() << endl;
}

感想

bitDP だとわかれば簡単（ほんまか？）

AOJ 1301 Malfatti Circles

競プロ AOJ

問題文

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=1301

解法

まず，一つの円を決め打ちします．
すると，残りの2円の位置は，二分探索などをして求めることが出来ます．
ここで，残りの2円が離れすぎている場合，最初に決め打ちした円が小さすぎるということになります．逆に，近すぎる場合（交わっている場合），最初に決め打ちした円が大きすぎるということになります．
なので，最初の円を二分探索でちょうどいいサイズにすると，答えが求まります．
二分探索の上限の円の中心位置は，三角形の重心にしておかないとめんどくさいことになりそうな気がします．

ソースコード

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

using ld = long double;

constexpr ld eps = 1e-8;

// ライブラリ部分は長いので略．

point rot(point p, ld ang) {
    ld x = real(p) * cos(ang) - imag(p) * sin(ang);
    ld y = real(p) * sin(ang) + imag(p) * cos(ang);
    return point(x, y);
}

circle calc(point cp, ld r, point const& limit, line l1, line l2) {
    point lb = l1.a, ub = limit;
    ld rad = 0;
    for(int i = 0; i < 60; ++i) {
        point m = 0.5l * (lb + ub);
        rad = dist_lp(l1, m);

        if(abs(cp - m) + eps < r + rad) {
            ub = m;
        } else {
            lb = m;
        }
    }
    return circle(lb, rad);
}

int main() {
    int x1, y1, x2, y2, x3, y3;
    while(cin >> x1 >> y1 >> x2 >> y2 >> x3 >> y3) {
        if(x1 == 0 && y1 == 0 && x2 == 0 && y2 == 0 && x3 == 0 && y3 == 0) {
            break;
        }
        point p1(x1, y1), p2(x2, y2), p3(x3, y3);
        point vec1 = rot(p2 - p1, (arg(p3 - p1) + arg(p2 - p1)) / 2 - arg(p2 - p1));
        point vec2 = rot(p1 - p3, (arg(p1 - p3) + arg(p2 - p3)) / 2 - arg(p1 - p3));
        point g = is_ll(line(p1, p1 + vec1), line(p3, p3 + vec2));

        point lb = p1, ub = g; 
        ld r1, r2, r3;
        for(int i = 0; i < 60; ++i) {
            point m = (ub + lb) * 0.5l;
            r1 = dist_lp(line(p1, p2), m);

            auto c2 = calc(m, r1, g, line(p2, p1), line(p2, p3));
            auto c3 = calc(m, r1, g, line(p3, p1), line(p3, p2));
            r2 = c2.r;
            r3 = c3.r;
            if(abs(c2.p - c3.p) + eps < c2.r + c3.r) {
                lb = m;
            } else {
                ub = m;
            }
        }

        cout << setprecision(10) << fixed;
        cout << r1 << ' ' << r2 << ' ' << r3 << endl;
    }
}

感想

多分かなり遅い実装をしていて，ループ回数をそれぞれ100ずつにするとTLEギリギリ通るぐらいでした（座標のサイズ的に100回も回す必要はないんですけどね…）．
まあ通ればなんでもいいよね．速度的には一方の二分探索をやめて，2次方程式で直接求めるほうが早いでしょう．

ARC067 E - Grouping

競プロ AtCoder

問題文

http://arc067.contest.atcoder.jp/tasks/arc067_c

解法

dp[i + 1][j] := i 人以下のグループのみを用いて，あと j 人グループ分けされてない人がいるときの場合の数
という dp で解けます．
遷移は，use == 0 または C <= use <= D のとき
$$\displaystyle dp[i + 1][j - use \times i] = dp[i][j] \times \frac{_jC_i \times _{j-i}C_i \times \ldots \times _{j-(use-1)\times i}C_i}{(use)!}$$
になります．これは，j 人から i 人を選ぶ，という操作を use 回繰り返すが，グループ自体の区別はつけないので，グループの順列で割るということです．

ソースコード

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

// mod は略

int main() {
    int N, A, B, C, D;
    cin >> N >> A >> B >> C >> D;
    vector<vector<mod>> dp(B + 2, vector<mod>(N + 1));
    dp[A][N] = 1;
    for(int i = A; i <= B; ++i) {
        for(int j = 0; j <= N; ++j) {
            if(dp[i][j] == 0) {
                continue;
            }
            mod t = dp[i][j];
            for(int use = 0; use * i <= j; ++use) {
                if(use == 0 || C <= use && use <= D) {
                    dp[i + 1][j - use * i] += t * fact<MOD>(use, false);
                }
                t *= comb<MOD>(j - use * i, i);
            }
        }
    }

    cout << (int)dp[B + 1][0] << endl;
}

感想

数え上げる問題が苦手すぎる．中学高校で訓練しておけばよかった…．

AOJ 2309 Vector Compression

競プロ AOJ

問題文

まず問題読解が難しかった．
http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2309

問題概要

N次元のベクトルが M 個与えられる．これを v[1], ..., v[M] とする．
M 個のベクトルを好きな順番に並べ替え，その順に記録していく．記録する際にコストがかかる．コストは，|v[i]|^2 か，それ以前に記録していたベクトルの中から好きなもの v[j] を選んで，|v[i] - r * v[j]|^2 のどちらかを選択できる．ただし，r は記録のたびに好きな値を適当にとって良いとする．
この時，すべてのベクトルを記録するのにかかるコストを最小化せよ．

・制約
1 <= N, M <= 100

解法

v[i] を記録するときに v[j] を使う，というのを有向辺 (j, i) （向きに注意）であらわし，そのコストを |v[i] - r * v[j]|^2 の最小値とする．
また，v[i] を直接使う場合は，有向辺 (M, i) で表し，そのコストを |v[i]|^2 とする．
すると，この問題はこうして出来たグラフの最小有向全域木を求めよ，と言い換えることができる．
r の求め方だが，幸い |v[i] - r * v[j]|^2 は明らかに下に凸であり，かつ r で微分するのも簡単なので，極値が簡単に求められる（平方完成してもいい）．

計算量は O(nm^2)．

ソースコード

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

using weight = double;
constexpr weight INF = 1e9;

constexpr weight eps = 1e-8;

struct edge {
    int from, to;
    weight cost;
    
    bool operator<(edge const& e) const {
        return cost < e.cost;
    }
};

using edges = std::vector<edge>;
using graph = std::vector<edges>;

void add_edge(graph& g, int from, int to, weight cost) {
    g[from].push_back(edge{from, to, cost});
}

int scc(std::vector<std::vector<int>>& G, std::vector<int>& cmp) {
    int V = G.size();
    std::vector<std::vector<int>> g(V), rg(V);
    std::vector<bool> used(V, false);
    std::vector<int> vs;
    cmp.resize(V);
    for(int i = 0; i < V; ++i) {
        for(auto to : G[i]) {
            g[i].push_back(to);
            rg[to].push_back(i);
        }
    }
    std::function<void(int)> dfs = [&g, &vs, &used, &dfs](int v) {
        used[v] = true;
        for(auto i : g[v]) {
            if(!used[i]) {
                dfs(i);
            }
        }
        vs.push_back(v);
    };
    std::function<void(int, int)> rdfs = [&rg, &cmp, &used, &rdfs](int v, int k) {
        used[v] = true;
        cmp[v] = k;
        for(int i : rg[v]) {
            if(!used[i]) {
                rdfs(i, k);
            }
        }
    };
    for(int v=0; v<V; ++v) {
        if(!used[v]) {
            dfs(v);
        }
    }
    std::fill(used.begin(), used.end(), false);
    int k = 0;
    for(int i=vs.size()-1; i>=0; --i) {
        if(!used[vs[i]]) {
            rdfs(vs[i], k++);
        }
    }
    return k;
}

std::vector<std::vector<int>> build_graph(std::vector<std::vector<int>> const& g, std::vector<int> const& cmp, int K) {
    int V = g.size();
    std::vector<std::set<int>> s(K);
    std::vector<std::vector<int>> res(K);
    for(int i = 0; i < V; ++i) {
        for(auto to : g[i]) {
            s[cmp[i]].insert(cmp[to]);
        }
    }
    for(int i = 0; i < K; ++i) {
        for(auto j : s[i]) {
            if(i != j) {
                res[i].push_back(j);
            }
        }
    }
    return res;
}


weight minimum_spanning_rooted_arborescence(graph const& g, int root, weight sum = 0) {
    int const n = g.size();

    std::vector<int> rev(n, -1);
    std::vector<weight> cost(n, INF);
    for(int i = 0; i < n; ++i) {
        for(auto& e : g[i]) {
            if(e.cost < cost[e.to]) {
                cost[e.to] = e.cost;
                rev[e.to] = i;
            }
        }
    }
    for(int i = 0; i < n; ++i) {
        if(i != root && rev[i] == -1) { // not exists
            return INF;
        }
    }

    std::vector<std::vector<int>> g2(n);
    for(int i = 0; i < n; ++i) {
        if(root == i) {
            continue;
        }
        g2[rev[i]].push_back(i);
        sum += cost[i];
    }
    std::vector<int> cmp(n);
    int const K = scc(g2, cmp);
    auto nxt = build_graph(g2, cmp, K);
    if(nxt.size() == n) {
        return sum;
    }

    graph ng(nxt.size());
    for(int i = 0; i < n; ++i) {
        for(auto& e : g[i]) {
            if(cmp[i] == cmp[e.to]) {
                continue;
            }
            ng[cmp[i]].push_back(edge{cmp[i], cmp[e.to], e.cost - cost[e.to]});
        }
    }

    return minimum_spanning_rooted_arborescence(ng, cmp[root], sum);
}

int main() {
    int N, M;
    cin >> N >> M;
    vector<vector<double>> v;
    for(int i = 0; i < M; ++i) {
        double t = 0;
        vector<weight> w(N);
        for(int j = 0; j < N; ++j) {
            cin >> w[j];
            t += w[j] * w[j];
        }
        if(t > eps) {
            v.push_back(w);
        }
    }
    M = v.size();
    graph g(M + 1);
    int const root = M;
    for(int i = 0; i < M; ++i) {
        for(int j = 0; j < M; ++j) {
            if(i == j) {
                continue;
            }
            double prod = 0, dom = 0;
            for(int k = 0; k < N; ++k) {
                prod += v[i][k] * v[j][k];
                dom += v[j][k] * v[j][k];
            }
            double r = prod / dom;
            double cost = 0;
            for(int k = 0; k < N; ++k) {
                cost += (v[i][k] - r * v[j][k]) * (v[i][k] - r * v[j][k]);
            }
            add_edge(g, j, i, cost);
        }
        double d = 0;
        for(int j = 0; j < N; ++j) {
            d += v[i][j] * v[i][j];
        }
        add_edge(g, root, i, d);
    }
    cout << fixed << setprecision(10) << minimum_spanning_rooted_arborescence(g, root) << endl;
}

感想

最初フローかなあと思っていて全然解けず困っていた．
最小有向全域木を書いたのははじめてなのでいい勉強になりました．

AOJ 2446 Enumeration

競プロ AOJ

問題文

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2446&lang=jp

解法

簡単のために，試しに2つの数字 a, b を用意して考えてみます．a, b が選ばれる確率をそれぞれ p, q とします．
まず，a のみを選んだ場合の期待値は，

  • (m / a) * p * (1 - q)

となります．
b のみを選んだ場合は，

  • (m / b) * (1 - p) * q

となります．
最後に a と b の両方を選んだ場合は，

  • (m / a + m / b - m / lcm(a, b)) * p * q

となります．
これらを足し合わせてみると，以下のようなきれいな式になります．

  • (m / a) * p + (m / b) * q - (m / lcm(a, b)) * p * q

これは一般の場合でも成り立つので，結局包除原理で簡単に求めることが出来ます． 2^n 通り試してやればよいです．
lcm のオーバーフローに注意．

ソースコード

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

constexpr ll INF = 2e18;
constexpr ll invalid = -1;

ll lcm(ll a, ll b) {
    ll g = __gcd(a, b);
    a /= g;
    if(a > INF / b) {
        return invalid;
    }
    return a * b;
}

int main() {
    int n;
    ll m;
    cin >> n >> m;
    vector<ll> a(n);
    vector<double> p(n);
    for(int i = 0; i < n; ++i) {
        cin >> a[i];
    }
    for(int i = 0; i < n; ++i) {
        cin >> p[i];
        p[i] /= 100;
    }

    double res = 0;
    for(int S = 1; S < 1 << n; ++S) {
        int cnt = 0;
        double q = 1;
        for(int i = 0; i < n; ++i) {
            if((S >> i) & 1) {
                cnt++;
                q *= p[i];
            }
        }

        ll t = 1;
        for(int i = 0; i < n; ++i) {
            if((S >> i) & 1) {
                t = lcm(t, a[i]);
            }
            if(t == invalid) {
                break;
            }
        }
        if(t != invalid) {
            if(cnt & 1) {
                res += q * (m / t);
            } else {
                res -= q * (m / t);
            }
        }
    }

    cout << fixed << setprecision(10) << res << endl;
}