CODE FESTIVAL 2017 qual B: E - Popping Balls

,

https://beta.atcoder.jp/contests/code-festival-2017-qualb/tasks/code_festival_2017_qualb_e

感想

  • 本番解けなかったD問題より圧倒的に簡単。悲しい
  • 簡単に$O(AB)$にできるが、$O(AB^2)$で間に合ってしまったのでそのまま

solution

赤色のボールの残っている数と青色のボールの残っている数をそれぞれ$y$軸$x$軸として、$2$次元の格子上の経路数を数えるような感じのDP。適切に書けば$O(N^2)$。

$s \lt t$としてよい。 赤色のボールは(残っているなら)常に取れるので、$s, t$は青色のボールを取るためにある。 $t$は、始めて青色のボールを取り出す時に青色のボールの先頭であるような位置にするのが最良。 同様に$s$は、$t$番目の位置にボールが存在しなくなってから始めて青色のボールを取り出す時の位置とするのが最良。

そのようにしたときの経路数は表を$2$段階で埋めるようにすれば計算できる。

implementation

#include <cstdio>
#include <vector>
#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define repeat_from(i, m, n) for (int i = (m); (i) < int(n); ++(i))
using ll = long long;
using namespace std;

constexpr int mod = 1e9+7;
vector<vector<int> > calc_choose(int n) { // O(n^2)
    vector<vector<int> > dp(n + 1);
    dp[0].assign(1, 1);
    repeat (i, n) {
        dp[i + 1].resize(i + 2);
        repeat (j, i + 2) {
            auto & it = dp[i + 1][j];
            if (j - 1 >= 0) {
                it += dp[i][j - 1];
                if (it >= mod) it -= mod;
            }
            if (j != i + 1) {
                dp[i + 1][j] += dp[i][j];
                if (it >= mod) it -= mod;
            }
        }
    }
    return dp;
}

int main() {
    // input
    int a, b; scanf("%d%d", &a, &b);
    // solve
    auto choose = calc_choose(a + b + 3);
    ll result = 0;
    result += a + 1;
    repeat_from (x, 1, b + 1) {
        ll acc = 0;
        repeat_from (y, b - x, a + 1) {
            acc += choose[b - 1][x - 1];
            if (acc >= mod) acc -= mod;
            int nb = b - x;
            repeat (z, nb) {  // NOTE: これ消せる
                int ny = y + (nb - 1 - z);
                if (ny <= a) result += acc * choose[nb - 1][z] % mod;
            }
        }
    }
    result %= mod;
    // output
    printf("%lld\n", result);
    return 0;
}

CODE FESTIVAL 2017 qual B: D - 101 to 010

,

https://beta.atcoder.jp/contests/code-festival-2017-qualb/tasks/code_festival_2017_qualb_d

感想

解けず。連続する1を数えて数列に変形するのは必須だと思ってその方向ばかり考えていたが、そうだとすると実装が面倒になりすぎるので完全に失敗だった。

solution

連鎖的に操作をするとき111...11101あるいは10111...111の形を000...00010などにできる。 そのようなものだけ見ればよいので、111...1110110111...111の形の部分文字列を関係としてDP。 愚直にやると$O(N^2)$だが、形が制限されているので前回の端を覚えておくなどすれば$O(N)$。

implementation

#include <algorithm>
#include <iostream>
#include <vector>
#define repeat_reverse(i, n) for (int i = (n)-1; (i) >= 0; --(i))
#define whole(x) begin(x), end(x)
using namespace std;
template <class T> inline void setmax(T & a, T const & b) { a = max(a, b); }

int solve(int n, string const & s) {
    vector<int> dp(n + 1);
    int r = n;
    repeat_reverse (l, n) {
        dp[l] = dp[l + 1];
        if (l + 2 < n and s[l + 2] == '0') {
            r = l + 2;
        }
        if (s[l] == '1') {
            if (s[l + 1] == '0') {
                if (s[r - 1] == '1') {
                    setmax(dp[l], dp[r] + r - l - 2);
                    if (s[r - 2] == '1') {
                        setmax(dp[l], dp[r - 1] + (r - 1) - l - 2);
                    }
                }
            } else {
                if (r + 1 < n and s[r + 1] == '1') {
                    setmax(dp[l], dp[r + 2] + (r + 2) - l - 2);
                }
            }
        }
    }
    return *max_element(whole(dp));
}

int main() {
    int n; string s; cin >> n >> s;
    cout << solve(n, s) << endl;
    return 0;
}

CODE FESTIVAL 2017 qual B: C - 3 Steps

,

https://beta.atcoder.jp/contests/code-festival-2017-qualb/tasks/code_festival_2017_qualb_c

感想

直前にICPCの練習で二部グラフに上手く落ちる問題を解いていたのですぐだった。

solution

$G$が二部グラフなら完全二部グラフに、そうでなければ完全グラフにできる。よって二部グラフかどうか判定すればよい。$O(N + M)$。

implementation

#include <algorithm>
#include <cstdio>
#include <functional>
#include <vector>
#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define whole(x) begin(x), end(x)
using ll = long long;
using namespace std;

int check_bipartite_graph(vector<vector<int> > const & g) {
    int n = g.size();
    vector<char> used(n, -1);
    function<bool (int, int)> dfs = [&](int i, int parent) {
        for (int j : g[i]) {
            if (used[j] != -1) {
                if (used[j] == used[i]) {
                    return false;
                }
            } else {
                used[j] = used[i] ^ 1;
                if (not dfs(j, i)) return false;
            }
        }
        return true;
    };
    used[0] = 0;
    if (not dfs(0, -1)) return -1;
    return count(whole(used), 0);
}

int main() {
    // input
    int n, m; scanf("%d%d", &n, &m);
    vector<vector<int> > g(n);
    repeat (i, m) {
        int a, b; scanf("%d%d", &a, &b); -- a; -- b;
        g[a].push_back(b);
        g[b].push_back(a);
    }
    // solve
    int k = check_bipartite_graph(g);
    ll result = k == -1 ?
        n *(ll) (n - 1) / 2 - m :
        k *(ll) (n - k) - m;
    // output
    printf("%lld\n", result);
    return 0;
}

CODE FESTIVAL 2017 qual B: B - Problem Set

,

https://beta.atcoder.jp/contests/code-festival-2017-qualb/tasks/code_festival_2017_qualb_b

solution

多重集合として$T \subseteq D$か見ればよい。sortして一緒に先頭から削っていく。$O(N \log N)$。

implementation

#!/usr/bin/env python3
n = int(input())
d = list(map(int, input().split()))
m = int(input())
t = list(map(int, input().split()))
d.sort()
t.sort()
j = 0
for t_i in t:
    while j < n and d[j] < t_i:
        j += 1
    if j == n or d[j] != t_i:
        result = False
        break
    j += 1
else:
    result = True
print(['NO', 'YES'][result])


DISCO presents ディスカバリーチャンネル コードコンテスト2017 予選: D - 石

,

http://ddcc2017-qual.contest.atcoder.jp/tasks/ddcc2017_qual_d

感想

$H, W$は偶数というのに気付かないまま通してしまった。ショック。

solution

南北方向に対称にしてから東西方向に対称にするのと、その逆をそれぞれ試す。$O(HW)$。 盤面を全部舐めて対応するものを持たなければ削る感じで書くと楽。

implementation

#include <iostream>
#include <vector>
#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
using ll = long long;
using namespace std;

ll solve(int h, int w, int a, int b, vector<string> m) {
    ll result = 0;
    bool modified = false;
    repeat (y, h) {
        repeat (x, w) {
            if (m[y][x] == 'S' and m[y][x] != m[h - y - 1][x]) {
                m[y][x] = '.';
                modified = true;
            }
        }
    }
    if (modified) {
        result += a;
    }
    repeat (y, h) {
        repeat (x, w) {
            if (m[y][x] == 'S' and m[y][x] != m[y][w - x - 1]) {
                m[y][x] = '.';
                m[h - y - 1][x] = '.';
                result += a;
                modified = true;
            }
        }
    }
    if (modified) {
        result += b;
    }
    repeat (y, h) {
        repeat (x, w) {
            if (m[y][x] == 'S') {
                if (2 * y + 1 == h and 2 * x + 1 == w) {
                    m[y][x] = '.';
                    result += a + b;
                } else if (2 * y + 1 == h) {
                    m[y][x] = '.';
                    m[y][w - x - 1] = '.';
                    result += 2 * a + b;
                } else if (2 * x + 1 == w) {
                    m[y][x] = '.';
                    m[h - y - 1][x] = '.';
                    result += a + 2 * b;
                } else {
                    m[y][x] = '.';
                    m[y][w - x - 1] = '.';
                    m[h - y - 1][x] = '.';
                    m[h - y - 1][w - x - 1] = '.';
                    result += max(a, b) + a + b;
                }
            }
        }
    }
    return result;
}

int main() {
    // input
    int h, w, a, b; cin >> h >> w >> a >> b;
    vector<string> m(h); repeat (y, h) cin >> m[y];
    // solve
    ll p = solve(h, w, a, b, m);
    vector<string> t(w, string(h, '\0'));
    repeat (y, h) repeat (x, w) t[x][y] = m[y][x];
    ll q = solve(w, h, b, a, t);
    // output
    printf("%lld\n", max(p, q));
    return 0;
}

DISCO presents ディスカバリーチャンネル コードコンテスト2017 予選: C - 収納

,

http://ddcc2017-qual.contest.atcoder.jp/tasks/ddcc2017_qual_c

solution

貪欲ぽく。 残っているなかで最も大きいのと最も小さいのを対にしていく。$O(N)$。

implementation

#include <algorithm>
#include <cstdio>
#include <vector>
#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define whole(x) begin(x), end(x)
using namespace std;

int main() {
    // input
    int n, c; scanf("%d%d", &n, &c);
    vector<int> l(n); repeat (i, n) scanf("%d", &l[i]);
    // solve
    sort(whole(l));
    int cnt = 0;
    int i = 0, j = n - 1;
    while (i < j) {
        while (i < j and l[i] + l[j] + 1 <= c) {
            ++ i;
            -- j;
            ++ cnt;
        }
        -- j;
    }
    // output
    int result = n - cnt;
    printf("%d\n", result);
    return 0;
}

DISCO presents ディスカバリーチャンネル コードコンテスト2017 予選: B - 鉛筆

,

http://ddcc2017-qual.contest.atcoder.jp/tasks/ddcc2017_qual_b

感想

$f(x) = a_nx^n + a_{n-1}x^{n-1} + \dots + a_1x + a_0$の$x = x_0$での値を計算するようなとき、$f(x) = (( \dots (a_nx + a_{n-1})x + \dots) x + a_1) x + a_0$とすると計算量が落ちるというあれはHorner法というらしいのですが、それだなあと思いながらコードを書きました。

implementation

#!/usr/bin/env python3
a, b, c, d = map(int, input().split())
print(((a * 12 + b) * 12 + c) * 12 + d)


HackerRank Performance Optimization: D. Skipping Subpath Sum

,

https://www.hackerrank.com/contests/optimization-oct17/challenges/skipping-subpath-sum

problem

頂点に重みの付いた木が与えられる。 次のクエリにたくさん答えよ:

  • 頂点$u, v$が与えられる。$u, v$間の唯一のpathについて、偶数番目の頂点の重みの総和と奇数番目の頂点の重みの総和のうち大きい方を答えよ。

solution

根に向かって$1$歩ずつLCAまで登っていき、その過程を貯め込んでいい感じにKadane’s algorithm。$O(QN)$。

知見

  • 数列から総和が最大となるような連続部分列を探すアルゴリズムとして、Kadane’s algorithm
    • 順に伸ばしていって総和が負になったらリセットする感じ
    • 名前を知らないだけでみんなやってるやつ

implementation

...

#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define whole(x) begin(x), end(x)
template <class T> inline void setmax(T & a, T const & b) { a = max(a, b); }

vector<int> skippingSubpathSum(int n, vector<int> const & c, vector<vector<int> > const & graph, vector<pair<int, int> > const & queries) {
    // prepare
    constexpr int root = 0;
    vector<int16_t> parent(n, -1);
    vector<int16_t> depth(n, -1); {
        function<void (int)> go = [&](int i) {
            for (int j : graph[i]) if (parent[j] == -1) {
                parent[j] = i;
                depth[j] = depth[i] + 1;
                go(j);
            }
        };
        parent[root] = root;
        depth[root] = 0;
        go(root);
    }
    // solve
    vector<int> answers;
    for (auto query : queries) {
        // prepare Kadane's algorithm
        int result[2] = {};
        int acc[2] = {};
        int path_size = 0;
        auto kadane = [&](int c_i) {
            int i = (path_size ++) & 1;
            acc[i] += c_i;
            setmax(acc[i], 0);
            setmax(result[i], acc[i]);
        };
        // run
        int u, v; tie(u, v) = query;
        if (depth[u] < depth[v]) {
            swap(u, v);
        }
        assert (depth[u] >= depth[v]);
        for (int depth_u = depth[u]; depth_u > depth[v]; ) {
            kadane(c[u]);
            u = parent[u];
            -- depth_u;
        }
        assert (depth[u] == depth[v]);
        vector<int> stk;
        while (u != v) {
            kadane(c[u]);
            stk.push_back(c[v]);
            u = parent[u];
            v = parent[v];
        }
        kadane(c[u]);
        while (not stk.empty()) {
            kadane(stk.back());
            stk.pop_back();
        }
        // answer
        int answer = max(result[0], result[1]);
        answers.push_back(answer);
    }
    return answers;
}

...

HackerRank Performance Optimization: C. Keywords

,

https://www.hackerrank.com/contests/optimization-oct17/challenges/keywords

problem

単語列$s$と単語の集合$\mathrm{keys}$が与えられる。$s$の連続部分列で$\mathrm{keys}$を単語として全て含むものを探し、そのようなものの文字列としての長さの最小値を答えよ。

implementation

...

#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define repeat_reverse(i, n) for (int i = (n)-1; (i) >= 0; --(i))
template <class T> inline void setmax(T & a, T const & b) { a = max(a, b); }
template <class T> inline void setmin(T & a, T const & b) { a = min(a, b); }

constexpr int inf = 1e9+7;
int minimumLength(string const & text, vector<string> const & keys) {
    int n = text.length();
    vector<uint16_t> found(n + 1);
    for (int l = 0; l < n; ) {
        if (text[l] == ' ') {
            ++ l;
        } else {
            int r = text.find(' ', l);
            if (r == string::npos) r = n;
            repeat (k, keys.size()) {
                if (text.compare(l, r - l, keys[k]) == 0) {
                    found[r] |= 1 << k;
                }
            }
            l = r;
        }
    }
    vector<vector<int> > next(keys.size(), vector<int>(n + 30, inf));
    repeat (k, keys.size()) {
        if (found[n] & (1 << k)) next[k][n] = n;
    }
    repeat_reverse (i, n) {
        repeat (k, keys.size()) {
            next[k][i] = next[k][i + 1];
            if (found[i] & (1 << k)) next[k][i] = i;
        }
    }
    int result = inf;
    repeat (i, n) {
        int acc = 0;
        repeat (k, keys.size()) {
            setmax(acc, next[k][i + keys[k].length()]);
        }
        setmin(result, acc - i);
    }
    if (result > n) result = -1;
    return result;
}

...

HackerRank Performance Optimization: B. Collatz Sequence Sum

,

https://www.hackerrank.com/contests/optimization-oct17/challenges/collatz-sequence-sum

problem

$g(K)$を$K$以下の自然数でCollatz数列が最長となるようなもの。複数あるなら最大。 $A, B$で生成される数列$N_1, \dots, N_T$に対し$\sum g(N_i)$を答えよ。

solution

いい感じにメモ化する。mapを入れてだいたい$O(N \log N)$とかなのでは。

implementation

...

#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
int collatzSequenceSum(int T, int A, int B) {
    constexpr int mod = 5003;
    vector<int> g(mod); {
        map<int, int> memo;
        memo[0] = 0;
        memo[1] = 1;
        function<int (int)> collatzSequenceLen = [&](int n) {
            if (not memo.count(n)) {
                memo[n] = n % 2 == 0 ?
                    1 + collatzSequenceLen(n / 2) :
                    1 + collatzSequenceLen(3 * n + 1);
            }
            return memo[n];
        };
        repeat (i, mod) {
            if (i >= 1) {
                g[i] = collatzSequenceLen(i) >= collatzSequenceLen(g[i - 1]) ?
                    i :
                    g[i - 1];
            }
        }
    }
    int n = 0;
    int result = 0;
    while (T --) {
        n = (A * n + B) % mod;
        result += g[n];
    }
    return result;
}

...

HackerRank Performance Optimization: A. Castle Towers

,

https://www.hackerrank.com/contests/optimization-oct17/challenges/castle-towers

感想

Performance Optimizationというから実行速度で点数あるいはペナルティが付くと思っていましたが、単にTLEコードがtextareaに始めから書いてあるというだけでした。 気付かずACした後に最適化して再提出みたいな真似をしてしまった。

problem

高さが最大のものの数を答えよ。

solution

sortで$O(N \log N)$だと遅いので、$O(N)$にする。

implementation

...

#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
int castleTowers(int n, vector<int> const & ar) {
    int maxi = 0;
    int cnt = 0;
    repeat (i, n) {
        if (maxi < ar[i]) {
            maxi = ar[i];
            cnt = 1;
        } else {
            cnt += (maxi == ar[i]);
        }
    }
    return cnt;
}

...

Tenka1 Programmer Contest: E - CARtesian Coodinate

,

http://tenka1-2017.contest.atcoder.jp/tasks/tenka1_2017_e

solution

$x$軸と$y$軸は独立。$x$座標$y$座標のそれぞれについて交点の座標の中央値を取ればよい。$O(N^2)$。定数倍最適化。

知見

  • #pragma clang loop vectorize(enable) は、付けなくてもvectorizeされているためあまり速度は変わらないだが、
  • vector<T>.push_back をloop内でするとvectorizeできないので、flagを作って外に出すテク

  • tanakhさんの提出をとても参考にした。ほとんど同じと言ってよい。
  • $23$WA $2$RE
  • 速度は足りてたが普通にバグ

implementation

#include <algorithm>
#include <cassert>
#include <cstdio>
#include <random>
#include <vector>
#define repeat(i, n) for (int i = 0; (i) < int(n); ++(i))
#define whole(x) begin(x), end(x)
using namespace std;

int main() {
    // input
    int n; scanf("%d", &n);
    vector<int16_t> a(n), b(n), c(n);
    repeat (i, n) scanf("%hd%hd%hd", &a[i], &b[i], &c[i]);

    // solve

    // // find the windows
    constexpr int sample_size = 500000;
    constexpr int window_radius = 1000;
    vector<double> samples_x(sample_size);
    vector<double> samples_y(sample_size);
    mt19937 gen;
    repeat (iteration, sample_size) {
        int i, j;
        do {
            i = uniform_int_distribution<int>(0, n - 1)(gen);
            j = uniform_int_distribution<int>(0, n - 1)(gen);
        } while (i == j);
        double den = a[i] * b[j] - a[j] * b[i];
        double x = (b[j] * c[i] - b[i] * c[j]) / den;
        double y = (a[j] * c[i] - a[i] * c[j]) / - den;
        samples_x[iteration] = x;
        samples_y[iteration] = y;
    }
    sort(whole(samples_x));
    sort(whole(samples_y));
    double dlx = samples_x[sample_size / 2 - window_radius];
    double drx = samples_x[sample_size / 2 + window_radius];
    double dly = samples_y[sample_size / 2 - window_radius];
    double dry = samples_y[sample_size / 2 + window_radius];
    bool determined_x = (dlx == drx);
    bool determined_y = (dly == dry);
    float lx = dlx;
    float rx = drx;
    float ly = dly;
    float ry = dry;

    // // get points in them
    int low_count_x = 0;
    int low_count_y = 0;
    vector<double> bucket_x;
    vector<double> bucket_y;
    vector<int> flag(n);
    repeat (j, n) {
        {
            int aj = a[j];
            int bj = b[j];
            int cj = c[j];
            repeat (i, j) {
                float den = a[i] * bj - aj * b[i];
                float x = (bj * c[i] - b[i] * cj) / den;
                float y = (aj * c[i] - a[i] * cj) / - den;
                low_count_x += (x < lx);
                low_count_y += (y < ly);
                flag[i] = (int(lx <= x and x <= rx) << 1) | int(ly <= y and y <= ry);
            }
        }
        {
            double aj = a[j];
            double bj = b[j];
            double cj = c[j];
            if (not determined_x) repeat (i, j) if (flag[i] & 2) {
                double x = (bj * c[i] - b[i] * cj) / (a[i] * bj - aj * b[i]);
                bucket_x.push_back(x);
            }
            if (not determined_y) repeat (i, j) if (flag[i] & 1) {
                double y = (aj * c[i] - a[i] * cj) / (aj * b[i] - a[i] * bj);
                bucket_y.push_back(y);
            }
        }
    }

    // // make results
    int half = (n * (n - 1) / 2 - 1) / 2;
    double x;
    if (determined_x) {
        x = dlx;
    } else {
        int ix = half - low_count_x;
        assert (0 <= ix and ix < bucket_x.size());
        nth_element(bucket_x.begin(), bucket_x.begin() + ix, bucket_x.end());
        x = bucket_x[ix];
    }
    double y;
    if (determined_y) {
        y = dly;
    } else {
        int iy = half - low_count_y;
        assert (0 <= iy and iy < bucket_y.size());
        nth_element(bucket_y.begin(), bucket_y.begin() + iy, bucket_y.end());
        y = bucket_y[iy];
    }

    // output
    printf("%.15lf %.15lf\n", x, y);
    return 0;
}

ISUCON7 オンライン予選: 参加記

,

http://isucon.net/archives/50949022.html

概要

やったこと

repoのcommit logを見て

反省

  • /icons/*と帯域の問題に関してまったく気付かなかった
    • せめてCPU利用率見てれば気付けたのでは
    • あるいは$3$台全部をweb鯖にするの試していればよかった
    • 改善を期待して行った施策が多数変化なしで原因不明なのだから、全体を律速する部分がどこかにあることを疑うべきだった
    • レギュレーションにもヒントたくさんだったのに
    • 飛んできているHTTP header見るのもしてない

感想

当初の目的は様子見(初参加なので)だったので達成してはいるが、悔しいね


AtCoder Regular Contest 055: B - せんべい

,

http://arc055.contest.atcoder.jp/tasks/arc055_b

hamkoさんが苦しんでて楽しそうだった: https://docs.google.com/document/d/19z-mgAde_D_Xtom_dShy4ggrhzSWxGnNnoZ9T4jT6tU/edit。なので解いた。 非独立である部分に気付かず私も苦しんだ。kmjpさんの解説を見た。

solution

確率DP。$\mathrm{dp} : (N + 1) \times (K + 1) \to [0, 1]$。$O(NK)$。

既に$i$個見て$j$個食べまだ$N$番目のせんべいは現れていないとし、その時点から最適に動いたときの$N$を食べられる確率を$\mathrm{dp}(i, j)$と定義する。 $i’ \gt i$について既知と仮定し$\mathrm{dp}(i, j)$を考えよう。 $i = N$なら確率は$0$。 そうでないと仮定し$i$番目のせんべいを確認したとする。 新たに確認したせんべいが過去最高でなければ、明らかに食べる必要はない。 過去最高であれば、食べた場合と食べなかった場合で良い方を選べばよい。 食べた場合は$N$である確率と$N$でなくそれ以降で$N$を食べる確率の和、食べなかった場合は$N$でなくそれ以降で$N$を食べる確率の和がそれ。 ここまでに嘘はないが罠はあり、$N$である確率は過去最高である確率とは独立でない。 事象は次に$3$種に分類される。これに注意して書けば通る。

  • $N$である (自明に過去最高)
  • $N$ではないが過去最高
  • 過去最高でない (自明に$N$でない)

implementation

#include <cmath>
#include <cstdio>
#include <functional>
#include <vector>
using namespace std;
template <typename X, typename T> auto vectors(X x, T a) { return vector<T>(x, a); }
template <typename X, typename Y, typename Z, typename... Zs> auto vectors(X x, Y y, Z z, Zs... zs) { auto cont = vectors(y, z, zs...); return vector<decltype(cont)>(x, cont); }

int main() {
    int n, k; scanf("%d%d", &n, &k);
    auto memo = vectors(n + 1, k + 1, NAN);
    function<double (int, int)> dp = [&](int i, int j) {
        if (not isnan(memo[i][j])) return memo[i][j];
        if (i == n or j == k) return memo[i][j] = 0;
        double is_n = 1.0 / (n - i);
        double fake = (1 - is_n) * (1.0 / (i + 1.0));
        double top = is_n + fake;
        return memo[i][j] =
            (1 - top) * dp(i + 1, j)
            + max(is_n + fake * dp(i + 1, j + 1),
                         fake * dp(i + 1, j));
    };
    printf("%.12lf\n", dp(0, 0));
    return 0;
}