Dinkelbach algorithm

D: 食塩水 - AtCoder Beginner Contest 034 | AtCoderのD問題「食塩水」で分数計画問題(分数計画問題 - ORWiki)が出ました.

分数計画問題のアルゴリズムについては次のブログに詳しく載っています.
Dinkelbach algorithm - 週刊 spaghetti_source - TopCoder部
また,上のブログの実装例の最小比全域木問題の詳しい説明は,同じ方の 様々な全域木問題 に載っています.

今回は,「D:食塩水」を上のブログで説明しているDinkelbach algorithmで解きました.コンテストの解説(AtCoder Beginner Contest 034 解説)では,上のブログの二分探索による解法について説明を行っています.これは,蟻本のpp.132-134に掲載されている「平均最大化」と同じです.

[D. 食塩水]

問題

日本語なので省略.

方針

まず与えられた問題を分数計画問題の形にします.
 p_i %の食塩水  w_i gに入っている食塩の量 c_i w_i \times p_i / 100 gです.i番目の容器とj番目の容器の食塩水を混ぜると濃度は, (c_i + c_j) / (w_i + w_j) となります.したがって,容器の集合Xを選んだ時の混ぜた後の濃度は,
   \sum_{i \in X}c_i / \sum_{j \in X}w_j
となります.
ここで,
   f(X) = \sum_{i \in X}c_i  g(X) = \sum_{j \in X}w_j
とすると,求める問題は,
  maximize:   f(X) / g(X)
  subjecto to:  X \subseteq 2^{\{1,2, \cdots, N\}}
という分数計画問題となります.ここで, 2^{\{1, 2, \cdots, N\}} は与えられた食塩水が入った容器の部分集合全体とします.

関数fは食塩量,関数gは食塩水の量への離散関数となっており,最小比全域木問題と同じく{0,1}整数線形計画となっています.したがって,実装は上のブログのDinkelbach algorithmの実装例とほとんど同じになっています.
異なる点は最大化(ソート関数の向きに注意)であることと,h(t)の終了判定の誤差を「1e-8」ではなく「1e-3」にしたところです.「1e-8」にすると時間内に処理が終わりませんでした. 100 \times t は求める濃度を表しているのですが,終了判定の誤差は求める値の誤差と関係がなさそうなのでどれぐらい小さくすれば良いのかが難しいです.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

struct Cost {
    double w, p;
    Cost(double w, double p) : w(w), p(p) {}
    Cost() {}
};
int n, k;
vector<Cost> c;

double Check(double t, double &fcost, double &gcost)
{
    fcost = gcost = 0;

    sort(c.begin(), c.end(), [t](Cost a1, Cost a2) {
            return a1.p - t * a1.w > a2.p - t * a2.w; });

    for (int i = 0; i < k; ++i) {
        fcost += c[i].p;
        gcost += c[i].w;
    }

    return fcost - t * gcost;
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);
    cout << setprecision(8) << setiosflags(ios::fixed); // printf("%.4f", x);

    cin >> n >> k;
    c.resize(n);

    for (int i = 0; i < n; ++i) {
        cin >> c[i].w >> c[i].p;
        c[i].p *= c[i].w * 0.01; // i番目の食塩の量をp[i]とする
    }

    double t = 0.0, fcost, gcost;
    while (true) {
        if (fabs(Check(t, fcost, gcost)) < 1e-3)
            break ;
        t = fcost / gcost;
    }

    cout << 100.0 * t << endl;

    return 0;
}

TopCoder Arenaでの文字化け対処法

長年悩まされたArenaでの日本語の文字化けを直す方法が分かったので書きます.

問題点

TopCoderのArenaを次のコマンドを実行して立ち上げていました.

$ javaws ContestAppletProd.jnlp

このとき,「Chat Arena」で日本語表示・入力をしようとすると次のように豆腐(□□□)が表示されて読めませんでした.
f:id:tatanaideyo:20160310012833p:plain

使用している環境は次の通りです.

解決策

原因はArenaで使用していたフォントが日本語に対応していなかったことです.気付けば当たり前なんですが.
解決方法は次の通りです.
「Options」→「Setup User Preferences」→「Chat」タブ→「Font:」で「IPA P□□□□」を選ぶ.

IPA」の他には「Takao」でも大丈夫なので好きなのを選ぶと良いと思います.「IPA P□□□□」で豆腐が現れていますが,この文字化けを直す方法は分かりません(直さなくても特に支障は無い).

POJ 1087

情報工学工房で「1087 -- A Plug for UNIX」という問題を解いたのですが,他の人の解法を聞いてなるほどと思ったので記事を書きます.

問題文は長いです.私の解法は最大流問題に帰着する方法です.なるほどと思った解法は二部グラフの最大マッチング問題に帰着する方法なのですがグラフの作成の部分に推移閉包を用いています.それぞれについて説明を行います.

[問題内容]

n個のコンセントのプラグ受けの形状とm個のデバイスの差込プラグの形状が与えられます.また,k個の変換プラグの種類が与えられます.それぞれの変換プラグは無限個あるとします.
デバイスをできる限り充電(差し込むコンセントを決める)したときの充電できなかったデバイスの数を答える問題です.この時,変換プラグの差込プラグとプラグ受けの形状が合えば連結して使用することができます.

制約:  1 \le m \le 100, 1 \le n \le 100, 1 \le k \le 100

[解法1(最大流問題)]

私の解法で最大流問題に帰着する方法です.
説明を書いている途中で思ったのですが弧の向きを逆にすればコンセントの差込プラグとプラグ受けの形状のイメージが合うと思って反省をしましたが,書いてしまったのでそのままで説明します.
次の図のようにグラフを作成してs-t最大流を求めます.
f:id:tatanaideyo:20151209141056p:plain

弧の付け方は次の通りで,すべての弧の容量は1です.
1. コンセント側とデバイス側で同じ形状
2. コンセント側と変換プラグの差込プラグが同じ形状
3. デバイス側と変換プラグのプラグ受けが同じ形状
4. 変換プラグのプラグ受けと変換プラグの差込プラグが同じ形状

mから求めた流量を引いた値が充電出来なかったデバイスの個数で答えとなります.この例では,phoneは上から3番目のcのコンセント,pagerは上から1番目のaのコンセント,clockは上から2番目のコンセント,combは上から4番目のコンセントに充電することが出来ます.充電出来なかったデバイスはlaptopだけで出力する答えは1になります.
f:id:tatanaideyo:20151209143025p:plain

[解法2(二部グラフの最大マッチング)]

次は聞いてなるほどと思った解法について紹介をします.二部グラフの最大マッチング問題に帰着しています.
まずは変換プラグのことを忘れて二部グラフを作成します.
n個のコンセントとm個のデバイスをそれぞれ二部グラフの部集合として,コンセントのプラグ受けとデバイスの差込プラグの形状が等しい場合に辺を張ります.
f:id:tatanaideyo:20151209144152p:plain

次に,この二部グラフに変換プラグを考慮して辺を張っていきます.
変換プラグのプラグ受けと差込プラグの形状を頂点として,二項関係(x,y)をプラグ受けの形状xと差込プラグの形状をyとする変換プラグが存在するとして弧を張ったグラフを考えます.
次に,「(x,z)かつ(z,y)ならば(x,y)」という推移関係が存在しなくなるまで弧を付け加えていきます.
この操作は,下の図で言うと(b,x)の変換プラグを(x,a)の変換プラグに差し込んでプラグ受けの形状をb,差込プラグの形状をaとする変換プラグを作ることになります.
f:id:tatanaideyo:20151209152042p:plain

このように,推移関係が存在しなくなるまで弧を付け加えたものを推移閉包と呼びます.すなわち,変換プラグを直列に接続して作ることの出来る可能な全ての変換プラグを作成したことに等しくなります.
推移閉包のきちんとした定義等は次の資料を参考にしてください.

離散数学 第13回 関係(4) 関係の閉包 岡本吉央」http://dopal.cs.uec.ac.jp/okamotoy/lect/2015/discretemath/handout13.pdf

この推移閉包はワーシャル–フロイド法 - Wikipediaで求めることが出来ます.

この推移閉包を用いて二部グラフの辺を加えていきます.
推移閉包の関係(x,y)に対して,二部グラフのデバイスの形状がxでコンセントの形状がyとなるデバイスとコンセントの間に辺を張ります.
そして,このようにして構成できた二部グラフに対して最大マッチングを求めて,mからその値を引いたものが答えになります.
f:id:tatanaideyo:20151209155906p:plain

[まとめ]

「ワーシャルフロイド法 = 全点対最短距離を求める」という感じで覚えており,ソースコードも短いから見ないでも書ける!!けどアルゴリズムの大事なところを覚えてないという残念な感じなので反省します.
推移閉包を用いる問題に出会ったことがなかったので良かったです.

[ソースコード]

解法1(最大流問題)のソースコードを載せます.

// Time: 788K   Memory: 63MS
#include <iostream>
#include <vector>
 
using namespace std;
 
typedef int Weight;
const Weight INF = 1e8;
 
struct Edge {
    int src, dst;
    Weight cap, rev;
    Edge(int f, int t, Weight c, Weight r) :
        src(f), dst(t), cap(c), rev(r) {}
};
 
typedef vector<Edge> Edges;
typedef vector<Edges> Graph;
 
void AddEdge(Graph &g, int from, int to, int cap) {
    g[from].push_back(Edge(from, to, cap, g[to].size()));
    g[to].push_back(Edge(to, from, 0, g[from].size() - 1));
}
 
Weight Dfs(Graph &g, vector<bool> &used, int v, int t, Weight f) {
    if (v == t)
        return f;
    used[v] = true;
    for (int deg = 0; deg < g[v].size(); ++deg) {
        Edge &e = g[v][deg];
        if (!used[e.dst] && e.cap > 0) {
            int d = Dfs(g, used, e.dst, t, min(f, e.cap));
            if (d > 0) {
                e.cap -= d;
                g[e.dst][e.rev].cap += d;
 
                return d;
            }
        }
    }
    return 0;
}
 
// 最大流問題:計算量 O(FE + V^2), F: 最大流量
Weight FordFulkerson(Graph g, int s, int t) {
    const int n = g.size();
    Weight flow = 0;
    while (true) {
        vector<bool> used(n, false);
        Weight f = Dfs(g, used, s, t, INF);
        if (f <= 0)
            break;
        flow += f;
    }
    return flow;
}
 
// グラフを作成
Graph MakeGraph(int &n, int &m, int &s, int &t) {
    int k;
    string tmp_s;
 
    // Input: n_name[1:n], m_name[n+1:n+m],k_name[n+m+1:n+m+k]
    cin >> n;
    vector<string> n_name(n);
    for (int i = 0; i < n; ++i)
        cin >> n_name[i];
 
    cin >> m;
    vector<string> m_name(m);
    for (int i = 0; i < m; ++i)
        cin >> tmp_s >> m_name[i];
 
    cin >> k;
    vector<pair<string, string> > k_name(k);
    for (int i = 0; i < k; ++i)
        cin >> k_name[i].second >> k_name[i].first;
 
    s = 0; // source
    t = n + m + k + 1; // sink
    Graph g(t + 1);
 
    // Make Graph
    // s -> n_name
    for (int to = 1; to <= n; ++to)
        AddEdge(g, s, to, 1);
    // m_name -> t
    for (int from = n + 1; from <= n + m; ++from)
        AddEdge(g, from, t, 1);
    // n_name -> m_name, k_name
    for (int from = 1; from <= n; ++from) {
        // n_name -> m_name
        for (int to = n + 1; to <= n + m; ++to)
            if (n_name[from - 1] == m_name[to - n - 1])
                AddEdge(g, from, to, 1);
        // n_name -> k_name
        for (int to = n + m + 1; to <= n + m + k; ++to)
            if (n_name[from - 1] == k_name[to - n - m - 1].first)
                AddEdge(g, from, to, 1);
    }
    // k_name -> k_name, m_name
    for (int from = n + m + 1; from <= n + m + k; ++from) {
        for (int to = n + m + 1; to <= n + m + k; ++to)
            if (from != to && k_name[from - n - m - 1].second == k_name[to - n - m - 1].first)
                AddEdge(g, from, to, 1);
        for (int to = n + 1; to <= n + m; ++to)
            if (k_name[from - n - m - 1].second == m_name[to - n - 1])
                AddEdge(g, from, to, 1);
    }
 
    return g;
}
 
int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);
 
    int n, m, s, t;
 
    Graph g = MakeGraph(n, m, s, t);
 
    cout << m - FordFulkerson(g, s, t) << endl;
 
    return 0;
}

双方向ダイクストラ

・(追記)2015年11月3日:Mi_Sawaさんとtmaeharaさんのソースコードへのリンクとtmaeharaさんのプログラムとの比較

[概要]

双方向ダイクストラを初めて実装したので備忘録としてまとめておきます.
時間がないのでざっくりとしか書いてません.
双方向ダイクストラ法を知らなかったのでMi_Sawaさんに教えていただきました.有り難うございます.

無向グラフ上で始点sと終点tが与えられたときのsからtへの最短路を求める問題を考えます.この問題は2頂点対最短経路問題と呼ばれています.グラフの上で最短経路を求める問題は一般的に最短経路問題と呼ばれておりいろいろな種類があります.
最短経路問題 - Wikipedia

特にsから各頂点への最短経路を求める問題は単一始点最短経路問題と呼ばれており,その問題を解く有名なアルゴリズムダイクストラ法(ダイクストラ法 - Wikipedia)があります.
ダイクストラ法の詳しい説明はしませんがイメージとしてはsを中心に近い頂点から最短路を確定していきます.
f:id:tatanaideyo:20151031190922p:plain

単一始点最短経路問題ではsからすべての頂点への最短経路を求める問題ですが,途中で頂点tへの最短経路が確定したらアルゴリズムを終了するとして2頂点対最短経路問題としても考えることが出来ます.
sから広がっていくダイクストラ法ではsを中心に近い方向へ優先的に広がっていくのでtに到達するまでに余分な頂点を見ている気がします.
そこで,sからのみではなくtからも同じくダイクストラ法をして,お互いからの最短路が確定した頂点が見つかったときに終了すると余分な頂点が少なく済みそうな気がします.
この双方向から始めるダイクストラ法を双方向ダイクストラ法(Bidirectional search - Wikipedia, the free encyclopedia)と呼びます.注意が必要なのですがダイクストラ法と双方向ダイクストラ法の最悪計算量は同じです.下の図はイメージ図です(下の参考のp.16の図が分かりやすかったのでそのまま同じです).
f:id:tatanaideyo:20151031192716p:plain

[双方向ダイクストラ法のざっくり説明]

ダイクストラ法を知っている前提で説明します.
sとtからダイクストラ法を行います.頂点vにおいて初めてsからの最短路距離とtからの最短距離が確定されたとします.このとき,2つのダイクストラ法を終了します.
そして,各頂点vに対して(sからvへの最短距離)+(tからvへの最短距離)が最小となるものが最短距離となります.
最後の部分が必要な理由があまり分かっていないので分かったら書きます.
(下のソースコードの方が伝わるかもしれません)

tmaeharaさんのツイートで終了条件を工夫すると速くなるそうです.

[実験]

実験内容

最悪計算量は同じなので実際のインスタンで時間計測をすることによって比較することにします.用いるデータはTSP(巡回セールスマン問題 - Wikipedia)のベンチマークとして知られているアメリカの都市のデータです.
  データ・セット:USA Traveling Salesman Problem

データは二次元平面上の点の座標として与えられています.都市の数は115475です.各々の都市を頂点とした完全グラフで頂点間の距離はユークリッド距離として無向グラフを作成します(実際には頂点数が多いので陽にグラフを持たない).

*注意 変更(2015/11/3)
二次元平面上の点の集合からグラフを作るときに完全グラフとしましたが,完全グラフだとsとtの間に辺があり三角不等式から他の頂点を経由するよりもsとtの間の辺のみを使ったほうが最短路となります.
完全グラフに変換した下の実験結果はいろいろおかしなことになっています.点の広がり具合は正しいのですがおかしいです.点の集合からグラフを作成する方法を考える必要がありそうです.

グラフの描画はOpenFrameworks(openFrameworksJp)を使っています.
f:id:tatanaideyo:20151101013508p:plain

実験結果

一様ランダムにsとtを選びs-t最短路を次の3つのアルゴリズムで求めます.
 (a) ダイクストラ
 (b) 双方向ダイクストラ
 (c) tmaeharaさんの上界を考慮した双方向ダイクストラ

10セットに対して各々4回時間を計測しました.それぞれのtimeは4回の平均時間です.解が正しいかどうかは各々のアルゴリズムが出力する経路長が等しいことで確認をしました.
(a)と(b)のアルゴリズムのそれぞれのデータ・セットに対して最短経路が確定した頂点を色付けしています.sとtから最短経路が確定している頂点をそれぞれ赤色と青色とします.sとtは少し大きく表示しています.


(1) Find 107983-37236 shortest path = 37341.5
               Dijkstra : time = 35.75秒
 Bidirectional Dijkstra : time = 50.75秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101024812p:plain:w270 f:id:tatanaideyo:20151101024824p:plain:w270



(2) Find 30352-64115 shortest path = 14454.2
               Dijkstra : time = 47.5秒
 Bidirectional Dijkstra : time = 27.75秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101024843p:plain:w270 f:id:tatanaideyo:20151101024858p:plain:w270



(3) Find 100118-98660 shortest path = 2203.29
               Dijkstra : time = 1秒
 Bidirectional Dijkstra : time = 0.25秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101024912p:plain:w270 f:id:tatanaideyo:20151101024921p:plain:w270



(4) Find 14236-84642 shortest path = 25333.2
               Dijkstra : time = 49.25秒
 Bidirectional Dijkstra : time = 59秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101024935p:plain:w270 f:id:tatanaideyo:20151101024944p:plain:w270



(5) Find 88274-17971 shortest path = 7347.34
               Dijkstra : time = 22.75秒
 Bidirectional Dijkstra : time = 5.25秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101024957p:plain:w270 f:id:tatanaideyo:20151101025005p:plain:w270



(6) Find 15312-9902 shortest path = 45132
               Dijkstra : time = 58.25秒
 Bidirectional Dijkstra : time = 64.5秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101025023p:plain:w270 f:id:tatanaideyo:20151101025030p:plain:w270



(7) Find 101765-80234 shortest path = 13710.6
               Dijkstra : time = 10.25秒
 Bidirectional Dijkstra : time = 5.25秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101025045p:plain:w270 f:id:tatanaideyo:20151101025054p:plain:w270



(8) Find 66463-50876 shortest path = 12381.7
               Dijkstra : time = 11.25秒
 Bidirectional Dijkstra : time = 9.5秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101025103p:plain:w270 f:id:tatanaideyo:20151101025112p:plain:w270



(9) Find 35090-62401 shortest path = 19442.1
               Dijkstra : time = 51.5秒
 Bidirectional Dijkstra : time = 31秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101025123p:plain:w270 f:id:tatanaideyo:20151101025131p:plain:w270



(10) Find 43497-89411 shortest path = 15385.5
               Dijkstra : time = 21.5秒
 Bidirectional Dijkstra : time = 27.5秒
  Maeharasan no Program : time = 0秒

f:id:tatanaideyo:20151101025143p:plain:w270 f:id:tatanaideyo:20151101025153p:plain:w270

まとめ

右側の都市たちは密なので広がるスピードが遅く,左側の都市たちは疎なので広がるスピードが速く感じます.

アメリカの都市のデータは完全グラフでtmaeharaさんのアルゴリズムでは始点と終点から隣接頂点を見るステップで終了してるので0秒となっています.なので,計算時間での正しい比較ができていないので対策を考えたいと思います.

時間があれば上の参考資料に書かれている他の実装もやってみたいです.

[ソースコード]

双方向ダイクストラ法の実装では2つのプライオリティキューを使っています.Mi_Sawaさんの実装でpriority_queueをvectorに入れて管理していてとてもスッキリしていたので真似ました.

Mi_Sawaさんのソースコードhttps://ideone.com/UayEyc
tmaeharaさんのソースコード: https://ideone.com/J9awus


#include <bits/stdc++.h>

using namespace std;

typedef pair<double, double>  P;

vector<P> p; // 二次元平面上の点の集合

// 二点間の距離を返す
inline double Weight(const pair<double, double> &a, const pair<double, double> &b)
{
    return sqrt((a.first - b.first) * (a.first - b.first) +
                (a.second - b.second) * (a.second - b.second));
}

// s-t 双方向ダイクストラ法
double BidirectionalDijkstra(int s, int t)
{
    const int n = p.size(); // 頂点数
    const double MAX_W = DBL_MAX / 2;

    // d[0][v] := sからvへの最短距離, d[1][v] := tからvへの最短距離
    vector<vector<double>> d(2, vector<double>(n, MAX_W));
    // visit[0][v] := sからvへ更新終了, visit[1][v] := tからvへ更新終了
    vector<vector<bool>> visit(2, vector<bool>(n, false));
    // pq[0] := sからの更新情報を管理,pq[1] := tからの更新情報を管理
    vector<priority_queue<P, vector<P>, greater<P>>> pq(2);

    pq[0].emplace(d[0][s] = 0, s);
    pq[1].emplace(d[1][t] = 0, t);

    bool loop_continue = true;
    while (loop_continue && !pq[0].empty() && !pq[1].empty()) {
        for (int i = 0; i < 2; ++i) {
            if (pq[i].empty())
                continue;

            int u;
            double w;

            tie(w, u) = pq[i].top();  pq[i].pop();
            if (w <= d[i][u]) {
                visit[i][u] = true;

                // 終了条件: 頂点uにおいてsとtからの最短経路が確定
                if (w != 0 && visit[0][u] && visit[1][u]) {
                    loop_continue = false;
                    break;
                }

                for (int v = 0; v < n; ++v) {
                    if (u == v)
                        continue;
                    double cost = Weight(p[u], p[v]) + w;
                    if (cost < d[i][v]) {
                        d[i][v] = cost;
                        pq[i].emplace(d[i][v] = cost, v);
                    }
                }
            }
        }
    }

    // 最短距離はすべての頂点vに対してs-v-tの最短距離の最小値
    double ans = MAX_W;
    for (int v = 0; v < n; ++v)
        ans = min(ans, d[0][v] + d[1][v]);

    return ans;
}

// s-t ダイクストラ法
double Dijkstra(int s, int t)
{
    const int n = p.size(); // 頂点数
    const double MAX_W = DBL_MAX / 2;

    // d[0][v] := sからvへの最短距離, d[1][v] := tからvへの最短距離
    vector<double> d(n, MAX_W);
    // sからの更新情報を管理
    priority_queue<P, vector<P>, greater<P> > que;

    d[s] = 0.0;
    que.push(P(0.0, s));

    while (!que.empty()) {
        P now = que.top();  que.pop();
        int u = now.second;

        if (now.first <= d[u]) {
            // 終了条件
            if (u == t)
                return d[u];

            for (int v = 0; v < n; ++v) {
                if (u == v)
                    continue;

                double cost = Weight(p[u], p[v]) + now.first;
                if (MAX_W <= d[v] || cost < d[v]) {
                    d[v] = cost;
                    que.push(P(cost, v));
                }
            }
        }
    }

    return d[t];
}


int main()
{
    int n, id;
    double x, y;

    // Input
    cin >> n;
    cout << "Vertex size = " << n << "\n\n";
    p.resize(n);

    for (int i = 0; i < n; ++i)
        cin >> id >> p[i].first >> p[i].second;

    // Output
    for (int cnt = 0; cnt < 10; ++cnt) {
        int s = rand() % n, t = rand() % n;

        cout << "(" << cnt + 1 << ") Find " << s << "-" << t << " shortest path\n";

        // Dijkstra
        auto start = std::chrono::high_resolution_clock::now();
        cout << "               Dijkstra : len = " << Dijkstra(s, t);
        auto end = std::chrono::high_resolution_clock::now();
        auto take_time = std::chrono::duration_cast<std::chrono::seconds>(end - start);
        cout << ", time = " << take_time.count() << "秒\n";

        // Bidirectional Dijkstra
        start = std::chrono::high_resolution_clock::now();
        cout << " Bidirectional Dijkstra : len = " << BidirectionalDijkstra(s, t);
        end = std::chrono::high_resolution_clock::now();
        take_time = std::chrono::duration_cast<std::chrono::seconds>(end - start);
        cout << ", time = " << take_time.count() << "秒\n\n";
    }

    return 0;
}

SRM 671 Div2 (練習)

[結果]

SRM 671 Div2: xo-
練習で解きました.easyが落ちて笑えない.

[Easy 250: BearPaints ]

問題

H行W列の格子上で面積がM以下で最大の矩形の面積を求めよ.

  •  1 \le W, H \le 10^6
  •  1 \le M \le 10^{12}
方針

矩形の縦または横の長さを固定すると,面積を最大にする片方の長さは簡単に求まります.例えば縦をhで固定した場合の横wは min(M/h, W) となります.
縦と横は最大で  10^6 以下なので,各々すべての長さの場合を計算して最大の面積を求めればよいです.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

class BearPaints {
public:
	long long maxArea(int W, int H, long long M) {
        ll res = 0;

        for (int w = 1; w <= W; ++w) {
            ll h = min((ll)M / w, (ll)H);
            res = max(res, w * h);
        }
        for (int h = 1; h <= H; ++h) {
            ll w = min((ll)M / h, (ll)W);
            res = max(res, w * h);
        }

        return res;
	}
};

[Med 500: BearDartsDiv2 ]

問題

N個の要素からなるwが与えられる.
a < b < c < d で w[a] * w[b] * w[c] = w[d] となる組合せ数を求めよ.

  •  4 \le w.size() \le 500
  •  1 \le w_i \le 10^6
方針

SRM671 Div1のMedと同様の方法で解きます.
式変形をして,
 w_a \times w_b = \frac{w_d}{w_c}
としたときにcを基準に全探索します.
a < b < c を満たすすべてのaとbに対して,a * bの値の出現回数をmapで管理して,c < d を満たすすべてのdに対して w[d]/w[c] となる値の出現回数を求める答えに足していきます.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

class BearDartsDiv2 {
public:
	long long count(vector <int> w) {
        const int n = w.size();
		map<ll, ll> cnt;

        ll ans = 0;
        for (int i = 1; i < n; ++i) {
            for (int j = 0; j < i; ++j)
                ++cnt[w[j] * w[i]];
            for (int j = i + 2; j < n; ++j)
                if (w[j] % w[i + 1] == 0)
                    ans += cnt[w[j] / w[i + 1]];
        }

        return ans;
	}
};

[Hard 1000: BearDestroysDiv2 ]

問題

H行W列の格子状の森を考える.始めそれぞれのセルには1本の木がある.
次に,それぞれのセルを右(East)に倒すか下(South)に倒すかの情報が与えられたとする.左上のセルから右に向かってセル内の木を倒すかどうかの操作をしていく.一番右のセルが終了したら次の行の左から同じ操作を続ける.すべてのセルの操作を終えたら終了.
- 木が倒されていたら次のセルへ.
- 情報通りに倒せるなら倒す.倒せないなら逆の操作で倒せれるなら倒す.どっちでも倒せないなら何もしない.終わったら次のセルへ.

情報通りに倒すとは,情報がEastなら現在のセルが右端でなく,右のセル内の木が倒されていない時に,現在のセルと右のセル内の木を倒す.
情報がSouthなら現在のセルが下端でなく,下のセル内の木が倒されていない時に,現在のセルと下のセル内の木を倒す.

すべての操作を終えた時に倒された木の総数 / 2を倒した木の数とする.
すべてのどの方向に倒すかの情報(2^(H*W))の場合の倒した木の総数を求めよ.

  •  1 \le W \le 7
  •  1 \le H \le 40
方針

分からなかったのでkmjpさんの解法を見て解きました.
TopCoder SRM 671 Div2 Hard BearDestroysDiv2 - kmjp's blog

kmjpさんの解法を見ての反省点を書きます.
全列挙をすると O(2^(H*W)*H*W) となり間に合いません.ここで,シュミレーション内で重複してカウントしている部分を考えます(dpで解けないか).

ある行に着目します.どの方向に切るのかの情報が与えられたときに,どの部分がこの行を切るシミュレーション結果に影響を与えるのかを考えます.
ある列の木を倒すことを考えると,その前の行の同じ列の結果によって倒し方が変わることに気付きます.すなわち,前の行の同じ列がSouthに倒した(現在着目している木が倒された)かによって変わります.
したがって,現在の行のすでに倒されている木の情報とどの方向に倒すのかの情報が同じならばシミュレーション結果は同じです(2つ前の行以前の倒され方やどの方向に倒すかの情報は必要ない).
今回は倒されたか倒されていないかの結果の場合の数ではなく倒された木の数が必要なので,現在まででいくつの木を倒したのかの状態も必要となります.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;
class BearDestroysDiv2 {

public:
	int sumUp(int W, int H, int mod) {
        const ll MOD = mod;
        const int MAX_NUM = H * W / 2 + 2;
        // dp[処理した行番号][倒れた木の総数][次の行に向けて倒れたbitmask]
        //  := その状態になる組合せ数
        vector<vector<vector<ll>>> dp(H + 1,
                                      vector<vector<ll>>(MAX_NUM + W / 2 + 2,
                                                         vector<ll>(1 << (W + 1), 0)));

        dp[0][0][0] = 1;
        for (int h = 0; h < H; ++h) {
            for (int num = 0; num < MAX_NUM; ++num) {
                for (int s = 0; s < 1 << W; ++s) {
                    if (dp[h][num][s] == 0)
                        continue;
                    for (int sw = 0; sw < 1 << W; ++sw) {
                        // swのiビット目が0 := South
                        // swのiビット目が1 := East
                        int nxt_num = num; // この行までで倒された木の総数
                        int nxt_s = 0; // 次の行に向けて倒された木のbitmask
                        int now = 0; // この行で倒された木のbitmask

                        for (int w = 0; w < W; ++w) {
                            if (s >> w & 1 || now >> w & 1) // すでに倒されている
                                continue;

                            if (sw >> w & 1) { // East
                                if (w < W - 1 && !(s >> (w + 1) & 1)) { // Right
                                    ++nxt_num;
                                    now |= (1 << w) | (1 << (w + 1));
                                }
                                else if (h < H - 1) { // Down
                                    ++nxt_num;
                                    now |= (1 << w);
                                    nxt_s |= (1 << w);
                                }
                            }
                            else { // South
                                if (h < H - 1) { // Down
                                    ++nxt_num;
                                    now |= (1 << w);
                                    nxt_s |= (1 << w);
                                }
                                else if (w < W - 1 && !(s >> (w + 1) & 1)) { // Right
                                    ++nxt_num;
                                    now |= (1 << w) | (1 << (w + 1));
                                }
                            }
                        }

                        dp[h + 1][nxt_num][nxt_s] = (dp[h + 1][nxt_num][nxt_s]
                                                     + dp[h][num][s]) % MOD;
                    }
                }
            }
        }

        ll ans = 0;
        for (int num = 0; num < MAX_NUM; ++num)
            ans = (ans + (num * dp[H][num][0]) % MOD) % MOD;

        return (int)ans;
	}
};

まとめ

Easyをアドホックに解いて間違った.計算機をもっと信用しようという気持ちになった.
Hardは蟻本に乗っているFliptileのdp風でいけるかな思ったけど上手くいかず.

SRM 671 Div1

結果

--- 1317 -> 1292
今回も1問も解けなかった.

[Easy 300: BearCries ]

問題

';'と'_'からなる長さNの文字列messageが与えられる.
messageにはいくつかの悲しいを表す顔文字が含まれる.
悲しい顔文字は';'が両端に1文字づつあり,その間に1文字以上の'_'が含まれる文字列のことを言う.
最小の悲しい顔文字は";_;"である(他は";___:"とか";_______;").
messageを正しい悲しい顔文字に分ける分け方の総数を求めよ.

1 <= N <= 200

方針

動的計画法で解きます.
状態は次のように持ちます.
  dp[i][j][k] := messageの初めからi文字みたときに,ペアを持っていない';'の数がj個(顔文字の左端か右端か決まっていない),";_・・・"と顔文字の左端が決まっているものの数がk個という状態の顔文字の総数

求める状態は dp[N][0][0](messageの全体を見た時にペアが決まってない';'の数が0個の時の顔文字の総数)となります.

dp[i][j][k] からの状態遷移はmessageのi文字目(0-index)によって変わります.配るdpと呼ばれているもの?になると思います.詳しくはソースコードを見てください.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

const int MAX_N = 210;
const ll MOD = 1000000007;

inline void Add(ll &x, ll y) {
    x = (x + y) % MOD;
}

class BearCries {

    ll dp[MAX_N][MAX_N][MAX_N];
public:
	int count(string message) {
        string s = message;
        const int n = s.length();

        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j)
                for (int k = 0; k < n; ++k)
                    dp[i][j][k] = 0;

        dp[0][0][0] = 1;
        for (int i = 0; i < n; ++i)
            for (int j = 0; j <= i; ++j)
                for (int k = 0; k <= i; ++k) {
                    if (dp[i][j][k] == 0)
                        continue;

                    if (s[i] == ';') {
                        Add(dp[i + 1][j + 1][k], dp[i][j][k]);
                        if (0 < k) // 左端候補';_'のk個のいづれかと顔文字を作る
                            Add(dp[i + 1][j][k - 1], dp[i][j][k] * k);
                    }
                    else {
                        Add(dp[i + 1][j][k], dp[i][j][k] * k);
                        if (0 < j) // 未決定候補';'のj個のいづれかを左端候補にする
                            Add(dp[i + 1][j - 1][k + 1], dp[i][j][k] * j);
                    }
                }

        return (int)dp[n][0][0];
	}
};

[Med 500: BearDarts ]

問題

N個の配列wが与えられる.
a < b < c < d で w[a] * w[c] = w[b] * w[d] となる組合せ数を求めよ.

1 <= w[i] <= 10^9, 4 <= N <= 10^3

方針

roitiさんの解法を見て解きました.
Topcoder SRM671 Div1 Med: BearDurts - roiti46's blog

w[a]/w[b] = w[d]/w[c] と式変形をして,両辺をそれぞれ分子分母の最大公約数で割った値をペアで持つ.
a' = w[a]/gcd(w[a],w[b]), b' = w[b]/gcd(w[a],w[b])
c' = w[c]/gcd(w[c],w[d]), d' = w[d]/gcd(w[c],w[d])
とすると,
(a',b') = (d',c')
となるペアを求める問題となる.ただし,a < b < c < d .
cを0からN-1まで見るときに,mapには a < b < c となる (a',b') のペアがすべて含まれるようにする.そして,c < d を満たす (c,d) のペアに対してmapで検索すると計算量は O(N^2 logN) となる.

(規約分数やmapに登録する探索順番など自分の注意力では解けなかった気がする)

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

class BearDarts {
public:
	long long count(vector <int> w) {
        const int n = w.size();

        ll ans = 0;
		map<pair<int, int>, ll> cnt;

        for (int i = 1; i < n; ++i) {
            for (int j = 0; j < i; ++j) {
                int gcd = __gcd(w[j], w[i]);
                ++cnt[make_pair(w[j] / gcd, w[i] / gcd)];
            }
            for (int j = i + 2; j < n; ++j) {
                int gcd = __gcd(w[j], w[i + 1]);
                ans += cnt[make_pair(w[j] / gcd, w[i + 1]/ gcd)];
            }
        }

        return ans;
	}
};

CODE FESTIVAL 2015 予選A

結果は3完の447位でした(もうダメすぎる).
問題文は日本語なので省略します.
Welcome to CODE FESTIVAL 2015 予選A - CODE FESTIVAL 2015 予選A | AtCoder

[A: CODE FESTIVAL 2015 ]

方針

与えられた文字列の末尾には必ず"2014"という文字列があるので最後の文字を'4'から'5'に変更.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

int main()
{
    string s;
    cin >> s;
    cout << s.substr(0, s.size() - 1) << "5\n";

    return 0;
}

[B: とても長い数列 ]

方針

観察するとi番目の数列A_iはpow(2, n - i - 1)回足される.
int型で計算結果は収まるけど念の為にlong long 型で計算.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    ll ans = 0, a;
    int n;
    
    cin >> n;
    for (int i = 0; i < n; ++i) {
        cin >> a;
        ans += pow(2, n - i - 1) * a;
    }

    cout << ans << endl;

    return 0;
}

[C: 8月31日 ]

方針

A_iの和をsum_aと置くと求める問題は次のように書ける.

minimize      x_1 + x_2 + \ldots + x_N
subjecto to     sum_a - (A_1 - B_1) \times x_1 - \cdots - (A_N - B_N) \times x_N \le T
          x_1, x_2, \ldots , x_N \in \{0, 1\}

変数x_iが1のときB_iを選択,0のときA_iを選択するとする.
最適解を求めるためには,初め初期解を(x_1, ..., x_N)=(0, ..., 0)として,条件式を満たすまで 係数 A_i - B_i > 0 の大きいものから変数X_iを1にしていく.

ソースコード
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

ll Solve()
{
    int n;
    ll t, sum_b = 0, sum_a = 0, a, b;

    cin >> n >> t;
    vector<int> d(n);

    for (int i = 0; i < n; ++i) {
        cin >> a >> b;
        sum_b += b;
        sum_a += a;
        d[i] = -(a - b);
    }

    if (t < sum_b)
        return -1;
    if (sum_a <= t)
        return 0;

    sort(d.begin(), d.end());

    int num = 0;
    for (int i = 0; i < n; ++i) {
        sum_a += d[i];
        ++num;
        if (sum_a <= t)
            return num;
    }
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    cout << Solve() << endl;

    return 0;
}

[D: 壊れた電車 ]

方針

点検作業時間kに関して二分探索をする.
点検作業時間がkのときすべての車両を点検できるかを判定する.

一番左端の車両は一番左の整備士が点検をする方が良い.
i番目の整備士が点検をしないといけない一番左端の車両をlとする.

・点検作業時間k内にlを点検できないならどのようにしてもk以内には点検できない.
・lを点検できるときは残りの時間でできるだけ右端の車両rを点検する.
 □ x_i <= l なら r = x_i + k
 □ l < x_i なら点検方法は次の2つの内rが最大となるものを選ぶ.
f:id:tatanaideyo:20150927031136p:plain

次のi+1番目の整備士はr+1をlとして同様に点検をする.
最後の整備士が点検を終えた時にrがN以上なら点検可能.それ以外は点検不可能となる.

ソースコード
// 反省点1:nとmをlong long型にしたら通った(オーバーフロー)
// 反省点2: 下のmaxで右辺しか考慮してなかった(よく考えるとkが大きければ左辺の方が得)
//                                         (kが小さいと右辺の方が得.式をいじると分かる)
#include <bits/stdc++.h>

using namespace std;

typedef long long  ll;

ll n, m;

bool C(ll k, const vector<int> &x)
{
    ll l = 1;
    for (int i = 0; i < m; ++i) {
        ll p = x[i];

        if (k < p - l)
            return false;

        if (p <= l)
            l = p + k + 1;
        else
            l = max(k + 2 * l - p, p + (k - (p - l)) / 2) + 1; // 反省点2
    }

    if (n <= l - 1)
        return true;
    else
        return false;
}

int Solve()
{
    cin >> n >> m;
    vector<int> x(m);
    for (int i = 0; i < m; ++i)
        cin >> x[i];

    ll lb = -1, ub = 3 * n;
    while (ub - lb > 1) {
        ll mid = (lb + ub) / 2;
        if (C(mid, x))
            ub = mid;
        else
            lb = mid;
    }

    return ub;
}

int main()
{
    cin.tie(0);
    ios::sync_with_stdio(false);

    cout << Solve() << endl;

    return 0;
}