プログラムの高速化について書きます。

UVaに現時点では2千題近いプログラミング問題が出題されています。プログラミングを楽しみながら稼ぎたいひとはTopCoderというサイトを利用するのもいいでしょう。

書いたプログラムの良し悪しは評価基準それぞれですが、プログラムの正確性、高速性、信頼性、可読性、メンテナンスのしやすさ等がよく評価基準に用いられると思います。

ここではもっとも比較しやすいプログラムの高速性を考えます。

2つのプログラムを同時に走らせて、先に終了するほうが速い。そういう観点からすると、高速性の比較がやりやすく、UVa でもランキングを決めるパラメータとして採用されています。

ちなみに、プログラムの正確性を判定するのは難しく、一般的には判定不能と証明されています。UVaでも言い方として「正解」というのではなく、「Accepted=受け入れられた、受理された」というレベルでわれわれの送ったプログラムを評価しています。(後にいくらでも実は前の判定が間違えていたとかの言い訳にするための逃げ道だというひともいますが。)

ほかに比較しやすいパラメータはメモリの使用量、ソースプログラムの長さ、等もあります。北京大学のPKU (acm.pku.edu.cn) ではソースの長さが競っていますが、あまり感心しません。ソースプログラムの可読性が著しく損なわれる恐れが大いにあると思います。1バイトでも短くしたいがために、変数名を1文字にしたり、スペースを全くソースに入れないことを平気でやってしまいますから。

UVaに参戦して2年間、つい3位になりました。

1位、2位を越す実力はないので、この辺で引退したいと思ったりします。あるいは、方針を変え、暇つぶしに付き合っていきたいと思います。

辛かった時期もありましたが、すべては努力の賜物だと感謝します。

1750問をそれぞれ2ページで解説することにしても、ものすごいボリュームになりそうです。いつかまとめてみます。

前の記事の続き.

<解けない盤面の存在>
 初期盤面は作り方のよって,ゴール状態に移動できないものも数多く(正確にはちょうど半分)存在する.15パズルに関し,懸賞問題も歴史的にあったらしいが,結局それが解けないことは後ほど証明されてしまった.

解答可能かどうかの判定に「転倒数」をまず数える.つまり,ゴール状態に比べ,盤面の1次元配列記述において,位置が前後に転倒してしまった数を累計する.

前の記事トップの盤面を例に取って「転倒数」を計算しよう.

数字の2が,ゴール状態では数字11の前に来るはずなので,転倒数は1となる.数字8の転倒数も1.数字5は数字11と8の前に来るべきなので,転倒数は2となる.以下は残りの数字すべてを列挙する.

 数字 11(に対応する転倒数は 0),数字 2(転倒数 1),数字 8(1),数字 5(2),13(0),14(0),10(3),12(2),4(7),3(8),1(10),7(6),9(5),15(0),6(9)

転倒数の合計は 0+1+1+2+0+0+3+2+7+8+10+6+5+9 = 54.

つぎに,その合計数と空白箇所(スペース)の行番号との和が偶数かどうかを見る.偶数であれば解答可能,奇数であれば解答不可能になる.

今の例では,転倒数の合計は54,空白位置の行番号は4なので,和は58,偶数である.従って解答可能な初期盤面といえよう.

また,歴史的に有名な解答不可能懸賞問題 14-15パズルでは,数字14だけが転倒しているので,その転倒数は1,転倒数の合計も1.空白位置の行番号4と合わせると和は5となり,奇数である.従って解答不可能である.いくら頑張っても,ゴール状態に決して辿りつけない訳で,懸賞金をもらうはずはない.

われわれのプログラムでも,まず解答可能かどうかを判定し,解答不可能な入力データをはじけておかないといけない.

以下はその判定部分に関わるプログラム.関数の引数 r0 は空白の行番号(C言語では習慣として 0からカウントするので,実際の行番号よりは1だけ数字が小さい).関数 solvable() のリターン値が1であれば解答可能,0であれば解答不可能のことを示す.

int solvable(int r0)
{
int i, j;
int n;
n = 0; for (i = 1; i < 16; i++) for (j = 0; j < 16; j++) { if (layout[j] == i) break; else if (layout[j] > i) n++; } return (n+r0) & 1; /* or (n+ r0+1) & 1 == 0 */ }

(プログラムの説明) 変数 n に転倒数の合計が入る.layout は盤面の1次元配列.r0 は空白位置の行番号(0~3).

<プログラムを組み立てる>
 すべての盤面やその盤面に関係するパラメータを管理しているのは move テーブル.その構造体は以下の通り.

#define MSIZE      50001  /* テーブルサイズ */
typedef struct { char t[16];  /* 盤面 */ char r0, c0;  /* 空白位置の行番号(r0: 0~3)と列番号(c0: 0~3)*/ int prev;   /* ひとつ前の盤面 ID */ char op;    /* ひとつ前の盤面からの移動の仕方.上下左右のどれか */ char step;   /* 初期盤面からの移動回数(ステップ数) */ } MOVE; MOVE move[MSIZE+1]; int msize;

以下は解答プログラム.キタナイ書き方をしていて,申し訳ない.とくに上下左右の移動に関する部分.テーブルにすれば4方向重複しなくて済むはずだが.

C言語解答プログラム: Solution for UVa 10181 15-Puzzle Problem (C program)
計算時間: 0.078秒

計算時間のランキング表に,0.000秒のような,極限にまで高速に計算できた解答も見受けられる.IDA*アルゴリズムをつかおうかな.あるいは,もっと凄いアルゴリズムが存在しているかも.

原問題はこちら

有名な 15パズル問題 を解くプログラムを公開する.ただし,断っておきたいが,本解答プログラムはあくまでも 上記のUVa問題を解くためのものであり,時間的・メモリ的制約から,このままでは一般的な15パズル問題を解けるものではないかも.

15パズルは子供の時から遊んできた知的ゲームのひとつ.空いているスペースの周りをスライドしていき,1から15までの順にきれいに揃えるとゲーム終了.

勿論,数字の代わりに絵柄を貼り付けることも可能だが,本質は変わらない.

われわれ人間が遊ぶ時に最短移動回数(最小手数)のことについてあまり意識しないのに対して,プログラミング問題では最小手数を問われることが多い.ただ,本問題では,出題側が明確に最小手数を求めてはいないようだが,50手を越えた解答を提出してもいけない.実際のテストデータを見たわけではないが,最小手数が45手のテストデータを複数用意されると,最大5手までの誤差しか許容されなくなる.そういう意味で,本問題の解答プログラムでも,最小手数を十分に意識したつくりでないといけないだろう.

15パズルを解くアルゴリズムはいくつもあると思うが,ここではA*アルゴリズムを使ってみた.

<盤面の記述>
 盤面(レイアウト)の記述には4X4の2次元配列でもいいが,ここでは計算の効率性から1次元配列を使う.順番の数え方は左上から右下まで1行ずつとする.

例えば,記事トップ(上の写真)の盤面では
   char layout[16] = { 11, 2, 8, 5, 13, 14, 10, 12, 4, 3, 1, 7, 9, 15, 0, 6 };
になるし,ゴール状態は
   char goal[16] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 0 };
になる.ただし,スペースは 0 で表す.

なお,2次元と1次元との位置対応は以下になろう.r は行番号(0~3),c は列番号(0~3),p は1次元での対応位置(0~15)とすると,

  p = 4*r + c
  r = p / 4, c = p % 4

効率向上のため,実際のプログラムではビットシフトやビット演算を使う.

#define RC2P(r,c)  (((r)<<2)+c)
#define P2R(p) ((p) >> 2) #define P2C(p) ((p) & 3)

<高速に解くことの困難性>
 15パズルは高々4X4なのだが,特定の盤面から完成のゴール状態までの最短移動回数(最小手数)は最長80手まであると証明されたそうだ.コンテスト制限時間 45秒内で解くことは大変困難であろう.そのためか,本問題のテストデータには 45手までのものしか使われていないようだ.

それでも、A*アルゴリズムや、IDA*アルゴリズム等高速探索法を使わないとタイムアウトになってしまうだろう.

<A*アルゴリズムと発見的関数>
 A*アルゴリズムは人工知能の研究によく使われる探索技法のひとつ.初期盤面からゴールまでの最短経路の探索に,ヒューリステック(発見的)関数を用いている.

つまり,各々の盤面について,以下のような発見的関数値を計算する.それらの値の中から,最も関数値の小さいもの(盤面)より,次の探索を始める.勿論,最小手数でなくても良ければ,s(t) を使うことはない,つまり,f(t) = e(t) とすることもできるが,本問題では最小手数に近い解答を出さないといけないことから,s(t) を全く使わなくて良い訳でもなさそう.

   f(t) = s(t) + e(t)
     f(t): 盤面 t での発見的関数の値
     s(t): 初期盤面から盤面 t までの手数(移動回数)
     e(t): 盤面 t からゴールまでの手数の推定値

ここでは、推定値 e(t) に,盤面 t での各数字の位置と,ゴール状態での各数字の位置とのズレ(マンハッタン距離)の合計を使う。

具体的に、例えば、記事トップの盤面では、数字の1がゴール状態では左上のところにあるはずなので、その位置ズレは横方向 2、縦方向 2、計 4 と計算する。その他の数字については以下の計算になる。

   (記事トップの盤面)
 数字1: 横ズレ2、縦ズレ2、計 4
    2: 横ズレ0、縦ズレ0、計 0
    3: 横ズレ1、縦ズレ2、計 3
    4: 横3、縦2、計 5
    5: 横3、縦1、計 4
    6: 横2、縦2、計 4
    7: 横1、縦1、計 2
    8: 横1、縦1、計 2
    9: 横0、縦1、計 1
    10: 横1、縦1、計 2
    11: 横2、縦2、計 4
    12: 横0、縦1、計 1
    13: 横0、縦2、計 2
    14: 横0、縦2、計 2
    15: 横1、縦0、計 1
従って,e(t) の値は最終的に,上記数字1から15までの位置ズレの合計値,すなわち 37 になる.

すぐに気づくことだが,e(t) はそこからの探索手数の下限を示すものでもある.例えば,記事トップの盤面では 36手以下で完成することはありえない.

以下はプログラムに使っている f(t) の計算 に関する部分.

int evaluate(char *t, int step)
{
int r, c, rr, cc;
int x, e;
e = 0; for (r = 0; r < 4; r++) for (c = 0; c < 4; c++) { if ((x = t[RC2P(r,c)] - 1) >= 0) { rr = P2R(x), cc = P2C(x); e += ABS(r-rr) + ABS(c-cc); } } return step + 2*e; }

(プログラムの説明) 変数 r, c は行番号、列番号を表す.RC2P(r,c) により,1次元配列での対応位置)を算出。x はゴール状態での数字の正しい位置(0~15)。それをまたゴール状態での行番号 rr および列番号 cc に戻し,マンハッタン距離を求める.なお,関数の引数,t は盤面配列,step は初期盤面からの手数(ステップ数)を表す.

何回ものテスト結果から,プログラムに実際に使われている f(t) は以下のようなものに変えた.

   f(t) = s(t) + 2*e(t)

e(t) の前の係数を 1 にすると計算時間および記憶領域が大幅に増えてしまい,2 を越える値にすると正しく推定できず,探索エラーになってしまう.

係数を実数(2.5 や 2.1 とか)にしたりして,多少のテストもしたが,整数との違いはいまいち判らない.

<優先度付きキューの出番>
 初期盤面からたどりつけられるすべての盤面から,発見的関数 f(t) で計算した最小値の盤面より次の探索を始める訳なので,すべての盤面を優先度付きキューに登録すると,最小値となる盤面を取り出すことは効率的にできる。

プログラムに,下の記事に説明した通りの優先度付きキュー(小さい順優先)を使用している.

<2重登録の条件付き禁止とハッシュテーブル>
 同一盤面を優先度付きキューに重複して登録すると大変効率が悪くなる.しかも,ゲームの性質上,重複盤面は沢山出てくるはず.時計周りあるいは反時計周りの移動を考えてもそのことにすぐ気づく.知的ゲームといわれたのも,われわれ人間が無限ループを如何に避け,ゴールまでの移動パターンを見つけ出すかを問うからであろう.

ということで,ここではハッシュテーブルを活用して2重登録を防ぐ.ただし,同一盤面にたどりつけたということは,初期盤面から違うルート(違う移動方法)で同一盤面にやってきたこともありうるので,初期盤面からの手数は小さければ、ハッシュテーブルに登録し直さないといけない.

以下はハッシュ関数に関する部分.引数 id は登録しようとする盤面の識別ID,lookup( ) 関数のリターン値は 登録拒否であれば 1,登録OKであれば 0 である.

#define HSIZE   39989
typedef struct { int id; short cno; } TBL; TBL tbl[HSIZE+1];
int lookup(int id) { int i; char *t; TBL *tp; unsigned long long d;
t = move[id].t; for (d = 0, i = 0; i < 16; i++) d = ((d << 4) + t[i]; i = d % HSIZE;
while ((tp = tbl+i)->cno == cno) { if (!memcmp(move[tp->id].t, t, 16)) { if (move[tp->id].step > move[id].step) { tp->id = id; return 0; } return 1; } if (++i == HSIZE) i = 0; } tp->cno = cno, tp->id = id; return 0; }

(プログラムの説明) 定数 HSIZE はハッシュテーブルのサイズ、素数にするといいだろう.ハッシュテーブルの構造体メンバー id は盤面の識別ID,cno はテストデータのシリアルナンバーを表す.

盤面の実体は move テーブル(詳細については次の記事を参照)にあるが,ハッシュテーブルにはその盤面に対応する識別ID が格納されている.

ハッシュ関数では,識別 id に対応する盤面 t を move テーブルから取り出し,盤面 t の1次元配列について,盤面の各数字を4ビット左シフトし,累計を取る.さらにそれをインデックスとしてハッシュテーブルに重複盤面ID の有無を照合する.盤面が重複していて,しかも,初期盤面からのステップ数(手数)が大きければ,2重登録だと認識し,登録を拒否する.

2重登録になるけれど,初期盤面からのステップ数が小さければ,そのまま登録を許可する.

——–

原問題はこちら

言語の文字列を認識する問題.つまり,長さ 100 までの文字列が与えられ,以上の表に従い,最終的に目標に帰納できるかどうかを判定する問題.

例1: 文字列 aaabbbba および,目標 a が与えられると,
   aaabbbba ⇒ babbbba ⇒ cbbbba ⇒ cbbba ⇒ cbba
      ⇒ cba ⇒ ca ⇒ a
   になるので,帰納できる.

例2: 文字列 bccbbbbb および目標 c が与えられると,どう組み合わせても,文字 c に帰納(収束)できない.つまり,
   bccbbbbb ⇒ c はありえない.

本問題で取り扱う文字は a, b, c の3種類のみ.ただし,与えられる入力文字列が最長 100文字のものもありうる.

1文字を2ビットで表現し,100文字を64ビット整数4つで表現しようと試みたが,シフト演算等がめんどくさくて,途中で諦めてしまった.最終的に,文字列を単に char s[101] で表現する.

また,解答アルゴリズムとして,掲示板ではDP(ダイナミックプログラミング)法が盛んに議論されているが,ここでは単に深さ優先探索+カット技法を用いる.計算結果からみて,DPよりも速いようで,自分では不思議に感じる.

つまり,解答プログラムでは,与えられた文字列に対し,その先頭から末尾まで,2文字ずつ結合させ,トライ(試みる)していく.最終的に目標となる文字に帰納できたら,そこで終了.帰納できなければ,バックトラックして再トライ.

同じ文字列を繰り返しトライしたら,時間的コストがかかりすぎ,所定の10秒以内ではプログラムは終了しない.それを防ぐため,ごく簡単なハッシュ表を利用する.トライするすべての文字列をハッシュ表に登録し,同じトライはやらないことにしている.

以下はここで使うハッシュ表関連の部分.

#define SIZE  100            /* length of string */
#define HSIZE 1999           /* size of hash table */
typedef struct { /* hash table */ char n[SIZE+1]; short cno; } TBL; TBL tbl[HSIZE+1]; int cno;
void init(void) { memset(tbl, 0, sizeof(tbl)); }
/* * return 1 if n exists in hash table, else return 0. */ int lookup(char *n, int len) { int i; TBL *tp;
i = (n[0]*101 + len*499) % HSIZE; while ((tp = tbl+i)->cno == cno) { if (strcmp(tp->n, n) == 0) return 1; if (++i == HSIZE) i = 0; } strcpy(tp->n, n); tp->cno = cno; return 0; }

入力データごとにハッシュ表をクリアすると時間がかかるので,入力データシーケンスに対応する変数 cno もハッシュ表に登録する.

ハッシュ表以外に特筆することはないので,他のことは省略する.

C言語解答プログラム: Solution for UVa 10981 String Morphing (C program)
計算時間: 0.012秒

DPよりも速く計算できて,なんとも不思議.入力データが甘く作られただけかもしれないな.

Golygon およびその数を求める問題.

原問題はこちら
Golygon についての解説は こちら

Golygon とは,始点から辺の長さが1,2,3...と1つずつ増やし,90度左折または右折しながら,水平または垂直に進み,最終的に始点に戻る多角形のこと.

一般的な認識としては,辺の数 n = 8k (つまり,Golygonの辺の数が8, 16, 24, 32, … )の時にしか Golygon が存在しないが,本問題では,始点に戻る最長辺は始点から出発する長さ1の辺と必ずしも直角でなくても良いらしく,n = 7, 8, 15, 16, … の時に本問題でいう Golygon がつくれる.

それぞれの n での異なる Golygon リストは以下の通り.

n = 7 では 8個:
  enwnwse, eswswne, neseswn, nwswsen,
  senenws, swnwnes, wnenesw, wsesenw

n = 8 では 8個:
  enwswsen, eswnwnes, neswswne, nwsesenw,
  senwnwse, swnenesw, wneseswn, wsenenws

n = 15 では 72個:
  enenwswnwswsene, enenwswswnwnese, eneswnwnwswnese,
  eneswswswswnene, enwneswnwseswne, enwseswswsenwne,
  enwswnenwsenesw, enwswseseswnwne, enwswseswsenenw,
  esenwnwnwnwsese, esenwswswnwsene, eseswnwnwswsene,
  eseswnwswnwnese, eswnenwnwneswse, eswnwnenenwswse,
  eswnwnenwnesesw, eswnwseswnesenw, eswsenwswnenwse,
  neneswseswswnen, neneswswsesenwn, nenwseseswsenwn,
  nenwswswswsenen, nesenwseswnwsen, neswnwswswnesen,
  neswseneswnenws, neswswnwnwsesen, neswswnwswnenes,
  nwneseseseswnwn, nwneswswseswnen, nwnwseseswswnen,
  nwnwseswsesenwn, nwsenesesenwswn, nwseseneneswswn,
  nwsesenesenwnws, nwseswnwsenwnes, nwswneswseneswn,
  seneswnenwswnes, senwnesenwseswn, senwnwswnwsesen,
  senwnwswswnenes, senwswnwnwsenes, sesenwnenwnwses,
  sesenwnwnenesws, seswnenenwnesws, seswnwnwnwneses,
  swneneseneswswn, swnenesesenwnws, swnenwswneswsen,
  swneseneneswnws, swnwsenwnesenws, swsenenenenwsws,
  swsenwnwnenwses, swswnenenwnwses, swswnenwnenesws,
  wnenwseneswsenw, wnesenwneswnwse, wneseswseswnwne,
  wneseswswsenenw, wneswseseswnenw, wnwneseneseswnw,
  wnwnesesenenwsw, wnwsenenesenwsw, wnwsesesesenwnw,
  wsenenwnenwswse, wsenenwnwnesesw, wseneswsenwswne,
  wsenwnenenwsesw, wseswnesenwnesw, wswneneneneswsw,
  wswneseseneswnw, wswseneneseswnw, wswsenesenenwsw

n = 16 では 112個 (数が多いので,詳細は次のC言語ソースプログラムを参照).

上の表現法では,東,北,南,西への進行方向をそれぞれ e, n, s, w で表現している.

例えば,n = 7 での最初の Golygon は enwnwse と書いてあるので,始点(原点)から,東へ1歩,北へ2歩,西へ3歩,北へ4歩,西へ5歩,南へ6歩,東へ7歩,と歩けば,最初の始点(原点)に戻れることを意味する.

なお,本問題の出力フォーマットに従い、上記のリストは e, n, s, w の辞書式順序に沿ってソートしてある.

通り道に障害物のない Golygon を表示し,その数を累計すれば,本問題はこれで解ける.

C言語解答プログラム: Solution for UVa 225 Golygons (C program)
計算時間: 0.047秒

数論では数少ない、きちんと整理された解法、平方剰余 (Legendre’s/Jacobi’s symbols) の応用問題があった。

ACM問題 #10831 Gerg’s Cake
意訳: パーティを披こうとした。p人の友達が到着したけれど、遅刻者はまだ a人もいる。長く待たせるのも悪いので、正方形のケーキを先着者にご馳走させようと考えた。遅刻者のa人分を残しながら、ケーキを p人で平等に分けたい。ただし、ケーキの一切れはまた正方形に切らないといけない。果たしてそんな切り方は可能だろうか。

<注意> 問題の記述文の、最初のパラグラフ(段落)のところに、p は合成数でないことが書いてある。言い換えると、p が1または素数だ。

例:先着者5人、遅刻者1人だとすると、ケーキを16等分の正方形にし、先着者はひとり3切れ食べられることになる。先着者3人、遅刻者1人なら、ケーキを4等分にすればOK。

上の話しはもちろんつぎの数式に合わせて書き上げられたものだろう。

  X2 = a + n*p

a と p はそれぞれ、遅刻者と先着者の数。n、X は非負の整数。

さて、式の両辺を p で割ると、

  X2 = a (mod p)

の形になり、平方剰余の計算に帰着される。

平方剰余のことは 既に他のところに書いた。つまり、奇素数 p と、 p で割りきれない a に対して、

    X2 ≡ a (mod p)

なる X が存在するときに、a を mod p での平方剰余といい、そんな解 X がないときは、平方非剰余という。解 X があるかどうかの判定 は効率的にできる。

ただ、上記の平方剰余の計算関数(解の存在を判定する関数)では、前提として、

  p は奇素数でなければならない。また、
  a は p で割れないものでなければならない。 

の2点があるので、ここでは、前提に合わない a と p をリストアップして、個別に考える。

  (1) p が偶素数(つまり2)、
  (2) p が 1、
  (3) a は p で割り切れるもの。

(1) については、もとの数式が X2 = a + 2n になるので、任意の奇数値の a に対しは、奇数の X を選び、任意の偶数値の a に対しては、偶数の X を選べば、整数 n が必ず存在することが判る。

(2) についても簡単に n が存在することが判る。なぜなら、 X2 = a + n という式になるので、n = X2 – a となる n が必ずあるからだ。

(3) については、a が p の倍数になるので、 a = mp とすると、X2 = (m+n) p になる。X を p の倍数、例えば、X = kp として選べば、 m+n = k2 p になるので、整数 n が必ず存在するわけだ。

前提に合わない a と p はすべて解をもつことが以上の分析で判っただろう。

最後に、参考のため、C言語プログラムをのせておく。
p10831.c

1566問もあるACM問題集に昨年11月末頃からチャレンジし、やっと50位ランキング入りを果たした。目標は今年にトップ10入りことだが、いまの進行ペースで考えると、その目標は通過点でしかないだろう。できればナンバー1になりたいが、無理だと思う。ほとんどの問題がその人に解けられたから。

ACMとはAssociation for Computing Machinery、情報工学に関するさまざまな情報交換を推進し、標準の普及等を通じて、情報工学の技術的、応用面での進歩発展をはかることを目的とした国際的な科学・教育組織のこと。現在会員は約80,000人。情報技術に関心を持つ実務者、開発者、研究者、教育者、管理者等から構成されている。日本では情報処理学会が有名だが、世界的にはACMだろう。

その手の問題集アーカイブサイトは世界中にあちこちあることだが、最強サイトは勿論そこ。英語が問題だと思うのだが、日本には残念ながらこの手のサイトは見当たらない。

ランキング表をみてビックリしたのがバングラデッシュという国の活躍。日本では全く知られない国だが、実力では世界有数のソフト大国なのかもしれない。一方、インドの方々の活躍は全くみられない。実情は果たしてどうだろう。

問題を解くことでどんなメリットがあるか、というと、まずは英語の勉強にはなる。問題の表現には最小限の英語しか使わず、図表入りの優しい説明にもなっていない。英語の微妙な言い回し等、自分の好きなプログラム作成を通して、勉強できるならばそれ以上のことはないだろう。

プログラミングやアルゴリズムのエッセンスも習得できる。大学に行かなくても、200問を自力で解けば、情報工学科とかは十分卒業できる実力が身につくはず。(日本の某情報系の大学教授も200問ぐらいしか解けなかった。それが世界の実力かもしれない。)そこにある問題はほとんどが世界各地で開催されたプログラミングコンテストで使われていたものなので、レベルは世界一。最先端のアルゴリズム研究も問題を解くことで問題意識やヒントとなるものは得られるだろう。

問題を解く環境は、Linuxを一台用意することだ。中古のPCでも構わない。LinuxにはGCC(
グニュC コンパイラ)が入ってるので、サーバと同じ実行環境になり、心強い。あとはネットに繋ぎ、サーバに登録。それだけ。

残りは気合と粘り精神。一日一問でも構わない。つねに考えることが大事。どうしても判らなければ、掲示板 http://online-judge.uva.es/board/index.php を利用するのもいいだろう。沢山解くよりも、一問一問完全に理解することを目標にしよう。解けたら、もっと速くするにはどうすればいいか、違う解き方はないか、を考えてみるのも悪くない。

本人の場合は、気合というよりもただ好き、それだけのことだ。飯よりもプログラムをつくるのは好きということ。アルゴリズムに魔法が存在する。今までに解けなかった問題、あっという間に解いてしまう技法、マジックというのはそういうものだろう。それを探求するのが楽しみなのだ。未解決の問題はまだ800問もあるので、まだまだ楽しめる。

ACM問題 #10465 Homer Simpson は一次不定式を解くもの。問題の記述に曖昧性はあるものの、ためになる良問のひとつ。

<問題の説明>
 2種類のハンバーガーをそれぞれM, N 分で食べ切ることができる。与えられた T 分ちょうどの時間に、最多何個のハンバーガーを食べ切ることができるか。また、どうしても時間が余ってしまうなら、余った時間を飲み物にすると良いが、飲む時間は最小にして欲しい。

 入力パラメータ  整数 M, N, T ( 0 < M, N, T < 10000 )
 出力値       K (食べたハンバーガーの個数、K >= 0)、t (飲み物の飲む時間。飲む時間がなければ出力不要)。

<例を考える>
 理解を助けるために、例を考えてみる。
 例えば、M=3、N=5、T=54の場合、18*M = 18*3 = 54 なので、K=18が解答になる。Nの値がMより大きいので、18が最大値であることが理解できよう。
 また、M=3、N=4、T=5の場合、ちょうど5分間でハンバーガーを食べ切ることはできないが、飲む時間をいれて調整すると解答が存在する。つまり、飲む時間を最小限の1分にすれば、ちょうどハンバーガー1個食べ切ることができるから、K=1、t=1 がその解答。

<問題の整理>
 ハンバーガー云々のストーリーが先にあって問題作られたのか、数学的モデルに無理やりストーリーを書きあげたのかは、出題者に聞かないと判らないが、 一次不定式の解き方が問われている。つまり、

  Mx + Ny = T
  x >= 0、y >= 0

を満たす整数解 x、y を求め、解があれば、x+y を最大にするということだ。
 もし、解 x、y がなければ、解があるまで、Tを減らすこと。減らした時間を飲む時間に当てる。勿論、減らす量を最小限に止めるべきだ。

<拡張剰余定理>
 整数 M, N の最大公約数を D=gcd(M, N) とすると、 1次不定式 Mx + Ny = D を満たす整数解 x, y は必ず存在する。詳細はここで省略するが、以下の関数は一組の X、Y 及び D を算出する。(関数では、入力引数の変数名に a, b を使うが、ここでいう M, N と考えてOK。)

int extended_gcd(int a, int b, int *x, int *y)
{
int d;
if (b == 0) { *x = 1; *y = 0; return a; } d = extended_gcd(b, a%b, y, x); *y -= a / b * (*x); return d; }

 また、証明は省略するが、問題の入力パラメータ T がDの倍数、つまり、T がMとNの最大公約数 D の倍数であれば、取りあえず整数解 x, y が見つかる(ただし、非負整数である保証はないので、我々の求める x, y でない可能性もある。)
 T がDの倍数でなければ、整数解が存在しない。飲む時間 t をいれて調整することが必要
 ⇒ t = T % D

 上記の拡張剰余定理で見つかった一組の整数解 x, y を元に、一般解を組み立てると以下になる。

<一次不定式の一般解>
 D = gcd(M, N) とし、 X, Y は不定式 MX + NY = D を満たすものとすると、

 Mx + Ny = T (ただし、T は Dの倍数でないときには、Dの倍数になるように値を減らす)の一般整数解 x, y は、

 x = (T/D)*X – (N/D)*S
 y = (T/D)*Y + (M/D)*S

で与えられる。 ただし、S は任意の整数。なお、(T/D) はC言語でいうところの割算。

<非負整数解にするには>
 x, y の値は、本問題では、ゼロ以上の非負整数でないといけないので、S の変域をここで求める。

 (N/D)*S <= (T/D)*X
 (M/D)*S >= -(T/D)*Y

 -(T/D)*Y / (M/D) <= S <= (T/D)*X / (N/D)

細かい記述は敢えてここで避けておくが、Sの変域の両端を整数値 Smin, Smax とすると、

 Smin = -(T/D)*Y / (M/D)
 Smax = (T/D)*X / (N/D)

となる。はずだが、慎重に計算して欲しい。
 1.4 <= S <= 4.5 では、Smin = 2 (切り上げ), Smax = 4 (正数の切り下げ=C言語の割算)
 -4.5 <= S <= -1.4 では、Smin = -4(負数の切り上げ=C言語の割算), Smax = -2 (切り下げ)

 また、x+y の値は x + y = (T/D)*X + (T/D)*Y + ((M/D – (N/D))*S で与えられるので、
 M >= N であれば、 S = Smax を取れば、 x + y は最大になり、
 M < N であれば、  S = Smin とすれば、 x + y は最大になる。

<話はまだ終わらない>
 一見すると解答は出たようだが、実は Smin と Smax は計算できない場合がある。つまり、Smin > Smax になってしまったりすることがある。その場合では、TをDの値分ずつ減らし、繰り返し計算すれば、必ず条件を満たす Smin, Smax が見つかる。

<解答プログラム> p10465.c

最後に、テストデータを掲載しておく。

<テスト用入力データ>
4 4 10
1 7 10
2 3 6
11 2 13
2 5 25
3 4 5
5 3 5
4 5 7
9 7 62
2 5 6
45 23 91
5 9 19
4 6 42
3 7 10
12 56 71
50 51 103
20 30 7
100 8 1
123 12 21

<対応する出力結果>
2 2
10
3
2
11
1 1
1
1 2
8
3
3
3
10
2
2 3
2 1
0 7
0 1
1 9

罠の多い問題。

問題 10093 An Easy Problem!

数字 R が与えられ、それを (N-1) で割り切れる最小の基数 N (つまりN進数の基底) を求めなさい。というのが問題の内容。

例えば、R = 2727 に対して、Rは9で割り切れるので、N = 10 がその答え。

10進数の場合、9で割り切れるかどうかの判定は極めて容易。数字の各桁の和を取り、それが9で割り切れるかどうかを調べるだけで必要十分条件になるわけだ。

例えば、2+7+2+7=18 は 9で割り切れるので、Rも 9で割り切れる。2725では各桁の数字の和は9の倍数ではないので、2725 も絶対に9の倍数ではないのだ。証明は簡単なので省略する。

さて、この事実はN進数にも成り立つ。N進数において、各桁の数字の和が N-1 で割り切れるならば、その数字も N-1で割り切れる。それを利用して、

各桁の数字が
  ’0′ ~ ‘9’ なら 0~9、
  ’A’ ~ ‘Z’ なら 10~35、
  ’a’ ~ ‘z’ なら 36~61
に置き換えて、素直に足していく。その合計を sum とする。

また、N進数において、各桁の表現できる最大数は N-1 (つまり2進数では1、10進数では 9 、16進数では15(便宜上Fと表現するが)のこと)なので、基数は最大数字の一つ上以上でないといけない。2727の例では基数は8からということだ。なお、1の基数 (1進数、死の世界)は認められていないので、最小基数は2となる。

その(最小基数-1) で 合計のsumを割って行き、割り切れたらそこで終了。1つずつ増やして62になっても割り切れないならば ”such number is impossible!” と出力。

例えば、 2727に対し、sum = 18。8では割り切れないが、つぎの9でOK。したがって、N = 10となる。

さて、これで問題が解けたのかと思ったら大間違い。苦戦苦闘の末判ったのが、

 1.入力の数字に符号+か ー が含まれている。
 2.入力数字の桁数は異常に長い。1000桁以上の数字が存在 している。2000桁の配列(入力バッファとして使っている)を取ったら、やっとACもらった。2000桁であっても、sumは整数で十分。各桁の大きさは61までなので、最大でも12万強。

最後にテストデータと結果。

Sample Input
0
1010101101
  +265
    -5464
z
Zz

Sample Output
2   ←最小基数は2だから
2   ←0と1からなる数字列は1で割り切れるから、すべては2進数
14
20  ←マイナス符号は関係ないので、読み飛ばしてOK。
62
such number is impossible!