mayoko’s diary

プロコンとかいろいろ。

AOJ 2397 Three-way Branch

AOJ
問題

Three-way Branch | Aizu Online Judge

W*H のグリッドが与えられる。これらのうち, N 個のセルは通れなくなっている。あるセル(x, y) からは (x-1, y+1) or (x, y+1) or (x+1, y+1) にのみ行ける場合, (0, 0) から (W-1, H-1) までの経路は何通りあるか。

解法

W がまぁまぁ小さくて H が超でかいので行列累乗な感じがします。

障害物は N 個しかないので,

  • 障害物のある行になるまでは行列累乗
  • 障害物ある行では 対応するところの場合の数を 0 にする

というようにやれば O(W^3 log H * N) でできます。…結構ギリギリですが。

typedef long long number;
typedef vector<number> vec;
typedef vector<vec> matrix;

const int MOD = 1e9+9;

// O( n )
matrix Identity(int n) {
    matrix A(n, vec(n));
    for (int i = 0; i < n; ++i) A[i][i] = 1;
    return A;
}
// O( n^3 )
matrix mul(const matrix &A, const matrix &B) {
    matrix C(A.size(), vec(B[0].size()));
    for (int i = 0; i < (int)C.size(); ++i)
        for (int j = 0; j < (int)C[i].size(); ++j)
            for (int k = 0; k < (int)A[i].size(); ++k) {
                C[i][j] += A[i][k] * B[k][j];
                C[i][j] %= MOD;
            }
    return C;
}
// O( n^3 log e )
matrix pow(const matrix &A, ll e) {
    if (e == 0) return Identity(A.size());
    if (e == 1) return A;
    if (e % 2 == 0) {
        matrix tmp = pow(A, e/2);
        return mul(tmp, tmp);
    } else {
        matrix tmp = pow(A, e-1);
        return mul(A, tmp);
    }
}
// O(n)
number tr(const matrix &A) {
    number ans = 0;
    for (int i = 0; i < (int)A.size(); ++i)
        ans += A[i][i];
    return ans;
}
// O( n )
number inner_product(const vec &a, const vec &b) {
    number ans = 0;
    for (int i = 0; i < (int)a.size(); ++i)
        ans += a[i] * b[i];
    return ans;
}
// O( n^2 )
vec mul(const matrix &A, const vec &x) {
    vec y(A.size());
    for (int i = 0; i < (int)A.size(); ++i)
        for (int j = 0; j < (int)A[0].size(); ++j)
            (y[i] += (A[i][j] * x[j] % MOD)) %= MOD;
    return y;
}


const int MAXN = 33;
int W, N;
ll H;

ll solve() {
	vector<pair<ll, int> > P(N);
	for (int i = 0; i < N; i++) {
		cin >> P[i].second >> P[i].first;
		P[i].first--; P[i].second--;
	}
	sort(P.begin(), P.end());
	vec v(W);
	v[0] = 1;
	ll now = 0;
	for (int i = 0; i < N; ) {
		matrix mat(W, vec(W));
		for (int j = 0; j < W; j++) {
			if (j-1 >= 0) mat[j][j-1] = 1;
			mat[j][j] = 1;
			if (j+1 < W) mat[j][j+1] = 1;
		}
		v = mul(pow(mat, P[i].first-now), v);
		now = P[i].first;
		int ni = i;
		while (ni < N && P[i].first == P[ni].first)
			ni++;
		for (int j = i; j < ni; j++)
			v[P[j].second] = 0;
		i = ni;
	}
	matrix mat(W, vec(W));
	for (int i = 0; i < W; i++) {
		if (i-1 >= 0) mat[i][i-1] = 1;
		mat[i][i] = 1;
		if (i+1 < W) mat[i][i+1] = 1;
	}
	v = mul(pow(mat, H-1-now), v);
	return v[W-1];
}

int main() {
	cin.tie(0);
	ios::sync_with_stdio(false);
	for (int t = 1; ; t++) {
		cin >> W >> H >> N;
		if (W == 0) break;
		cout << "Case " << t << ": " << solve() << endl;
	}
	return 0;
}

Educational Codeforces Round 16 D. Two Arithmetic Progressions

CodeForces
問題

codeforces.com

x = a1 * k + b1 = a2 * l + b2, L <= x <= R を満たす数 x がいくつあるか求めよ。

解法

k, l の組を一つ求められれば, 後は x が lcm(a1, a2) の周期で解になるので簡単です(実装が楽とは言っていない)。

a1 * k + b1 = a2 * l + b2 を変形すると
a1 * k - a2 * l = b2 - b1 となります。

これの解は, 拡張ユークリッドの互除法を使えば求められます。これを使うと
a1 * x + a2 * y = g (g は a1, a2 の最大公約数) という形の x, y, g を求めることができ, 全体を (b2-b1)/g 倍すれば k = x, -l = y となっている, というわけです。

後はめんどくさいですがゴニョゴニョします。
下のコードでは, この後

  • k, l を 0 以上に調整する
  • x が L 以上になるように調整する
  • L <= x <= R となるものの数を求める

というようにやりました。全部二分探索です。

なんか素直に実装すると 10^18 まわりがめんどくさそうなので python で実装しました。
python ってなんかカッコよくかかないといけない感じがしますが下のコードはダサダサです。

def extgcd(a, b):
    x, y = 0, 0
    d = a;
    if b != 0:
        d, y, x = extgcd(b, a%b)
        y -= (a//b) * x
    else:
        x, y = 1, 0
    return (d, x, y)

def main():
    a1, b1, a2, b2, L, R = map(int, input().split())
    # とりあえず k, l の候補を出す
    g, k, l = extgcd(a1, a2);
    b = b2-b1;
    if (b%g != 0):
        print (0)
        return
    k *= b//g
    l *= -b//g
    # k >= 0, l >= 0 を満たすようにする
    low = -2**100
    high = 2**100
    while high-low > 1:
        med = (low+high)//2
        tk = k+med*a2//g
        tl = l+med*a1//g
        if (tk >= 0 and tl >= 0):
            high = med
        else:
            low = med
    k = k+high*a2//g
    # x >= L を満たすようにする
    x = a1*k+b1
    low = -1
    high = 2**100
    lcm = a1*a2//g
    while high - low > 1:
        med = (low+high)//2
        tx = x+med*lcm
        if tx >= L:
            high = med
        else:
            low = med
    x = x+high*lcm
    # x <= R を満たすものの数を求める
    low = 0
    high = 2**100
    while high-low > 1:
        med = (low+high)//2
        tx = x+med*lcm
        if (tx <= R):
            low = med
        else:
            high = med
    if low == 0 and x > R:
        print (0)
        return
    print (low+1)
    return

if __name__=="__main__":
    main()

Codeforces Round #369 (Div. 2) E. ZS and The Birthday Paradox

CodeForces
問題

codeforces.com

整数 n, k が与えられる。

p = 2^n として, 1 - (p! / p^k * (p-k)!) の分母と分子を求めよ。ただし分母と分子は非常に大きな値になることがあるので, 10^6+3 で割ったあまりを求める。
(本当は p < k の場合は場合分けしないといけない)

解法

上では p! / p^k * (p-k)! と書きましたが, p*(p-1)*...*(p-k+1) / p^k と書く方が考えやすいです。

まず分母のほうから考えましょう。p^k を考えるだけなら簡単に出来ます。ただ, 分子で割られるので, そこがちょっと注意です。p が必ず偶数であることに注意すると, 2 で割られる回数は, p の分の n 回と, (p-1) * ... * (p-k+1) の分(これは k-1 にいくつの 2 の成分が含まれているかで求められる)です。

次に分子です。最終的に 10^6+3 で割ったあまりを求めるわけなので, k > 10^6+3 である場合は, 上記の掛け算の結果は 0 です。そうでない場合は, 直接計算すれば良いでしょう。2^m の成分がありますが上と同じように k-1 に含まれる 2 の成分の数分だけ割り算すれば良いです。

const int MOD = 1e6+3;

ll mod_pow(ll x, ll p, ll MOD) {
    ll a = 1;
    while (p) {
        if (p%2) a = a*x%MOD;
        x = x*x%MOD;
        p/=2;
    }
    return a;
}

// mod_inverse
ll mod_inverse(ll a, ll m) {
    return mod_pow(a, m-2, m);
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    ll n, k;
    cin >> n >> k;
    if (n <= 61 && 1ll<<n < k) {
        cout << "1 1" << endl;
        return 0;
    }
    // 分母を求める
    ll mother = 1;
    // 減らされる分母の数
    ll cnt = n;
    {
        ll K = k-1;
        while (K) {
            cnt += K/2;
            K /= 2;
        }
        mother = mod_pow(2, n, MOD);
        mother = mod_pow(mother, k, MOD);
        mother *= mod_pow(mod_inverse(2, MOD), cnt, MOD);
        mother %= MOD;
    }
    // 分子を求める
    ll child = 1;
    {
        if (k >= MOD) child = 0;
        else {
            ll first = mod_pow(2, n, MOD);
            for (int i = 0; i < k; i++) {
                (child *= first-i) %= MOD;
            }
            if (child < MOD) child += MOD;
        }
        // 2 で割りまくる作業をする
        (child *= mod_pow(mod_inverse(2, MOD), cnt, MOD)) %= MOD;
    }
    child = (mother-child+MOD)%MOD;
    cout << child << " " << mother << endl;
    return 0;
}

Codeforces Round #369 (Div. 2) Directed Roads

問題

codeforces.com

n 個の頂点からなるグラフが与えられる。各頂点 i から, A[i] への有向グラフが一本ずつ生えている。

このグラフの辺の向きをいくつか入れ替えることにより, このグラフにサイクルができないようにしたい(つまり x0 -> x1 -> ... -> xm -> x0 みたいな辺の向きができないようにしたい)。このような辺の向きの場合の数はいくつあるか。

解法

n 頂点, n 辺なので, グラフが連結ならこの中に閉路は一つだけあります。この閉路のなかで辺の向きがすべて同じ(このようなものは 2 通りある)時上記のサイクルができてしまうので, サイクルを構成する頂点数が cnt であるとして, 2^n - 2*2^(n-cnt) 通りが答えになります。

グラフが連結でない場合は これらを掛け算すれば良いです。

const int MAXN = 200200;
const int MOD = 1e9+7;
int A[MAXN];

vi G[MAXN];
int degree[MAXN];
bool done[MAXN];

void dfs(int v, int& num, int& cnt) {
    done[v] = true;
    num++;
    if (degree[v]) cnt++;
    for (int ch : G[v]) {
        if (!done[ch]) dfs(ch, num, cnt);
    }
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    int n;
    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> A[i];
        A[i]--;
        G[i].push_back(A[i]);
        G[A[i]].push_back(i);
    }
    queue<int> que;
    for (int i = 0; i < n; i++) {
        degree[i] = G[i].size();
        if (degree[i] == 1) {
            que.push(i);
            degree[i] = 0;
        }
    }
    while (!que.empty()) {
        int now = que.front();
        que.pop();
        for (int v : G[now]) {
            if (degree[v] == 0) continue;
            if (--degree[v] == 1) {
                degree[v] = 0;
                que.push(v);
            }
        }
    }
    ll ans = 1;
    for (int i = 0; i < n; i++) {
        if (!done[i]) {
            int num = 0, cnt = 0;
            dfs(i, num, cnt);
            //cout << i << " " << num << " " << cnt << endl;
            ll tmp = 1;
            for (int j = 0; j < num; j++) {
                tmp = (2*tmp)%MOD;
            }
            ll minus = 1;
            for (int j = 0; j < num-cnt+1; j++) {
                minus = (2*minus)%MOD;
            }
            tmp = (tmp-minus+MOD)%MOD;
            (ans *= tmp) %= MOD;
        }
    }
    cout << ans << endl;
    return 0;
}

AtCoder Regular Contest 060 F - 最良表現 / Best Representation

AtCoder
問題

arc060.contest.atcoder.jp

解法

n = |w| とします。

全て同じ文字で構成されている場合答えが n 1 となること, および w が良い文字列である場合答えが 1 1 になることは明らかです。

その他の場合は最良表現が 2 になることが分かるようです。

よって, 以下のようなアルゴリズムが考えられます：

  • 答えが n 1 になる場合は自明なので先に処理する。
  • 答えが 1 1 になるかどうかを判定する。
  • そうでない場合, [0, i) と [i, n) に分けてそれぞれが良い文字列かどうかを判定し, その数を求める。

ここで問題になるのは「良い文字列になるかの判定をどうやるか」ということです。文字列の長さが m であるとすると, 周期的になる文字の長さ d は, m の約数しかありません。

よって, 各 d に対してその周期の繰り返しでないことが分かれば良いわけですが, これは [0, n-d) の範囲の文字列と [d, n) の文字列が一致しているかどうかで判定することができます(最初の d 文字が一致していると 最初の 2 周期が同じであることがわかって, 次の d 文字が一致していると 3 周期一致していることがわかって, … というような感じ)。
文字が同じかどうかはローリングハッシュを使えば簡単に判定できます。ローリングハッシュって神。

また, [1, n] の約数の数を数えようとするアルゴリズムを考えると, 以下のプログラムより 約数の数の総和は O(N log N) であることがわかります。

for (i = 1; i <= n; i++) {
    for (j = i; j <= n; j += i) {
        cnt++;
    }
}

よって, 最良表現が 2 の場合でも, 全体で O(N log N) で判定できます。

std::random_device rnd;
const ull mul[2] = {10000 + (rnd()%10000), 10000 + (rnd()%10000)};
const ull mod[2] = {1000000007, 1000000009};

class RollingHash {
public:
    RollingHash() {}
    template<class T>
    RollingHash(const T& s) {
        init(s);
    }
    template<class T>
    void init(const T& s) {
        n = s.size();
        for (int i = 0; i < 2; i++) {
            pow[i].resize(n+1);
            h[i].resize(n+1);
            pow[i][0] = 1; h[i][0] = 0;
        }
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < n; j++) {
                pow[i][j+1] = (mul[i]*pow[i][j]) % mod[i];
                h[i][j+1] = (s[j] + h[i][j]*mul[i]) % mod[i];
            }
        }
    }
    inline pair<ull, ull> hash(int i) const {return make_pair(h[0][i], h[1][i]);}
    // [l, r)
    inline pair<ull, ull> hash(int l, int r) const {
        ull p0 = (h[0][r] - (h[0][l]*pow[0][r-l]%mod[0]) + mod[0]) % mod[0];
        ull p1 = (h[1][r] - (h[1][l]*pow[1][r-l]%mod[1]) + mod[1]) % mod[1];
        return make_pair(p0, p1);
    }
private:
    int n;
    vector<ull> pow[2], h[2];
};

const int MAXN = 500100;
vector<int> divs[MAXN];

// [l, r)
bool isGood(const RollingHash& rh, int l, int r) {
    int n = r-l;
    for (int d : divs[n]) {
        if (rh.hash(l, r-d) == rh.hash(l+d, r)) return false;
    }
    return true;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    string s;
    cin >> s;
    int n = s.size();
    {
        vector<int> cnt(26);
        for (char c : s)
            cnt[c-'a']++;
        for (int i = 0; i < 26; i++) {
            if (cnt[i] == n) {
                cout << n << endl;
                cout << 1 << endl;
                return 0;
            }
        }
    }
    for (int i = 1; i <= n; i++) {
        for (int j = 2*i; j <= n; j += i) {
            divs[j].push_back(i);
        }
    }
    RollingHash rh(s);
    if (isGood(rh, 0, n)) {
        cout << 1 << endl;
        cout << 1 << endl;
        return 0;
    }
    int ans = 0;
    for (int i = 1; i < n; i++) {
        if (isGood(rh, 0, i) && isGood(rh, i, n)) ans++;
    }
    cout << 2 << endl;
    cout << ans << endl;
    return 0;
}

良い文字列の判定の仕方が思いつかなくて死んでいた。言われてみれば当たり前なので悔しい。

AtCoder Regular Contest 060 D - 桁和 / Digit Sum

問題

arc060.contest.atcoder.jp

解法

本番はよくわからない解法で通したのですが, そちらも紹介したいと思います。

ある数 n に対して b を決めると, 次のような規則性があることに気付きます。

n/2+1 から n までは f(n/2+1, n) から始まって 1 ずつ f(b, n) の値が減っていきます。
n/3+1 から n/2 までは f(n/3+1, n) から始まって 2 ずつ f(b, n) の値が減っていきます。

...

n/(i+1)+1 から n/i までは f(n/(i+1), n) から始まって i ずつ f(b, n) の値が減っていきます。

この規則性通り調べることを考えると,

  • \sqrt{n} までは普通に調べる
  • それ以降は, 上記の規則性に従って考えると, n/i の i の候補は高々 \sqrt{n} 個しかないのでそれをすべて調べる

というようにやれば良いです。

想定解もやはり \sqrt{n} で区切って場合分けをします。 \sqrt{n} までは普通に調べ, そのあとは再帰の深さが 2 にしかならないことを利用します。b > \sqrt{n} の時, n = p * b + q とおけます。一方で, s = p + q です。ゆえに, n-s = p*(b-1) となり, b の候補は n-s の約数に絞られます。これらを列挙して, 実際に解になっているかを確かめれば良いです。

想定解法っぽいコード

ll calc(ll n, ll b) {
    if (n < b) return n;
    return calc(n/b, b) + n%b;
}

ll solve(ll n, ll s) {
    for (ll i = 2; i*i <= n; i++) {
        if (calc(n, i) == s) {
            return i;
        }
    }
    vector<ll> cand;
    for (ll i = 1; i*i <= n-s; i++) {
        if ((n-s)%i == 0) cand.push_back((n-s)/i+1);
    }
    sort(cand.begin(), cand.end());
    for (ll el : cand) {
        if (calc(n, el) == s) return el;
    }
    if (n == s) {
        return n+1;
    }
    return -1;
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    ll n, s;
    cin >> n >> s;
    cout << solve(n, s) << endl;
    return 0;
}

本番通したコード

ll calc(ll N, ll div) {
    if (N < div) return N;
    return N%div + calc(N/div, div);
}

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    ll N, S;
    cin >> N >> S;
    //for (int i = 2; i <= N; i++) {
    //    cout << i << " " << calc(N, i) << endl;
    //}
    ll last = -1;
    for (ll i = 2; i*i <= N; i++) {
        last = i;
        if (calc(N, i) == S) {
            cout << i << endl;
            return 0;
        }
    }
    for (ll i = last; i > 1; i--) {
        ll start = N/i+1;
        ll last = N/(i-1);
        ll tmp = calc(N, start);
        if (tmp >= S && (tmp-S)%(i-1) == 0) {
            ll o = (tmp-S)/(i-1);
            if (start+o > last) continue;
            cout << start+o << endl;
            return 0;
        }
    }
    if (S == N) cout << N+1 << endl;
    else if (S == 1) cout << N << endl;
    else cout << -1 << endl;
    return 0;
}

AtCoder Regular Contest 060 E - 高橋君とホテル / Tak and Hotels

AtCoder

なんだかんだ 3 完できてうれぴー。

問題

arc060.contest.atcoder.jp

解法

クエリでは a < b として良いです。

戦略としては, 「その日にたどり着けるところまでなるべく遠くへ行く」というのが明らかに最適です。ということで, nxt[0][i] = (頂点 i からなるべく遠くへ行こうとしたときたどり着く頂点) というのをあらかじめメモしておきます。これは二分探索を使えば各 i について O(log N) でできます。

で, nxt[j+1][i] = nxt[j][nxt[j][i]] というのを繰り返すことにより, i から 2^j 日貪欲に移動した際にたどり着ける最も遠い頂点を得ます。これがあれば各クエリに対して二分探索をして O(log N) で答えを求められます。

クエリでの二分探索の書き方は蟻本の LCA と同じです。b の直前になるところまで行って, 後 +1 みたいな感じ。

const int MAXN = 100100;
int N, L, Q;
int X[MAXN];

int nxt[20][MAXN];

int main() {
    cin.tie(0);
    ios::sync_with_stdio(false);
    cin >> N; 
    for (int i = 0; i < N; i++)
        cin >> X[i];
    cin >> L;
    memset(nxt, -1, sizeof(nxt));
    for (int i = 0; i < N; i++) {
        if (i == N-1) {
            nxt[0][i] = N-1;
            continue;
        }
        int low = i+1, high = N;
        while (high-low > 1) {
            int med = (high+low)/2;
            if (X[med] - X[i] <= L) low = med;
            else high = med;
        }
        nxt[0][i] = low;
    }
    for (int i = 0; i < 20-1; i++) {
        for (int j = 0; j < N; j++) {
            if (nxt[i][j] == -1) nxt[i+1][j] = -1;
            else nxt[i+1][j] = nxt[i][nxt[i][j]];
            assert(nxt[i][j] != -1);
        }
    }
    cin >> Q;
    for (int t = 0; t < Q; t++) {
        int a, b;
        cin >> a >> b;
        a--; b--;
        if (a > b) swap(a, b);
        if (a < b) {
            int ans = 0;
            for (int i = 20-1; i >= 0; i--) {
                if (nxt[i][a] < b) {
                    a = nxt[i][a];
                    ans |= 1<<i;
                }
            }
            cout << ans+1 << endl;
        }
    }
    return 0;
}