ベスパリブ

プログラミングを主とした日記・備忘録です。ベスパ持ってないです。

2022年ですけど数独ソルバーを作りました。

作ったもの

taketakeyyy.github.io

動機

WEBで適当に「数独 ソルバー」と検索すると先人たちのソルバーがヒットするのですが、ブラウザ上で動作するものをいくつか使ってみると(いわゆる)「仮置き」に対応していなくて途中で計算を打ち切ったり、複数解持つ数独は解かない仕様となっていました(人が楽しむ上で数独の答えが複数解あるものはあまりよろしくないと思いますが、数独の定義に「唯一解である」はない……と思う)。 インストール形式の数独ソルバーもいくつかあったのですが、ちょっとした数独の確認のためにインストールするのもうーん……だったので、以下の要望を満たすソルバーを作りました。

  • 仮置きに対応してちゃんと最後まで解く
  • 複数解ある場合でも最後まで解いてひとつの解を出力する
  • ブラウザ上で動く

数独ソルバーの作り方(雑)

DFS(深さ優先探索)が書けるならソルバーを書けます。

以下、雑なアルゴリズムです。

  • 【ステップ1】:マスに入る数字の候補が1つしかない場合、その数字で確定させる
  • 【ステップ2】:マスに入る数字の候補が複数ある場合でも、どの縦、横、3x3マス内にもない、自分しか持っていない候補の数字を持っているなら、その数字で確定させる(たとえば自分のマスの候補が「1,2,3」で、縦、横、3x3のどのマスも3を候補に持っていないなら、自分は3で確定)
  • 【ステップ3】:確定するマスがなくなるまで【ステップ1】に戻る。
  • 【ステップ4】:確定するマスがなくなり、かつまだ全マス埋まっていないなら、最も候補数字の少ないマスをひとつ見つける。そのマスの候補数字の1つを「仮置き」して、【ステップ1】へDFSで再帰する。仮置きした数字で全てのマスが埋まらないなら、他の候補数字を仮置きして同様に再帰させる。
  • 【終了条件】:【ステップ3】または【ステップ4】の過程で全てのマスに数字が埋まったら数独が解けた。埋まらなかったら数独は解けなかった。

DFSで深い探索になるとStack Overflowになる可能性があるので、「仮置き」する前に確定できるマスがあるなら埋める、というのが良いソルバーのコツでしょう。調べた感じX-wingなど様々な方法がありましたが、今の所基本的な【ステップ1】と【ステップ2】で解けています。

一番難しかったのはJavaScriptでした

TypeScriptで数独ソルバーを書いたのですが、数独ソルバーのアルゴリズム云々よりも、ブラウザで動作させる以上非同期処理を逐次処理に置き換える作業が必要になってきます。VSCodeデバッグのやり方もようわからんし、配列をコピーしたいだけなのにJSON.parse(JSON.stringify())しろとか言われるし、Setのdeepcopyはそれに対応してないから自分で作らないといけなかったり、いわんやArray<Array<Set<number>>>の場合をや……など、JavaScriptはDOM操作をする言語であってデータ構造をいじる言語ではないと思いました。次にソルバーを書くならRustで書いてwasmで処理させると思います。

おわりに

Arto Inkala氏が作成した世界で最も難しい数独や、Hiroshi Watanabeスパコンで作成した数独や、全空白の数独なども解ける数独ソルバーを作りました。

「この問題解けるはずなのに解けてないよ~」等の不具合があれば、数独ソルバーサイトの「不具合報告」からお願いします。

ABC233: E - Σ[k=0..10^100]floor(X/10^k) を解いた

問題:AtCoder Beginner Contest 233: E - Σ[k=0..10100]floor(X/10k)

解説

 \sum_{k=0}^{10^{100}} \lfloor \frac{X}{10^k} \rfloor

を計算せよという問題。

Xの制約が

 1 <= X < 10^{500000}

と大きいし、k=0から10100まで計算は普通にしていたら制限時間に間に合わない。

Xの桁数が最大500000なので、問題をエスパーすれば計算量O(桁数)で解くのかなぁと当たりはつくが、まだどうすればいいのかわからない。


とりあえず紙に書いてみる。

とりえあず書いてみる

するとなんか法則性が見えてくる。

特にX=98764のとき、98764 + 9876 + 987 + 98 + 9 + 0 + ... と一桁ずつ減っていくのを足し算する形になっていることに気づけば、他のパターンも同様になっていることに気づく。


他に思いつかないし、じゃあこれを愚直にプログラムに落とし込んで書いたコードが以下の形。(TLEする)

配列ansは、答えの各桁の数字を格納する。

void __solve() {
    // TLE lol
    string X; cin >> X;
    const ll MAX_ANS = 500010;
    vector<ll> ans(MAX_ANS);

    while(X.size()!=0) {  // 最大 (桁数) 回
        string tmpX = X;
        ll i = 0;
        while(tmpX.size()!=0) {  // 平均 (桁数+1)/2 回
            ans[i] += tmpX[tmpX.size()-1] - '0';
            if (ans[i] >= 10) {
                ans[i+1] += ans[i]/10;
                ans[i] = ans[i] % 10;
            }
            tmpX = tmpX.substr(0, tmpX.size()-1);
            i++;
        }
        X = X.substr(0, X.size()-1);
    }
    for(ll i=0; i<MAX_ANS; i++) {
        if (ans[i] >= 10) {
            ans[i+1] += ans[i]/10;
            ans[i] = ans[i] % 10;
        }
    }


    // 出力
    bool is_still_zero = true;
    for(ll i=MAX_ANS-1; i>=0; i--) {
        if (is_still_zero && ans[i]==0) continue;
        if (ans[i] != 0) is_still_zero = false;
        printf("%lld", ans[i]);
    }
    printf("\n");
}

書いたあとから気づくのだが、最初のwhileループが最大 桁数 回、二個目のwhileループが平均 (桁数+1)/2 回実行されるので、だいたい計算量O((桁数)2)程度となり、これはTLEする。

コンテスト中はTLEして終了。


寝て起きたら解法が思いつく。こんなことってあるんですね。

X=987654のときのものを、筆算の形で書いてみたのが以下の形。

筆算の形で書いてみる

各桁が累積和になっていることに気づく。

累積和はO(桁数)で計算できる。

各桁の数字の計算も、累積和を使えばO(桁数)で計算できる。

よって線形時間で解けた。


この「筆算の形で考えると解法が見えてくる」というのは、ABC231: E - Minimal paymentsでも見たので、結構よくあるテクニックかもしれない。

コーナーケース

  • 特になし

実装

C++で実装した。

#define _USE_MATH_DEFINES  // M_PI等のフラグ
#include <bits/stdc++.h>
#define MOD 1000000007
#define COUNTOF(array) (sizeof(array)/sizeof(array[0]))
#define rep(i,n) for (int i = 0; i < (n); ++i)
#define intceil(a,b) ((a+(b-1))/b)
using namespace std;
using ll = long long;
using pii = pair<int,int>;
using pll = pair<long,long>;
const long long INF = LONG_LONG_MAX - 1001001001001001;
void chmax(int& x, int y) { x = max(x,y); }
void chmin(int& x, int y) { x = min(x,y); }


void solve() {
    string X; cin >> X;
    const ll MAX_ANS = 500010;
    vector<ll> ans(MAX_ANS);

    // 累積和作成
    vector<ll> r(X.size()+1, 0);
    for(ll i=X.size()-1; i>=0; i--) {
        r[i] = r[i+1] + (X[X.size()-1-i]-'0');
    }

    // ans作成
    for(ll i=0; i<X.size(); i++) {
        ans[i] += r[i];

        ans[i+1] += ans[i]/10;
        ans[i] = ans[i]%10;
    }

    // 出力
    bool is_first_zero = true;
    for(ll i=MAX_ANS; i>=0; i--) {
        if (ans[i]==0 && is_first_zero) continue;
        if (ans[i]!=0) is_first_zero = false;
        if (!is_first_zero) {
            printf("%lld", ans[i]);
        }
    }
    printf("\n");
}

int main() {
    solve();
    return 0;
}

AHC007参戦記

THIRD プログラミングコンテスト 2021 (AtCoder Heuristic Contest 007)の参加記録です。

ヒューリスティックコンテストには参加したこと無かったのですが、今回は4時間と短めの時間設定だったので参加してみました。

問題概要

800x800の二次元座標上に400個の頂点と1995個の辺がランダムに与えられるので、すべての頂点が連結になるように辺を採択せよ。ただし、採用した辺のコストの総和が小さいほうがより良い得点になる(理想は最小全域木になる)。

得点:6,600,871,568(何も考えずUnionFind)

問題文を誤読していて、入力で各辺の実際のコストlが与えられるのに気づかずに提出して、意図しない形でACしました。

順番に与えられた辺に対して、「連結じゃないなら併合する」というUnionFindをし続けるとこの点数になります。

その後、実際の辺のコストの決まり方が一様分布と仮定して95%予言的中区間的なものを計算しようとしたり、「UnionFindで2点間の辺を取り除いたときにまだ連結しているか」みたいなのを実装してごちゃごちゃしてましたが、結局、最後まで時間内にこの得点を上回ることができないという悲しい結果になりました。

以下はUnionFindの実装。

#include <bits/stdc++.h>
#define _USE_MATH_DEFINES  // M_PI等のフラグ
#define MOD 1000000007
#define COUNTOF(array) (sizeof(array)/sizeof(array[0]))
#define rep(i,n) for (int i = 0; i < (n); ++i)
#define intceil(a,b) ((a+(b-1))/b)
using namespace std;
using ll = long long;
using pii = pair<int,int>;
using pll = pair<long,long>;
const long long INF = LONG_LONG_MAX - 1001001001001001;
void chmax(int& x, int y) { x = max(x,y); }
void chmin(int& x, int y) { x = min(x,y); }
 
template <typename T>
struct UnionFindVerSize {
    private:
        vector<T> sizes;   // グループのサイズ
        vector<T> parents;  // 親の番号を格納する。自分が親の場合は自分の番号になり、それがそのグループの番号になる
        T group_num;  // 連結成分の個数
 
    public:
        vector<set<ll>> graph;  // graph[i] := 頂点iと隣接している頂点リスト
        UnionFindVerSize(T N=0): sizes(N,1), parents(N) {
            for (T i=0; i<N; i++) {
                parents[i] = i;
            }
            this->group_num = N;
            graph.resize(N);
        }
 
        T find_root(T x) {
            /* ノード番号xの木の根(xがどのグループか)を求める */
            if (this->parents[x] == x) return x;
            this->parents[x] = this->find_root(this->parents[x]);  // 縮約
            return this->parents[x];
        }
 
        void unite(T x, T y) {
            /* ノード番号xとyが属する木を併合する(グループの併合) */
            T gx = this->find_root(x); T gy = this->find_root(y);
            if (gx == gy) return;
 
            // 深い方が浅い方を併合する(木の偏りが減るので)
            if (this->sizes[gx] < this->sizes[gy]) {
                this->parents[gx] = gy;
                this->sizes[gy] += this->sizes[gx];
            }
            else {
                this->parents[gy] = gx;
                this->sizes[gx] += this->sizes[gy];
            }
            this->group_num--;
 
            graph[x].insert(y);
            graph[y].insert(x);
        }
 
        void disunite(T x, T y) {
            graph[x].erase(y);
            graph[y].erase(x);
        }
 
        /* xとyの辺を取り除いても連結でいられるか?*/
        bool is_cycle_if_disunite(T x, T y) {
            deque<T> deq;
            set<T> visited;
            deq.push_back(x);
            while(!deq.empty()) {
                T u = deq.front(); deq.pop_front();
                if (visited.count(u)) continue;
                visited.insert(u);
 
                for(T v: graph[u]) {
                    if (u==x && v==y) continue;
                    if (visited.count(v)) continue;
                    if (v == y) {
                        return true;
                    }
                    deq.push_back(v);
                }
            }
            return false;
        }
 
        T get_size(T x) {
            /* ノード番号xが属するグループの深さを返す */
            return this->sizes[this->find_root(x)];
        }
 
        bool is_same_group(T x, T y) {
            /* ノード番号xとyが同じ集合に属するか否か */
            return this->find_root(x) == this->find_root(y);
        }
 
        T get_group_num() {
            /* 連結成分の個数を返す */
            return this->group_num;
        }
 
        void print_parents() {
            /* デバッグ用parentsの中身を出力する */
            for (T i=0; i<this->parents.size(); i++) {
                cout << this->parents[i] << endl;
            }
        }
 
        void print_sizes() {
            /* デバッグ用sizesの中身を出力する */
            for (T i=0; i<this->sizes.size(); i++) {
                cout << this->sizes[i] << endl;
            }
        }
};

以下はこの手法の実装。誤読して実際のコストlの入力を受け取ってない。うーん。

/* 2点間のユークリッド距離を整数に四捨五入する計算 */
ll distance(ll x1, ll y1, ll x2, ll y2) {
    return (ll)round(sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2)));
}
 
void solve() {
    ll N = 400, M = 1995;
    vector<ll> X(N), Y(N);
    for(ll i=0; i<N; i++) {
        cin >> X[i] >> Y[i];
    }
 
    UnionFindVerSize uf = UnionFindVerSize<ll>(N);
 
    for(ll i=0; i<M; i++) {
        ll u, v;
        cin >> u >> v;
        if (uf.get_group_num() == 1) {
            printf("0\n"); cout << flush;
            continue;
        }
        if (uf.is_same_group(u, v)) {
            printf("0\n"); cout << flush;
            continue;
        }
        printf("1\n"); cout << flush;
 
        uf.unite(u, v);
    }
}
 
 
int main() {
    solve();
    return 0;
}

得点:1,338,761,965(バグらせ)

すべての出力を1にする(すべての辺を採択する)とこの点数になります。

後述するクラスカル法をバグらせてこれをやりました。

得点:12,599,835,049(クラスカル法)

ここからはコンテスト後の話になります。

与えられる辺の順番はわかっているので、2点間の距離dは求めることができます。しかしこの問題のミソは、実際の辺のコストはdじゃなくて[d, 3d]の区間でランダムに決まるというものです。

でもまあ辺のコストを全て2dと仮定してクラスカル法で最小全域木を求めてもそれなりに良いスコアがでるだろう(そしてそれを改良していけばいいだろう)という仮説のもと、普通にクラスカル法を実装しました。するとこの点数になります。

以下はクラスカル法の実装。

struct edge {long long u, v, cost; };  // 頂点uと頂点vを結ぶ無向辺(コスト:cost)
 
bool comp(const edge& e1, const edge& e2) {
    return e1.cost < e2.cost;
}
 
pair<UnionFindVerSize<long long>, long long> kruskal(vector<edge> es, long long V, long long E) {
    sort(es.begin(), es.end(), comp);  // edge.costが小さい順にソートする
    UnionFindVerSize uf = UnionFindVerSize<long long>(V);
 
    long long res = 0;  // 最小全域木のコスト
    for (long long i=0; i<E; i++) {
        edge e = es[i];
        if (!uf.is_same_group(e.u, e.v)) {
            uf.unite(e.u, e.v);
            res += e.cost;
        }
    }
    return {uf, res};
}

以下はこの手法の実装。

void solve() {
    ll N = 400, M = 1995;
    ll MAX_DISTANCE = distance(0,0,800,800);
 
    vector<ll> X(N), Y(N);
    for(ll i=0; i<N; i++) {
        cin >> X[i] >> Y[i];
    }
 
    vector<edge> es(M);
    for(ll i=0; i<M; i++) {
        ll u, v;
        cin >> u >> v;
        edge e = {u, v, 2*distance(X[u], Y[u], X[v], Y[v])};
        es[i] = e;
    }
 
    auto [kuf, kuf_cost] = kruskal(es, N, M);
    for(ll i=0; i<M; i++) {
        ll l; cin >> l;
        auto [u, v, d] = es[i];
        if (kuf.graph[u].count(v)) {
            cout << "1" << endl << flush;
        }
        else {
            cout << "0" << endl << flush;
        }
    }
}

int main() {
    solve();
    return 0;
}

得点:13,901,664,333(毎回クラスカル法)

頂点数Nが400で辺数Mが1995なので、O(M2)くらいでもいけることに(ようやく)気づきます。

クラスカル法の計算量がO(|M| log|N|)なので(実際はソートやUnionFindを使うので違うが)、実際のコストlが与えられたときに毎回クラスカル法で最小全域木を構築しても計算量はO(M2 log|N|)で間に合います。

具体的には以下のような実装をしました。

  • 最初に辺(u, v)が与えられるとき、2点間の距離dを求めて、この辺のコストを2dと仮定しておく。
  • 実際の辺(u, v)のコストlが与えられるとき、辺のコストを2dからlに更新し、クラスカル法で最小全域木を求める。求めた最小全域木で辺(u, v)が採択されていれば、この辺を採択する。そうでなければこの辺は採択しない。
    • 採択した辺は覚えておく。これ以降のクエリで最小全域木を構築するときに、必ずこの辺を使わなければならないから。
    • 採択しなかった辺も覚えておく。これ以降のクエリで最小全域木を構築するときに、必ずこの辺を除外しなければならないから。

以下がこの手法の実装になります。

pair<UnionFindVerSize<long long>, long long> kruskal2(vector<edge> es, long long V, long long E, const vector<edge> &connected, const set<pair<ll,ll>> &disconnected) {
    sort(es.begin(), es.end(), comp);  // edge.costが小さい順にソートする
    UnionFindVerSize uf = UnionFindVerSize<long long>(V);
    long long res = 0;  // 最小全域木のコスト
 
    // すでに接続してるものがあるので、それは必ず使う
    for(edge e: connected) {
        uf.unite(e.u, e.v);
        res += e.cost;
    }
 
    for (long long i=0; i<E; i++) {
        edge e = es[i];
        if (disconnected.count({e.u, e.v})) continue;
 
        if (!uf.is_same_group(e.u, e.v)) {
            uf.unite(e.u, e.v);
            res += e.cost;
        }
    }
    return {uf, res};
}


void solve() {
    ll N = 400, M = 1995;
    ll MAX_DISTANCE = distance(0,0,800,800);
 
    vector<ll> X(N), Y(N);
    for(ll i=0; i<N; i++) {
        cin >> X[i] >> Y[i];
    }
 
    vector<edge> es(M);
    for(ll i=0; i<M; i++) {
        ll u, v;
        cin >> u >> v;
        edge e = {u, v, 2*distance(X[u], Y[u], X[v], Y[v])};
        es[i] = e;
    }
 
    vector<edge> connected;   // 接続することが決まった辺
    set<pair<ll,ll>> disconnected;  // 接続しないことが決まった辺
    for(ll i=0; i<M; i++) {  // O(M)
        ll l; cin >> l;
        auto [u, v, d] = es[i];
        es[i] = {u, v, l};
        auto [kruskal_uf, kruskal_cost] = kruskal2(es, N, M, connected, disconnected); // O(M logN)?
 
        if (kruskal_uf.graph[u].count(v)) {
            cout << "1" << endl << flush;
            // connected.insert(es[i]);
            connected.push_back(es[i]);
        }
        else {
            cout << "0" << endl << flush;
            disconnected.insert({es[i].u, es[i].v});
        }
    }
}

最初の辺のコストをすべて2dと仮定しましたが、うまい値に設定すれば上振れ狙えると思います。(モンテカルロ法を使えばいい?)

まとめ

  • コンテスト中の得点は意図せずACした6,600,871,568で、順位も539/632と残念な結果でした。最終的には得点を倍以上の13,901,664,333まで上げる実装を見つけることができました。
  • 解説はプリム法と書いてあるけどプリム法を知らなかった。
  • 4時間というコンテスト時間は個人的にはちょうど良かったです。
  • ビジュアライザがうまく実装されていてきれいで、動かすのが楽しかった。
  • C++のstd::setに構造体を入れることができないことを知りました。

得点:13,921,417,428【追記2021/12/15】

Twitter見てたら上のように書いている方がいて、なるほど多数決を採る手法があるのかと思いました。

というわけで以下の実装をしてみました。

  • 並行世界world[WORLD_NUM]を作成する。並行世界は各辺の情報をそれぞれ持つ。
  • 最初に辺(u, v)が与えられるとき、2点間の距離dを求めて、この辺のコストを[d, 3d]の区間でランダムに決める。これをすべての並行世界にやる。
  • 実際の辺(u, v)のコストlが与えられるとき、辺のコストをlに更新し、クラスカル法で最小全域木を求める。求めた最小全域木で辺(u, v)が採択されていれば、この辺を採択する。そうでなければこの辺は採択しない。これをすべての並行世界でやる。
    • 各並行世界で「辺を採択した/しなかった」のどちらが多かったかを多数決を採り、多かったほうを採用する。

それを実装したのが以下になるのですが、WORLD_NUM = 3程度でもTLEしました。

void solve() {
    random_device rnd;
    mt19937 mt(rnd());
    ll N = 400, M = 1995;
    // ll MAX_DISTANCE = distance(0,0,800,800);
 
    vector<ll> X(N), Y(N);
    for(ll i=0; i<N; i++) {
        cin >> X[i] >> Y[i];
    }
 
    ll WORLD_NUM = 19;  // 19
    vector<vector<edge>> world(WORLD_NUM, vector<edge>(M));  // 並行世界を作る
    for(ll i=0; i<M; i++) {
        ll u, v;
        cin >> u >> v;
        ll d = distance(X[u], Y[u], X[v], Y[v]);
        uniform_int_distribution<> myrand(d, 3*d);
 
        for(ll w=0; w<WORLD_NUM; w++) {
            // edge e = {u, v, (ll)myrand(mt)};  // ランダムに[d, 3d]に設定
            edge e = {u, v, 2*d};
            world[w][i] = e;
        }
    }
 
    // 毎回クラスカルする
    vector<edge> connected;   // 接続することが決まった辺
    set<pair<ll,ll>> disconnected;  // 接続しないことが決まった辺
    for(ll i=0; i<M; i++) {  // O(M)
        ll l; cin >> l;
        ll u = world[0][i].u;
        ll v = world[0][i].v;
 
        // 並行世界で多数決を採る
        ll vote = 0;  // 投票数
        for(ll w=0; w<WORLD_NUM; w++) {
            world[w][i] = {u, v, l};
            auto [kruskal_uf, kruskal_cost] = kruskal2(world[w], N, M, connected, disconnected); // O(M logN)?
 
            if (kruskal_uf.graph[u].count(v)) {
                vote++;
            }
        }
 
        // 過半数を越えてたら辺を採択
        if (vote > WORLD_NUM/2) {
            cout << "1" << endl << flush;
            connected.push_back({u, v, l});
        }
        else {
            cout << "0" << endl << flush;
            disconnected.insert({u, v});
        }
    }
}

TLEの主な原因は、毎回最初からクラスカル法をするのが重いからでしょう。

よくよく考えてみれば毎回最初からUnionFind木を構築する必要はなく、実際の辺のコストが与えられてその辺を採択する/しないを決めたあとは、そのときのUnionFind木を保存しといて次回はそこから再開すればいいです。

また辺のソートについても毎回最初からする必要はなく、確定した辺についてはソートしなくてよいので、sort(es.begin()+i, es.end(), comp);のように必要のある辺だけソートすればいいはずです(本当?)。

そのように改良したのが以下になり、これをするとWORLD_NUM = 5でも実行可能になり、得点は13,918,225,210~13,921,417,428くらいになりました。

pair<UnionFindVerSize<long long>, long long> kruskal3(ll start, vector<edge> es, long long V, long long E, const set<pair<ll,ll>> &disconnected, UnionFindVerSize<ll> uf, ll target_u, ll target_v) {
    // TODO 毎回ソートせんでもいいのでは?
    sort(es.begin()+start, es.end(), comp);  // edge.costが小さい順にソートする

    // 最小全域木は最後まで作る必要はなく、辺(target_u, target_v)が接続するかどうかわかったら打ち切ってよい
    for (long long i=start; i<E; i++) {
        edge e = es[i];
        if (disconnected.count({e.u, e.v})) continue;

        if (!uf.is_same_group(e.u, e.v)) {
            uf.unite(e.u, e.v);
        }

        if (e.u==target_u && e.v==target_v) break;
    }
    return {uf, -1};
}


// 並行世界を作る戦略
void solve() {
    random_device rnd;
    mt19937 mt(rnd());
    ll N = 400, M = 1995;
    // ll MAX_DISTANCE = distance(0,0,800,800);

    vector<ll> X(N), Y(N);
    for(ll i=0; i<N; i++) {
        cin >> X[i] >> Y[i];
    }

    ll WORLD_NUM = 5;  // 今の所5が限界…
    vector<vector<edge>> world(WORLD_NUM, vector<edge>(M));  // 並行世界を作る
    for(ll i=0; i<M; i++) {
        ll u, v;
        cin >> u >> v;
        ll d = distance(X[u], Y[u], X[v], Y[v]);
        uniform_int_distribution<> myrand(d, 3*d);

        for(ll w=0; w<WORLD_NUM; w++) {
            // edge e = {u, v, rand_distance(d)};
            edge e = {u, v, (ll)myrand(mt)};  // ランダムに[d, 3d]に設定
            // edge e = {u, v, 2*d};
            world[w][i] = e;
        }
    }

    // 毎回クラスカルする
    set<pair<ll,ll>> disconnected;  // 接続しないことが決まった辺
    UnionFindVerSize uf = UnionFindVerSize<ll>(N);  // クラスカル法を途中から再開できるように、予め確定した辺までを保存しておく
    for(ll i=0; i<M; i++) {  // O(M)
        ll l; cin >> l;
        ll u = world[0][i].u;
        ll v = world[0][i].v;

        // 並行世界で多数決を採る
        ll vote = 0;  // 投票数
        for(ll w=0; w<WORLD_NUM; w++) {
            world[w][i] = {u, v, l};
            auto [kruskal_uf, kruskal_cost] = kruskal3(i, world[w], N, M, disconnected, uf, u, v);

            if (kruskal_uf.graph[u].count(v)) {
                vote++;
            }
        }

        // 過半数を越えてたら辺を採択
        if (vote > WORLD_NUM/2) {
            cout << "1" << endl << flush;
            uf.unite(u, v);
        }
        else {
            cout << "0" << endl << flush;
            disconnected.insert({u, v});
        }
    }
}

元ツイートの方は並行世界を19個作れているのでまだ改良の余地はありますが、キリがないのでこのくらいにしようと思います。

073 - We Need Both a and b(★5)解いた

問題:競プロ典型 90 問 073 - We Need Both a and b(★5)

解説

以下の木DPを構築して解きたい。

dp[u][j] := (頂点uを根とする部分木において、)頂点uを含む連結成分の状態がjのときの場合の数。ただしj=0は'a'しかない状態、j=1は'b'しかない状態、j=2は'ab'両方ある状態

詳しい解説は公式解説こちらの記事に任せるとして、dpの更新式の勘所のみを解説したい。

dp[u][0]の更新式の考え方

  • 頂点uが'a'の場合、
    • (1)子の連結成分が'a'の状態かつ、連結している
    • (2)子の連結成分が'ab'の状態かつ、連結していない

f:id:takeg:20211203110534p:plain

「(1)子の連結成分が'a'の状態かつ、連結している」のパターンは明らかにOKである。なぜなら、子の連結成分が'a'の状態で頂点uと連結していれば、頂点uは'a'の状態のままだからである。

「(2)子の連結成分が'ab'の状態かつ、連結していない」のパターンも明らかにOKである。なぜなら、子の連結成分が'ab'の状態なら問題の制約上連結しなくてもOKであり、頂点uは'a'の状態のままだからである。

他のパターンも考えてみる。たとえば、「子の連結成分が'b'の状態かつ、連結していない」パターンはなぜ足しては駄目なのか。これだと「子の連結成分が'b'の状態のまま」になってしまい、これは問題の制約上「どの連結成分も'ab'の状態になっていなければならない」を満たせないからである。

dp[u][2]の更新式の考え方

  • 頂点uが'a'の場合、
    • (1)子の連結成分が'ab'の状態かつ、連結している
    • (2)子の連結成分が'b'の状態かつ、連結している
    • (3)子の連結成分が'a'の状態かつ、連結している
    • (4)子の連結成分が'ab'の状態かつ、連結していない

f:id:takeg:20211203100651p:plain

「(1)子の連結成分が'ab'の状態かつ、連結している」のパターンは明らかにOKである。なぜなら、子の連結成分が'ab'の状態で頂点uと連結していれば、頂点uは状態'ab'になるからだ。

「(2)子の連結成分が'b'の状態かつ、連結している」のパターンも明らかにOKである。 なぜなら、子の連結成分が'b'の状態で頂点uと連結していれば、頂点uは状態'ab'になるからだ。

「(3)子の連結成分が'a'の状態かつ、連結している」のパターンはなぜOKなのか。子の連結成分が'a'の状態で頂点uと連結すると頂点uは状態'a'のままだ。しかし、他の子の部分木がひとつでも(1)(2)のパターンであれば、頂点uは状態'ab'になれるので、(3)のパターンを含んでもよい。 ただし、他の子の連結成分も状態'a'の場合や(4)のパターンの場合、頂点uは'a'の状態のままである。なのでこの状態はあとから引かなければならない(★1)。

「(4)子の連結成分が'ab'の状態かつ、連結していない」のパターンはなぜOKなのか。問題の制約上、子の連結成分が'ab'の状態であるなら、連結していなくてもOKであり、頂点uは他の子の連結成分と連結して'ab'状態になれればOKだからだ。ただし、他の子の連結成分も(4)のパターンの場合、頂点uは'a'の状態のままである。なのでこの状態はあとから引かなければならない(★2)。

あとから引く(★1)+(★2)の場合の数とは、「頂点uが'a'の状態」のことであり、それはdp[u][0]のことに他ならない。

コード

提出コード(C++)

コールバック関数の処理も込みでawaitしたい

コールバック関数の処理込みでawaitしたいときは、以下のようにreturn new Promiseで囲むようにすればOK。

// 大雑把な説明
const f = () => {  // (3) 全体を関数化
    return new Promise((resolve, reject) => {  // (2) return new Promiseで囲む
        SomeFunction(callback);  // (1) コールバック関数込みでawaitしたい関数
    });
};

以下は、setTimeout()のコールバック関数の処理込みでawaitしたい場合の例。

const f1 = () => {
  return new Promise((resolve, reject) => {
    setTimeout(() => {
      return resolve("f1 ends.");  // コールバック関数の返り値
    }, 1000);
  });
};

const f2 = () => {
  return new Promise((resolve, reject) => {
    setTimeout(() => {
      return resolve("f2 ends.");
    }, 2000);
  });
};

const main = async () => {
  console.log(await f2());  // コールバック関数の返り値を受け取れるようにもなる
  console.log(await f1());
};

main();

// 出力結果は以下の順番になる。
// f2 ends.
// f1 ends.

CodePenの動作確認:https://codepen.io/taketakeyyy/pen/rNGBLeL

他にも、chrome.storage.sync.get()のようなものも同様に書ける。

/* TypeScriptで記述 */
// chrome.storage.sync.get(読み込みたいkey, callback関数)
const get_chrome_storage = (): Promise<string> => {
    return new Promise((resolve, reject) => {
        chrome.storage.sync.get({
            key: "default"
        }, (result) => {
            return resolve(result.key);  // これが返り値
        });
    });
};


const main = () => {
  const key = await get_chrome_storage();  // コールバック関数の処理が終了するまで待ち、返り値を受け取る
  console.log(key);
}

main();

IBM Quantum Challenge Fall 2021 体験記

challenges.quantum-computing.ibm.com

10月27日 ~ 11月06日の間に開催されて、量子コンピュータライブラリであるQiskitを使った4つのチャレンジにパスするとバッジがもらえるよ(何のバッジだ)、というハンズオンチャレンジでした。

バッジを貰える14問はクリアできましたが、最後のchallenge4-cはクリアできなかったです。

f:id:takeg:20211110010549p:plain
クリア時ステータス

今はCloseされたので、解法や感想を書こうと思います。といっても書くこと少ないのですが。

基本的にわからないことがあればQiskitのSlack(https://ibm.co/joinqiskitslack)に質問すれば誰かが答えてくれます。また、よくある質問や詰まるところのヒントもSlack内で検索すれば出てくるので、そういうのを駆使して問題を解いていきました。

このChallengeはあともう少しくらいは公開するらしく(バッジやスコアには反映されない)、それを過ぎたら削除される予定らしいので、やりたい方は早めに。やった方はバックアップ推奨。

私の提出コードはここにあります。

Challenge 1 Qiskit Finance - Optimizing your portfolio with quantum computers

金融においてはリスクを最小化しリターンを最大化することが大事で、そんなポートフォリオを作ることはポートフォリオ最適化問題と呼ばれています。それをQiskitを使って量子コンピュータで問題を解こうというチャレンジでした。

金融のこと何もわからないですが、書いてあるとおりに書き下していくと問題にパスできます。ライブラリの使い方を調べるゲームでした。

Challenge 2 Qiskit Nature - Band gap calculation of OLED molecules

有機EL分子のエネルギーバンドギャップの計算?

とたんに化学の知識を要求されるようで精神的にきつかったです。実際は書いてあるとおりにやって参考ビデオや文献を読めばやるべきことが見えてくるので、原子軌道がなんなのかとかは理解せずともクリアできます。

とはいえ唐突に出てきた分子軌道の計算はさっぱりだったので、適当にYouTubeで「分子軌道」と検索して出てきた動画(https://youtu.be/ZICW_IJRk4o?t=44)によると分子軌道=原子軌道ということが判明したのでこれでどうじゃい、をするとACしました。

あとはライブラリの使い方を調べて引数にこれとこれを与えたらこれが返ってくるので……をひたすらします。

ライブラリの使い方を検索したら古いバージョンのが最初に引っかかったりとしたので苦戦しました。

そういう試行錯誤をするゲームでした。

hallenge 3 Qiskit Machine Learning - Image classification by QSVM

量子機械学習による画像分類

Quantum SVMを体験しようというチャレンジです。

という感じのことをします。

ここのChallenge 3cが難関で、N_DIMrepsを単に大きくすればスコアが増えるわけではないので注意。ZZFeatureMapを使って試行錯誤でクリアしました。

Challenge 4 Qiskit Optimization - Battery Revenue Optimization

電池による収益の最適化

組合せ最適化をQAOAで解こうというチャレンジです。というかナップサック問題をQAOAで解こうというものです。

ナップサック問題競技プログラミング頻出なので、競プロやっててよかったと思いました。競技プログラミングIBM Quantum Challenge Fall 2021の役に立つ。

古典的な解法では、

# dp[i][C] := i日目のときの、C_max以内のコストでの、収益の合計の最大値
dp = [[0 for _ in range(C_max+1)] for __ in range(8)]
for i in range(1, 8):
    for c in range(0, C_max+1):
        dp[i][c] = dp[i-1][c]
        # 市場M1を選んだ場合
        if c-C1[i-1] >= 0:
            dp[i][c] = max(dp[i][c], dp[i-1][c-C1[i-1]]+L1[i-1])
        # 市場M2を選んだ場合
        if c-C2[i-1] >= 0:
            dp[i][c] = max(dp[i][c], dp[i-1][c-C2[i-1]]+L2[i-1])
print(dp[7][C_max])

程度のナップサック問題です。とはいえ古典的なナップサック問題の解法を書くわけでなく、QiskitのライブラリKnapsakで使える形式に問題を変形せよ、という問題でした。

参考文献の[2]にあるhttps://arxiv.org/pdf/1908.02210.pdfのペーパーの1.2節と1.3節を読めば大体わかります。

配列の長さt = 7
市場M1の収益L1 = [5,3,3,6,9,7,1]
市場M2の収益L2 = [8,4,5,12,10,11,2]
市場M1のコストC1 = [1,1,2,1,1,1,2]
市場M2のコストC2 = [3,2,3,2,4,3,3]
コストはC_max = 16 以下になるようにする

という2つの市場からうまく選んで収益を最大化するナップサック問題は、

newL = [L2[i]-L1[i] for i in range(len(L2))]
newC = [C2[i]-C1[i] for i in range(len(C2))]
newC_max = C_max - sum(C1)

という風に1次元に変換できる、というのが問題の主題です。

つまらないTypoで小一時間詰まりましたが、なんとかACできました。

Challenge 4cは49名の人がパスしたそうです。すごい!私は時間が足りずに無理でした。

バッジ

f:id:takeg:20211110005942p:plain
バッジ

Credlyというサイトでバッジがもらえました。Credlyを知らなかったのですが、これはオフィシャルバッジを作成・管理・配布するためのバッジプラットフォームらしいです。へぁ~ニッチ産業。と私は思いました。

まとめ

最大の敵はIBM Quantum LabのKernelフリーズでした。なんか全然パスできないと思っても、Kernelを再起動したらパスしたということがありました。

全体を通した学びとしては量子機械学習のお気持ちだけわかったかなといった感じです。

量子コンピュータは各Big Tech、スタートアップが力を入れて開発していますが、実用的なレベルにはまだまだ時間がかかりそうという話もよく聞きます。三体読んで夢が広がっているので、生きているうちに量子コンピュータで制御された恒星間航行用宇宙機器を作ってみたいものです。

以上。

0が重要な条件なら「i > 1」は「i-1 > 0」と書いて欲しいという提案

まずは数学の話から

高校生の頃は数学は好きな科目でしたが、必ずしも得意とは言えない科目でした。特に苦手だったのが教科書を読んでいて、

 a \ne -1のとき 」

のような条件を読み解くことが苦手でした。

たとえば、以下の数式を見てください。

 a \ne -1 のとき、 \lbrace f(x) \rbrace^{a} f'(x) = \frac{1}{a+1}( \lbrace f(x) \rbrace^{a+1})' が成り立つ。 」

この「 a \ne -1」 という条件はなぜ付いているのでしょうか? 数学に慣れている人は一瞬で気づくと思いますが、 \frac{1}{a+1} のゼロ除算を避けるためにあります。

こういうとき私は、「 a+1 \ne 0のとき」と書いてくれたほうがゼロ除算を避けるためなことが一目瞭然でわかりやすいのになあ、と思います。

ところでここで a \ne -1と書くモチベーションは何かというと、

  • 【動機1】こちらのほうが短くすっきりと書けるから
  • 【動機2】これ以降の展開で必要だから予め移項させている

だと思います。

これらのモチベーションは正しいと思います。そして私は数学を専門としているわけではなくこれ以上話を展開させる気はないので、数学の話はここで終わりにします。

プログラミングの話

ここからが本題なのですが、プログラミングにおいて「リーダブルな条件式の書き方とは何か?」を考えたとき、表題にあるように『0が重要な条件なら i > 1i-1 > 0 と書いたほうが読みやすい』という主張です。

例えば以下のコードを見てください。

for (int i=0; i<n; i++) {
  for (int w=0; w<=W; w++) {
    if (w >= weight[i]) {  // (★)
      dp[i+1][w] = max(dp[i][w], dp[i][w-weight[i]] + value[i]);
    }
    else {
      dp[i+1][w] = dp[i][w];
    }
  }
}

上記の if (w >= weight[i]) { // (★) という条件式はなんのためにあるのでしょうか?

上記のコードは典型的な動的計画法のコードなので、慣れている人は一瞬で気づくと思います。

もし読み解くことが難しいと思った人は、以下のようなコードならどうでしょうか? (★)の行しかコードを変えていません。

for (int i=0; i<n; i++) {
  for (int w=0; w<=W; w++) {
    if (w-weight[i] >= 0) {  // (★)
      dp[i+1][w] = max(dp[i][w], dp[i][w-weight[i]] + value[i]);
    }
    else {
      dp[i+1][w] = dp[i][w];
    }
  }
}

答えを書くと、(★)の直下の dp[i][w-weight[i]] の範囲外アクセスを避けるための条件式になります(配列の添字は0以上でなければならない)。

つまり上記の場合、条件式は if (w-weight[i] >= 0) と書いたほうが、直下の dp[i][w-weight[i]] と形が一致しているのでよりリーダブルになる、という主張です。

さて、「プログラミングにおけるリーダブル」という観点に重きを置いたとき、

  • 【動機1】こちらのほうが短くすっきりと書けるから

という動機においては、どちらの書き方も大差ないでしょう。

一方、

  • 【動機2】これ以降の展開で必要だから予め移項させている

に関しては、if (w-weight[i] >= 0) の形の方が、これ以降のコードで w-weight[i] の形が一致している箇所があるので軍配が上がります。

コードミスもこっちのほうが気づきやすいと思いませんか?