AOJ 1338 Clock Hands

解法

まず,それぞれの針がなす角を求めましょう.短針が1周するのに$H$時間かかるとします.
今時刻が $h:m:s$ である時,針が上を向いている状態を 0 度とすると,長針,短針,秒針はそれぞれ
$\displaystyle ang_h = \frac{360}{H} \times h + \frac{6}{H} \times m + \frac{s}{10H}$
$\displaystyle ang_m = 6m + \frac{s}{10}$
$ang_s = 6s$
の角度にあります.
簡単にするために,$ang_h = 0$ となるように時計を回転させます.
$ang_h = 0$
$\displaystyle ang_m = 6m + \frac{s}{10} - \left(\frac{360}{H} \times h + \frac{6}{H} \times m + \frac{s}{10H}\right)$
$\displaystyle ang_s = 6s- \left(\frac{360}{H} \times h + \frac{6}{H} \times m + \frac{s}{10H}\right)$
これらの値が負になっていれば,360 を加えて辻褄を合わせます.

秒針と長針の位置関係は二通りありますが,とりあえず $ang_s < ang_m $の場合を考えましょう.
このとき,秒針との短針および長針の間の角はそれぞれ
$\displaystyle ang_s - ang_h = 6s- \left(\frac{360}{H} \times h + \frac{6}{H} \times m + \frac{s}{10H}\right)$
$\displaystyle ang_m - ang_s = 6m - \frac{59s}{10}$
これらが仮に等しいとして等式を立て,式を整理すると,
$\displaystyle s = \frac{n}{d} = \frac{60m(H + 1) + 3600h}{119H - 1}$
が成り立ちます.秒針と長針の位置関係が異なる場合も,似た式になります.

つまり,目的の時刻において,$d$ は $119H - 1$ の約数であるということが分かります.
ここから,$d = 119H - 1$ とおいてしまって,その上で n を少しずつ増やしていく全探索が考えられます.実際,秒針は高々2周程度する間に,目的の時刻は必ず存在します.(最後に $n$ と $d$ を gcd で割れば良い)

浮動小数点を扱ったり,分数をきちんと持つのは面倒なので,すべての式に $10 \times H \times d$ をかけてしまって,一周が $360 \times 10 \times H \times d$ 度であるということにすれば実装が楽です.

計算量も,$H <= 100$ という制約から,十分間に合います.

ソースコード

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

using ll = long long;

int main() {
    int H, h0, m0, s0;
    while(cin >> H >> h0 >> m0 >> s0, H) {
        const ll d = 119 * H - 1;
        ll n = s0 * d;
        ll ho = h0, mo = m0;
        while(true) {
            ll ang_h = 3600 * ho * d + 60 * mo * d + n;
            ll ang_m = 60 * H * mo * d + n * H;
            ll ang_s = 60 * H * n;
            ang_m = (ang_m - ang_h) % (360 * 10 * H * d);
            ang_s = (ang_s - ang_h) % (360 * 10 * H * d);
            if(ang_m < 0) {
                ang_m += 360 * 10 * H * d;
            }
            if(ang_s < 0) {
                ang_s += 360 * 10 * H * d;
            }
            ang_h = 0;

            if(ang_h != ang_m) {
                ll ang1, ang2;
                if(ang_m < ang_s) {
                    ang1 = 360 * 10 * H * d - ang_s;
                    ang2 = ang_s - ang_m;
                } else {
                    ang1 = ang_s;
                    ang2 = ang_m - ang_s;
                }
                if(ang1 == ang2) {
                    break;
                }
            }

            n++;
            if(n == 60 * d) {
                n = 0;
                mo++;
            }
            if(mo == 60) {
                mo = 0;
                ho++;
            }
            if(ho == H) {
                ho = 0;
            }
        }
        ll g = __gcd(n, d);
        cout << ho << ' ' << mo << ' ' << n / g << ' ' << d / g << endl;
    }
}

感想

THE ICPC って感じの問題だった.

ARC057 D - 全域木

解法

メモ化再帰で解きます.
$rec(C) := $連結成分の状態が $C$ であるような状態から始めたときに,条件を満たすような辺の張り方の総数

$A[i]$ を小さい方から張って順に張っていくとき,$A[i]$に含まれない辺は,それまでにクラスカル法の要領で作られていた連結成分同士を結ばないように辺をはる必要があります.
逆に,$A[i]$ に含まれる辺は,連結成分同士を結ぶように辺を張ることになります.

今 $A[i]$ の辺を張ろうとした時,仮に連結成分の集合が $C = \{C_1, C_2, \ldots, C_k\}$ であったとしましょう.この時,重みが $A[i - 1] + 1$ から $A[i] - 1$ である辺は,$\displaystyle \sum_{i=1}^{k} \frac{|C_i|(|C_i| - 1)}{2} - (N - |C|)$ 本のまだ辺が張られていないところから適当に選ぶ必要があるので,これは単純に順列の総数を求めれば良いです.($N - |C|$ はMSTにすでに使われている辺の数)
次に $A[i]$ を張るときは,連結成分から適当に2個選んでそれらを併合して,再帰的に解きます.例えば $C_i$ と $C_j$ を併合するときは,どの頂点同士を結ぶかで $|C_i| \times |C_j|$ 通りあります.
最後にすべての併合を試したものの総和に,$A[i - 1]$ から $A[i]$ までの辺の張り方をかけると答えになります.

$C$ の状態数は $N = 30$ の時高々5000程度(分割数という)なので,十分間に合います.

実装面では,番兵として $A[i]$ に 0 と $\frac{N(N - 1)}{2}$ を最初と最後に入れておくと楽かも?

ソースコード

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

using pii = pair<int, int>;
using ll = long long;

constexpr ll M = 1e9 + 7;

map<vector<int>, ll> memo;

ll rec(vector<int> C, vector<int> const& A) {
    sort(C.begin(), C.end());
    if(memo.count(C) == 1) {
        return memo[C];
    }

    int const N = A.size() - 1;

    ll t = 0; // 連結成分内で,使うと閉路が出来てしまう辺
    for(auto sz : C) {
        t += sz * (sz - 1) / 2;
    }
    int now = N - C.size(); // 次に加える辺
    t -= A[now];
    ll unuse = 1; // MSTで使わない辺の張り方
    for(int i = 0; i < A[now + 1] - A[now] - 1; ++i, --t) {
        unuse = unuse * t % M;
    }
    if(unuse == 0) { // どうしても A[i] 以外の辺を使うハメになる
        return memo[C] = 0;
    }

    // 連結成分同士の結び方
    if(C.size() == 1) {
        return memo[C] = unuse;
    } else {
        ll res = 0;
        int const m = C.size();
        for(int i = 0; i < m; ++i) {
            for(int j = i + 1; j < m; ++j) {
                vector<int> nC;
                nC.push_back(C[i] + C[j]);
                for(int k = 0; k < m; ++k) {
                    if(k != i && k != j) {
                        nC.push_back(C[k]);
                    }
                }
                res = (res + rec(nC, A) * C[i] * C[j]) % M;
            }
        }
        return memo[C] = (res * unuse) % M;
    }
}

int main() {
    int N;
    cin >> N;
    vector<int> A(N + 1);
    for(int i = 1; i < N; ++i) {
        cin >> A[i];
    }
    A[N] = N * (N - 1) / 2 + 1; // 番兵
    vector<int> C(N, 1);
    cout << rec(C, A) << endl;
}

AOJ 2436 Card

解法

カードの数字同士を足すわけではなくそのまま並べるので,桁に注目したいなーという気持ちになります.
しかも $a[i]$ は高々4桁なので嬉しいです.

こういう問題は,すべての数字をまとめてやるのは大変な気がするので,ある $a[i]$に着目してやってみましょう.
つまり,$a[i]$ を $j$ 桁目に並べるような並べ方が何通りあるかさえわかれば,それを $C$ 通りだとすると,$a[i] \times 10^j \times C$ を答えに加えるということをすればよいことになります.

さて,例えば $X$ が $j$ 桁めにあるようなものをカウントしましょう.
それは以下のような感じになります.
.......($X$ の前に数字が並ぶ) $X$ ($X$ のあとに $j$ 桁並ぶ)........
当たり前ですが,前と後ろで独立に考えられます.まず後ろを考えましょう.
後ろは以下のようなdpで求めることが出来ます.
$dp[i][j] := X$ 以外の $i$ 枚のカードを並べて $j$ 桁にする場合の数(ゼロリーディングを含める)
遷移の式は簡単ですね.
しかしこのままでは,注目する数 $X$,今追加して並べようとする数 $a[j]$,すでに並べている枚数 $k$,すでに出来ている桁数 $l$,の4つでループが回るので間に合いません.
しかしよく考えると,同じ桁数である $a[i]$ と $a[j]$ について dp の値は全く同じであることがわかります.なので,この計算はループの外に逃がすことが出来ます.以降,dp のキーに桁数を追加して説明することとします.

さて,dp テーブルが求まったら,前に並ぶ数字も考えて最終的な答えを求めにいきます.
今注目している数字が $a[i]$ で,$a[i]$ が $d[i]$ 桁であり,後ろで $j$ 枚のカードを使っていて,それが $k$ 桁であり,前に $l$ 枚のカードを並べるとき,$\displaystyle \sum_{l = 0}^{n - 1 - j}(_{n-1-j}C_{l} \times l! \times a[i] \times 10^j \times dp[d[i]][j][k])$
残念ながら,これも4重ループになっていますが,先と同様,同じ桁である $a[i]$ をまとめて処理してやるとオーダーが落とせます.$sum[x] := x$桁である$a[i]$の総和とすれば,$a[i]$ の部分を $sum[i]$ にしてやればOKです.

最後にゼロリーディングの処理ですが,これは最初に0である$a[i]$を1つだけ取り除いたカード列について同じ計算をして,最後に0を先頭にくっつけたものを引く,という処理で実現できます.先頭の0をどの0にするか,を書け忘れないように.

計算量は O(4 * N^3) になります.

ソースコード

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

using pii = pair<int, int>;
using ll = long long;

constexpr ll M = 1e9 + 7;

ll ten_pow[1024];
ll fact[256];
ll comb[256][256];

void init() {
    ten_pow[0] = 1;
    for(int i = 1; i < 1024; ++i) {
        ten_pow[i] = (ten_pow[i - 1] * 10) % M;
    }
    fact[0] = 1;
    for(int i = 1; i < 256; ++i) {
        fact[i] = (fact[i - 1] * i) % M;
    }
    for(int i = 0; i < 256; ++i) {
        for(int j = 0; j < 256; ++j) {
            if(i == j || j == 0) {
                comb[i][j] = 1;
            } else {
                comb[i][j] = (comb[i - 1][j - 1] + comb[i - 1][j]) % M;
            }
        }
    }
}

ll solve(vector<int> const& a, vector<int> const& d) {
    int const n = a.size();

    vector<vector<vector<ll>>> dp(5, vector<vector<ll>>(n, vector<ll>(n * 4 + 1)));
    ll res = 0;
    vector<ll> sum(5);
    for(int i = 0; i < n; ++i) {
        sum[d[i]] = (sum[d[i]] + a[i]) % M;
    }
    for(int i = 0; i < n; ++i) {
        if(dp[d[i]][0][0] == 0) {
            dp[d[i]][0][0] = 1;
            for(int j = 0; j < n; ++j) {
                if(i == j) {
                    continue;
                }
                for(int k = n - 1; k >= 0; --k) {
                    for(int l = 0; l < n * 4 + 1; ++l) {
                        if(dp[d[i]][k][l] == 0) {
                            continue;
                        }
                        dp[d[i]][k + 1][l + d[j]] = (dp[d[i]][k + 1][l + d[j]] + dp[d[i]][k][l] * (k + 1)) % M;
                    }
                }
            }
        }
    }
    for(int i = 0; i < 5; ++i) {
        for(int j = 0; j < n; ++j) {
            for(int k = 0; k < 4 * n + 1; ++k) {
                if(dp[i][j][k] == 0) {
                    continue;
                }
                for(int l = 0; l <= n - 1 - j; ++l) {
                    ll t = (sum[i] * dp[i][j][k]) % M;
                    t = (t * ten_pow[k]) % M;
                    t = (t * comb[n - 1 - j][l]) % M;
                    t = (t * fact[l]) % M;
                    res = (res + t) % M;
                }
            }
        }
    }

    return res;
}

int main() {
    init();

    int n;
    cin >> n;
    int zero_cnt = 0;
    vector<int> a(n);
    vector<int> d(n);
    for(int i = 0; i < n; ++i) {
        cin >> a[i];
        d[i] = to_string(a[i]).size();
        zero_cnt += (a[i] == 0);
    }

    ll res = solve(a, d);
    if(zero_cnt != 0) {
        int idx = find(begin(a), end(a), 0) - begin(a);
        a.erase(a.begin() + idx);
        d.erase(d.begin() + idx);
        ll leading_zero = (zero_cnt * solve(a, d)) % M;
        res = (res - leading_zero + M) % M;
    }

    cout << res << endl;
}

感想

自明なDPから考えて,そこから少しずつオーダー落としていくと楽だった.
でもまあ結構めんどくさかった.

AOJ 2725 Live Programming

問題概要

催し物を時間 $T$ だけ開催することになった.
そこで歌う歌の候補が $N$ 個あり,それぞれの歌はパラメーター $t[i], p[i], f[i]$ を持つ.$t[i]$ は歌の長さ,$p[i]$ は歌の満足度,$f[i]$ は歌の特徴量を表します.これらの歌からいくつか選んで,歌うことにする.
この時,以下の制約を満たすような歌の列 $\{s_1, ..., s_k\}$ で,満足度が最大を取るときの満足度を求めよ.

  • 歌う歌の時間の総和が $T$ 以下である.
  • 最初に歌う歌の満足度は $p[s_0]$ とする.
  • 以降に歌う歌の満足度は,$p[s_i] - (f[s_{i-1}] - f[s_i])^2$ である.$s_{i-1}$ は,直前に歌った歌である.

・制約
1 <= N <= 4000
1 <= T <= 4000

解法

DPっぽい感じがするのでDPでやるんですが,その前に考えるべきことがあります.

今回問題になるのは,「どの歌の集合を選ぶか」と「歌をどの順番で並べるか」の2つあり,このままではしんどいです.
では,仮に歌の集合を決め打ちした場合,それらをどう並べるのがいいか考えてみましょう.これは明らかに,f[i] の大きさでソートされた順番に歌うのが最適です.

つまり,あらかじめ f[i] の値で歌をソートしてしまえば,以下のDPでできることは明白です.
$dp[i][j] := i$ 番目の歌を最後に歌うとして,総時間が $j$ であるときの最大満足度
このDPの遷移は以下のとおりです.
$\displaystyle dp[i][j] = \max_{k < i}\{dp[k][j - t[i]] + p[i] - (f[k] - f[i])^2\}$
式を整理して,さらに(別に置き換えなくてもできますが)max を min に置き換えると以下の式になります.
$\displaystyle dp[i][j] = -\min_{k < i}\{-dp[k][j - t[i]] + (f[k])^2 - 2f[k] \times f[i]\} + p[i] - (f[i])^2$
これはConvexHullTrickチャンスなので高速化でき,計算量 O(NT) になります.

注意すべき点は j のループは大きい方から小さい方に回すことです.これは普通のナップザックと同じです.
こういう細かいミスをなくすには,先に自明なDPを書いてから,あとで高速化するとよいですね.(慣れているなら別)

ソースコード

#include <bits/stdc++.h>
using namespace std;
 
using ll = long long;
 
class convex_hull_trick {
    using T = long long;
 
    bool check(int f1, int f2, T aa, T bb) {
        return (a[f2] - a[f1]) * (bb - b[f2]) >= (b[f2] - b[f1]) * (aa - a[f2]);
    }
    bool check(int f1, int f2, int f3) {
        return (a[f2] - a[f1]) * (b[f3] - b[f2]) >= (b[f2] - b[f1]) * (a[f3] - a[f2]);
    }
 
    T f(int idx, T x) {
        return a[idx] * x + b[idx];
    }
 
public:
    convex_hull_trick() = default;
    // a : sorted sequence
    convex_hull_trick(std::vector<T> const& a_, std::vector<T> const& b_) {
        assert(a_.size() == b_.size());
        for(int i = 0; i < a_.size(); ++i) {
            add_f(a_[i], b_[i]);
        }
    }
 
    void add_f(T aa, T bb) {
        while(a.size() >= 2 && check(a.size() - 2, a.size() - 1, aa, bb)) {
            a.pop_back();
            b.pop_back();
        }
        a.push_back(aa);
        b.push_back(bb);
    }
    T min(T x) {
        assert(a.size() > 0);
        while(a.size() >= 2 && f(0, x) >= f(1, x)) {
            a.pop_front();
            b.pop_front();
        }
        return a[0] * x + b[0];
    }
 
    bool empty() const {
        return a.size() == 0;
    }
 
private:
    std::deque<T> a, b; // functions
};
 
int main() {
    int N, T;
    cin >> N >> T;
    // (f, t, p)
    vector<tuple<ll, ll, ll>> songs;
    for(int i = 0; i < N; ++i) {
        ll t, p, f;
        cin >> t >> p >> f;
        if(t > T) {
            continue;
        }
        songs.emplace_back(f, t, p);
    }
    N = songs.size();
    sort(songs.begin(), songs.end());
 
    vector<vector<ll>> dp(N + 1, vector<ll>(T + 1));
    vector<convex_hull_trick> cht(T + 1);
    ll res = 0;
    for(int i = 0; i < N; ++i) {
        ll fi, ti, pi;
        tie(fi, ti, pi) = songs[i];
        for(int j = T; j >= ti; --j) {
            if(cht[j - ti].empty()) {
                continue;
            }
            dp[i][j] = -cht[j - ti].min(fi) - fi * fi + pi;
            res = max(res, dp[i][j]);
            cht[j].add_f(-2 * fi, -dp[i][j] + fi * fi);
        }
        dp[i][ti] = max(dp[i][ti], pi);
        res = max(res, dp[i][ti]);
        cht[ti].add_f(-2 * fi, -dp[i][ti] + fi * fi);
    }
    cout << res << endl;
}

Codeforces Round #185 (Div. 1) B. Cats Transport

問題概要

N個の地点があり,順に 0 から N-1 とする.地点 i と i + 1 の間の距離は d[i] である.
また,猫が m 匹いて,猫 j は時刻 t[j] に地点 h[j] にあらわれるものとする.それぞれの猫は,餌をもらえるまで出現してからその地点で待機し続ける.
今,猫に餌を上げる人が p 人いて,それぞれ地点 0 から N-1 の方向に単位時間あたり1だけ移動する.(ただし初期位置は自由.つまり負の時刻に地点 0 を出発するような人がいても良い.これは制約に書いてないが…)
移動中に餌を待っている猫と遭遇した場合,その猫に必ず餌をやる.ただし,餌をやる時間は無視できるものとする.
では,全ての猫に餌をやるように人を配置した時に,それぞれの猫が餌を待つ時間の総和を最小化せよ.

・制約
2 <= N <= 2 * 10^5
1 <= M <= 10^5
1 <= p <= 100

解法

まず,今のままでは位置と時間という軸が2つあり面倒なので,これをひとつにまとめます.人は単位時間に1だけすすむので,時刻 t[j] に地点 h[j] に現れるという情報を a[j] := t[j] - pos[j] とまとめてしまいます(pos[j] は地点 h[j] の座標).a[j] は,時刻 a[j] に人を地点 0 に配置すれば,その人は猫 j が現れるよりあとに通過する(つまり餌をあげることができる)ことを意味します.
こうして得られた a[j] を昇順にソートしておきます.
時刻 a[j] に人を地点 0 に配置した場合,それ以前の直近で人が配置された時刻を a[k] をすると,k + 1 <= j となるような全ての猫について餌をやることになります.

ここからDPパートです.
p が小さいので,以下のような dp が自然に思い浮かびます.
dp[i][j] := i 人目まで使い,かつ i 人目を時刻 a[j] に出発させた時の,j 以前の猫の待ち時間の総和の最小値
このdpの更新は以下のようになります.
$\displaystyle dp[i][j] = \min_{k < j} \{dp[i - 1][k] + \sum_{l=k+1}^{j}(a[j] - a[l])\}$
式を変形すると,以下のようにできます.
$\displaystyle dp[i][j] = \min_{k < j} \{dp[i - 1][k] - k \times a[j] - \sum_{l=k+1}^{j}a[l]\} + j \times a[j]$
さらに累積和 $\displaystyle S[j] = \sum_{l=0}^{j-1} a[l]$を用いて整理すると
$\displaystyle dp[i][j] = \min_{k < j} \{dp[i - 1][k] - k \times a[j] + S[k+1] \} + j \times a[j] - S[j+1]$
このままだと O(pM^2) ですが,この形のDPは ConvexHullTrick という方法で高速化できます.(この方法自体の解説は省略します.)
よって O(pM) でこの問題を解くことができます.

ソースコード

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

class convex_hull_trick {
    using T = long long;

    bool check(int f1, int f2, T aa, T bb) {
        return (a[f2] - a[f1]) * (bb - b[f2]) >= (b[f2] - b[f1]) * (aa - a[f2]);
    }
    bool check(int f1, int f2, int f3) {
        return (a[f2] - a[f1]) * (b[f3] - b[f2]) >= (b[f2] - b[f1]) * (a[f3] - a[f2]);
    }

    T f(int idx, T x) {
        return a[idx] * x + b[idx];
    }

public:
    convex_hull_trick() = default;
    // a : sorted sequence
    convex_hull_trick(std::vector<T> const& a_, std::vector<T> const& b_) {
        assert(a_.size() == b_.size());
        for(int i = 0; i < a_.size(); ++i) {
            add_f(a_[i], b_[i]);
        }
    }

    void add_f(T aa, T bb) {
        while(a.size() >= 2 && check(a.size() - 2, a.size() - 1, aa, bb)) {
            a.pop_back();
            b.pop_back();
        }
        a.push_back(aa);
        b.push_back(bb);
    }
    T min(T x) {
        assert(a.size() > 0);
        while(a.size() >= 2 && f(0, x) >= f(1, x)) {
            a.pop_front();
            b.pop_front();
        }
        return a[0] * x + b[0];
    }

private:
    std::deque<T> a, b; // functions
};

using namespace std;

using ll = long long;

constexpr ll INF = 2e18;

int main() {
    ll n, m, p;
    cin >> n >> m >> p;
    vector<ll> d(n - 1);
    for(int i = 0; i < n - 1; ++i) {
        cin >> d[i];
    }
    vector<ll> pos(n);
    for(int i = 0; i < n - 1; ++i) {
        pos[i + 1] = pos[i] + d[i];
    }
    vector<ll> x(m);
    for(int i = 0; i < m; ++i) {
        ll h, t;
        cin >> h >> t;
        x[i] = t - pos[h - 1];
    }
    sort(x.begin(), x.end());
    vector<ll> sum(m + 1);
    for(int i = 0; i < m; ++i) {
        sum[i + 1] = sum[i] + x[i];
    }

    vector<vector<ll>> dp(p + 1, vector<ll>(m, INF));
    vector<convex_hull_trick> cht(p + 1);
    for(int i = 0; i < m; ++i) {
        dp[1][i] = i * x[i] - sum[i];
        cht[1].add_f(-i, dp[1][i] + sum[i + 1]);
    }
    for(int i = 2; i <= p; ++i) {
        for(int j = 0; j < m; ++j) {
            if(j == 0) {
                dp[i][j] = 0;
            } else {
                dp[i][j] = cht[i - 1].min(x[j]) + j * x[j] - sum[j + 1];
            }
            cht[i].add_f(-j, dp[i][j] + sum[j + 1]);
        }
    }
    cout << dp[p][m - 1] << endl;
}

感想

ConvexHullTrickを初めて書いたので大変だった.かなり典型問題なので勉強になりました.

AOJ 2606 Perm Query

解法

まず,制約を読むと,ある項についてたかだか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

解法

この問題を解く上でポイントとなるのは,大気を循環させたときに,(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もあるかと言われると微妙….