成績

POH第2弾は平面上の長方形領域の探索です。背景や課題の詳細は
https://paiza.jp/poh/paizen
を参照してください。#paizahack_02

最後のテストケースで0.02秒の壁がなかなか破れず苦戦しましたが、アルゴリズム上のちょっとしたブレークスルーにより0.01秒達成。
http://paiza.jp/poh/paizen/result/8bd4048d0f2d21fb9512df541404b45b

テストケースの作成

前回の経験から、最大規模($H=W=300, N=HW$)のテストケースを用意します。性格の異なる次の6種類を使いました。最初から先を見越して全部用意したのではなく、C以降は後から作りました。

A. 高さ1で幅がランダムな0と1を交互に配置
B. ランダムなサイズの0の長方形をいくつもランダムに配置
C. 高さ1の横ストライプ150本
D. 幅1の縦ストライプ150本
E. 格子（CとDの論理積）
F. すべて0

初期アルゴリズム

とりあえず各テストケースの正解出力を用意しておかないとデバッグできないので、愚直なアルゴリズムでプログラムを作成します。とはいえ、与えられた長方形サイズごとに、画面内の空き領域を全探索していたのでは、全体として計算量は$O(H^3W^3)$となり、たぶん終わるのを待っていられないので、少しだけ工夫しました。

入力時点で、各箇所から上方向に何個連続して隙間があるかを調べてその座標に保存しておきます。続いて、各長方形サイズについて、高さTの隙間が連続してS個以上横に並んでいる回数を各行でカウントすることにしました。計算量は$O(H^2W^2)$です。

TとSの記号の意味を、問題文と逆に使ってしまいました。
X,Yと同じ順序でS,Tと使うものと決めつけていたもので。
というわけで、以下、そのつもりで読み替えてください。

また、これも前回の経験からですが、入力処理ではどうせscanf()は使いものにならないので、fread()だけして、空白や改行も決め打ちし、解析はすべて手作業で行うことにしました。strtol()すらも呼びません。

実行時間は、テストケースにより差があり、手元のマシンで（以下同様）3.9〜6.0秒。最初としてはまずまずですが、これではpaizaに提出しても時間内の通過は無理でしょう。しかし、この段階で制限時間3秒が目前となると、問題規模がちょっと小さすぎるような気がします。

以下に肝部分のコードを示します。定数やグローバルな変数は最終版と同じです。t, sは長方形サイズ、y, xは座標です。配列インデックスはわかりやすいように1オリジンで使っていますが、これには番兵を置きやすいという効果もあります。

void
countcandidates(void) {
  int y, x, t, s;
  int count;
  int ncontinuouscols;
  for (t = 1; t <= MAX_T; t++) {
    for (s = 1; s <= MAX_S; s++) {
      count = 0;
      for (y = t; y <= H; y++) {
        ncontinuouscols = 0;
        for (x = 1; x <= W; x++) {
          if (nzerosabove[y][x] >= t) {
            ncontinuouscols++;
          } else {
            ncontinuouscols = 0;
          }
          if (ncontinuouscols >= s) count++;
        }
      }
      counts[t][s] = count;
    }
  }
}

出力処理の改善

実際はもっと後の段階で実施したことですが、各処理のプロファイリングをすると、90000件をprintf()する出力処理だけで0.02秒かかっていることがわかりました。後々、0.01秒の攻防になるかもしれないので、入力と同様に手作業で出力バッファを完全に作成し、fwrite() 1回で出力してしまうことにしました。

次回からは、最初から、入力も出力も今回のコードを使い回すつもりです。

アルゴリズムの改善

カウント対象の長方形幅を領域サイズで限定

長方形のサイズを決めてから、それが何箇所に入るかをカウントしていては、ループ回数を減らせる気配がありません。

高さTの長方形を考えたとき、ある箇所から左に高さT以上の隙間がS個並んでいるなら、その箇所には幅1からSまでの長方形だけが配置できます。Sより大きい長方形は考慮する必要がないわけです。Tが小さければSは大きいこともあるでしょうが、Tが大きいときはSはほぼ0であることが期待されます。

最悪計算量のオーダーは変わりませんが、テストケースFの場合でも、すべての箇所ですべての長方形の幅を調べるわけではないので、ループ回数が半減します。

結果、0.02秒から0.03秒、大きな隙間のあるテストケースBが0.11秒、最悪ケースのテストケースFでは2.15秒でした。これなら提出すれば100点取れるレベルでしょう。

コードは以下のとおりです。

void
countcandidates(void) {
  int y, x, t, s;
  int ncontinuouscols;
  for (t = 1; t <= MAX_T; t++) {
    for (y = t; y <= H; y++) {
      ncontinuouscols = 0;
      for (x = 1; x <= W; x++) {
        if (nzerosabove[y][x] >= t) {
          ncontinuouscols++;
          for (s = 1; s <= ncontinuouscols; s++) {
            counts[t][s]++;
          }
        } else {
          ncontinuouscols = 0;
        }
      }
    }
  }
}

配置可能最大幅が確定するまでカウントを延期

高さT以上の隙間が幅Sだけ続きそれ以上の幅にできないとき、その隙間領域には、高さT幅Sの長方形は1箇所、幅S-1の長方形は2箇所、…、幅1の長方形はS箇所に配置できます。このことから、高さT以上の隙間が続かないことが確定した後、これらの配置箇所数を加算していけば、1ずつちまちまカウントアップするより回数が減ります。

コードは以下のようになります。見かけ上、forループが4重になっていますが、最深部のループはxが1からW+1まで変化する間に、最大で合計W回しか回らないので、計算量のオーダーは$O(H^2W)$となります。実行時間は、どのテストケースも0.02〜0.03秒となりました。

void
countcandidates(void) {
  int y, x, t, s;
  int delta;
  int ncontinuouscols;
  for (t = 1; t <= MAX_T; t++) {
    for (y = t; y <= H; y++) {
      ncontinuouscols = 0;
      for (x = 1; x <= W+1; x++) {
        if (nzerosabove[y][x] >= t) { // use sentinel nzerosabove[y][W+1]==0
          ncontinuouscols++;
        } else {
          delta = 1;
          for (s = ncontinuouscols; s > 0; s--) {
            counts[t][s] += delta;
            delta++;
          }
          ncontinuouscols = 0;
        }
      }
    }
  }
}

ここまではわりとすぐに到達したのですが、ここから先に手間取りました。

ベクトル加算問題への転換

ループカウンタを逆に回して0と比較できるようにしたり、番兵的なものを置いたり、手動でループアンローリングしてみたりと、様々な小手先技を駆使してもいっこうに改善しませんでしたが、一つ光明が見えたのが、上のコードの本質が三角形をいくつも加算していることに気付いたときでした。

つまり、高さtの長方形について、高さt以上最大幅3の隙間領域が一つあれば、[3 2 1 0 0...] という三角形状のベクトルをcount[t]に加算しています。最大幅2の隙間領域が一つあれば、[2 1 0 0 0..] を加算します。ということは、高さT以上で最大幅Sの隙間領域が何個あるかだけあらかじめカウントしておいて、後からこれらの三角形状のベクトルを再現して足していくこともできます。計算量のオーダーは変わりませんが、とりあえず最深部のループを追い出せれば何かできるかもしれない、あわよくば実行時間半減もありうると期待して実装しました。

以下がそのコードです。このカウントを入れる配列は別途用意した方がリーダブルになりますが、キャッシュの効果などを考えて結果を入れる配列をそのまま使いました。実行時間は内部ループがなくなり、候補箇所のカウントも長方形サイズごとに1回しか行わなくなったので多少減り、どのテストケースも0.02秒で実行できました。

void
countcandidates(void) {
  int y, x, t, s;
  int delta;
  int ncontinuouscols;
  for (t = 1; t <= MAX_T; t++) {
    for (y = t; y <= H; y++) {
      ncontinuouscols = 0;
      for (x = 1; x <= W+1; x++) {
        if (nzerosabove[y][x] >= t) { // use sentinel nzerosabove[y][W+1]==0
          ncontinuouscols++;
        } else {
          if (ncontinuouscols > 0) {
            counts[t][ncontinuouscols]++;
            ncontinuouscols = 0;
          }
        }
      }
    }
    delta = 0;
    for (s = W; s > 0; s--) { // must be "downto"
      delta += counts[t][s];
      counts[t][s] = counts[t][s+1] + delta; // use sentinel counts[t][W+1]==0
    }
  }
}

ブレークスルー直前

ここからは完全に手詰まりでどうしようもなくなり、なかばやけぎみに、ループの手順を入れ替えたりしてみました。つまり、tによるループを内側に入れてみただけです。そのためには、ある高さTの列が何回連続したかを表すカウンタが、Tの個数だけいるので、その配列アクセスの手間の分だけ遅くなる可能性があります。

以下がそのコードですが、カウンタへのアクセス回数の多いテストケースDとEで0.05秒になってしまいました。

void
countcandidates(void) {
  int y, x, t, s;
  int delta;
  int ncontinuouscols;
  for (y = 1; y <= H; y++) {
    for (x = 1; x <= W+1; x++) {
      for (t = 1; t <= y; t++) {
        if (nzerosabove[y][x] >= t) { // use sentinel nzerosabove[y][W+1]==0
          runcounter[t]++;
        } else {
          if (runcounter[t] > 0) {
            counts[t][runcounter[t]]++;
            runcounter[t] = 0;
          }
        }
      }
    }
  }
  for (t = 1; t <= MAX_T; t++) {
    delta = 0;
    for (s = W; s > 0; s--) { // must be "downto"
      delta += counts[t][s];
      counts[t][s] = counts[t][s+1] + delta; // use sentinel counts[t][W+1]==0
    }
  }
}

最終版

ところがこの改造が奏功します。

すべてのTについて毎回カウンタを増やすのではなく、その箇所の隙間の高さT分だけカウンタを増やすことにし、Tが直前の隙間の高さT'より減ってしまったときは、T'〜T-1までの高さの隙間領域がそこで終わったことになるので、その範囲の高さ分のみ一気に処理できるようになります。

このアルゴリズムでも最悪計算量のオーダーは変わりません。隙間領域の高さが頻繁に増減し、かつ減少時の幅が大きい場合が最悪で、テストケースDとEが該当します。しかし現実にはそのような隙間の状況はまずないし、paizaのテストケースにそのような状況が含まれている可能性は低そうです。通常の状況では$O(H^2W)$に対してかなり小さい定数倍の計算量になっていると思われます。

コードは最後にまとめて示します。実行時間はDとEだけ0.04秒、他はみごとに0.00秒となりました。（FreeBSDの/usr/bin/timeは小数点以下を tv_usec/10000 で計算しているので0.01未満切り捨て）

これをpaizaに提出したところ、全テストケース0.01秒という結果が出ました。

まとめ

論点がぼやけるので書きませんでしたが、コンパイラの差異や、2次元配列の構成方法の差異、2次元配列アクセス方向によるキャッシュ容量と連想度の影響の有無など、技術的に興味深い検討がいろいろできました。これは大きな収穫だったと思います。

結局のところ、計算量のオーダーは$O(H^2W)$から下げられませんでした。わりと満足している反面、まったく違うアプローチでより高速なアルゴリズムがありそうな気もします。もっと大きなテストケースで他の解答と比べてみたいところです。

コード

以下に最終コードを示します。

# include <stdio.h>

# define MAX_H 300
# define MAX_W 300
# define MAX_T 300
# define MAX_S 300
# define NCOLS 512 // > MAX_W, > MAX_S, > MAX_H, 2^n
# define MAX_N (MAX_H*MAX_W)

# define MAX_INPUT_BYTES (3+1+3+1+(MAX_W+1)*MAX_H+5+1+(3+1+3+1)*MAX_N)
# define MAX_OUTPUT_BYTES ((5+1)*MAX_N)

int H;
int W;
int N;
int T[MAX_N];
int S[MAX_N];

// nzerosabove[y][x]: holds how many zeros have continued vertically
// at the position.
int nzerosabove[MAX_H+1][NCOLS];

// runcounter[t]: counts how many vertical spaces with height t or less
// continues horizontally.
int runcounter[MAX_T+1];

// counts[t][s]: holds how many continuous maximal runs of width s
// that have height t.
// also converted to the final results.
int counts[MAX_T+1][NCOLS];

void
parseinput(void) {
  char buffer[MAX_INPUT_BYTES+1];
  int bufindex, y, x, i;
  int nread;

  nread = fread(buffer, 1, MAX_INPUT_BYTES, stdin);
  buffer[nread] = '\n'; // make sure input is terminated by LF.

  for (bufindex = 0; buffer[bufindex] != ' '; bufindex++)
    H = H*10 + buffer[bufindex] - '0';

  for (bufindex++; buffer[bufindex] != '\n'; bufindex++)
    W = W*10 + buffer[bufindex] - '0';

  for (y = 1; y <= H; y++) {
    for (bufindex++, x = 1; x <= W; bufindex++, x++) {
      if (buffer[bufindex] == '0') {
        nzerosabove[y][x] = nzerosabove[y-1][x] + 1;
      }
    }
  }

  for (bufindex++; buffer[bufindex] != '\n'; bufindex++)
    N = N*10 + buffer[bufindex] - '0';

  for (i = 0; i < N; i++) {
    for (bufindex++; buffer[bufindex] != ' '; bufindex++)
      T[i] = T[i]*10 + buffer[bufindex] - '0';
    for (bufindex++; buffer[bufindex] != '\n'; bufindex++)
      S[i] = S[i]*10 + buffer[bufindex] - '0';
  }
}

void
countcandidates(void) {
  int y, x, t, s;
  int delta;
  int lastnz, nz;

  lastnz = 0;
  for (y = 1; y <= H; y++) {
    for (x = 1; x <= W+1; x++) {
      nz = nzerosabove[y][x]; // use sentinel nzerosabove[y][W+1]==0
      for (t = lastnz; t > nz; t--) { // must be "downto"
        counts[t][runcounter[t]]++;
        runcounter[t-1] += runcounter[t]; // runcounter[0] is dummy
        runcounter[t] = 0;
      }
      runcounter[nz]++;
      lastnz = nz;
    }
  }

  for (t = 1; t <= MAX_T; t++) {
    delta = 0;
    for (s = W; s > 0; s--) { // must be "downto"
      delta += counts[t][s];
      counts[t][s] = counts[t][s+1] + delta; // use sentinel counts[t][W+1]==0
    }
  }
}

void
buildoutput(void) {
  char buffer[MAX_OUTPUT_BYTES+1];
  int bufindex, i, digitindex, number;
  char digits[20];

  bufindex = 0;

  for (i = 0; i < N; i++) {
    number = counts[T[i]][S[i]];
    digitindex = 0;
    do {
      digits[digitindex] = number % 10 + '0';
      number /= 10;
      digitindex++;
    } while (number > 0);
    for (digitindex--; digitindex >= 0; digitindex--) {
      buffer[bufindex++] = digits[digitindex];
    }

    buffer[bufindex++] = '\n';
  }

  fwrite(buffer, 1, bufindex, stdout);
}

int
main(void) {
  parseinput();
  countcandidates();
  buildoutput();
  return 0;
}