311の昨日に、AOJとの愛称で呼ばれ、日本国内では有名なプログラミングオンラインジャッジシステム(Aize Online Judge)で解く問題数ランキングで3位になった。昨年の7月末に参戦し、7ヶ月半で達成した。つまり、平均して一日5問解いた計算になる。

UVaでの教訓として、今回は問題をただ解くだけでなく、なんとか各問題の計算スピードが1位になるように頑張った。0:00 秒となった解答は提出された時間順で並べられ、それ以上の順位に食い込むことはできないが、それ以外の問題についてはチャレンジチャンスはいくらでもあるから。

ほぼC言語一本で解いてきたので、C++やJAVAに比べて不利な面は否定できない。ただ、C言語の自由度、スピード性は21世紀の今日になってもその特徴は生きているところも事実だと再認識した。

1位になるか。未解答問題の多くは難しく、ここからさらに1位との差である116問を解くのに大きな山はいくつもあると思うが、根性という点では誰にも負けないことを期待しよう。

昔のUVa時代と違って、競技プログラミングの重要性は広く認知されているし、アルゴリズムの素晴らしさに感動したひとも多くなった。小中学校生のプログラミング必修化は競技プログラミング愛好家にとっても喜ばしい状況ではある。

そういえば、UVaも最高3位(2006年11月21日に達成)だった。その成績は日本からの挑戦者記録としていまも生きている。

計算速度が至上命令か 。本問題のプログラム実行制限時間はついに0.666秒に設定されてしまった。自分としては嫌なんだね、はっきりと。

アルゴリズムの優劣性については、実行速度がとても大事だが、可読性、そこからの発展性、理解しやすさ、など場面や用途によって様々なものがあっていいはず。ソートアルゴリズムについても、O(N)のものがあっても、ソートの本質を理解する上では、様々なものが用意されている。実行制限時間をきつく設定したら、アルゴリズムの強制採用にほかならない。大反対だね。

1. 問題文
n (1<=n<=10000) 個のデータ (整数、10進8桁以内、負数もありえる)が与えられている。さらに、様々な値m (1<=m<=n) が与えられ、m個のLISシーケンス(Longest Increasing Subsequence、最長単調増加シーケンス)をデータからピックアップして表示せよ、というのが問題の意味。なお、LISシーケンスはなるべくデータの先頭にあるものを選択しないといけない。

2. 実例
例えば、以下のn=10個のデータがあったとする。m=1,2,3,...10のLISシーケンスはそれぞれ以下の通り。

-3 19 12 -7 18 23 -5 0 11 22

m=1 LIS=-3

m=2 LIS=-3 19

m=3 LIS=-3 19 23

m=4 LIS=-3 12 18 23

m=5 LIS=-7 -5 0 11 22

m>5 LISは存在しない

3. 計算法
O(nlogn)バージョンのLISアルゴリズムを使う。ブログに紹介があるので、参照されたい。

つまり、n個のデータに対するそれぞれのLIS長をまず算出する。上の実例では、次になる。

data (-3 19 12 -7 18 23 -5 0 11 22)

LIS (4 2 3 5 2 1 4 3 2 1)

最長LISは5という結果も上のLISテーブルから分かる。

つぎに、それぞれのmに対して、そのLISシーケンスをピックアップする。

m=1 LIS1=4 >= 1 なので、data1=-3がその結果。
m=2 LIS1=4 >= 2、LIS2=2 >= 1 しかも data2=19 > data1なので、-3 19が結果になる。
m=3 LIS1=4 >=3、LIS2=2 >=2 しかも data2 > data1、data2=19がLISシーケンス2番目のものになる。さらに、LIS3=3 は 1以上だが、値はdata2以下なので、却下。LIS4, 5も同様の理由で却下。LIS6=1が1以上だし、値data6=23は19よりも大きいので、3番目はdata6になる。
m=4については、data2=19のLIS2は2なので、3よりも小さいことから、2番目にはならない。隣のdata3=12が2番目になる。つぎの12よりも大きく、LISの値が2以上のものは18なので、3番目になる。4番目は18よりも大きく、LISの値が1以上のもの、つまり、23となる。
m=5については省略するが、先頭からまずLISの値が5よりも大きいものを探し、-7を取る。つぎに-7よりも値が大きく、LISの値が4よりも大きいものを取る。

全体の実行時間は、LISテーブルの算出にO(nlogn)、LISシーケンスの抽出にO(n)でできるので、全体でもO(nlogn)になる。

4. 結果
入力部分を工夫して、最終的にランク1位の実行速度と同じにした。出力部分はprintf()をつかったままなので、まだ余力は残っているが、誤差範囲内の時間(0.01秒)を縮めても意味ないので、そのままにしておいた。

ともかく、スピードが命という風潮に強く警告したい気持ちでいっぱい。いくらなんでもやり過ぎだな。

p11271.jpg

XY平面に無限に広がる抵抗ネットワークに対し、任意の両点(格子点)間の抵抗値を求める問題。自分は大学では電子工学を勉強したのに解けない。情けないことかな。

1. 問題
任意の両点だが、無限に広がっているので、原点をどこに取っても同じはず。そこで、両点の片方を原点にし、原点からの抵抗値を求めるというのが本問題。なお、1Ωとは隣接両点間の線分の抵抗値。

ネット上を調べたら、(0,0)と(1,0)間の抵抗値は0.5Ωとなる解説は多いが、ほとんどが事後解釈の領域を出ない。(0,0)と(1,1)間の抵抗値が分かるのであれば、はじめて本物と認めてあげたい気持ちになる。

2. 計算法
自分の知識では全く歯が立たないので、一生懸命ネットを調べた。すごいことにそれを計算した論文をついにみつけた。

(PDF) Infinite resistive lattices

細かい中身はご自身でみてもらうことにして、必要な計算式と計算結果の一部を紹介しておく。

p11271-1.jpg
p11271-2.jpg

変数nとpはそれぞれy座標、x座標と考えてよい。対称性があるから、入れ替えても計算結果は同じになるはず。

p11271-3.jpg

上のテーブルには (5,5)までの抵抗値がリストアップされている。例えば、(0,0)-(1,1)間の抵抗値は2/Phi=0.637Ωになる。円周率Phiがここに出てくる。

ということで、論文を見る限り、数値積分をしないといけなさそう。実際にプログラムを組んで計算してみたら、上記のテーブルに近い値が出てきた。ただ、積分部分のプログラム(シンプソン法)は精度が良くないせいか、xとyの対称性に若干の誤差があったりする。UVaの自動判定システムにもAcceptedされていない。

単純そうに見えた計算問題に大変複雑な定積分が出てくるとはびっくる。漸化式では無理かな。

でも、問題文のヒントがとても気になる。もっとうまい計算法はあるような言い方。
Hint: There is a surprising Dynamic Programming solution, but how do you get it to fit under the memory requirement? 🙂 .

最後に、プログラム内に使っている、計算式に関連するコードをピックアップしておく。

int x, y;  /* 計算したい (x,y) 座標 */
double alpha(double beta)
{
    double t = cos(beta);
    return log(2-t+sqrt(3+t*(t-4)));
}
double func(double beta)
{
    double a = alpha(beta);
    val = (1.0-1/exp(x*a)*cos(y*beta))/sinh(a);
    return val;
}

上のやり方では問題ないように思われるが、確信はまだない。Acceptedされていないし。ちなみに、現時点では問題を解いたひとは二人だけ、出題者本人とUVaシステムの関係者?

新しいUVaサイトがよく落ちるね。広告を入れるのはしょうがないが、もう少しセンスがあっていいじゃないかな。無駄なスペースを取りすぎ、とくに垂直方向。肝心な問題文をスクロールしてみるという有様。IEもfirefoxも開けない。なんとかしてほしいよ。

今日から使えたGoogle CromeはSafariと同様、ちゃんとアクセスできてよかった。

#11281番の問題はよく勘違いされるようで、自分も最初では間違っていた。やはりちゃんと図形を描いてみることだね。

1. 問題文
丸い円をいくつも与えられて、さらに様々な形の三角形が与えられている。円に入る三角形を教えろ、というのが問題の意味。実際の問題では、円の代わりに孔を使ったり、三角形の代わりに立体の三角プリズムを使っているが。

2. 考え方
三角形の外接円の半径を求め、円の半径と比較すれば簡単に問題が解けそうに思われるが、大きな落とし穴が待ち構えている。

たとえば、下の三角形をよくみてみよう。ブルーの円は外接円だけど、赤の円よりは明らかに半径が大きい。つまり、鈍角三角形に対しては、その最長辺が最小円の直径になる。外接円では全然大きい。鋭角三角形に対しては、外接円が最小円になるけど。

p11281.jpg

上のことを理解すれば、問題は簡単に解けるはず。

1. 問題
高速性が要求される問題のひとつ。N (3<=N<=5000) 個の入力データから、任意の3つ (x, y, z) を選んで、x+y=z となるTripleの数を求めよ。
なお、入力データはすべて正の整数、32ビット範囲内に収まる。プログラムの実行制限時間が8秒。

2. 実例
次の入力データをみてみる。
1 2 3 4 5 6

1+2=3, 1+3=4, 1+4=5, 1+5=6, 2+3=5, 2+4=6 になるので、答えは6になる。

3. 計算法
入力データの個数が5000までなので、一見簡単に解けそうだが、やってみて、大変さが分かった。出題者の意図は O(n^2) アルゴリズムの提出にある。

自分の順位は現時点では18位なので、まだまだ改善するところを考えないといけないが、つぎのような方法でやってきた。

入力したデータに対し、ソートを行い、グループ化を行う。グループ化というのは、同じ値の入力データをまとめることだ。たとえば、1 1 1 2 2 3 3 という入力データを (1,3) (2,2) (3,2) にクループ化する。ソートにはO(nlogn)、グループ化には O(n) で済むが、ソートの高速化は問題として残っている。

次に、ソートされ、グループ化されたデータに対し、Tripleの数を算出する。

言葉で説明するのが面倒なので、ソースプログラムを使う。勿論分かりやすくするため、一部を書き換えているし、もっと高速化できるところもある。この部分の実行時間はO(n^2)。なお、ans変数は64ビット整数(long long int)でないとオーバーフローになってしまう。

typedef long long llong;
#define SIZE  5000
int N;
unsigned data[SIZE+5];  /* 入力データ(ソート済み) */
int freq[SIZE+5];  /* 同じ値を持つ入力データの個数 */
ans = 0;
if (N == 1) { puts("0"); exit(1); }
for (k = 1; k < N; k++) {
    i = 0; j = k-1;
    while (i < j) {
        if (data[i]+data[i] == data[k]) {
            ans += (llong)freq[i]*(freq[i]-1)/2*freq[k];
            break;
        }
        if (data[j]+data[j] == data[k]) {
            ans += (llong)freq[j]*(freq[j]-1)/2*freq[k];
            break;
        }
        if (data[i]+data[j] == data[k]) {
            ans += ((llong)freq[i])*freq[j]*freq[k];
            i++, j--;
        } else if (data[i]+data[j] &lt; data[k]) i++;
        else j--;
    }
}
printf("140733424702768\n", ans);

4. スピードの改善
スピードを上げるために、データのソートとグループ化に、2分探索木を使ってみた。順位は14に上がったが、劇的な改善にまだ繋がっていない。最後のtriple数の計算がボトルネックになっているから。

2分検索(バイナリサーチ)のような手法が使えるのではないか、そんな気がしている。ランク1位の時間は0.540秒なので、自分としてはせめて1秒を切りたい。

元問題

余計なことは書かない。一生懸命導いた数式の結果だけをメモしておく。

p11250.jpg

あとは、入力データのm, nに対し、それぞれ多倍長整数(本ブログ参照)を使って、上記の通り、分子、分母に分けて計算し、最後に、ユークリッド互除法で分子、分母を約分し、結果を得る。

数式のネット検索は、Googleといえども、まともにできない。入力の仕方にも、使う変数の文字にも、可変要素が多すぎるかな。

平方数とは平方根が整数の数のことだ。小さい順で並べていくと、0, 1, 4, 9, 16, 25, 36, ... となる。あるいは、素因数分解して、すべての素因数の指数が偶数である場合も平方数になる。たとえば、36 = 22x32、素因数の指数がそれぞれ2であるので、36は平方数といえる。平方根を計算することは結構時間がかかるので、素因数が分かるのであれば、素因数の指数を調べるほうがはるかに高速。

1. 問題
 n (0<n<=200000) 個の整数(絶対値は1018よりも小さいが、すべての素因数は30以下。)が与えられた時に、n個の中から任意の2つを選んで、その積が平方数となる(組合せの)数Xを、n個の中から任意の3つを選んで、その積が平方数となる(組合せの)数Yを求め。

2. 実例
 例をもって説明する。
 たとえば、3つの整数 2, -2, 2 が与えられた時、2x2 = 4 になるので、X = 1。また明らかに Y = 0。
 簡単だな、と思うかもしれないが、データの個数が20万、個々のデータの値の大きさを考えると、そう容易く解けないはずだ。しかも、実行時間は3秒に制限されている。
 20万個のデータすべてをここに書き出すことは非現実的だが、20万個のデータが何百組も与えられても、3秒以内に計算を終了させ、正解を出すプログラムを作ることが簡単ではない。すべては高速性にかかってくる。

下は多少複雑な例。それでもデータの個数は100でしかない。

答えは、X=4950 Y=161700。

3. 計算法
 まず、1~2500000の間の、素因数が30以下のすべての値に対し、エラトステネスのふるいを応用して、素因数分解をしておく。2500000を超えた入力値に対しては、単純に1つ1つ素因数分解をしておく。
 ただ、本問題では素因数の指数の奇偶性が問題になるので、個々の素因数に対応して、以下のビットパターンを採用する。

 最下位ビットから bitpat[11] = { マイナス, 2の数&1, 3の個数&1, 5の個数&1, 7の個数&1, 11の個数&1, 13の個数&1, 17の個数&1, 19の個数&1, 23の個数&1, 29の個数&1 };

 「マイナス」はデータが負数を表すビット、
 「2の数&1」は素因数2の指数が奇数なら1、偶数なら0になる。
 「3の数&1」は素因数3の指数が奇数なら1、偶数なら0になる。
などを意味する。

 例えば、6のビットパターンは110、-6のビットパターンは111、10のビットパターンは1010、平方数のビットパターンは0。

 データの積を取ることは、すなわち、ビットパターンの XOR 計算にほかならない。そのXORの結果がゼロになれば平方数、ゼロ以外であれば非平方数となる。

 また、入力データを値がゼロのものとそうでないものに分ける。値がゼロのデータの数がzだとすると、XとYに相当する部分は以下のように算出する。

  X = nC2 - n-zC2 = z(2n-1-z)/2
  Y = nC3 - n-zC3

あとは、値がゼロ以外のデータについてのみ、XとYの値を算出し、上の値に足し加えておけばOK。

4. 結果
 書くのは嫌になるほど、多くのことをトライして、最終的に実行時間を0.480秒に押さえ、現時点でのランク1位にできた。

細かいことは省略するが、ヒントになるものをリストアップしておく。

  1. 入力の値に対するビットパターンは2047以下。つまり、異なる入力データの数は2048通りしかない。
  2. XOR計算 a^b = 0 となるa, bは、a = bに限る。a^b^c = 0となるa, b, cは a^b = c に限る。従って、異なる2048個の値から、任意の2つを取り出してXORがゼロになる計算は2048回(O(n))で済むし、任意の3つを取り出してXORがゼロになる計算は2048x2048回(O(n2))で済む。
  3. 入力データを平方数と非平方数に分けた。平方数同士では、XORがゼロになる。平方数と非平方数とのXORはゼロにならない。非平方数と非平方数とのXORがゼロになるのは、その2つの非平方数が等しい場合に限る。
  4. 入力データの処理にscanf()ではなく、read()関数を使った。getchar()を使ってランク2にできたが、1位になるには、read()を使ってしまった。
  5. 計算のオーバーフローに気をつけよう。整数同士の積に64ビット整数を使おう。

5. ソースの一部

typedef long long llong;
typedef unsigned short ushot
#define MAX    2500000
ushot bitpat[MAX+5];
ushot pp[10] = { 2,3,5,7,11,13,17,19,23,29 };
/* 1~MAXまでの各整数に対するビットパターン */
void sieve(void)
{
    register int i, j;
    int k, p;
    ushot b;
    for (k = 0; k < 10; k++) {
        p = pp[k], b = 1 << (k+1);
        for (j = p; j <= MAX; j *= p) {
            for (i = j; i <= MAX; i += j) {
                if (bitpat[i] & b) bitpat[i] &= ~b;
                else               bitpat[i] |= b;
            }
        }
    }
}
/* MAXを超えた値に対するビットパターン */
ushot factor(llong n)
{
    int k, p, m;
    ushot ans;
    ans = 0;
    if (!(n & 1)) {
        m = 0;
        do n >>= 1, m++;
        while (!(n & 1));
        if (m & 1) ans |= 2;
    }
    for (k = 1; n > 1 && k < 10; k++) {
        p = pp[k];
        if (!(n  0x7fff0dc99530)) {
            m = 0;
            do { m++; n /= p; }
            while (n  0x7fff0dc99530 == 0);
            if (m & 1) ans |= 1 << (k+1);
        }
    }
    return ans;
}

XとYを算出する部分(分かりやすくするため、ソースの一部を書き直した。ソースそのものとは異なる)。なお、freq[a]はビットパターンがaである入力データの個数を示す。

X = Y = 0;
for (i = 1; i < 2048; i++) {
    if (freq[i]) {
        X += ((llong)freq[i]*(freq[i]-1))/2;
        for (j = i+1; j < 2048; j++) {
            if (freq[j]) {
                if ((i^j) > j && freq[i^j]) Y += (llong)freq[i]*freq[j]*freq[i^j];
            }
        }
    }
}

1. 問題

p11200.jpg

上記の通り、各セルには必ず、右下あるいは左下斜めの壁が置かれる。迷路の入口は上、出口は下。左右両側は通行禁止。こういう条件の下で、入口から出口までの異なる通路は果たして何本あるか。ちなみに、右上の迷路では通路が3本、左下の迷路では通路が1本という結果になっている。

2. 入力データの例
 上記の左下迷路に対応する入力データは以下の通り。大きくなると、目視での確認は相当きつい。¥は本来はバックスラッシュ\、右下斜めを表しているはずなのに。

3. 計算法

p11200-2.jpg

セル(r,c)に位置することは、そのセルの下の端の中央にいることと考える。そうすると、周りの壁の形によって、上記の8通りの進み方がある。

証明はまだしていないが、入口1箇所に対応する出口は高々1箇所。従って、1つの入口に入って探索していき、出口ひとつに出くわしたら、それ以降の探索は打ち切り、つぎの入口から探索し続ける、という戦略を取って問題なさそう。また、異なる入口から入り、最終的に同じ出口にいく通路もないと推測する。

1. 問題
 様々なサイズの長方形(または正方形)の紙Aが与えられている。1枚の紙Aを、折り紙するための4枚の正方形Bに切り分けたいけど、もっとも大きな折り紙Bが得られるような、最適なAがどれかを選べ。

2. 実例
 以下の実例を見てみる。

1枚目の紙10x20では最大5x5の折り紙しか得られない。2枚目の紙40x8では最大8x8の折り紙が得られる。3枚目の紙12x12では最大6x6の折り紙しか得られない。ということで、2枚目の紙を選べば、最大サイズの折り紙を手にすることができるわけだ。

3. 計算法
 縦幅に比べ、横幅の短くない長方形(または正方形。そうでなければ90°回転すればOK)に対して、最大サイズの折り紙(正方形、同一サイズ)4枚は3通りの並び方にしかならない。

p11207.jpg

長方形の縦幅サイズをY、横幅サイズをXとすると、3パターンによって得られる折り紙の最大サイズはそれぞれ、

 パターン1: X/4 と Y の小さいほうの値
 パターン2: X/2 と Y/2 の小さいほうの値
 パターン3: X/3 と Y/2 の小さいほうの値

になる。3パターンに対応する3つのサイズから、最大サイズを得、それが長方形 X, Y に対して入手できる折り紙の最大サイズになる。

さらに、与えられたすべての紙の中から、最大サイズのものをピックアップすれば、答えが得られる。

ということで、本問題は簡単に解ける。残りは実行時間の短縮ということだ。

自分の書いたプログラムは最初0.020秒だったけど、最適化を施し、最終的に0.000秒に到達した。つまり、乗算、除算しないで、ビットシフトでの計算を心がけた。

例えば、元のサイズX, Y を最初から12倍すれば、2, 3, 4で割り切れるので、実数ではなく、整数だけで計算可能になる。

となると、3つのパターンはそれぞれ

 パターン1: min(3X, 12Y)
 パターン2: min(6X, 6Y)
 パターン3: min(4X, 6Y)

に変わり、ビットシフトの性質を考えて、

 3X = (X << 1) + X  4X = X << 2  6X = 3X << 1
 3Y = (Y << 1) + Y  12Y = 3Y << 2

などと工夫すれば、除算も乗算も必要なくなる。

まあ、そんなテクニックは実行スピードの本質的な改善に繋がらないから、分からなくても構わないが、少しでも他人よりもランクを上げたい気持ちがあれば、覚えて損することはないだろう。

1. 問題
 使えるカセットテープの長さ(録音時間)と、録音したい曲のリストが与えられ、すべての曲が入るだけの、最も短いカセットテープを探せ。

2. 実例
 実例を見たほうが分かりやすい。

出力結果を見てみよう。

3. 考え方
 詳しいことは問題文に述べてないので、常識的に考えて問題を解いてみた。
 一番録音時間の長いカセットテープは120分だと思う。それはA面+B面での録音時間。片面ではその半分しか録音できない。また、曲の数は最大30とした。

本問題にBFS(幅優先探索)を使った。

構造体
 typedef struct {
  char no; /* 曲の通し番号 */
  char side; /* A面かB面か */
  int a, b; /* A面とB面の録音時間 */
  int prev; /* 前へのポインタ */
 } QUE;
と定義し、1曲目はA面録音として、2曲目以降はA面に、またはB面に録音するように探索キューを展開していく。

つまり、n曲目においては、検索キューに存在しているn-1曲目までの組合せすべてに対し、次の操作を経て、再びキューに入れる。

すべての曲の探索キューを作成した後、短いテープの順に、A面もB面も録音時間に収まるものが探索キューにあるかどうかを調べる。

4. 感想
 DPによる方法も考えてみたが、配列はとてつもなく大きいので断念した。今回のBFS探索法は強引的な面もあったが、実行時間がゼロ秒になっているので、本問題に使えることは間違いなさそう。