2016-09-24

凸包の頂点数の話

まとめ

$m\times m$ のグリッドに収まる格子点を頂点とする凸多角形の頂点数は $O(m^{2/3})$ 。
詳しい話は論文を読んでください。
On the maximal number of edges of convex digital polygons included into an m × m-grid

はじめに

凸包を取ると考慮すべき点の数が減って嬉しい、という問題が競プロでたまに出題されます。
この類の問題の解説で頂点数のオーダーが $O(\sqrt{m})$ と書いてあることがある*1のですが、実際は $O(m^{2/3})$ だよという話です。*2

方針

頂点数を大きくするにはマンハッタン距離が短い辺から使うのが最適です。
マンハッタン距離が $n$ 以下の辺を使うと、頂点数が $\Theta(n^2)$ の凸多角形が $\Theta(n^3)$ のグリッドに収まります。

計算

めんどうなので、成分が正の部分だけ考えます。*3
そうすると使う辺は(1,1),(1,2),(2,1),(1,3),(2,2),(3,1),...となります。
これを角度でソートして並べると凸多角形の右下部分になります。
本当は成分が互いに素な辺だけを使う必要がありますが、とりあえず無視します。

マンハッタン距離が $n$ 以下の辺を全て使う場合、各座標の総和は
$\sum_{i=1}^{n-1} i(n-i)=(n-1)n(n+1)/6$
となります。

辺の数は
$\sum_{i=1}^{n-1}i=(n-1)n/2$
となります。

グリッドのサイズが $\Theta(n^3)$ 、頂点数が $\Theta(n^2)$ なので、 $m\times m$ のグリッドに $\Theta(m^{2/3})$ 頂点の凸多角形を収めることができています。

互いに素

上の計算では成分が互いに素でない辺まで使っていましたが、互いに素である確率は定数なので、これを使わないようにしてもオーダーは変わりません。
互いに素 - Wikipedia

*1:どうやら蟻本の誤りに起因するらしい

*2:Codeforcesの記事で見たはずなんですがURLがわからなくなった

*3:辺を反時計回りにベクトルで持つとすると、一番下から一番右の部分

2016-03-20

ナップサック問題でマラソンマッチ入門

マラソンマッチって？

競技プログラミングのうち、「より良い解を求める」ことを競うコンテストをマラソン形式と呼びます。
例えば厳密解を求めることができない問題について、近似解のスコアを競ったりします。

マラソンでよく使われるアルゴリズム

マラソンで頻出なのは「ビームサーチ」と「焼きなまし法」です。
この記事ではナップサック問題を例にしてこの２つのアルゴリズムを解説します。
もちろん、この２つのアルゴリズムはどちらも近似アルゴリズムなので最適解は求められません。

ビームサーチ

まずは順番にナップサックに入るだけ入れるコードを書いてみます。

#include <iostream>

using namespace std;

int main() {
  // 個数
  const int N = 100000;
  // ナップサックの大きさ
  const int W = 100000;

  // 重さ・価値
  int w[N], v[N];
  for (int i = 0; i < N; ++i) cin >> w[i] >> v[i];

  // 今までに入れた物の重さの総和
  int s = 0;
  // 今までに入れた物の価値の総和
  long long t = 0;

  // 入るだけ入れる
  for (int i = 0; i < N; ++i) {
    if (s + w[i] <= W) {
      s += w[i];
      t += v[i];
    }
  }

  // 結果
  cout << t << endl;
}

これをビームサーチに書き換えるとこうなります。

#include <algorithm>
#include <iostream>
#include <vector>

using namespace std;

int main() {
  // 個数
  const int N = 100000;
  // ナップサックの大きさ
  const int W = 100000;
  // ビーム幅
  const int BEAM = 100;

  // 重さ・価値
  int w[N], v[N];
  for (int i = 0; i < N; ++i) cin >> w[i] >> v[i];

  // 現在の状態（価値・重さ）
  vector<pair<long long, int>> curr;
  curr.emplace_back(0, 0);

  // ビームサーチ
  for (int i = 0; i < N; ++i) {
    // 次の状態（価値・重さ）
    vector<pair<long long, int>> next;
    for (const auto& s : curr) {
      next.emplace_back(s);
      if (s.second + w[i] <= W) {
        next.emplace_back(s.first + v[i], s.second + w[i]);
      }
    }
    // 価値が大きい順に並び替え
    sort(next.rbegin(), next.rend());
    // 状態をビーム幅に減らす
    if (next.size() > BEAM) next.resize(BEAM);
    curr = move(next);
  }

  // 結果
  cout << curr[0].first << endl;
}

それぞれの状態について物を入れるか入れないかの遷移を行って、スコアの良い上位１００個の状態を残しています。
この１００個のことをビーム幅と呼び、基本的にはこれが大きければ大きいほど良いので、ビームサーチを使うときは高速化が効果的である場合が多いです。

どれくらい改善するかを試した所、以下のようになりました。

アルゴリズム	スコア
愚直	829374
ビームサーチ	3252218
厳密解	26383646

実験に使った入力はこれです。
愚直解に比べるとだいぶ改善していますが、厳密解はかなり遠いです。
これは、問題の性質を上手く使えていないのが原因です（最後に書きます）。

焼きなまし

今度はまだナップサックに入れていない荷物を選んで、選んだ荷物が入るように入っている荷物をランダムに捨てるとしてスコアが改善するなら実際に荷物を入れる、という方法を考えます。

#include <algorithm>
#include <iostream>

using namespace std;

// 乱数生成器
inline uint32_t rnd(void) {
  static uint32_t y = 2463534242;
  y = y ^ (y << 13); y = y ^ (y >> 17);
  return y = y ^ (y << 5);
}

int main() {
  // 個数
  const int N = 100000;
  // ナップサックの大きさ
  const int W = 100000;
  // ループ回数
  const int T = 1000000;

  // 重さ・価値
  int w[N], v[N];
  for (int i = 0; i < N; ++i) cin >> w[i] >> v[i];

  // 今までに入れた物の重さの総和
  int s = 0;
  // 今までに入れた物の価値の総和
  long long t = 0;
  // 今までに入れた物・入れていない物
  vector<int> used, unused(N);
  // はじめは全て使っていない
  iota(unused.begin(), unused.end(), 0);

  for (int i = 0; i < T; ++i) {
    // 次に入れる物をランダムに選ぶ
    int k = rnd() % unused.size();

    // ナップサックに収まるようにランダムに捨てる
    int ss = 0;
    long long tt = 0;
    vector<int> del;
    while (tt < v[k] && s + w[unused[k]] - ss > W) {
      int kk = rnd() % used.size();
      if (find(del.begin(), del.end(), kk) != del.end()) continue;
      ss += w[used[kk]];
      tt += v[used[kk]];
      del.emplace_back(kk);
    }

    // 価値が大きくなるなら更新
    if (tt < v[k] && s + w[unused[k]] - ss <= W) {
      s += w[k] - ss;
      t += v[k] - tt;
      unused[k] = unused.back();
      unused.pop_back();
      used.emplace_back(k);
      sort(del.rbegin(), del.rend());
      for (int kk : del) {
        used[kk] = used.back();
        used.pop_back();
        unused.emplace_back(kk);
      }
    }
  }

  // 結果
  cout << t << endl;
}

これを焼きなまし法に書き換えるとこうなります。

#include <algorithm>
#include <iostream>

using namespace std;

// 乱数生成器
inline uint32_t rnd(void) {
  static uint32_t y = 2463534242;
  y = y ^ (y << 13); y = y ^ (y >> 17);
  return y = y ^ (y << 5);
}

bool next(long long b, int t) {
  if (b > 0) return true;
  return exp(b / 10000.0 / pow(0.999999, t)) > double(rnd()) / numeric_limits<uint32_t>::max();
}

int main() {
  // 個数
  const int N = 100000;
  // ナップサックの大きさ
  const int W = 100000;
  // ループ回数
  const int T = 1000000;

  // 重さ・価値
  int w[N], v[N];
  for (int i = 0; i < N; ++i) cin >> w[i] >> v[i];

  // 今までに入れた物の重さの総和
  int s = 0;
  // 今までに入れた物の価値の総和
  long long t = 0;
  // 今までに入れた物・入れていない物
  vector<int> used, unused(N);
  // はじめは全て使っていない
  iota(unused.begin(), unused.end(), 0);

  for (int i = 0; i < T; ++i) {
    // 次に入れる物をランダムに選ぶ
    int k = rnd() % unused.size();

    // ナップサックに収まるようにランダムに捨てる
    int ss = 0;
    long long tt = 0;
    vector<int> del;
    while (tt < v[k] + 10000 && s + w[unused[k]] - ss > W) {
      int kk = rnd() % used.size();
      if (find(del.begin(), del.end(), kk) != del.end()) continue;
      ss += w[used[kk]];
      tt += v[used[kk]];
      del.emplace_back(kk);
    }

    // 焼きなまし
    if (next(v[k] - tt, i) && s + w[unused[k]] - ss <= W) {
      s += w[k] - ss;
      t += v[k] - tt;
      unused[k] = unused.back();
      unused.pop_back();
      used.emplace_back(k);
      sort(del.rbegin(), del.rend());
      for (int kk : del) {
        used[kk] = used.back();
        used.pop_back();
        unused.emplace_back(kk);
      }
    }
  }

  // 結果
  cout << t << endl;
}

スコアが改善しない時も、選んだ荷物をたまには使うようにします。
選んだ荷物を使う確率は、スコアの下落幅が小さい・ループの始めのほうで大きくなるようにします。
ループの後半になるに従って荷物を使う確率を下げていきますが、これに用いるパラメータを温度と呼びます。
このような温度管理によって局所解に陥ることをある程度回避することができます。

スコアは以下のようになりました。

アルゴリズム	スコア
愚直	18241125
焼きなまし法	21275856
厳密解	26383646

終わりに

ビームサーチと焼きなまし法を紹介しましたが、これを使えばマラソンマッチで簡単に勝てるかというとそんなことはありません。
一番大事なのは問題の性質を理解することです。
例えば今回の問題の場合、価値÷重さでソートした後に順番に詰めていくだけでかなり良いスコアが出ます。

#include <algorithm>
#include <iostream>

using namespace std;

struct Item {
  int w, v;
  bool operator>(const Item& i) const {return (double)v / w > (double)i.v / i.w;}
};

int main() {
  // 個数
  const int N = 100000;
  // ナップサックの大きさ
  const int W = 100000;

  // 荷物
  Item item[N];
  for (int i = 0; i < N; ++i) cin >> item[i].w >> item[i].v;
  sort(item, item + N, greater<Item>());

  // 今までに入れた物の重さの総和
  int s = 0;
  // 今までに入れた物の価値の総和
  long long t = 0;

  // 入るだけ入れる
  for (int i = 0; i < N; ++i) {
    if (s + item[i].w <= W) {
      s += item[i].w;
      t += item[i].v;
    }
  }

  // 結果
  cout << t << endl;
}

アルゴリズム	スコア
ビームサーチ	3252218
焼きなまし法	21275856
ソート	26382284
厳密解	26383646

この解法を元にしてビームサーチや焼きなまし法をすることはできますが、そのようなアルゴリズムの知識よりも問題の考察が重要になることが多いです。
この辺りの話はtakaptさんの記事によくまとまっています。
もちろん、知識は多いほど有利になるのは間違いありません。
ぜひマラソンマッチに参加して経験を詰んだり他の参加者から知識を吸収したりして下さい。
そう言えば今日、Chokudai Contest 001ってのがAtCoderで開催されるらしいので参加してみてはどうでしょうか（ステマ）。

2015-12-07

短縮王「割り算と足し算」解説（C言語部門）

これは Competitive Programming （その2） Advent Calendar 2015 の7日目の記事です。

CODE FESTIVAL 2015の短縮王で出題した「割り算と足し算」のC言語部門の解説です。

この問題のC言語部門で1位を獲得した ehaさんのコードを元に解説します。

n;f(x,i){n=++i>x?:fmax(f(x%i?0:x/i,1)+x/10+x%10,f(x,i));}
main(){n=scanf("%d",&n)/printf("%d\n",99/n?f(n,1):19);}

(112byte)
(表示の都合上、改行を入れています。以下同様。)

ehaさんのコードの特徴は i についてのループを回す代わりに再帰で処理している、という点です。
この発想がないと、なかなかコードが縮みません。
ジャッジ側が用意していた解答もループを回すものだったので、ehaさんの解答はなかなか衝撃的でした。

さて、この解は x/10+x%10 の部分が余分なので1byte削ることが出来ます。
また、入力が100の場合のために処理を分ける必要もありません。

n;f(x,i){n=++i>x?:fmax(f(x%i?0:x/i,1)+x-x/10*9-x/100*9,f(x,i));}
main(){n=scanf("%d",&n)/printf("%d\n",f(n,1));}

(111byte)

CODE FESTIVAL中にはここまでしか縮みませんでした。
後日、さらに頑張るとまだまだ無駄があることがわかりました。
まずは、この問題の全言語部門で2位を獲得した hirokazuさんのコードを参考にすると3byte縮みます。

n;f(x,i){n=++i>=x?:fmax(f(i,x%i*x)+x-x/10*9-x/100*9,f(x,i));}
main(){n=scanf("%d",&n)/printf("%d\n",f(n,0));}

(108byte)

ずいぶんスッキリしましたが、x-x/10*9-x/100*9 はいかにも縮みそうな気配がします。
この式は x/10*9+x/100*9 → (x/10+x/100)*9 → x/10*11/10*9 という変形をすると2byte縮めることが出来ます。

n;f(x,i){n=++i>=x?:fmax(f(i,x%i*x)+x-x/10*11/10*9,f(x,i));}
main(){n=scanf("%d",&n)/printf("%d\n",f(n,0));}

(106byte)

これが現在見つかっている最短コードです。
まだまだ縮みそうな気がしてきますね。
そう思った人は頑張って考えてみましょう。

さて、明日は

動作を停止しました。'8日目の担当' は決まっていません。

2015-11-16

CODE FESTIVAL 2015 参加記

N日目(N<0)

本戦・短縮王・リレーの問題準備をする。短縮王をテストプレイしたら徹夜だった。

短縮王の順位表を宇宙ツイッタラーさんに依頼する。

0日目

スタッフで前夜祭をする。

1日目

朝

5時起床に失敗して若干TLEする。

リハをする。徐々に目が覚めてくる。

選手が到着したので雑用をする。

ネットワークトラブルが発生する。できることがないので早く解消してくれ〜と祈る。

本戦

解説スライドのチェックをしたり短縮王ランキングのチェックをしたり書道コーディングしたり雑談したりする。

個別指導塾

あんまり人が来ないので、雑談しに来た人に強引に解説したりする。

エキシビション

ボドゲしてた。

夜

僕「リレーの説明の打ち合わせしたいんでそっちの部屋に行っていいですか」

y3「ok」

〜2分後〜

部屋を尋ねる。

y3「(˘ω˘)ｽﾔｧ」

僕「!?」

しょうがないので打ち合わせは明日することにして短縮王の提出を眺めたりする。

2日目

朝

4時起床。短縮王のC言語部門のジャッジ解を整理する。 C問題のehaさんの解を1バイト縮める。D問題のジャッジ解を2バイト縮める。

短縮王の順位表の修正を依頼する。

あさプロ

リレーの説明の原稿を作ってた。

昼

shinhさんがいらしたので最短コードとか表彰とかの打ち合わせをする。

リレー準備

今回の最大の佳境。

問題文が印刷されていないことが発覚する。 JAG夏合宿2日目のCSSを流用してその場で印刷用問題文を作る。

設営の人手が足りないと言われたのでディスプレイの設置とかをする。

コンテストの設定あってるかな〜とぼんやり眺めてたら、待機スペースのPCが管理者権限でコンテストにログインしてたのであわててけんしょーさんを呼ぶ。正しいIDでログインしなおす。

選手が入場してからチーム番号の1〜20とA〜Tが対応していないことに気づき大移動をしてもらう。そういえば打ち合わせで断片的に聞いた気がする…。

リレーの説明をする。それなりに笑いをとれたのでまぁ良し。

印刷された問題文が届かない。自分がやるべきだなぁと思ったのでスタッフ間の連絡役をやる。インタビューで時間が持たなくなったら何しようかなぁと考える。

リレー

問題文が届いたので15分遅れで開始。

システムがちゃんと動いてそうで安心する。

みんながわいわい問題を解いてるのを眺める。

あまりに解かれる速度が速いので、頼む〜詰まってくれ〜と祈る。

祈りが通じず40分すぎに全完がでる。マジかよ。

MVTを決めるとかいう無茶ぶりな仕事をこなす。

バタバタしたけど無事終了。表彰式をやる。

表彰式

大量のクリーナーを押し付ける。

きゅうりは天才。

天才情報です。 #code_festival pic.twitter.com/Zs9jDwmips
— すめけ (@snuke_) 2015, 11月 15

打ち上げ

頼んでもないのに突然運ばれてきてびっくりする。

pic.twitter.com/orwnFacBfc
— not (@not_522) 2015, 11月 15

感想

めっちゃ楽しかったです。

リレーでめっちゃトラブったのは申し訳なかったけど、トラブルを収束させるの好きなので楽しかった。

短縮王はshinhさんがガッツリ解説してくれるなら、もっとマニアックな問題にすればよかったかなぁと思った。ショートコーディング初めての人でも、がんがん縮むから楽しいかなぁと思って問題作ったけど、どうなんだろう。

おまけ

【焼肉チャレンジ】リレーの問題を4問作りました。僕が作った問題を一発で4問当てた人に焼肉おごります。 #code_festival
— not (@not_522) 2015, 11月 15