2023-12-01

ICPC2023 Asia Yokohama Regional 参加記（pointN）

競プロ旅行

pointNです。

ICPC2023アジア横浜地区大会に筑波大学からBig O of N cubedとして参加したのでそのことを書きます。

チームについて

チームメイトはniuezさんとruthenさんです（以下敬称略）。参加時のAtCoderのレートは3人とも2000ぐらいでした。全員前の年までは別のチームで活動しており、チームを組むのは今年が初めてでした。チーム名は要するに $\mathcal{O}({N}^{3})$ であり、由来は3人のハンネの最長共通部分列がnだからです。pointNとruthenは年齢の制約上今年がラストイヤーでした。ちなみに筑波にはもっと強いGoodBye2023というチームがおり、つまりBig O of N cubed（われわれ）は二軍です。

データ構造やアルゴリズムに対する知識量が他より多く重実装にも強いniuezが実装の中心になり、それを感覚寄りのpointNと理論寄りのruthenがいい感じに支えたりそうでなかったりというスタイルのチームでした。

金曜日（前々日）

大会後に持ち越すとまずい宿題があったのでそれをやり、残りはそわそわしていました。チームメイトはライブラリを準備してくれていたりしました。この辺、任せっきりになってたのは申し訳ないと思っています。恥ずかしながら蟻本すら持っておらず、螺旋本ぐらいしかないのでまあそれは持って行くねーとだけ言いました。

土曜日（前日）

11時に秋葉原で集合ということになっていたので、頑張って起きてTXに乗りました。途中で競プロイベントで使えそうな名札を忘れたことに気づきましたが、もう遅いし、どうせ知り合いも少ないので気にしないことにしました。niuezは先に秋葉原に行ってゲーセン行ってたっぽいです。改札前に待機していると上半身オレンジでめちゃくちゃ目立つ顧問（YesNoおじさんとして知られるアランニャ・クラウス先生）が来たので合流、まもなくniuez、ruthenの順に到着しました。GoodByeの人たちは別行動するらしいので、4人で京浜東北線に乗って横浜に向かいました。

僕は横浜に行くのは初めてだったのですが、練習で解いた問題の話をしていると思ったよりすぐ目的地の関内に着きました。駅を出たところに野球場があり、しかもこの日はファンミーティングがあってユニフォームを着ている人が多かった（このせいでホテルが取りにくかったそうな）ので、道中は野球の話をしていました。とはいえ、チーム全員あまり野球に明るくなく、どの球団がセ・リーグでパ・リーグなのか、DeNAとベイスターズって同一の球団だっけ、などと、かなり程度の低い会話をしていました。広島出身のpointNはもうちょっとプロ野球に精通しているべきだと思うのですが。

そこから中華街でお昼を食べました。金香楼というお店で、壺料理を推しているお店でした。店内には水が流れていました。僕がAランチ（スペアリブ壺スープ）、ruthenとniuezがCランチ（牛バラ壺煮）、Y/NおじさんがDランチ（麻婆豆腐）を選択しました。らっきょうを食べたことない人（niuez）とか杏仁豆腐を食べたことない人（niuez）とかがいて、世界は広いなあと思いました。

集合時間までには余裕があったので中華街を歩いて産貿ホールに向かいました。ずっと産貿ホールのことを産賀ホールだと思っていたのは内緒です。お昼時で人が多く、それなりに荷物が多い自分は歩くのが大変でした。鳥の揚げたでかいやつがおいしそうでした。

産貿ホールに着いてすぐ受付が始まりました。英語なのでちょっとだけ異世界感がありました。自分たちのスペースは入口にかなり近い位置でした。机上にはリハーサルで使う触れてはいけない封筒、名札、Tシャツなどが置かれており、名札はつけて、Tシャツは着てねーというアナウンスがモニターに映っていました。Tシャツを着るのは日曜の本番だけだと思っていたので、パーカーの上から半袖Tシャツという、かなり奇抜なファッションになってしまいました。

リハーサルは4問で行われました。システムを使うのは初めてだったのですが、なんとなく理解できました。Dで伝達ミスがありniuezに違う問題を解かせてしまいました。ごめんなさい。すぐ修正してくれるので心強かったです。ちなみに、リハーサルに置かれていた問題は横浜が初めての自分は知らないものでしたが、niuezが正しいほうのDを実装しながらなんかこれ書いたことある、などと言っていました。Clar（質問）の練習～とか言って、中華街でおすすめの料理店を聞いたりしました。関係ないことを聞くのはやめたほうがいい気がします（ちなみにYaku-Shoronpoという答えが返ってきたのですが、終わってから検索してもいまいちわかりませんでした）。

その後はチーム紹介でした。順に壇上に上がって30秒で簡単に紹介していくものでしたが、マイクを通した英語は全然わからないし、チーム数多いし、知り合い全然いないしでほとんど記憶に残りませんでした。インパクトのあるスライドと簡単なフレーズを意識するのが大事そうです。右後ろにSpeed Star（強者揃いの東大チーム、翌日ぶっちぎりで優勝する）がいることに気づき、はえーと思っていました。

終わってから一度ホテル（東横INN横浜スタジアム前2）に戻って荷物を置き、夜の中華街に3人で繰り出しました。最初に中華っぽい食品を扱うスーパーを覗きました。僕以外の二人は何か買っていました。

おなかがすいたので、お昼を食べた金香楼の隣にある東光飯店というお店に入りました。pointNが名物っぽい東光チャーハンを、チャーハンに目がないniuezも何らかのチャーハンを、ruthenはエビの入ったラーメンを頼み、別途、回鍋肉と小籠包を注文しました。スープが飛び散って美味しくなるなどと競プロerにしか通じないネタで盛り上がりました。3人で割って1800円ぐらいだったと思います。

その後は中華街をゆるっと歩きました。さっきの中華スーパーに戻り、pointNは5個入りの胡麻団子を、niuezは肉まんを買いました。niuezが肉まんの大きさに驚いていました。胡麻団子は僕が独り占めしてすべて平らげたのですが、あれは2個目がピークです。好きなので5個なら平気ですが。その後、pointNとruthenはパンダまんを買いました。率直な感想として、僕は「中のあんこは美味しいけど外の皮が観光客の味がする」と述べました。

ホテルに戻り、一人で食べるには多いと、ruthenが中華スーパーで購入したバナナチップスを分けてくれました。お腹に余裕がまだあったので受け取り、明日7時にホテルで朝食という約束を再確認して解散しました。バナナチップスは自室でポリポリ食べました。

甘いものが連続して少し辛い物が食べたくなり、一人で近くのコンビニに行きました。レモンの炭酸と、カラムーチョを選択しました。お風呂に入ってから食べたのですが、カラムーチョは意外と量が多く（しかもちょっと増量キャンペーン中）、後半食べるのに苦戦しました。ABCもあったのですが、普通に無理なので寝ました。

寝たはいいものの、3時、5時、6時などに目が覚めてしまいました。途中、昔の同級生が出てくる夢を見たのですが、残念ながら全然楽しくないストーリーでした。トコジラミがツイッターのトレンドに入っていたりして若干不安だったことなども眠りが浅かった要因かなと思っています。

日曜日（当日）

無事全員起床に成功し、ホテルのバイキングに並びました。麻婆豆腐があって珍しいと言ったらruthenから「中華街だからでは」という真っ当な指摘が入り、なるほどと思いました。自分たちが席に着いた頃から急に客の列が長くなり、早めに行って正解でした。

その後はスマホの充電器を忘れそうになりながらも荷物をまとめました。起きたらちょっとだけ残しておいた炭酸が凍っており、その画像をツイッターにアップし、飲まずに部屋に忘れました。3人そろってチェックアウトして集合場所（隣のホテルのロビー）に行きました。GoodByeの人から、近くの台湾チームに貰ったというパイナップルケーキをもらいました。いくら競プロerとはいえ、"闇討ち"というフレーズがすぐに浮かんでしまったのは反省です。コンテスト後にでもありがたくいただくことにして出発しました。ギリギリ傘がなくても耐えそうな雨が降っており、寒かったです。

会場で受付を済ませ、自分たちのスペースに入って準備をしました。荷物は半透明の袋に詰める必要があったのですがそこまで大きくなく、詰め放題の主婦の要領で袋を引き延ばすことでどうにか入れ込んだのですが、気づくと袋が破けてしまっていました。スタッフの人に袋を指さしながら Excuse me, my plastic bag is broken ... などと怪しい英語で事情を説明すると一回り大きな新しい袋を持ってきてくれました。Thank you.

ほどなくして、コンテストが始まりました。初動は、ログイン情報の入った小封筒の開封とPCのセットアップをruthen、Aをniuez、問題文の入った大封筒の開封とBをpointNという分担で動きました。Bの問題文が若干難しく、僕がモタモタしている裏でAはあっという間に通っていた印象です。後で問題を見返したら普通に重実装だったので、niuezがうまくはまったという印象です。

いっぽうBを読んでも全然解法がわからない僕はniuezと協力してBを、ruthenが以降を眺める方針をとりました。niuezがLazy Segment Tree（有名ながらおよそBで使うとは思えないデータ構造）を使った解法に気づいたので、完全に信頼して僕も問題読みに加わりました。

Fが驚くほど簡単で、ruthenに軽く説明をし、Bを通したniuezにruthenから説明を行い（僕の説明があまりにも雑で窘められた、リハーサルのことがあったので完全にruthenが正しい）、後は実装マシーンniuezに任せました。Dを読んで区間DPっぽいことを確認したり、Kのインタラクティブなどを読みました。インタラクティブを優先的に読んだのは、形式の違う問題なら多くのチームに解ける可能性のある難易度調整にしているだろうというメタ読みをしていたためです。他にもいろいろ読んだのですが、僕の力では解法に至れる問題はありませんでした。

ruthenとniuezがFとDをやっている間にKを考えました。円の位置と大きさを特定する問題なのですが、とりあえずの大まかな位置の特定にかなりのクエリ数を使うため、残りの部分がなかなか詰められませんでした。その間にFDが通り、他は解法すら浮かばないのでKをとりあえず書いていくことになりました。乱択でクエリ数の期待値を減らすなど迷走し、二分探索で行ける、いやこれ三分探索にならないか、いややっぱり二分探索か、にしてもこれ場合分け面倒すぎませんか、などの議論を経て、算数を頑張ることで比較的シンプルに解ける方針を思いついたので、それを実装しました。しかしバグらせてしまい、いったん印刷してPCを離れたりしました。

その後はniuez-ruthenチームとpointNチームで唸る時間があったような気がします。なんやかんやあってKのバグに気づいたので修正し、サンプルが通ることを確認しました。他のサンプルがなくて不安などとほざいていたらniuezがインタラクティブツールの何かしらをいじってうまいことやってくれました。大感謝。それも通ったので問題ないと判断し、投げるとACが返ってきました。

その後は全員でGとTLEしているEを考えました。途中でトイレに行きたくなったのでスタッフに Can I go to the restroom? などと言い、連れて行ってもらいました。戻ってきて僕もruthenにEの解法をなんとなく聞いて高速化の検討をしたのですが、あまり活躍できませんでした。niuezがEのゼータメビウス変換？とかいうやつを思いつく裏で定数倍高速化を頑張ったruthenがそのEを通し6完、終了間際にniuezが投げた愚直のGはWAが返ってきました。

システムトラブルがあったチームがあるらしく（GoodBye？）、延長したりしなかったり、変な感じでコンテストは終了しました。途中、なんかやばそうなアナウンスがあった気がしますが、大丈夫だったようでよかったです。

コンテストが終わり、いつの間にか配られていたまい泉のすずらん弁当を食べました。配られたドーナツも食べました。僕はポン・デ・リングの亜種（ショコラ？）をもらいました。スポンサーの紹介を聞き、解説を聞き、我らがYesNoおじさんによるYes/Noがありました。Big O of N cubedは27位でした。3の3乗なのは何の因果でしょうか。個人的には、全員がACに貢献できたのは平和でよかったと思っています。

懇親会では大量の食糧が置いてあったのですが、昨晩のカラムーチョのせいかあまり食べられず、また知り合いが少ないのでやや心細い思いをしました。niuez周辺に人だかりができており、ゲームとかそういうので繋がりがあるのかな～と思いました。企業ブースも何ヶ所か回り、いくらかの戦利品を入手しました。Fixstarsのクイズに正解してポロシャツをもらったのが特に大きかったです。その他、先輩であるTsumoiYorozuさんにご挨拶をしたり、ITF.PCでお世話になったtatyamさんにお礼を言ったりしました。

友達をチャーハンに連れていくというniuezとはここでお別れし、ruthenと一緒に（幸せそうな人であふれる）赤レンガ倉庫やハンマーヘッドを散策しながら競プロ引退老人トークをしました。みなとみらい駅まで歩き、横浜でお土産を買って秋葉原へ、そしてTXに乗って家に帰りました。休日の格安回数券を使い忘れたことを今も少し後悔しています。

ICPC5回目、最後の年に横浜に行けて本当によかったです。ruthenさん、niuezさん、ありがとうございました。

おわりに

競プロと出会ってからそれなりに長い時間が経ちました。未だに青と黄色を反復してしまっていますが、競プロerとしてまあまあ成長できたのではないかと思います。これまでにチームを組んだ人たち、サークルの人たち、その他競プロを通じて出会った人たち、chokudaiさん、Mike Mirzayanov、コンテストを運営してくれている皆さん、ありがとうございました。やり残したことだらけなので、今後も競プロとはゆるく付き合っていければなと思っています。

2022-12-11

車輪を最発明しようの会リード・ソロモン符号

プログラミング

この記事は呉高専アドベントカレンダー2022 11日目の記事です。

ごあいさつ
- 注意点
はじめに
ふんわり理解するリード・ソロモン符号
- 符号化・復号の流れ
  - 符号化
  - 復号
忙しくない人と僕のためのRS符号の理論
実装
おわりに
おまけ
参考文献

ごあいさつ

こんにちは、点Nです。電気情報工学科卒です。高専在学中は高専プロコンに出たり、呉高専史上初めてICPC（大学対抗競プロ合戦）に参加して予選敗退したりといった感じです。

アドベントカレンダーというと「（最新の何かすごい技術）を使って（こんな面白いもの）作ってみた！」みたいな（勝手な）イメージを持っているのでそういうキラキラした記事を書きたいところですが、あいにくそういうのは性に合わないので黒い画面に白い文字が出るタイプのプログラミングをやります。タイトルにもある通り車輪の最発明をやりますので、みなさんは僕が車輪を最発明しているようすを見ていてください。目次の長さや画面右のスクロールバーの大きさから察している方も多いと思いますが、この記事はめちゃくちゃ長いです。せっかくの日曜日ですが（なので）、どうか最後までお付き合いください。

注意点

素人プログラミング
ふわふわ理論
ガタガタ構成
エンドレス文章

それでもOKな人はどうぞ

はじめに

（呉高専アドベントカレンダーなので、高専を絡めた導入をどうにか考えました。）

高専といえばロボコン。世の中そういうイメージの人が大半ですが、高専にはそれ以外のコンテストもあります。代表的なところでは、ごあいさつで述べた高専プロコンというプログラミングの大会が有名です。その中の競技部門という部門では与えられた課題（対戦ゲームや一人用ゲームなど）を解くプログラムの強さを競います。しかしこの部門、稀に「プログラミングするより人間がやったほうが強い」という現象が起こることがあります。例えば、第23回有明大会の競技課題は「サイコロの山を見てその数を数える」というものでしたが、上位の成績を収めた多くのチームが（画像処理ではなく）人の目でサイコロを数えるという手段をとったそうです。

第27回鳥羽大会でも似たようなことが起きました。競技内容は木製のジグソーパズルが（現物で！）与えられるので解け、というものでした。パソコンにピースの形状を読み込ませ、正しい配置を探索するというのが正攻法なのですが、それでも解くのが難しく、また外枠内に収めたピースの数で評価したために「最大のピースを諦めて残りを人が急いで詰め込む」という計算機完全無視の手段がかなり強くなってしまいました。実際、翌年の募集要項に「昨年の大会では，手でパズルを並べると上位に来ることができました」という記述があります（https://www.procon.gr.jp/wp-content/uploads/2017/04/8db30280f6638520be40ad86b4bf1b37.pdf）。

勘のいい人はお分かりかと思いますが、募集要項にこのようなことが書かれた第28回大島大会の競技課題は、第27回の競技課題をアレンジしたものでした。競技者はペナルティを支払うことでパズルを解くヒントを利用することができるようになり、ヒントを使ってでも「パズルを完成させる」ことを重視する評価基準になりました。そして、ヒントは紙に印刷されたQRコードで配布されました。仕組みを理解していれば不可能ではないのですが、よほどの天才でない限りQRコードは人力では読めないので、さすがにプログラミングをせざるを得ません。

高専プロコンの話はこれぐらいにして、本題に入ります。このQRコード、多少汚れるなどしていても読み取れるのをご存知でしょうか。これは本来のデータだけでなく、いっしょに「データが壊れた時に気付いて、元に戻すための手がかり」も記録しておくことにより、汚れなどで欠けた部分を補って正しいデータを復元できるようにしているためです。

これを実現するための具体的なやり方はいろいろありますが、QRコードではリード・ソロモン符号というものを用いています。今回はこの符号を最発明します。といっても理論を思いつくのは不可能なので、意図的に壊したデータを頑張って復旧させるプログラムを書きたいと思います。

余談：人力という戦略について

人の力に頼ること自体は悪いことではありません。人もシステムの一部とみなして有効に使うことで非常に優れた問題解決能力を実現できる場合があります。例えば、第30回都城大会の競技課題は陣取りゲーム形式の通信対戦ゲームでした。1ターンの時間は数秒から十数秒と短く、人力で最大8人もの味方の動きを制御するのは難しいため、基本的にはコンピュータが考えることになります。しかし、味方全員は不可能であっても、その中の一人か二人程度であれば、その動きを人が決定することが可能です。「シンプルな戦略に則ってコンピュータが個々の味方の行動を計算、人間は必要があればそれを修正し、大局的に見てより良い手を柔軟に選択できるようにする」という方針をとった東京高専はこの年の競技部門で優勝しました（というかこの年は東京高専が3部門全てで優勝したのですが）。

当時の東京高専の学生によるツイート：

shibh308 on Twitter: "#procon30 深夜になってしまいましたが、東京高専競技部門の方針について説明します" / Twitter

計算能力ではパソコンに遥かに劣りますが、人の力も決して弱くはありません。問題を人力で解くことが十分強いのであればそれも立派な戦略だと思います。プログラミングの大会としてコンピュータを一切使わないのが最適ってのはいかがなものかと思いますが。

ふんわり理解するリード・ソロモン符号

リード・ソロモン符号（Reed-Solomon Coding, RS符号）はリードさんとソロモンさんが考えた誤り訂正符号（データにくっつけておくことでデータに発生した間違いを発見し、直せる符号）です。QRコードをはじめ、Blu-rayディスク、地デジ放送、ADSLその他諸々に応用されています。この世界はリード・ソロモン符号なしには成り立たないと言っても過言ではありません。ちなみに背景にある数学的な理論はちゃんと難しいです。

符号化・復号の流れ

いったん、細かい話は抜きでざっくり書きます。かなり説明不足で厳密性に欠けます。

符号化

線形代数の気持ちになります。小文字がベクトルで大文字が行列です。まず元の情報 $m$ があります。これに生成行列 $G$ をかけて符号化します。それを $w$ とします。 $w=mG$ です。 $w$ は $m$ よりちょっと長くなります。 $m$ は壊れるとおしまいですが、これを多少壊れてもいい符号語 $w$ に変換しました。終わりです。

復号

壊れている可能性のある、 $w$ と思しきデータ $w'$ を $m$ に復元しましょう。まず、 $w'$ が本当に壊れているか確認します。これには検査行列 $H$ を使います。 $s=w'H ^ {T}$ とします。 $s$ はシンドロームと呼ばれます。 $s \neq 0$ ならデータが壊れています。

$w' \neq w$ の場合を考えます。これは混入したエラー $e$ があって、 $w'=w+e$ な状態です。 $w$ を復元するために $e$ を特定します。方程式 $eH ^ {T}=s$ を解けば $e$ が特定できます（ざっくり過ぎてこの辺まではだいたいの実用的な符号化に当てはまる話です）。ですが、今は未知の変数が多すぎて解けません。そこで、先にエラーのある場所を調べ、エラーでないと分かった変数を消します。

そのために、 $s$ からシンドローム行列 $S$ なるものを作ります。なんと $rank(S)$ がエラーの個数です。さらに、 $S$ から誤り訂正多項式 $\sigma (x)$ を得ます。その根（ $=0$ となる解）は $w'$ のどこが $w$ と異なるかを示します。無事エラーの位置が分かりました。変数が減ったので方程式が解け、 $e$ がわかります。

$w'$ から $e$ を除去し、 $w$ がわかります。 $mG=w$ なる $m$ を求めれば、復号完了です。

忙しくない人と僕のためのRS符号の理論

実装するだけなら理解しなくても大丈夫ですが、RS符号がなんなのか全く分かっていないのもまずいので、ここでRS符号とはつまり何なのかについて知っておきたいと思います。適宜飛ばしてください。僕の理解が足りないので、間違った記述も含まれていると思います。

そもそもRS符号は何がすごいのか

とにかく訂正能力が高いです（そのぶん処理はやや複雑です）。正しいかは素人なので判断できませんが、僕は「理論値じゃん」という印象を受けました。

まず、符号どうしの距離を考えます。符号の最小距離 $d$ は、異なる符号どうしのハミング距離の最小値です（ハミング距離は違っている場所の数です。たとえば、"abcde"と"abcxy"という符号のハミング距離は2です）。距離の小さい、つまり似たような符号があると、ちょっとしたエラーが発生するだけでそれらを間違えてしまいそうです。なので、 $d$ はできるだけ大きくなるように符号化したい気持ちになります。この最小距離は符号の誤り訂正能力に直結し、最小距離 $d$ の符号が訂正できるのは $d/2$ 未満の誤りだけです。直感的には、もとの符号から出発して、他の符号とのちょうど中間のところまで来てしまうと、どこに帰ればよいかわからなくなってしまう、というイメージです。

しかし当然ながら $d$ はどこまでも大きくできるわけではありません。 $k$ 文字を $n$ 文字に符号化する符号において、 $d$ は $n-k+1$ 以下になります。これをシングルトン限界と呼びます。

ここがすごいところなのですが、RS符号は、最小距離がシングルトン限界 $n-k+1$ に等しい「最大距離分離符号」です。だから訂正能力が高いです。

RS符号は $(n-k+1)/2$ 未満のエラーであれば訂正できます。言い換えれば、訂正できるエラーの最大数は $\lfloor (n-k)/2 \rfloor$ 個です。

RS符号とは

ずばり、こういうものです。

$x _ {1}, ..., x _ {n}$ を体 $\mathbb{F} _ {q}$ の異なる元とする． $k\le n$ に対し， $\mathbb{F} _ {q} [x]$ における次数が $k$ より小さいすべての多項式の集合 $\mathbb{P} _ {k}$ を考える．リード・ソロモン符号 (RS)符号 (Reed-Solomon code)は，以下の符号語の全体からなる．

$(f(x _ {1}), f(x _ {2}), ..., f(x _ {n}) )$ 　（ただし $f\in \mathbb{P} _ k$ ）

（参考書籍[1] p.63, 定義 4.1.1）

長さ $k$ のメッセージを長さ $n$ の符号語に変換する話をしています。体とは（ざっくり言うと）四則演算ができる世界のことです。その世界における多項式のうち、次数が $k$ 未満のもの（例えば $k=3$ なら $x ^ {3}$ 以上の項がないもの、 $ax ^ {2}+bx+c$ や $bx+c,ax ^ {2}+c, c$ など）をひとつ用意します。また、その世界から $n$ 個の元（要素、≒数）を持ってきて、その多項式に代入していった結果を並べて $(f(?), f(!),...)$ のような $n$ 次元ベクトルを作ります。これがひとつの符号語になります。考えうるすべての多項式に同じことをして得られる符号語の集合がRS符号ということです。

あるメッセージと対応する符号語の関係をざっくり言うと、符号語は「メッセージを多項式関数に置き換え、1から順に数を代入した結果の一覧」です。符号化はこの一覧表の作成をやるだけです。復号は結果一覧から多項式関数を特定する作業になります。多項式が違えば結果はある程度異なるものになるので、結果一覧に誤植があっても、少しなら問題ありません（誤り訂正）。

$q$ は素数のべき乗なら何でもよいのですが、実用上は $q$ を2のべき乗、 $n=q-1$ 、 $x _ {i}=\alpha ^ {i-1}\ (1\le i\le q-1)$ とすることが多いです。 $\alpha$ は体 $\mathbb{F} _ {q}$ の原始元です。今回は全体的にこれにならうこととし、 $q=2 ^ {8}=256$ とします。ASCIIコードでは1bitが1文字に対応し、1bitは256種類の数が表現できるため、体の元と文字の対応付けも簡単です。

なぜ体が出てくるのか

なぜ文字をそのまま0から255までの整数だけの世界（256で割った余りの世界/mod 256の世界）で考えようとしないのかですが、それだと困ることがあります。加減算と乗算はよい（結果を256で割った余りに置き換えれば0~255の範囲に収め直せる）のですが、割り算が困ります。なんと割り算ができないことがあります（例えばmod 256の世界に2x=1を満たすxは存在しません、これはmod 非素数の世界を考えた時に起こります）。この先で割り算を使うので、四則演算ができる体、特に要素が有限で閉じている（演算結果も必ずその体の要素な）もの（=有限体/ガロア体）の世界で考える必要があります。

原始元とは

体の元であって、 $\mathbb{F} _ {q}$ の原始 $(q-1)$ 乗根（= $(q-1)$ 乗してはじめて1になるもの）です。これにより、体の0以外の元は原始元 $\alpha$ を用いて $\alpha ^ {i}$ の形で表せます。

RS符号の世界（ガロア体）について

$\mathbb{F} _ {q}$ は要素数が $q$ 個と有限で、四則演算について閉じている（演算結果もまた体の元）という性質をもつ体で、これをガロア体とか有限体と呼びます。

ガロア体の例

最も簡単なガロア体は、ずばり、1bitの世界、0と1の世界です。 $\mathbb{F} _ {2}$ です。加算や乗算は整数の時とほぼ一緒ですが、2で割った余りを取ります。実際には、1+1=0であることだけが違いです。減算と除算もその逆を考えればよいです。注意すべきは0-1=1なことぐらいです。あと、この世界でもゼロ除算はダメです。

この加減算は排他的論理和（xor）という演算で、0どうし1どうしでは0に、0と1では1になります。プログラミング言語でもこの演算が（大抵）サポートされていて、（大抵）^という演算子を使います。

同様に、素数 $p$ に対してガロア体 $\mathbb{F} _ {p}$ を作ることができます。加算と乗算は $p$ で割った余りになります。

ガロア拡大体（今回扱う体の話）

ガロア体 $\mathbb{F} _ {p}$ を拡大し、 $\mathbb{F} _ {p}$ の $m$ 次元ベクトルを元（全 $p ^ {m}$ 個）とするガロア拡大体 $\mathbb{F} _ {p ^ {m} }$ を作ることができます。面倒なのでここからは $p=2$ で進めます。

$\mathbb{F} _ {2}$ を拡大し、 $\mathbb{F} _ {2 ^ {m} }$ を作ります。 $m$ 次元ベクトルの要素を $\mathbb{F} _ {2}$ 上の $(m-1)$ 次多項式の係数だと思うことにします。 $\{0,1\}=1, \{1,0\}=x,\{1,1\}=x+1$ のような感じです。加算はベクトルの要素ごとに和をとればよいです。つまりはbitごとにxorを計算するということです。乗算は少し大変で、 $a\times b\ (a,b\in \mathbb{F} _ {2 ^ {m} })$ は、 $a,b$ が多項式 $a(x), b(x)$ に対応するとして、 $a(x)b(x)\ mod f(x)$ に対応する元（ $f(x)$ は $\mathbb{F} _ {2}$ 上の $m$ 次の既約な（=因数分解できない）多項式）とします。

（ここから今回の $m=8$ での話になります） $x ^ {8}+x ^ {4}+x ^ {3}+ x ^ {2}+1$ は原始多項式という既約多項式で、この根（多項式に代入すると0となる解）は原始元 $\alpha$ になり、体の元は0を除いて $\alpha$ のべき乗で表せます。 $\alpha$ は原始多項式の根なので $\alpha ^ {8} = \alpha ^ {4}+\alpha ^ {3} + \alpha ^ {2} +1$ です。多項式的には、乗算のmodをこの多項式にします。 $x ^ {8}=x ^ {4}+x ^ {3}+ x ^ {2}+1$ です。

結局、こんな感じになります。

$\mathbb{F} _ {2 ^ {8}}$ の元は $0$ と $\alpha ^ {i}(0\le i \le 254, \alpha ^ {0} = 1)$ の合計256個です。
$\alpha ^ {255} = \alpha ^ {0} = 1$ です。
$GF(2 ^ {8})$ における四則演算は以下のようになります。乗除算については上の性質を使います。
- 加算/減算：bitwise xor（排他的論理和）
- 乗算： $\alpha ^ {a} \times \alpha ^ {b} = \alpha ^ {(a+b)\ mod 255}$
- 除算： $\alpha ^ {a} / \alpha ^ {b} = \alpha ^ {(a-b)\ mod 255}$

体の元とbit列、多項式は次の表のような対応をします。

元	対応するbit列	多項式としての解釈
$0$	00000000	0
$1=\alpha ^ {0}$	00000001	1
$\alpha$	00000010	$x$
$\alpha ^ {2}$	00000100	$x ^ {2}$
$\alpha ^ {3}$	00001000	$x ^ {3}$

$\alpha ^ {7}$	10000000	$x ^ {7}$
$\alpha ^ {8}=\alpha ^ {4}+\alpha ^ {3} + \alpha ^ {2} +1$	00011101	$x ^ {4}+x ^ {3}+x ^ {2}+1$
$\alpha ^ {9}=\alpha ^ {8} \times \alpha = \alpha ^ {5}+\alpha ^ {4}+\alpha ^ {3}+\alpha$	00111010	$x ^ {5}+x ^ {4}+x ^ {3}+x$

$\alpha ^ {i} =\alpha ^ {i-4}+\alpha ^ {i-3} + \alpha ^ {i-2} +\alpha ^ {i-8}\ (8\le i\le 254)$

符号化

前提を再掲します： $q=2 ^ {8}=256,n=q-1=255,x _ {i}=\alpha ^ {i-1}\ (1\le i\le n)$

符号化は、長さ $k$ のメッセージ $m$ を $\mathbb{F} _ {2 ^ {8} }$ 上の $(k-1)$ 次多項式 $m(x)=\sum _ {i=1} ^ {k}m _ {k-i} x ^ {k-i}$ とみて、 $(m(x _ {1}), m(x _ {2}), ..., m(x _ {n}) )$ とすればよいです。実際には、生成行列 $G$ をかけるという操作によって符号化します。

生成行列と検査行列

生成行列 $G$ は $k\times n$ の行列で、これを用いてメッセージの符号化を行います。符号の定義より、 $G _ {i,j}=x _ {j} ^ {i-1}\ (1\le i\le k, 1\le j\le n)$ です。今回の場合は $x _ {j} =\alpha ^ {j-1}$ ですので、 $G _ {i,j} = \alpha ^ {(i-1)(j-1)}\ (1\le i\le k, 1\le j\le n)$ となります。DFT行列に似ています。というか実際フーリエ変換そのものっぽいです（参考[6]）。

検査行列 $H$ は $Gh ^ {T}=0$ をみたす（線形独立な）行ベクトル $h$ を $n-k$ 個並べた $(n-k)\times n$ の行列です。上の生成行列 $G$ に対応する検査行列は $H _ {i,j} =\alpha ^ {i(j-1)}\ (1 \le i\le n-k,1\le j\le n)$ です。

理由

実際に計算してみるとわかります。 $M=GH ^ {T}$ とします。 $M _ {i,j} = \sum _ {k} G _ {i,k}H^{T} _ {k,j}=\sum _ {k} \alpha ^ {(i-1)(k-1)} \alpha ^ {j(k-1)}=\sum _ {k=1} ^ {n} (\alpha ^ {i+j-1} ) ^ {k-1}$ です。 $1\le i+j-1 \le n-1$ より $\alpha ^ {i+j-1}\neq 1$ なので、等比数列の和の公式より $M _ {i,j}=( (\alpha^ {i+j-1}) ^ {n}-1) / (\alpha ^ {i+j-1}-1)$ ですが、分子の $(\alpha^ {i+j-1}) ^ {n}$ に注目すると $(\alpha^ {i+j-1}) ^ {n}=(\alpha ^ {n}) ^ {i+j-1}=1 ^ {i+j-1}=1$ であるため、 $M _ {i,j}=0$ となります。

$w'$ に誤りがなく $w'=w(=mG)$ ならシンドロームは $s=wH ^ {T}=mGH ^ {T}=0$ になるはずですが、そうならない場合はエラーがあることになります。そのときのシンドロームはどうなっているかですが、あるエラー $e$ があって $w'=w+e$ と表せるので、ここに検査行列をかけた値 $w'H ^ {T}=wH ^ {T}+eH ^ {T}=eH ^ {T}$ になっています。

エラー検出：ピーターソン法

$eH ^ {T}=s$ でした。これを解けば $e$ が分かるのですが、 $n=|e|>|s|=n-k$ より未知数に対して方程式の数が足りません。これを何とかします。

復号のアイデアは以下の通りだそうです。正直僕はあまり理解できてないです。

$r=c+e$ を受信語とし， $w(e)\le t=\lfloor (n-k)/2 \rfloor$ を仮定する．復号法のアイデアは，2変数多項式

$Q(x,y) = Q _ {0}(x)+yQ _ {1}(x)\in \mathbb{F} _ {q}[x,y] \backslash \{0\}$

であって，以下の条件を満たすものを決定することである．

$Q(x _ {i}, r _ {i})=0, i=1,...,n$

$\deg (Q _ {0})\le n-1-t$

$\deg (Q _ {1})\le n-1-t-(k-1)$

（参考書籍[1], p.65）

本文中で使ってきた $w',w$ はそれぞれ $r,c$ に対応しています。条件1より方程式が $n$ 個できますが、条件2,3より $Q _ {0}, Q _ {1}$ の未知の係数はそれより多くなるので、このような多項式があることがいえます。

また、この方程式に関して次の定理が成り立ちます。

送信語が $g(x)$ で生成され，そして，誤りの個数が $d/2$ より小さければ， $g(x)=-Q _ {0}(x)/Q _ {1}(x)$ が成り立つ．

（参考書籍[1], p.65, 定理4.2.2）

なぜ？

参考書籍[1], p.65, 定理4.2.2の証明に基づく説明です。

恒等的に $Q(x,g(x) )=0$ であることを示せばあとは単純な式変形で定理を導けるので、それを示します。

「送信語が $g(x)$ で生成され」るとは、 $c _ {i}=g(x _ {i})\ (1\le i\le n)$ ということです。条件1より $Q(x _ {i}, g(x _ {i}) + e _ {i})=0$ です。エラー個数は $t=\lfloor (n-k)/2\rfloor$ 以下です。ということは $e$ のうち $n-t$ 個以上の要素がゼロなので、それらに関しては $Q(x _ {i}, g(x _ {i}) )=0$ が成立しています。これにより $Q(x,g(x) )$ は少なくとも $n-t$ 個の根をもつことになります。ところがこの多項式の次数は条件2,3より高々 $n-t-1$ で、それより多い $n-t$ 個の根をもっています。ということは、 $Q(x,g(x) )=0$ です。

エラーの位置は $Q _ {1}(x)$ の根の中にあります（ $Q(x,y)=Q _ {1}(x)(y-g(x) )$ と書けます。エラー位置 $x _ {i}$ が根でないと仮定した場合、 $x=x _ {i}, y=r _ {i}$ を代入すると $g(x _ {i})=c _ {i}\neq r _ {i}$ より $y-g(x)\neq 0$ 、仮定より $Q _ {1}(x)\neq 0$ で $Q(x,y)\neq 0$ になってしまいます）。この性質から $Q _ {1}(x)$ は誤り位置多項式と呼ばれます。

以後 $Q _ {0}, Q _ {1}$ の最大次数を $l _ {0} = n-1-t, l _ {1} = n-1-t-(k-1)$ と表します。

さて、ここから復号までには何通りかの方法があります。ここではピーターソン法と呼ばれる方法を取り上げます。この方法ではまずエラー位置特定（ $Q _ {1}(x)$ の特定）を行い、次に $e$ を求めます。

第一段階、エラー位置特定です。次の定理を用います。（ $n-k$ が奇数であることを仮定しています。偶数なら $l _ {1}$ は $l _ {1}-1$ になります）

誤り位置多項式 $Q _ {1}$ の係数は，次の連立一次方程式の解である．
（参考書籍[1], p.72, 定理4.4.1）

本当に？

参考[1] pp.72-74の証明にもとづいて説明をします。

$Q$ のみたす条件1から、次のような関係が成立します。

両辺に左から次の行列をかけると、（生成行列に検査行列をかけたのと同様の原理で）左辺第一項が消えます。

残った左辺第二項は、次のように変形できます。

あとは、

であることを確認します。行列を分解したことで計算が簡単になって、右辺の計算結果は

になることが確認できます。いっぽう、左辺については $S=rH ^ {T}=Hr ^ {T}$ で計算したことを思い出すと、

であることから一致します。

証明は省きますが、 $Q _ {1}$ がこの方程式の解ならば最初の方程式の解 $Q _ {0}$ も存在します。

$Q _ {1}$ がわかればその根を総当たりで求めることができ、誤り位置が判明します。

第二段階です。 $eH ^ {T}=s$ でした。エラーの数は $\lfloor (n-k)/2\rfloor$ 以下なので、 $e$ はそれらを残してあとは全部ゼロです。未知数が減ったので方程式を解くことができ、 $e$ が定まります。

最後に方程式 $mG=w'-e$ を解いて、元の $m$ が復元できます。

実装

さて、だいぶ理論の話が長くなりましたが、ここから実装していきます。実装してから理論を書いたので若干記号が違ったりします。

何をするか

ユーザの入力：ASCII文字列 $m$ とデータ破壊数 $t$

プログラムから見た入力： $w'$ 、つまり破壊された $w$

出力：もとのメッセージ $m$

手抜きPoint

今回はエラーがひどすぎて復元できない場合を想定しません。めんどくさいし、復号できないことがわかっても嬉しくないし、直せるものがきちんと直せると分かればいいので。
ヒントとして、元のメッセージの文字数は復号プログラムに教えます。固定長になるように分割したり付け足したりするか、特定の書式を用意したりすれば解決するのですが、そこは面白くない場所なので。
計算資源に関する制約はあまり意識しません。最良の方法であることよりも実装の手軽さなどを優先し、一般的な家庭用のパソコンで実行したときに即座に実行が完了する程度に動けばよいことにします。

入力

入力された一行のメッセージ $m$ をgetline()で受け取ります。cinは空白で区切られてしまうので不便です。符号化したメッセージ $w$ の長さを $N = 2 ^ {8} -1 = 255$ とするため、メッセージ長 $K \lt N$ である必要があります。もっというと破壊したいので $K\le N-2$ がよいですが、ここではとりあえず200文字以内のメッセージなら受理ということにします。自分が使うプログラムなので基本的には性善説に基づいた実装ですが、いちおう長すぎる入力や空の入力は弾くようにします。僕の環境ではどうにかなるようですが、マルチバイト文字は入れないでほしいです。

ついでに $w$ の破壊文字数 $t$ も入力してもらいます。最大で $t _ {0} = \lfloor (N-K)/2 \rfloor$ 文字破壊できます。それより上の要求は撥ね退けます。わけの分からん入力に関しては、文字列をint型の整数に変換する関数stoi()様の解釈を正義とします。エラーの場合は再入力です。データを破壊しないこと、つまり $t=0$ は許容します。

const string input = Input::getRangeString("Input the string (<="+to_string(maxK)+" letters)",1,maxK+1);
const int K = input.size();
const int maxChange = (N-K)/2;
const int change = Input::getRangeNum("Input the number of data to be corrupted (<="+to_string(maxChange)+")",0,maxChange+1);

Input::xxx

namespace Input{
  bool inRange(const int n, const int l, const int r){
    return l<=n&&n<r;
  }
  bool isRangeNumStr(const string& s, const int l, const int r){
    int a;
    try{ a = stoi(s); }
    catch(...){ return false; }
    return inRange(a,l,r);
  }
  string getRangeString(const string& message, const int l, const int r){
    string S;
    while(true){
      cout << message << endl;
      getline(cin,S);
      if(inRange(S.size(),l,r)) return S;
    }
  }
  int getRangeNum(const string& message, const int l, const int r){
    string numString;
    while(true){
      cout << message << endl;
      getline(cin,numString);
      if(isRangeNumStr(numString,l,r)) return stoi(numString);
    }
  }
}

拡大体

$\mathbb{F} _ {2 ^ {8} }$ 上の元を表す型を用意しましょう。これを構造体alphaとして定義します。中身は8bitならいいのでunsigned char型にします。signedにするとろくなことにならない（筆者は論理シフトと算術シフトの違いで痛い目を見たことがあります）のでunsignedです。ここに加減乗除を定義します。演算子オーバーロードというやつです。これにより（ソースコード上で）普通の整数と同様に扱うことができます。

配列A[256]を用意し、Aを参照することで $\alpha ^ {i}$ が即座に求められるようにします。 $A[i]=\alpha ^ {i}\ (0\le i\le 254), A[255]=0$ のような対応にしておきます。プログラムを起動したときに初期化関数を呼ぶことにし、配列Aは初期化関数内部で漸化式的に求めます。

加減算はbitwise xorなので^演算子を用いてシンプルに書けます。

乗除算も頑張ればbit演算でできる気がしますが、簡単な方法でやります。まず、0が絡む場合は場合分けして個別に処理します。それ以外は肩の数字のmod 255での加減算になります。各元は自分が $\alpha$ の何乗なのかをわかっている必要がありますが、これも初期化の際にカンペを作っておいてそれを参照することにしましょう。

alpha

剰余をとるのは負の数を正の数に直してから。一部省略していますが、実際は複合代入演算子（+=など）も定義しています。

struct alpha{
  unsigned char v;
  alpha():v(0){}
  alpha(unsigned char a):v(a){}
  const alpha operator+ (const alpha& r) const{ return alpha(v^r.v); }
  const alpha operator- (const alpha& r) const{ return alpha(v^r.v); }
  const alpha operator* (const alpha& r) const{
    if (*this==0 || r==0) return 0;
    int expL = revA[v], expR = revA[r.v];
    return A[(expL+expR)%255];
  }
  const alpha operator/ (const alpha& r) const{
    assert(r!=0);
    if(*this==0) return 0;
    int expL = revA[v], expR = revA[r.v];
    return A[(expL+255-expR)%255];
  }
  static alpha pow(int x){
    if(x<0) x+=((-x)/255+1)*255;
    return A[x%255];
  }
  static array<alpha,256> A;
  static array<unsigned char,256> revA;
  static void initialize(){
    A[255]=0;
    revA[0]=255;
    for(unsigned char i=0;i<8;i++){
      A[i]=1u<<i;
      revA[A[i].v]=i;
    }
    for(unsigned char i=8;i<255;i++){
      A[i]=A[i-4]+A[i-5]+A[i-6]+A[i-8];
      revA[A[i].v]=i;
    }
  }
}

あと、後々面倒なのでalphaのベクトルや行列、すなわちvector<alpha>やvector<vector<alpha>>は略記ができるようにしておきます。

using Vec = vector<alpha>;
using Matrix = vector<vector<alpha>>;

符号化

生成行列

$G _ {i,j}= \alpha ^ { (i-1)(j-1) } \ (1\le i\le K, 1\le j\le N)$

$K\times N$ の行列です。これを用いてメッセージの符号化を行います。行列（2次元vector）を作るだけなのでプログラミング的な面白さはないです。

const Matrix G = RS::MatrixGen::Generate(K,N);

RS::MatrixGen::Generate

Matrix RS::MatrixGen::Generate(const int K, const int N){
  Matrix G(K,Vec(N));
  for(int k=0;k<K;k++) for(int n=0;n<N;n++) G[k][n] = alpha::pow(k*n);
  return G;
}

符号語生成

$w=mG$ を計算します。 $|w|=N$ です。行ベクトルと行列を受け取ってその積（行ベクトル）を返す関数を作り、それを呼び出してあげます。

const Vec message = RS::MatrixGen::message(K,input);
const Vec word = Math::mul(message,G);

RS::MatrixGen::message, Math::mul

Vec RS::MatrixGen::message(const int K, const string& S){
  assert((int)S.size()==K);
  Vec word(K);
  for(int i=0;i<K;i++) word[i] = S[i];
  return word;
}

Vec Math::mul(const Vec& a, const Matrix& M){
  int H = M.size(), W = M[0].size();
  Vec res(W);
  for(int i=0;i<W;i++) res[i] = 0;
  for(int y=0;y<H;y++) for(int x=0;x<W;x++) res[x] += a[y]*M[y][x];
  return res;
}

データ破壊

$w$ の成分を指定の数だけランダムな別の値に書き換えた $w'$ を用意します。面倒なので書き換える場所と内容はプログラムが乱数で決定します。

const Vec brokenWord = RS::MatrixGen::errorMessage(word,change);

RS::MatrixGen::errorMessage

そうする意味は特にないのですが、アルゴリズムの癖として、連続した場所を破壊する傾向にあります（既に破損させた場所を再選択した場合に右方向で最も近い未破損データを壊すため）。

mt19937 mt(chrono::system_clock::to_time_t(chrono::system_clock::now()));
uniform_int_distribution<int> distN(0,N);
uniform_int_distribution<unsigned int> dist255(0,256);

Vec RS::MatrixGen::errorMessage(const Vec& word, const int change){
  Vec brokenWord(word);
  vector<bool> changeFlag(N,false);
  for(int i=0;i<change;i++){
    int pos = distN(mt);
    while(changeFlag[pos]) pos=(pos+1)%N;
    changeFlag[pos] = true;
    alpha after;
    do{ after = dist255(mt); } while(after==brokenWord[pos]);
    brokenWord[pos] = after;
  }
  return brokenWord;
}

符号語比較

人間が見てわかるようにここで $w$ と $w'$ の比較を行いたいと思います。

符号語表示

alpha型の中身はunsigned char型の変数ですが、文字としてそのまま出力することはできません（ASCIIが対応するのは0-127までなので範囲外になることがあり、かつ範囲内にも改行やヌル文字などの表示できない文字=制御文字がある）。unsigned charが扱う0~255の値を2桁の16進数（適宜0埋め）に変換して表示します。

表示

ostream& operator << (ostream& ost, const alpha& a){
  return ost << hex << setw(2) << setfill('0') << (int)a.v << setfill(' ') << dec;
}
ostream& operator << (ostream& ost, const Vec& v){
  for(auto itr=v.begin();itr!=v.end();itr++) ost << *itr;
  return ost;
}

比較

上の要領で符号語を人に読める文字にできますが、ふたつの符号語を上下に並べて表示したとして、それを比較するのはおよそ人間のやることではありません。違っているところを目立たせるために、文字と背景の色を反転させます。coutで出力するとき、\033[7mを噛ませることで反転、\033[0m=\033[mで戻すことができます（参考：もう一度基礎からC言語第47回特殊な画面制御～コンソール入出力関数とエスケープシーケンスエスケープシーケンスによる画面制御）。

compare(word,brokenWord,"Original","Broken");

compare

せっかくなので符号語以外のvectorも対応できるようにしておきます。

template <typename T>
void compare(const vector<T>& A, const vector<T>& B, const string& s1, const string& s2){
  assert(A.size()==B.size());
  const string flip = "\033[7m";
  const string endFlip = "\033[m";
  cout << endl;
  cout << "------compare------" << endl;
  int N = A.size();
  vector<int> dif;
  cout << s1 << ':' << endl;
  cout << A << endl;
  cout << s2 << ':' << endl;
  for(int i=0;i<N;i++){
    if(A[i]==B[i]) cout << B[i];
    else{
      dif.push_back(i);
      cout << flip << B[i] << endFlip;
    }
  }
  cout << endl;
  if(dif.empty()) cout << "No error." << endl;
  else cout << "Error at:" << dif << endl;
  cout << "-------------------" << endl << endl;
}

検証

検査行列

$H _ {i,j} = \alpha ^ {(i-1)j} \ (1\le i\le N,1\le j\le N-K)$

誤り検出に使用します。転置した状態で使うので最初から転置した状態の $N\times (N-K)$ の行列にします。式にしたがって生成するだけなので特にいうことはありません。

const Matrix H = RS::MatrixGen::Check(N,N-K);

RS::MatrixGen::Check

Matrix RS::MatrixGen::Check(const int N, const int NK){
  Matrix H(N,Vec(NK));
  for(int n=0;n<N;n++) for(int x=1;x<=NK;x++) H[n][x-1] = alpha::pow(n*x);
  return H;
}

シンドローム

$s = w' H$ でシンドロームを計算します。符号語を生成したときと同じ要領なので同じ関数を使いまわします。結果が $s=0$ の場合はそのままメッセージを復号、そうでない場合は復号の前にエラー処理を行います。

const Vec syndrome = Math::mul(brokenWord,H);
Vec restoredWord(word);
if(!Math::isZeroVec(syndrome)){
  const vector<int> errorPos = RS::specifyError(maxChange,syndrome);
  restoredWord = RS::restoreWord(brokenWord,errorPos,H,syndrome);
}
else { cout << "Data is not broken." << endl; }

エラー処理

$eH=s$ を解けば $e$ がわかりますが、このままでは解けないので先に $e$ から零の要素を除去し、未知の変数を減らします。

RS::specifyError

vector<int> RS::specifyError(const int maxChange, const Vec& syndrome){
  Matrix S = RS::MatrixGen::Syndrome(maxChange,maxChange+1,syndrome);
  const int rnk = Math::gauss_jordan(S);
  const Vec sigma = MatrixGen::sigma(rnk,S);
  return RS::findErrorPos(sigma);
}

シンドローム行列と誤り位置多項式

$S _ {i} = s [ i;\ t _ {0}+i ]$ の $t _ {0} \times (t _ {0}+1)$ 行列になります。生成行列や検査行列に比べると少しややこしいです。

RS::MatrixGen::Syndrome

Matrix RS::MatrixGen::Syndrome(const int H, const int W, const Vec& syndrome){
  Matrix S(H,Vec(W));
  for(int y=0;y<H;y++) for(int x=0;x<W;x++) S[y][x] = syndrome[y+x];
  return S;
}

$S$ に掃き出し法を使い、対角成分に $1$ が $t$ 個並ぶように変形します（const int rnk = Math::gauss_jordan(S);）。下の $t _ {0} - t$ 行はゼロになります。エラー個数 $t=rank(S)$ が特定できます。掃き出し法についてはもう少し後で。

掃き出し法を使った後の $S$ から誤り位置多項式 $\sigma (x) = x ^ {t} + \sigma _ {1} x ^ {t-1} + \sigma _ {2} x ^ {t-2} + ... + \sigma _ {t-1} x + \sigma _ {t}$ を得ます。 $S$ の $t+1$ 列目の成分上から $t$ 個が、 $\sigma _ {t}, \sigma _ {t-1}, ..., \sigma _ {1}$ となります（const Vec sigma = MatrixGen::sigma(rnk,S);）。

RS::MatrixGen::sigma

Vec RS::MatrixGen::sigma(const int rnk, const Matrix& S){
  Vec sig(rnk+1);
  for(int y=0;y<rnk;y++) sig[y]=S[y][rnk];
  sig[rnk] = 1;
  return sig;
}

理論と言ってることが微妙に違いませんか

文献[4]を参考に、エラー個数と誤り位置多項式を一気に求める方法をとっています。

ここから適当なことを言うのですが、根を求めるだけなので、誤り位置多項式には適当な定数をかけてよく、最高次の係数を1にすることができます。その状態から変形すると掃き出し法に持っていける感じがしてきます。

で、見つかるまで次数をあげていきながら多項式を求めている、という感じの処理になっている（のだと思っていますが、厳密な原理は参考元をあたってください）。

$\sigma (x) = 0$ をみたす $x$ を求めます。実際に $x$ に元を代入して全探索で求めます（RS::findErrorPos(sigma);）。

RS::findErrorPos

vector<int> RS::findErrorPos(const Vec& sigma){
  vector<int> errorPos;
  for(int j=0;j<N;j++){
    if(Math::polynomial(sigma,alpha::pow(j))==0) errorPos.push_back(j);
  }
  return errorPos;
}

alpha Math::polynomial(const Vec& coef, const alpha a){
  alpha x(1), sum(0);
  for(const alpha& c:coef){
    sum += c*x;
    x *= a;
  }
  return sum;
}

今回は排除しているので心配いりませんが、一般に $t=t _ {0}$ の場合、エラーが $t _ {0}$ 個とは限りません（もっと多い可能性があります）。 $t _ {0}>t$ の場合は、この後出てくる誤り位置多項式の根が見つからなくなるようです（参考：文献[4]）。

掃き出し法

掃き出し法の実装がプログラム全体で最も難しいと思います。とはいっても、今回は誤差やオーバーフローの心配がなく、その点ではやや簡単ともいえます。基本的には人間がやるような感じをプログラムに書き直すだけです。

具体的な手続き

大まかな流れを記述します。行列 $A$ に対して掃き出し法を行います。左の列から（ $j=0$ から）順番に、繰り返し処理していきます。

その列にゼロでない要素 $A _ {i,j}\neq 0$ をもつ（、固定されていない）行 $i$ をひとつ探し、 $\alpha _ {x} := A _ {i,j}$ とします。
見つかったら、行 $i$ の要素をすべて $\alpha _ {x}$ で割ります。
行 $i$ を除くすべての行について、 $\alpha _ {y _ {k}}:=A _ {k,j}$ とし、各要素 $A _ {k,l}$ に対して、 $A _ {k,l} :=A _ {k,l}-\alpha _ {y _ {k} } \times A _ {i,l}$ とします。
行 $i$ を（固定されていない）いちばん上の行とswapし、その位置に固定します。

プログラムは雑に書くとそれ以降も参照する必要のある変数をうっかり書き換えてしまうバグを埋め込みやすいです。

Math::gauss_jordan

一部の処理の順番を入れ替えてあります。

int Math::gauss_jordan(Matrix& A){
  const int H = A.size(), W = A[0].size();
  int rnk=0;
  for(int x=0;x<W;x++){
    for(int y=rnk;y<H;y++){
      if(A[y][x]!=0){
        if(y!=rnk){
          swap(A[y],A[rnk]);
          y=rnk;
        }
        const alpha ayx = A[y][x];
        for(int xx=0;xx<W;xx++){
          A[y][xx] /= ayx;
        }
        for(int yy=0;yy<H;yy++){
          if(yy==y) continue;
          const alpha f = A[yy][x];
          for(int xx=0;xx<W;xx++){
            A[yy][xx] -= A[y][xx]*f;
          }
        }
        rnk++;
        break;
      }
    }
  }
  return rnk;
}

エラー除去

未知数が減ったので、方程式 $eH=s$ が解けるようになりました。 $e,H,s$ から必要な部分（ $e$ の非零要素、それとの掛け算が発生する $H$ の行）を抜き出して、掃き出し法で方程式を解けば $e$ がわかります。 $w=w'-e$ として、破損前のデータ $w$ を得ることができます。

RS::restoreWord

Vec RS::restoreWord(const Vec& brokenWord, const vector<int>& errorPos, const Matrix& H, const Vec& syndrome){
  Matrix ErrorEquation = RS::MatrixGen::ErrorEq(errorPos,H,syndrome);
  Math::gauss_jordan(ErrorEquation);
  const Vec error = RS::MatrixGen::error(errorPos,ErrorEquation);
  return Math::sub(brokenWord,error);
}

Matrix RS::MatrixGen::ErrorEq(const vector<int>& errorPos, const Matrix& H, const Vec& syndrome){
  int error = errorPos.size();
  Matrix E(error+1,Vec(H[0].size()));
  for(int i=0;i<error;i++) E[i] = H[errorPos[i]];
  E[error] = syndrome;
  return Math::transpose(E);
}

Vec RS::MatrixGen::error(const vector<int>& errorPos, const Matrix& E){
  vector<alpha> err(N,0);
  int e = errorPos.size();
  for(int i=0;i<e;i++) err[errorPos[i]] = E[i].back();
  return err;
}

メッセージ復号

遥か昔のことなので忘れそうですが、 $w=mG$ というふうに符号化したのでした。 $w,G$ はわかっているので連立方程式を解きます。これも掃き出し法で解けます。最後に復号したメッセージを出力して比較すれば無事return 0;です。

const string output = RS::decode(G,restoredWord,K);
cout << "Restored your input: " << output << endl;
compare(vector<char>(input.begin(),input.end()),vector<char>(output.begin(),output.end()),"Original message","Recovered message");
return 0;

RS::decode

string RS::decode(const Matrix& G, const Vec& word, const int K){
  Matrix Gw = Math::AugCoef(G,word);
  Math::gauss_jordan(Gw);
  string message = "";
  for(int i=0;i<K;i++) message += Gw[i].back().toChar();
  return message;
}

ソースコード

全体実装はそれなりに長いのでリンクを貼ります。

reedsolomon.cpp - Google ドライブ

動作確認

プログラムを動かしてみて、誤り訂正できていることを確かめました。僕の環境はC++14ですがC++11でも動くと思います。実行する場合は画面処理（文字色と背景色の交換）があるので対応できる環境で。Windowsのコマンドプロンプトではうまくいきません。どうにかできるとは思いますが、おとなしくPowerShellとかWSLを使ったほうが手っ取り早いと思います。

おわりに

リード・ソロモン符号（のおいしいところだけ、ですが）を実装してみました。何かの役に立つわけでもないし、画面に文字が並ぶだけで見栄えもしませんが、我々の生活を陰で支えるアルゴリズムに目を向けて、実装してみるのもけっこう面白いと思います。それでは。

XXXXXX_XXX_XXXXXXX_XX_XXX_XXXX_XXXXX

209be895e02f0bfdfc911a984fc51f516687cb55ce343c3c38b806c0a72b0a371f012706e6adc468c2a1ad20891f107d647b3b5935e9137aaace5db2fcd1164b9df34c3102c10dbd317af006ec43ea45e5c5c0ef28c7dfb606071457244d31c345a83df373d6aac8f0139b05bb84dd4aa3fb21c3742d5fe3dcefec442a1f01e2f5ba4d11f909ffe5ebdbfef722ad5fc15e0ed2160d5753ae08c9de9fa58ebd9c46fc8753d2cb5c92cf15f638d0796efddfd60b3fd89e5f6bdcc25940777b32343fdc677fa4bb98d42cefe74c18f15c190846f6fb32dd0489d296154aedd499c140936a8768e04470e70a5a0252413c6855df9604b895ccf44fee14fa103b5a

おまけ

当初は高専要素皆無だった導入部分（こっちのほうが自然でいいと思う）

情報は、時に欠落や破損をしてしまうものです。手の甲のメモが消えて読めなくなってしまったり、「ニンニク増しで」って言ったのに「ニンニク無し」が出てきたり、（この先は血で汚れていて読めない）

しかし、頑張ればなんとかなるケースも多々あります。なんか文字化けし縺ｦ繧けどギリギリ読める、誤字脱字がるけど陽はこういう意味だよな、と人間の脳は情報を補完します。

同じことをコンピュータの世界でもやっています。QRコードは多少汚れていても読めるし、ディスクに多少傷がついていてもCDは再生できます。必要なデータだけでなく、いっしょに「データが壊れた時に気付いて、元に戻すための手がかり」も記録しておくことで欠けた部分を補っているのです。

具体的なやり方はいろいろあり、その一つにリード・ソロモン符号というものがあります。今回はこれを最発明します。といっても理論を思いつくのは不可能なので、意図的に壊したデータを頑張って復旧させるプログラムを書きたいと思います。

参考文献

[1] 誤り訂正符号入門（第2版）, イエルン・ユステセン, トム・ホーホルト原著, 阪田省二郎, 栗原正純, 松井一, 藤沢匡哉訳, 森北出版(2019)

[2] 有限体（ガロア体）の基本的な話 | 高校数学の美しい物語

[3] 有限体でのガロア拡大（その３） - zuruyasumineko2002’s blog

[4] QRコードの符号化・復号アルゴリズム解説：学習読み物としての試み, 佐藤創, 専修大学情報科学研究所所報,76,37-52 (2011-06-30)

[5] Reed-Solomon符号を15分で理解する【M3 Tech Talk 第187回】 - YouTube

[6] リード・ソロモン符号とその復号法, 松嶋敏泰, 映像情報メディア学会誌 Vol. 70, No. 4, pp. 571～575（2016） https://www.jstage.jst.go.jp/article/itej/70/7/70_571/_pdf

2022-12-05

Cookie Clicker（ブラウザ版）の配色が気に食わない

プログラミングゲーム

OBが1人で3枠は取りすぎな気がしますが、記事が余ってた&カレンダーが23:30時点でまだ空いてたので呉高専アドベントカレンダー 5日目に放流します。呉高専なーんも関係ない短編です。

ネタバレを多く含みますのでご注意ください。

Cookie Clickerの説明は省きます。知らない人はまずここにアクセスして遊んでください：https://orteil.dashnet.org/cookieclicker/

僕は攻略wikiは見るけどゲームを有利に進めるmodは入れない立場です。なのでmodには基本消極的なのですが…

ゲームを進めていくと遊べるstock marketというミニゲーム、その配色がどうも気に食わないのです。

valueが安いときに買って高いときに売るってだけなのですが、どれが安くてどれが高いのかマジで分かりません。

この配色を、modで変更します。こんな感じ。

背景色とかをいい感じにするほうがいい気がしますが、とりあえずこれでひとつ。

0: Cookie Clickerのmod

ブラウザのブックマーク機能を使い、ブックマークレットを実行することで実現します。

1: 値段に応じてpriceの色を変化させる

stock marketの各アイテムの情報はGame.Objects.Bank.minigame.goods[?]を見ます。wikiにある計算式に基づいて標準価格を計算し、その0.5倍から2倍にかけてpriceの色を赤から緑に変化させていきます。プレイしている感じ、あまりよい色の付け方ではないです。安いものは5倍ぐらい振れるし高いのは1.5倍にすらならないので。あとfff4f4とかfafff4とか人の目じゃ分からないです。

// 赤<-> 緑
let rgAry = ["#ff0000","#ff2b2b","#ff5555","#ff8080","#ffaaaa","#ffd5d5","#ffeaea","#fff4f4","#ffffff","#fafff4","#f4ffea","#eaffd5","#d5ffaa","#bfff80","#aaff55","#95ff2b","#80ff00"];
for(let [key,val] of Object.entries(Game.Objects.Bank.minigame.goods)){
  let price=val.val; // 現在の価格
  let center=10+val.id*10+Game.Objects.Bank.level-1; // 標準価格
  let left=center/2; // の半分
  let right=center*2; // の倍
  if(center < price){ // 高い
    let d=Math.min(8,Math.floor(8*(price-center)/(right-center))); // 1倍~2倍を8段階に分けた時、どれぐらいか
    val.valL.style.color=rgAry[8+d]; // 高いほど緑
  }
  else if(price < center){ // 安い
    let d=Math.min(8,Math.floor(8*(center-price)/(center-left))); // 0.5倍~1倍を段階づける
    val.valL.style.color=rgAry[8-d]; // 安いほど赤く
  }
}

2: 在庫があるときの色は緑から黄色に

価格表示に緑を使う関係から、在庫があるときのstockの色が緑だとまずいです。黄色にします。在庫が0になったら元の薄い白に戻すのも忘れずに。

if(val.stock!==0){
  val.stockBoxL.style.color = "#ffff00"; // 黄色
}
else{
  val.stockBoxL.style.color = "#ffffffb3"; // デフォルト
}

3: 隠し実績

Third-partyという隠し実績があり、特定のmodを入れると解除されます（通常プレイで使われない実績解除コマンドがmod内で叩かれてるってだけで、やろうと思えば開発者モードからの解除も可能なのですが。）。せっかくなので解除します。

if(Game.Achievements['Third-party'].won==0){
  Game.Win('Third-party');
}

4: 以上をまとめて毎秒更新

stock marketの情報は時々刻々変化するので一定時間で更新します。1秒に1回で十分です。

setInterval(
  function(){
    // ここに1-3の内容を書くことで毎秒実行される
  }
  ,1000
);

全体

setInterval(
  function(){
    if(Game.Achievements['Third-party'].won==0){
      Game.Win('Third-party');
    }
    if(Game.Objects.Bank.minigameLoaded){
      let rgAry = ["#ff0000","#ff2b2b","#ff5555","#ff8080","#ffaaaa","#ffd5d5","#ffeaea","#fff4f4","#ffffff","#fafff4","#f4ffea","#eaffd5","#d5ffaa","#bfff80","#aaff55","#95ff2b","#80ff00"];
      for(let [key,val] of Object.entries(Game.Objects.Bank.minigame.goods)){
        let price=val.val;
        let center=10+val.id*10+Game.Objects.Bank.level-1;
        let left=center/2;
        let right=center*2;
        if(center < price){
          let d=Math.min(8,Math.floor(8*(price-center)/(right-center)));
          val.valL.style.color=rgAry[8+d];
        }
        else if(price < center){
          let d=Math.min(8,Math.floor(8*(center-price)/(center-left)));
          val.valL.style.color=rgAry[8-d];
        }

        if(val.stock!==0){
          val.stockBoxL.style.color = "#ffff00";
        }
        else{
          val.stockBoxL.style.color = "#ffffffb3";
        }
      }
    }
  },
  1000
);

5: 使う

ネットにアップして、ブックマークレットを作成します

javascript:(function(){Game.LoadMod('httpsから始まるmodソースファイルの場所.js');}());

これをブックマークしてゲームを開いているタブで開けばok

javascript:(function() {Game.LoadMod('https://pointn.github.io/cookiepricecheck/cookiepricecheck.js');}());

2022-12-04

変なソートを作りたい

プログラミング競プロ

この記事は呉高専アドベントカレンダー2022に11月30日まで空きがあったら適当に埋めるための記事です。→ 4日目の記事になりました ~~/ 空きがないので放流します~~。

ごあいさつ
はじめに
既存の有名ソートアルゴリズムたち
変なソートを作ろう
- 平方分割バブルソート
  - 追記
- 最長増減部分列分割ソート
  - 最長増加部分列
実験
- ソースコード
おわりに

ごあいさつ

こちらは11日の記事の息抜きに書いていたものです。さっそく始めます。

はじめに

ソート。配列を昇順（あるいは降順）に並び替える処理です。単純ながらあらゆる場面で頻繁に使用するアルゴリズムであり、その実現方法には様々なものがあります。呉高専OBである私も、情報処理の授業で各種ソートを実装した覚えがあります。（余談：当時はC言語でしたが、シラバスを見たところ今の子たちはPythonでこれらを学ぶっぽいです。またシラバスにはグラフ理論や動的計画法（ナップサック問題）、乱択アルゴリズムなどに関する記述もあり、時代の変化を感じます）しかし、教科書に書いてあることを写してソート関数完成やったーというのは（必要かもしれませんが）あまり楽しいものではありません。自分でソートのアルゴリズムを考えて実装したほうが楽しいはずです。なのでそれを今からやりたいと思います。まずアルゴリズムの教科書によく書いてある有名なソートをいくつか振り返ります（知っている人は飛ばしてください）。次に、変なソートを作ってみて、その性能について考えてみたいと思います。途中でちょっとだけ競技プログラミング（競プロ）の話もします。

お約束

ソートしたい配列をA、その長さをNとします。0-index（AはA[0]からA[N-1]まで）です。ソートは常に昇順です。
計算量はたいてい平均時間計算量の意味です。計算量の表現にオーダ記法を用います。使い方が雑ですが目をつぶってください。知らない人は、中に書かれた関数をプロットしたときに急に増えていかないほうが嬉しい、と思ってください。
Aの連続した区間は半開区間でA[a:b)のように表します。aを含みbを含まない、A[a]からA[b-1]まで、の意味です。

既存の有名ソートアルゴリズムたち

バブルソート

Aを走査しながら、A[i]>A[i+1]となっている部分をswapします。走査回数が $N$ 回、走査する部分の長さが平均 $N/2$ なので、計算量は $O(N ^ {2})$ です。単純かつ実装難易度が低いですが、計算量は（実用に耐えうる）ソートアルゴリズムの中では悪く、優れたアルゴリズムではありません。

選択ソート

A[0:x)がソートされているとします。未ソートの範囲A[x:N)を走査します。最小の値がA[min]にあるとして、これをA[x]とswapするとA[0:x+1)がソート済みになります。これを繰り返します。計算量はバブルソートと同じ理屈で $O(N ^ {2})$ です。

クイックソート

正直この記事とあまり関係ないので紹介しなくてもいいけどこれ紹介しないのは「アンパンマンワールドのゆかいななかまたちを紹介します！」って言ってアンパンマンを紹介しないのと同じようなものじゃんという感じなので仕方なく紹介します。

A[0:N)をソートすることを考えます。適当な値をピボットとして設定し、それ未満の値をA[0:x)、それ以上の値をA[x:N)の範囲に持ってきます。A[0:x)とA[x:N)がソートされれば全体をソートしたことになるので、これで小さな2つの問題に分割されました。この2つの問題も、同じようにそれぞれ小さな2つの問題に分割して解きます。問題の最小単位は長さ2以下の配列に対するソート（単純な大小比較とswap）です。

検索すればいくらでも出てくるので理屈は他の文献に委ねますが、計算量は平均 $O(N\log N)$ 、かつ同程度に速い他の手法と比較しても優れていると言われる優秀なアルゴリズムです。しかし、最悪計算量 $O(N ^ {2})$ という問題児な一面も兼ね備えています。

マージソート

分割統治の考え方に基づくアルゴリズムです。A[0:N)をソートすることを考えます。これをA[0:N/2)をソートする問題とA[N/2:N)をソートする問題に分割します。これを解いたら、この2つの数列をマージします。これは「A[0:N/2)とA[N/2:N)に対し、両者の先頭の値を比較して小さいほうを空の配列に加えて、元の配列から除去する」という操作を繰り返すことで可能です。言葉で聞くより図を見るほうがわかりやすいと思います。再帰の深さがおよそ $\log _ {2} N$ 、マージの計算量が各再帰深さについて合計 $O(N)$ なので、計算量は $O(N\log N)$ です。

このほか、挿入ソート： $O(N ^ {2})$ やヒープソート： $O(N\log N)$ などが有名だと思います。また、C++11以降のstd::sort()はイントロソートというアルゴリズムを使っています（規格として定まってはいませんが、そのように実装されていることが多いそうです。参考: sort - cpprefjp C++日本語リファレンス）。イントロソートは、クイックソートをベースに、再帰が深くなるとヒープソートに切り替えることによって最悪計算量が $O(N\log N)$ に改善されるというものです。PythonやJavaやRustではティムソートなんてのを使っているらしいです。

変なソートを作ろう

思い付きで変なソートを作ってみたいと思います。目標は「バブルソートよりはマシ」、つまり $O(N ^ {2})$ よりもよい計算量です。

平方分割バブルソート

上で述べた通り、バブルソートは計算量が $O(N ^ {2})$ と悪いのですが、これにマージソートの要素を加えて計算量を改善してみたいと思います。

Aをいくつかのブロックに分割することを考えます。何個に分割するかは全体を俯瞰してから適切なものを決めたいので、ここではB個のブロックに分割するものとします。すると、Aは長さN/Bの配列B個に分割されます。

次に各ブロックに対し、ブロック内の範囲でバブルソートします。この操作の計算量は全体で $O((N/B) ^ {2}\times B)=O(N ^ {2} / B)$ です。これによりソートされた長さN/Bの配列がB個できあがります。

最後にこれらの配列を1つの配列にマージします。B個の配列の先頭を走査し、最小の値をソート済み配列に追加し、元の配列から削除する、という一連の操作を繰り返します。配列の要素数の総和は当然Nですから、B要素の走査がN回行われることになり、この処理の計算量は $O(NB)$ です。

全体の計算量は $O(N ^{2}/B + NB)$ となりました。 $B$ はいくらにするのが最適でしょうか。大きくすると $NB$ が、小さくすると $N ^ {2}/B$ が支配的になってしまうので、その中間の $N ^ {2}/B=NB$ になるあたりで手を打つのがよさそうです。これを変形すると $B=\sqrt{N}$ になります。よって配列を長さ $\sqrt{N}$ の配列 $\sqrt{N}$ 個に分割して上の処理をすれば、ソート部分とマージ部分の計算量がどちらも $O(N\sqrt{N})$ になり、 $O(N ^ {2})$ から改善されます（もちろん通常の $O(N \log {N})$ のマージソートには劣ります）。

このように大きさ $\sqrt{N}$ のブロックに分割して計算量を落とす手法を平方分割と呼んだりします。競プロでもよく使う手法です。

例えばこんな問題が解けます

参考にしたのはこちらです。平方分割の練習をしようにも難しい問題ばかり、そんなお悩みに狙いを決めて手取り足取りのレクチャーです！ - CADDi Tech Blog

長さNの配列Aが与えられる。以下の2種の命令を順に処理せよ。

A[i]をxに変更せよ。
A[l:r)の最大値を答えよ。

上の命令は簡単で、即座にできます。下の命令は単純に計算すると $O(N)$ です。重たい処理に引っ張られるので、命令あたりの計算量は $O(N)$ です。

ここで配列を長さ $\sqrt{N}$ のブロック $\sqrt{N}$ 個に分割します。また事前に各ブロックごとに要素の最大値を計算し、メモしておきます。こうすると、範囲[l:r)はいくつかのブロックといくつかの端数に分割されます。ブロックの数は $\sqrt{N}$ 個以下、端数は右と左で各 $\sqrt{N}$ 未満ずつなので、それらの最大値を求めるための計算量は $O(\sqrt{N})$ に高速化されます（もちろん、変更命令の際にA[i]だけでなくA[i]が属するブロックの最大値も変更する必要があります。素直にブロック内を走査すればよく、 $O(\sqrt{N})$ で計算できます）。なお、この問題はsegment treeというデータ構造（配列の一点変更と区間最大値をどちらも $O(\log N)$ で処理できます、考え方はマージソートで使ったものに似ています）などを用いてもっと効率的に解くことができます。

ところで、このソートにおいて $\sqrt{N}$ 個の配列に分割した後、それらに対してこのソートを再帰的に適用すると、計算量は改善されるでしょうか。答えは否です。これは最終的に長さ $\sqrt{N}$ の配列 $\sqrt{N}$ 個をマージするパートの計算量が $O(N\sqrt{N})$ から落ちないためです。再帰的に処理する元気があるなら普通のマージソートを使ってください。

追記

先行研究（バブルソートか選択ソートかの違いでほぼ同じもの）が見つかったので紹介します、"ソート平方分割"などのワードで検索すると他にも出てくると思います：

平方分割でマージソートするやつ（O(N√N)なので実用性は無）https://t.co/NeZtktnHlB
— 🐟 (@0x3b800001) 2021年11月16日

最長増減部分列分割ソート

ソートの話をしすぎて、ソートをサボりたい気分になりました。配列からすでにソートされている部分を抜き出して、それらをマージするというのはどうでしょうか。例として

A={1,4,7,3,5,6,2,8}

を考えます。分割の際、先頭は必ず選ぶとして、それと同じか大きいものになっていくように選んでいきます。

A={1,4,7,3,5,6,2,8}

A={3,5,6,2}

A={2}

これで3個の配列{1,4,7,8} {3,5,6} {2}に分割されたのでこれを先ほどと同じ要領でマージします。

しかし、問題があります。こんな配列はどうでしょう。

A={8,7,6,5,4,3,2,1}

これは先頭が常に最大になるので8個の数列（{8}{7}{6}{5}{4}{3}{2}{1}）に分割されてしまいます。逆順ソートされた配列に対してはその要素数ぶんだけ長さ1の配列が作られるので、マージの計算量が $O(N^ {2})$ になってしまいます。処理が複雑なぶん、バブルソートにも負けてしまいそうです。

では、ソートされている部分を抜き出すのと逆順ソートされている部分を抜き出すのを交互にやったらどうでしょうか。こうすれば、先ほどの配列は{8}{7,6,5,4,3,2,1}に分割され、逆順ソートされているものをもとに戻してマージすればいい感じになります。最初の{1,4,7,3,5,6,2,8}もできる配列は3つのままです。

本当にこれでいいのでしょうか。答えは否です。例えば、

A={8,1,7,2,6,3,5,4}

は8つの配列に分割されてしまいます。（逆順ソートを先に抜き出すと若干改善しますが、ほぼ変わりません）

では、この考えかたでは $O(N ^ {2})$ から改善することはできないのでしょうか。もう少し頑張ってみます。

なぜ先ほどのやり方ではうまくいかないのでしょうか。それは、「目先のものに後先考えず飛びついているから」です。{8,1,7,2,6,3,5,4}がいい例で、先頭の8に飛びつかず、次の1が来るのを待って、7を我慢して、2を取って、という感じで、できれば{1,2,3,4}{8,7,6,5}に分割してもらいたいのです。ようは、できるだけ長い部分列を取ってきてほしいわけですが、これを実現するアルゴリズムがあります。

最長増加部分列

上の問題は最長増加部分列（Longest Increasing Subsequence、LIS）といい、競プロでも頻出です。これを解くための計算量は $O(N\log N)$ で、{8,1,7,2,6,3,5,4}から部分列{1,2,3,4}を抜き出すことができます。アルゴリズムの中身は動的計画法（Dynamic Programming）と呼ばれるものです。配列dpを用意して、Aを走査しながら

$dp [ i ]$ = 長さ $i$ の増加列の末尾の最小値

となるように更新していきます（今回はLISを配列から抜き出す必要があるので、末尾の最小値ではなく末尾の最小値が格納されているインデックス、としています）。

LISについてもう少し詳しく

例としてA={1,6,2,5,8,3,4,7}のLISを考えます。最初、dp={inf, inf, inf, inf, inf, inf, inf, inf}です（infは無限大の意）。また、今だけ1-indexにします（そのほうが理解が容易です）。

A[1]=1に注目します。これを使うと長さ1の増加列{1}ができます（増加列というのは違和感ですが）。dp[1]:=1と更新します。dp={1, inf, ..., inf}です。
A[2]=6です。先ほどの1と合わせて長さ2の増加列{1, 6}ができます。dp[2]:=6です。dp={1, 6, inf, ..., inf}です。
A[3]=2です。残念ながら長さ3の列{1,6,2}は増加列ではありません（6→2で減ってしまう）。ですが、最初の1と合わせて長さ2の増加列{1, 2}を作ることができます。これは{1, 6}よりも嬉しいです。なぜなら、6より大きい数よりも2より大きい数のほうが多いため、増加列を伸ばしやすくなるからです。dp[2]:=2、dp = {1, 2, inf, ..., inf} です。
A[4]=5です。さっきのdp[2]=2が効いて、長さ3の増加列{1, 2, 5}が作れます。dp[3]:=5、dp = {1, 2, 5, inf, ..., inf}です。
A[5]=8です。dp[4]:=8、dp = {1, 2, 5, 8, inf, ..., inf}です。
A[6]=3です。今度の更新はdp[3]:=3になります。5を3に書き換えることで後ろに続くことができる数の範囲を広げます。dp = {1, 2, 3, 8, inf, ..., inf}です。
A[7]=4です。dp[4]:=4とできます。dp = {1, 2, 3, 4, inf, ..., inf}です。
A[8]=7です。5. の終わりではdp={1, 2, 5, 8, ...}でしたがこれが6. と7. で{1, 2, 3, 4, ...}に変わったのが効いて、dp[5]:=7とできます。dp = {1, 2, 3, 4, 7, inf, ..., inf}でフィニッシュです。dpからinfでない部分を抜き出して、LISは{1, 2, 3, 4, 7}となります。

さて、dp[?]:=A[i]という更新をたくさんやってきましたが、この?の位置はどうやって決めればよいでしょうか。よく観察すると、この位置は「A[i]より大きい数が入っている場所のうち、最も左」です。これは二分探索で効率的に求められます。

（LISについてはこちらも合わせてご覧ください： LIS でも大活躍！ DP の配列使いまわしテクニックを特集 - Qiita）

二分探索とは

二分探索は単調増加/減少なものから何かを探索するときに使うアルゴリズムです。

今回は、「dp配列のうちA[i]より大きい数値が入っている場所の左端」を求めます。dp配列が以下のように単調増加になっているので二分探索が使えます。

$dp[1]\le dp[2]\le ... \le dp[N]$

まず、配列全体のど真ん中を見ます。もしその値がA[i]より大きいなら、更新すべき場所は配列のど真ん中それ自身か、それより左側のどこかです。逆に小さければ、更新すべき場所は配列のど真ん中より右側のはずです。これだけで考えるべき範囲が一気に半分になりました。次に、半分になった探索範囲のど真ん中を見て、それよりも右か左かを判定します。すると、全体の半分だった探索範囲をさらに半分にできます。これを繰り返すことで探索範囲を毎回半分に減らしていくことができ、いずれ1か所に特定できます。繰り返し回数はおよそ $\log _ {2} N$ 回で、Nに対してかなり小さい値になります。

このように、探索範囲を半分に狭めながら行う探索が二分探索です。

これを用いて配列の最長増加部分列と最長減少部分列のうち長いほうを順次抜き出していき、最後それらをマージする、これならいい感じに計算量が減ってくれそうです。

さて、最長増加部分列または最長減少部分列に分解していくと、結局いくつの配列ができあがるのでしょうか。これについて、興味深い定理があります。それがErdös Szekersの定理です。

a,b を正の整数とする。各項が相異なる長さ ab+1 の数列があるとき，以下の二つのうち少なくとも一つは必ず存在する。

長さ a+1 の部分列で，単調増加なもの

長さ b+1 の部分列で，単調減少なもの

(Erdös Szekersの定理とその証明 | 高校数学の美しい物語, https://manabitimes.jp/math/909)

ここでは証明には触れませんが（鳩の巣原理という競プロでもよく使う原理を使います。そこまで難しくないと思うのでぜひ引用元のページを読んでみてください）、要はどんな場合も現時点の配列の長さの平方根以上の長さの単調増加or減少部分列を取ることができると解釈できます。

この性質を知ったうえで、ある長さの配列が何個ぐらいの部分列に分割されるかを実際にプログラムで試してみます。具体的には漸化式

$f(0)=0,f(n) = f(n- \lfloor \sqrt{n} \rfloor ) +1$

を $n$ を変えて計算します。すると、

$f(1000)=63, f(10000)=199, f(10 ^ {6}) = 1999, f(10 ^ {10})=199999$

と、 $f(n)\approx 2\sqrt{n}$ になりそうです。

ざっくりとした理屈を考えてみましょう。まず $\lfloor\sqrt{n}\rfloor$ が気持ち悪いので $n=x ^ {2}$ として式を書き直すと、

$f(x ^ {2})=f(x ^ {2}-x)+1=f( (x-1/2) ^ {2}-1/4)+1\approx f( (x-1/2) ^ {2})+1$

となります。 $g(x)=f(x ^ {2})$ とおくと $g(x)=g(x-1/2)+1$ となり、これは一次関数 $g(x)=2x$ です。結果的に、

$f(n)=g(x)=2x=2\sqrt{n}=O(\sqrt{n} )$

になりそうです。

配列数が $O(\sqrt{N})$ なので、マージにかかる計算量は $O(N\sqrt{N})$ です。問題はLISを求めるパートの更新計算にかかる $\log {N}$ です。二分探索を行う回数が $O(N)$ とかになってくれないかなーと思ったのですが残念ながら $O(N\sqrt{N})$ になるようなので、こちらがボトルネックで $O(N\sqrt{N}\log{N})$ っぽいです。

二分探索回数の解析

配列がいくつに分割されるか考えた時とほぼ同様の手法です。その時点での配列長 $n$ における二分探索の実行回数 $S(n)$ は、

$S(0)=0,S(n)=n+S(n-\lfloor\sqrt{n}\rfloor)$

であると考えられます。 $n=x ^ {2}$ として書き直すと

$S(x ^ {2})\approx x ^ {2}+S( (x-1/2) ^ {2}) =x ^ {2}+(x-1/2) ^ {2}+S( (x-1) ^ {2})=...=\sum _ {i=1} ^ {2x} (i/2) ^ {2}=(\sum _ {i=1} ^ {2x} i ^ {2}) /4$

で、自然数の二乗の和の公式

$\sum _ {i=1} ^ {x} i ^ {2}=x(x-1)(2x-1)/6$

より、

$S(n)\approx (16x ^ {3}/6) / 4=2x ^ {3}/3=O(x ^ {3})=O(n\sqrt{n})$ です。実際にプログラムで検証してみると、これに近い結果が得られます。

実験

手元のノートパソコンで、要素数N={10000, 50000, 100000, 500000}、0以上 2³⁰ 未満のランダムな値の配列各100個に対して各種ソートを実行し、実行時間を計測してみます。せっかく計算量を評価したので、変なソートたちがバブルソートより速いこと（そして標準のソートより数段遅いこと）、あとできれば変なソート達の計算量の差が確認できると嬉しいです。

というわけでプログラムを書きました。200行近いソースコードをここに貼ると大変なことになるのでもう少し下にリンクを貼っています。

下表に平均処理時間[ms]を示します。xと書いてある部分のデータは膨大な時間がかかるのでとっていません。処理時間がきわめて短いものは範囲で示しています。

N	標準ソート	バブルソート	平方分割	部分列分解
1e4	0~1	118	1~4	12
5e4	1~4	3915	27	158
1e5	4~8	x	74	463
5e5	34	x	740	5732

この表から分かる通り、C++標準のソートは速いしバブルソートは使い物になりません。平方分割によってバブルソートは格段に速くなりますが標準のソートよりだいぶ遅いです。部分列分解は平方分割よりさらに数倍遅く、 $\log$ や実装の重さが響いている感じでしょうか。

ソースコード

リンクはこちら（N=100や1000のデータもとるようになっているのですが、軽すぎて無意味でした）。自然な実装にしてあるので、「おいおい何やこの無駄な処理は」と思う方は我慢してください。Intel Core i7-7500U(3.5GHz)、gcc8.1.0、C++17、最適化オプションO2です。

strangesort.cpp - Google ドライブ

おわりに

使い物にならないソートを作り、その性能を確かめました。皆さんも面白いソートを作ってみてはいかがでしょうか。

2022-08-31

ひとりで行きにくい場所にふたりで行った話メイド喫茶とスイーツバイキング編

日記

男は思った。「メイド喫茶でメイドさんを眺めながらお茶を楽しみたい」と。

男は思った。「スイーツバイキングで甘いものをお腹いっぱい食べたい」と。

男たちは思った。「だが、ひとりでは行きにくい」と。

人
計画
集合
ワンダーパーラー
スイーツパラダイス
解散
その後

人

点N：大学生。スイーツバイキングに行きたいと思っているが、ひとりでは行きにくい。僕。

グレタ：社会人。点Nの高専時代の同期。メイド喫茶に行きたいと思っているが、ひとりでは行きにくい。いやこいつは僕と違ってその気になればひとりで行ってたと思うけど

計画

しょっぱなから申し訳ないのだが、ここだけ書くのをさぼっていたので全然思い出せない（3か月以上前の話なので）。お互いこういう場所に行きたいという話をしていて、何かの拍子にじゃあふたりで行こうということになった。お互いがお互いの心理的な盾になれるし、僕もグレタも相手の行きたい場所に多少の興味があった。実行に関してはグレタのほうが積極的だったような気がする。まあ、スイパラよりもメイド喫茶のほうが「人生で一度は行ってみたい度」は高いだろう。

多忙なグレタと暇な僕が相談し、最終的に立てたプランはこちら。ワンダーパーラーにはグレタが30分ほど心の準備をしたのち電話で予約を入れ、スイパラには僕がネットでサクッと予約を入れた。予約状況の関係もあったが、スイパラはあえてアウェイ中のアウェイと思われる原宿店を選択した。

5月28日（土）

12:30 池袋駅29番出口

13:00 池袋ワンダーパーラー

16:00 スイーツパラダイスSoLaDo原宿店

集合

遅寝遅起き生活を送っていた僕にとって12時半とはすなわち普段の起床時刻を少し過ぎたぐらいの時間帯である。つくばから東京までは何をどうやっても1時間以上かかるので、この日ばかりはきちんと目覚ましをかけて起きなければならない。5時からアラームの波状攻撃を仕掛け、無事6時に起きることに成功した。起きたらとりあえず昨夜サボったお風呂に入って洗濯機を回しつつ朝ご飯のロールパンを野菜ジュースで流し、洗濯が終わりしだい干して家を出る。グレタと前回会った時お土産にヘルエスタ名物うなぎパイを貰ったので何か買っていこうと思い、つくば駅で適当なお土産（芋のなんか、安めのやつ）を買って、帰りのTXの格安券を買ってバス乗り場へ。この格安券を使うと通常1210円が970円になるので買わないわけにはいかない。

高速バスの出発は10時半。1時間と少しで東京駅に着く。土曜ということもあっていつもより座席は混雑していた。僕の隣にはおばあちゃんが座り、ペンを片手にナンプレを解いていた。バスに揺られながら、寝て、起きた。おばあちゃんのナンプレはそこまで進んでいないようだった。どれぐらい難しい問題を解いていたのかはわからないが、こういうのはじっくりゆっくり楽しめる難易度がちょうどいいのだろう。

東京駅の日本橋口に降り立ったら、地下鉄丸ノ内線に乗り換える。山手線しか知らなかったころの僕とは違う。こっちのほうが速い。でも丸ノ内線の乗り場までが遠いんだよな、どっちがよかったんだろう。

地下鉄に揺られること十数分、池袋駅に到着したのだがここからが長かった。29番出口を探して行ったり来たり出たり入ったり、スマホを見て駅の案内を見て、29番出口に到着したころにはほぼ定刻。駅本体から道路を挟んで向かい側だとは思わなかった。グレタと合流し、今日このあとの予定に対する期待と不安を共有した。

ワンダーパーラー

まずは先攻・グレタのターンである。池袋駅から歩いてだいたい10分、帝京平成大学のすぐ近くの路地を少し入ったところに目的の喫茶店はある。メイド喫茶というと真っ先に思い浮かぶのはピンクの衣装に身を包んだメイドさんがオムライスにおいしくなる魔法をかけてくれるやつだが、こちらのお店はそうじゃないほうのメイド喫茶。落ち着いた雰囲気の空間でお茶を楽しむ場所である。

おそらくグレタは移動に不慣れな僕が10分程度遅刻することも想定したスケジュールを立ててくれていたのだが、僕が無事定刻に到着したのもあり、まっすぐ行くと暇を持て余すことになりそうだったので、近況の話をしながらグレタの勘を頼りになんとなくで歩いてみることにした。だいたいこの辺だろうというところまで来たが店までは辿り着けず、万が一あさっての方向に進んでいたら取返しがつかなくなりそうになったところでスマホを確認した。結果かなりいいところまで来ており、右に曲がるべきところを左に曲がってしまっただけということが分かった。後ろを向くと電柱にワンダーパーラーの広告が。ふたりで惜しかったねと笑いあった。

思ったより惜しい位置にいたことが判明したためお店の前を通り過ぎて一旦撤退し、帝京平成大学を経由して公園でひと休み。グレタは大学の写真を誰かに送り付けていた。帝京平成大学出身の知り合いがいるらしい。公園にはなんか典型地雷女子みたいな人がいた。

予約の時刻が迫ったのでお店に向かう。外観を撮ろうとする僕にグレタは最初、「後からでもいいだろ」的なことを言っていたのだが、急に何かを思い出したらしく、手のひらを返して「今撮っておけ」と言い始めた。なぜかグレタは理由を言わなかったが、ここがメイド喫茶であることを考えれば、きっとメイドさんがお見送りしてくれるから撮りにくくなるんだろうなと思った。

着いた、お店の写真も撮った、あとは入るだけである。グレタがドアを開け、僕は後に続いた。緊張していたせいか席に着くまでの記憶があまりないのだが、予約していたグレタです、お待ちしておりましたこちらへ、みたいな会話がグレタとメイドさんの間で行われていたような気がする。椅子には手触りのいい良さげな布が敷いてあった。

当店のご来店は初めてなのでメイドさんから懇切丁寧なご説明を賜り、グレタと顔を合わせる。正直僕はかなり場の空気に圧倒されてしまっていたのだが、とにかく自分が注文するメニューを考えることにした。

僕にとってはこのあとに控えるスイパラがメインの用事であり、スイパラとはすなわちバイキングであるので、ここでお腹いっぱいになってしまうのはよくない。ということで紅茶とデザートのアフタヌーンセットにした。デザートは何種類か用意されていて、店内のボードに書いてあるものから選ぶ形。あまおう苺のモンブランに決めた。紅茶はもっと種類が多いので、紅茶のあれこれがよくわからない僕はメニューに書いてあるチャートを参考に決めた。人気の品種とかも書いてあるので何もわからん陰キャでもなんとかなると思う。ちなみに、これを書いている今の時点でもう何を頼んだか忘れてしまった。

一方のグレタはワンダーパーラー目的で参加しているので、悩む僕と対照的に事前対策ばっちりという感じだった。僕が注文を決めたことを確認すると、グレタはベルを鳴らした。あとで聞いたのだが、どれぐらい鳴らせばいいのかわからなくて不安だったらしい。僕もわからなかったので、鳴らしてもらって少し安心していた。

メイドさんがやってきた。グレタの注文はフードセットメニューで、料理は鶏のクリームシチュー、紅茶はカーニバル、デザートはごまと豆乳のモンブラン。鶏のクリームシチューはもとは限定のメニューだったようだが人気があって常設メニューになったらしく、グレタはこれを食べに来たとのこと。またカーニバルはお米さんのメイドティーで、たぶんこれもグレタのお目当て。お米さんのツイートめっちゃリツイートしてたし。メイドティーはそのメイドさんがいるときしか提供されないのだが、シフトは基本的に公表されていないので、運が良かったようだ。

僕の注文が先に運ばれてきた。一杯目の紅茶はメイドさんがカップに注いでくれる。ポットが冷めないためのカバーをして、熱いのでお気を付けくださいとの忠告をし、レモンと砂糖を置いてメイドさんは戻っていった。何回も見てるけど、戻っていくときにひざを曲げる（カーテシー、という挨拶らしい）のがかわいい。置き場所を定めるために利き手を聞かれることにグレタは感心していた。

グレタなんか待っていたらせっかくメイドさんがいれてくれたお茶が冷めてしまうので、先にいただくことにした。カップの下の皿って持つんだろうかとか余計なことを考えつつ、まったく手馴れていない動作でひと口。おいしい。紅茶のことなんも知らんけど。あと皿はべつに持たなくていいっぽいとグレタが教えてくれたのでもう持たないことにした。デザートのケーキもおいしいし、お皿にチョコで描かれたイラストがかわいい。このイラストはメイドさんのサービスである。

少ししてグレタのクリームシチューが到着したのだが、グレタの様子がおかしい。話を聞いてみるとスプーンがないらしい。少し待ってみても運ばれてくる様子がないのでおそらくお店側のミスだろうと結論付け、通りかかったメイドさんにグレタがスプーンが欲しいと伝えたところ（ベル鳴らせばよかったと後から悔やんでいた）、すぐに持ってきてくれた。なんかそういう変化球ぎみのサービスの一環というよりはシンプルミスなんだろうと思うが、そういうのも含めて楽しむのがいいんだと思う。

僕は二杯目の紅茶を注ぎ（メイドさんの言う通りポットが本当に熱かったので皆様もお気を付けください）、お皿のチョコ絵もケーキも崩さぬよう集中しながら、ゆっくり食べ進めていった。

クリームシチューを食べ終えたグレタにも紅茶とデザートが届き、だんだん場の空気に慣れてきた僕はグレタと互いの近況を話しあった。引きこもり大学生の近況なんて面白いものではないが、それぐらいしか話題がないのでやむを得ない。グレタが紅茶に蜂蜜を入れるときの目つきがやけに真剣なのでそう言ったら、仕事の癖が出てしまっているのかもしれないと言っていた。どうやら彼は機械に慎重に油を注ぐ仕事をしているらしい。

なんやかんやで滞在可能時間の90分近くが経過したのでお店を出ることにした。お金はグレタが払ってくれた。ポイントカードを（せっかくなので僕のぶんも）作ってもらい、メイドさんのブロマイドを受け取った。

入店前の予想通りメイドさんにお見送りしてもらってお店をあとにした。グレタによると我々が見えなくなるまでお辞儀していたらしい。ちゃんとしてる。

この日のためにグレタは入念な下調べをしていたのだが、ブロマイドをもらうことを忘れていたらしく、持ち運ぶうちにへしゃげてしまわないよう、近くの百均に行ってそういうのを保存するための透明で固くて薄いやつを買った。百均が地下にあるのが都会だなぁと思った。道路には知らないアイドルが中で歌い踊るバスが走っていた。

スイーツパラダイス

後攻・点Nのターン。スイーツパラダイス、通称スイパラ。都内と各地方の主要都市に分布する、スイーツ食べ放題のお店。女性客が圧倒的に多く（事実）、成人男性が単独で入るには強い心が必要である（個人の見解）。

ワンダーパーラーをニコニコで退店してから30分、時刻は15時。ここから1時間後にスイーツバイキングというのは若干予定を詰めすぎてしまったような気もしたが、予約したからには向かうしかない。

というわけで、まずは池袋から原宿への移動である（池袋にもスイパラはあるのだが、2日前の時点で予約がいっぱいだったので避けた）。集合時、1時間半あるし歩くのもいいかもねなどと言っていたが、池袋から原宿は思いのほか遠く、シンプルに不可能だったのでおとなしく山手線に乗ることにした。道中池袋のスイパラの前を通ったのだが、案の定列ができていた。

休日午後のやや混みぐらいの山手線に乗って数分、原宿駅に到着した。目指すスイパラSoLaDo原宿店は竹下通り沿いに歩いて行くだけで着く、のだが。

レインボーカラーの綿菓子や、ちょっとした武器になりそうな長さのトルネードポテトなど、食べやすさを犠牲に"映え"に全振りしたような食べ物たちと、それに負けない派手な人たちが視界に入った。

やがて右手にスイパラのあるSoLaDoが見えてきた。予約の時間まで少し時間があったので通りを抜けたところまでいったん避難して、ふたりで竹下通りを歩いてみた感想を話し合った。内容は到底ここには書けない。

5分前になったのでスイパラに向かった。奥の階段から上に行く必要があるのだが、建物の構造がよくわからなかったので1階を突っ切ることにした。どこを見てもガールズファッションの店で、どう見ても用のない陰気な男ふたりが高専の真逆みたいな男女比の空間を歩くのはかなり気まずかったので、急いでフロアを突っ切った。

外の階段で3階に上るとそこはいよいよスイパラ。入り口から続く行列が階段まで伸びており、列を構成するのはやはり若い女の子であった。足元のこれ（下画像）を見て、一瞬、スイパラって男性のみのご利用不可だったっけと思ってしまった。そんなことはないし、これはスイパラ横にあるトイレの待機列表示である。

幸運なことに、並び始めてすぐ我々の後ろにやや年配の夫婦が並んでくれたため、非常に心強かった。自分たちはここにいていいのだと、心から思えた瞬間であった。

勝手がわからないので予約していたが、特にそれで優先的に入れるというようなこともなく、普通に順番通り案内された。おすすめバイキングのチケットを2枚購入したのち、およそ10分遅れで店内に入り、店員さんからの簡単な説明を右から左へ受け流すと、いよいよ80分間のバイキングが始まった。開始時刻は16時15分だったので、終了は17時35分。

記念すべき1皿目。16時24分。仕方のないことだが、取りに行くときに手袋しなきゃならないのがちょっとめんどくさい。とりあえず王道っぽいケーキを中心に取った。ムースの類は酸味があって、甘いばっかりにならないのがちょっと助かった。

ケーキ右からふわふわショート・チョコレートケーキ・ミルクレープ。小さいのは手前からティラミス・プチ抹茶ムース・プチベイクドチーズ・プチブルーベリームース

ドリンクはアイスコーヒーを選択した。甘いものはどちらかというと得意だと思っているが、さすがに甘味を打ち消さないとケーキバイキングはやってられない。

ちなみにこのミルクレープ、ホールの最後の1カットだったのだが、以後補充された様子がなく、結果として最後の1個だったらしい。にじさんじのB級バラエティ（仮）第23回によればミルクレープは手間のかかるケーキだそうで、大量に生産するのが難しいとのこと。食べることができてよかった。

www.youtube.com

バイキングが始まってしまえばこちらのもので、皿の上のものを食べていれば周囲の客は気にならない。飛ばし気味だったのでグレタに若干引かれつつ、1皿目を完食。

続けて2皿目。16時41分。まだまだ余裕だったので多めに取ってきた。めっちゃ食うやんみたいなことをグレタが言ってきた気がする。めっちゃ食ったほうがいいだろバイキングなんだから。

スイパラのケーキは一切れがかなり小さいため、たくさんの種類が選べるのがよい。その半面、取りにくいし皿に立てにくいから盛り付けが大変なことになるのだが。

3皿目。16時57分。さすがに甘さと油できつくなる。食べたいものは2皿目までに取っておくべきである。

そういえば、スイパラは通常のバイキングのほか、事あるごとに何かしらとコラボして限定のメニューを販売しており、それを目的に来店する客も多い。僕たちの周囲でもそれっぽいものを食べている人がいて、そういう人たちどうしの交流のようなものを少しだけ垣間見ることができた。

4皿目、17時11分。甘いものがきつくなったのでスパイシースイパラカレーでいったん休憩。今求めているのはルーの辛味であって米に用はないのだが、いちおう気持ち程度に米も食べた。ルーだけ食べるのはなんというか、RTA走者みたいな感じになっちゃうし。

5皿目、17時24分。4皿目の時点でかなり満足だったのだが、もう少しだけ食べておきたくなってしまった。

カレーを挟んでお腹も復活したかに思えたが、体はしっかり甘いものを受け付けなくなっていた。正直ひと口食べた時点で取ってきたことを後悔したが、コーヒーで脳を騙しながらどうにか時間内に食べ終わった。せっかくアイスがあるんだからアイスにすればよかった。

初めてのスイパラだったが、僕としてはおおむね満足であった。一部の食べ物の補充が間に合っていなかったり、洗い方が雑に思える皿があったり（女の子多い店でそれはちょっとどうでしょうという気持ちにはなった。ご時世的な心配もある）と、店を回し切れていない様子は見られたが、まあいいとしよう。食べ放題時間の80分は思ったより短かったが、時間を持て余すでもなく食べ足りないでもない、かなりいい感じの時間設定に思えた。

ちなみにグレタは昼のクリームシチューが効いていたのか、最初にケーキを一皿食べたのと途中でカレーを一皿食べたぐらい。さすがにメイド喫茶のそれのほうが美味しいと思うが、ここでも紅茶を飲んでいた。後で聞いた話だが、グレタは別に甘いものが得意なわけでもないらしい。なんで来たんだ。

解散

YouTuberを名乗る調子こいた怪しい若者に声をかけられつつも完全に無視してなんとか竹下通りを脱出し、（主に僕ひとりが）かなり苦しかったので日比谷公園でひと休み。公園内には外国人が沢山いてワイワイしていた。上空を飛行機がたくさん飛んでいた。

日比谷公園には長いベンチがあり、競技プログラマーなら誰もが知っている（解けるとは言ってない）あの問題を思い出さずにはいられない。腹の中のケーキに苦しみながら、グレタに一方的に競プロの話をしつつ園内を軽くぐるっと回った。公園の出口のところに駆け出しの大学生シンガーがいた。ミュージシャンとして大成するかは知らないけど、そのメンタルがあるならこの先どうなろうときっとお前ならなんとかなる、などとやや悪口寄りのコメントをふたりで言いながら原宿駅に向かった。道中、ボンゴとコンガってごっちゃになるよねという話をした。大きいほうがコンガ。

池袋駅のコインロッカーに荷物を預けたグレタについて戻り、構内で迷いまくりながらもどうにか荷物を取り出し、くたびれた足で少しパルコをうろついてから解散した。グレタは何かに感化されたのか、どこかでティーポットと紅茶の茶葉を買ったらしい。なお、翌日さらに同期二名を加えて遊びに行ったりもしたのだが、ここに書くことでもないので割愛。

その後

あれ以来グレタとは会っていないのだが、ツイッター見る限りあいつめっちゃワンダーパーラー行ってる。働いてない僕が言うことでもないのだが、破産だけはしないでほしい。

2022-08-06

南京錠マインスイーパソルバを書く

ゲームプログラミング

ストーリー

「このネジを外せば、この部分がとれるのかな？」

フタを外しました。

1*	2	1
2	3	1*
1*	2	1

「なんだろうね……この数字は？」

フタを外しました。

0	1	1*
1	2	1
1*	1	0

「あっ、また数字の羅列！」

「この数字が穴の位置を示している、ということかな？」

フタを外しました。

2	3	2	2	2	4	3	2	1	1	2	2	3	2	2	2	2	1	2	2	3	3	3	3	1
2	4	3	3	3	5	4	2	2	2	3	2	4	3	3	3	4	4	5	5	6	6	5	4	1
2	4	3	3	3	4	4	2	4	3	5	4	5	4	2	4	4	5	5	6	8	7	6	5	3
1	4	4	5	5	4	4	1	3	2	3	2	4	4	3	3	4	6	7	7	7	6	5	5	3
2	4	3	4	5	4	5	2	4	2	3	2	4	4	4	4	5	6	7	6	6	4	4	4	3
3	4	3	5	6	4	4	2	3	1	1	1	3	3	4	4	6	7	7	5	4	4	4	4	2
3	4	2	3	3	2	3	3	4	2	1	2	3	3	4	5	7	7	6	5	5	6	5	4	2
3	4	3	4	3	2	3	4	3	1	0	3	4	5	4	5	5	5	3	4	4	7	6	7	4
3	5	4	4	2	2	3	5	4	3	2	3	5	5	5	5	6	6	4	4	3	5	5	6	4
3	4	3	3	2	2	4	5	4	2	3	4	6	5	5	5	6	6	5	4	2	2	3	5	4
3	5	3	4	2	3	3	5	5	5	6	5	5	2	4	4	7	6	6	3	2	1	2	2	2
3	5	3	4	2	3	3	5	4	4	5	5	4	1	3	3	6	5	6	3	3	2	2	1	1
3	5	3	4	2	4	4	6	4	5	5	5	4	2	3	2	5	4	5	2	3	3	4	3	2
4	5	3	3	2	3	3	4	2	3	3	4	3	2	2	1	4	4	5	2	2	2	3	4	3
3	3	1	2	3	5	4	3	2	4	5	5	4	3	2	0	2	3	4	3	3	3	3	4	3
3	3	1	2	4	5	3	2	2	5	5	4	2	1	1	0	2	3	3	3	4	4	2	2	2
2	3	1	2	4	6	5	3	3	6	7	5	4	3	3	1	2	4	5	5	5	5	4	4	3
3	4	2	2	3	5	4	3	3	6	6	5	3	3	3	3	5	6	6	5	5	6	6	7	5
2	3	3	3	3	4	4	4	4	5	5	4	3	3	3	3	5	7	8	6	5	5	6	7	5
2	3	5	5	4	3	2	4	4	4	3	4	3	3	3	5	6	6	5	4	5	6	7	7	5
3	4	6	5	4	1	1	3	3	2	1	3	4	5	4	5	5	5	4	3	4	5	6	6	4
4	6	8	6	5	2	1	2	2	3	2	4	4	6	6	8	7	5	2	2	3	5	5	6	4
4	5	6	5	6	3	2	2	2	4	2	4	3	6	5	6	6	6	4	4	4	5	4	5	4
2	3	4	5	7	5	3	2	2	5	3	4	1	3	3	5	5	5	4	5	5	4	2	3	3
1	1	1	3	5	4	2	1	1	3	2	3	1	2	1	2	2	3	3	4	4	3	2	2	2

「これで３つ目だけど……やっぱり内装は同じなのね」

「同じ……？」

「鍵が開いてた錠を元にして、法則を解くしかなさそうよ」

「これじゃ法則がわかっても不可能だろ……」

南京錠マインスイーパソルバを書く。法則はもう知っているので、ざっくり、次のような問題が解きたいということになる。

問題文1

$H$ 行 $W$ 列の二次元配列 $A$ が与えられます。

以下の条件を満たす $H$ 行の文字列 $S_{i} (1\le i\le H)$ を出力してください。

すべてのテストケースについて、解はただひとつ存在することが保証されます。

条件

$S_{i}$ は. と#からなる $W$ 文字の文字列である
$S_{ij} (1\le i\le H, 1\le j\le W)$ およびその8近傍に#はちょうど $A_{ij}$ 個存在する

制約

$1\le H,W\le 25$

入力例

出力例

#.#
.#.
..#

解法を考える

bit全探索は $O(HW2^{HW})$ で無理。深さ優先探索も厳しそう。

よく考えると、連立一次方程式が見えてくる。

簡単のために $3\times 3$ の場合を考える。 $A_{ij}$ は結局 $S_{ij}$ とその8近傍の#の和なので、

a	b	c
d	e	f
g	h	i

と記号をふってやると、

$a+b+d+e=A_{11}$

$a+b+c+d+e+f=A_{12}$

...

$e+f+h+i=A_{33}$

みたいな感じで式が9個できるので、これを解けばいい。一般に式の数と変数の数は $HW$ 個。この連立方程式は $O((HW)^{3})$ で解けるので、 $H$ や $W$ が25ぐらいなら現実的な時間で動く。係数行列を正しく構築してやれば、あとはPythonかなんかにやらせればいい。

問題文2

正整数 $H,W$ が与えられます。上の問題1において、 $H$ 行 $W$ 列の $A$ のうち矛盾のないもの全てに対して解が一意に定まるならYes、そうでないならNoと出力してください。

制約

$1\le H,W \le 10^{18}$

解が一意に定まる条件

係数行列 $A$ （入力の $A$ ではない）と $S_{ij}$ を#にするかの $\{0,1\}$ ベクトル $x$ 、あと入力を1列に並べ直したベクトル $b$ でもって $Ax=b$ として、 $x=A^{-1}b$ で解が求まる（最初はまじめにLU分解でも書こうかと思ったけどめんどくてPythonにやらせることにしたので細かい実装は省略）予定なのだが、気になるのは「どんな盤面サイズに対しても係数行列は逆行列をもつか、解は一意に定まるか」ということである（上の問題文2は要はこれの判定問題である）。ソルバは書いただけではだめでテストしなきゃならないので、こういうことを考えておこうという趣旨（もちろん行列式を求めて確認すればいいが、その辺の理論をちょっとぐらい考えたほうがいいだろう）。

たとえば $2\times 2$ の盤面ではどうだろうか。

2 2
1 1
1 1

1か所だけ#にすればどこでも同じになるので、答えは一意に定まらない。ほかにも

2 5
2 2 1 0 0  
2 2 1 0 0

に対する答えは以下のどっちでも成立してしまう。

#....  .#...
.#...  #....

はっきりとした答えがだせなかったので小さめの係数行列を構築して行列式を求め、それが非零であるかをチェックしてみる。

import numpy as np

table = [[0] * 10 for i in range(10)]

for h in range(1,11):
    for w in range(1,11):
        A = [[0]*(w*h) for i in range(w*h)];
        for a in range(w*h):
            ay = a//w
            ax = a%w
            for b in range(w*h):
                by = b//w
                bx = b%w
                if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
                    A[a][b] += 1
        A = np.array(A)
        d = np.linalg.det(A)
        table[h-1][w-1] = d
        print('%2d x %2d det A = %d' % (h,w,d))

print('h/w ', end='')
for i in range(1,11):
    print('%2d ' % i, end='')
print('')
for h in range(10):
    print(' %2d ' % (h+1), end='')
    for w in range(10):
        print('%2d ' % table[h][w], end='')
    print('')

必要な部分だけ出力を抽出するとこんな感じ。行列式は $\pm 1$ か $0$ で、 $H\equiv 2(mod\ 3)$ または $W\equiv 2(mod\ 3)$ のいずれかが成立していれば $0$ になるっぽい。そうじゃない場合は（ある程度パターンは見えるが）いまいち規則性が分からない。

h/w  1  2  3  4  5  6  7  8  9 10
  1  1  0 -1 -1  0  1  1  0 -1 -1
  2  0  0  0  0  0  0  0  0  0  0
  3 -1  0  1 -1  0  1 -1  0  1 -1
  4 -1  0 -1  1  0  1 -1  0 -1  1
  5  0  0  0  0  0  0  0  0  0  0
  6  1  0  1  1  0  1  1  0  1  1
  7  1  0 -1 -1  0  1  1  0 -1 -1
  8  0  0  0  0  0  0  0  0  0  0
  9 -1  0  1 -1  0  1 -1  0  1 -1
 10 -1  0 -1  1  0  1 -1  0 -1  1

とりあえず $H\equiv 2(mod\ 3)$ と $W\equiv 2(mod\ 3)$ の場合を処理したい。どっちでも一緒なので紙面の都合を考えて $W$ が該当した場合で進める。 $A_{1}=\{1,1,...,1\},A_{2}=\{1,1,...,1\}$ の場合、

$S_{1}=$ #.[.#.][.#.]...[.#.]、 $S_{2},S_{3},...=$ ... ... ...
$S_{1}=$ [.#.][.#.]...[.#.].#、 $S_{2},S_{3},...=$ ... ... ...

の2パターンが考えられる。[]の中が3文字、外に2文字で、[]の数を変えることで、3で割った余りが2であるようないかなる $W$ に対しても解が一意に定まらない $A$ が存在するとわかる。

ということで、残る問題は「 $H mod\ 3 \neq 2$ かつ $W mod\ 3\neq 2$ なら係数行列の行列式は必ず非零か」となる。

実際に係数行列がどうなっているのか見ないことには始まらないので見てみる。

import numpy as np

h,w = map(int,input().split())

A = [[0]*(w*h) for i in range(w*h)];
for a in range(w*h):
    ay = a//w
    ax = a%w
    for b in range(w*h):
        by = b//w
        bx = b%w
        if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
            A[a][b] += 1
for y in range(h*w):
    print(A[y])

A = np.array(A)
d = np.linalg.det(A)
print('det A = ',d)

$H=3,W=4$ のときの係数行列はこんな感じ。

[1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0]
[1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0]
[0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0]
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0]
[1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0]
[1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0]
[0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1]
[0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1]
[0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0]
[0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0]
[0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1]

$H=4,W=3$ はこんな感じ。

[1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0]
[1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0]
[0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0]
[1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0]
[1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]
[0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0]
[0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0]
[0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1]
[0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0]
[0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1]

目が疲れそうだが、要は「サイズWの三重対角行列の、サイズHの三重対角行列（特に非零成分はすべて1）」のような感じになっている。主対角線とその上下の成分がすべて1でサイズ $n$ の三重対角行列を $T_{n}$ と書くことにすれば、

$H=3,W=4$ の係数行列は

$\displaystyle \begin{pmatrix} T_{4}&T_{4}&0 \\ T_{4}&T_{4}&T_{4} \\ 0&T_{4}&T_{4} \\ \end{pmatrix}$

$H=4,W=3$ の係数行列は

$\displaystyle \begin{pmatrix} T_{3}&T_{3}&0&0 \\ T_{3}&T_{3}&T_{3}&0 \\ 0&T_{3}&T_{3}&T_{3} \\ 0&0&T_{3}&T_{3} \\ \end{pmatrix}$

といった感じである。

さて、三重対角行列の行列式は簡単な漸化式で求めることができる。英語版Wikipediaにまとめられているので参考にしつつ、今回は特殊なケースゆえもっと簡単になって、結局 $|T_{n}|=|T_{n-1}|-|T_{n-2}|, (|T_{0}|=1,|T_{-1}|=0)$ となる。 $|T_{n}|$ は周期6で循環し、以下のようになる。

$n$	1	2	3	4	5	6
$\|Tn\|$	1	0	-1	-1	0	1

で、今問題なのはこれの三重対角行列の行列式が欲しいというところだが、こちらもほぼ同様の性質が現れる。

$T_{n}$ の三重対角行列 $U_{m}$ の行列式は、実際に計算してみると

$|U_{1}|=|T_{n}|$

$|U_{2}|=0$

$|U_{3}|=(-1)^{n}|T_{n}|^{3}$

$|U_{4}|=(-1)^{n}|T_{n}|^{4}$

$|U_{5}|=0$

$|U_{6}|=|T_{n}|^{6}$

以下想像の通りである（と願いたい）。結局先ほどプログラムを書いて確認した行列式の値はこれに基づいていたということになる。さらに、ここから $H mod\ 3 \neq 2$ かつ $W mod\ 3 \neq 2$ なら係数行列に逆行列が存在することがわか（った気にな）る。証明でも何でもない。

ソルバ実装

あとは心ゆくままにソルバと問題ジェネレータを作るだけ。というわけで完成。三重対角行列の性質がいい感じにしてくれるらしく逆行列の成分もどうやら-1か0か1のようで、精度に注意する必要はまずないはず（？？？）。あまり変な入力をされても壊れちゃうだけなので乱暴はしないでもらえるとありがたい。効率面での改善の余地に関しては、正直その道の人が見ればオーダーがひとつふたつ下がるとか、それぐらいはありそうな気もする。逆行列はどうも高速に計算できそう。

ソルバ。最大 $25\times 25$ まで対応できる。それなりな環境で実行すれば $25\times 25$ でも1秒いらないぐらいだと思う。解が複数ある場合にその一例を示す、はさすがに諦めた。

import numpy as np
import sys

print('')
print('input problem size (if HxW -> H W)')
print('')
print('---matrix size---')
H,W = map(int,input().split())

if H%3==2 or W%3==2:
    print('cannot calc')
    sys.exit()

N = H*W

if N>25*25:
    print('too big!')
    sys.exit()

print('')
print('input problem matrix like this')
print('a b c')
print('d e f')
print('g h i')
print('')
print('---matrix---')
B = [list(map(int,input().split())) for i in range(H)]

vb = []
for y in range(H):
    for x in range(W):
        vb.append(B[y][x])
vb = np.array(vb).reshape(N,1)

A = [[0]*N for i in range(N)];
for a in range(N):
    ay = a//W
    ax = a%W
    for b in range(N):
        by = b//W
        bx = b%W
        if abs(ay-by) <= 1 and abs(ax-bx) <= 1:
            A[a][b] += 1

d = np.linalg.det(A)
if -0.5 < d < 0.5:
    print('cannot calc')
    sys.exit()

Ainv = np.linalg.inv(A)

X = np.matmul(Ainv,vb).reshape(H,W)
ans = [''] * H
Bcheck = [[0] * W for i in range(H)]
for y in range(H):
    for x in range(W):
        if X[y][x] > 0.5:
            ans[y] += '#'
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        Bcheck[y+dy][x+dx] += 1
        else:
            ans[y] += '.'

if B!=Bcheck:
    print('solution not found')
    sys.exit()

print('')
print('---answer---')
for y in range(H):
    print(ans[y])

問題ジェネレータ。ネジ存在確率指定、ネジ個数指定、手動ネジ座標設定の3つのモードで簡単問題生成。ほんとは機械的に生成した盤面に人の手を加えられたりするといい気がするけど、それはやりたい人がやって。

import numpy as np
import sys

print('')
print('input problem size (if HxW -> H W)')
print('')
print('---matrix size---')
H,W = map(int,input().split())

if H%3==2 or W%3==2:
    print('cannot calc')
    sys.exit()

N = H*W

if N>25*25:
    print('too big!')
    sys.exit()

print('')
print('select mode')
print('1:probability')
print('2:mine number')
print('3:handmade')

mode = int(input())

A = [[0]*W for i in range(H)]
used  = [[0]*W for i in range(H)]

print('')
if mode==1:
    print('input probability, range:[0.0 ~ 1.0]')
    print('probability:',end='')
    p = float(input())
    arr = np.random.rand(H,W)
    for y in range(H):
        for x in range(W):
            if arr[y][x] < p:
                used[y][x] = 1
                for dy in range(-1,2):
                    for dx in range(-1,2):
                        if 0 <= y+dy < H and 0 <= x+dx < W:
                            A[y+dy][x+dx] += 1
elif mode==2:
    print('input number of mines, range:','[%d ~ %d]' % (0,H*W))
    print('number:',end='')
    n = int(input())
    n = min(H*W,max(0,n))
    while n>0:
        y = np.random.randint(H)
        x = np.random.randint(W)
        if used[y][x]==0:
            used[y][x] = 1
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        A[y+dy][x+dx] += 1
            n -= 1
elif mode==3:
    print('input position, format:y x range:y...[0 ~ %d], x...[0 ~ %d]' % (H-1,W-1))
    print('you can set mines as many as you like by entering multiple lines')
    print('to finish setting, input \'-1 -1\'')
    print('---your input---')
    while True:
        y,x = map(int,input().split())
        if y<0 and x<0:
            break
        y = min(y,H-1)
        x = min(x,W-1)
        if used[y][x]==0:
            used[y][x] = 1
            for dy in range(-1,2):
                for dx in range(-1,2):
                    if 0 <= y+dy < H and 0 <= x+dx < W:
                        A[y+dy][x+dx] += 1
else:
    print('bad input')
    sys.exit()

print('---mines---')
ans = [''] * H
for y in range(H):
    for x in range(W):
        if used[y][x]==1:
            ans[y] += '#'
        else:
            ans[y] += '.'

for y in range(H):
    print(ans[y])

print('---problem---')
print(H,W)
for y in range(H):
    print(*A[y])

問題例（冒頭の問題）

25 25
2 3 2 2 2 4 3 2 1 1 2 2 3 2 2 2 2 1 2 2 3 3 3 3 1
2 4 3 3 3 5 4 2 2 2 3 2 4 3 3 3 4 4 5 5 6 6 5 4 1
2 4 3 3 3 4 4 2 4 3 5 4 5 4 2 4 4 5 5 6 8 7 6 5 3
1 4 4 5 5 4 4 1 3 2 3 2 4 4 3 3 4 6 7 7 7 6 5 5 3
2 4 3 4 5 4 5 2 4 2 3 2 4 4 4 4 5 6 7 6 6 4 4 4 3
3 4 3 5 6 4 4 2 3 1 1 1 3 3 4 4 6 7 7 5 4 4 4 4 2
3 4 2 3 3 2 3 3 4 2 1 2 3 3 4 5 7 7 6 5 5 6 5 4 2
3 4 3 4 3 2 3 4 3 1 0 3 4 5 4 5 5 5 3 4 4 7 6 7 4
3 5 4 4 2 2 3 5 4 3 2 3 5 5 5 5 6 6 4 4 3 5 5 6 4
3 4 3 3 2 2 4 5 4 2 3 4 6 5 5 5 6 6 5 4 2 2 3 5 4
3 5 3 4 2 3 3 5 5 5 6 5 5 2 4 4 7 6 6 3 2 1 2 2 2
3 5 3 4 2 3 3 5 4 4 5 5 4 1 3 3 6 5 6 3 3 2 2 1 1
3 5 3 4 2 4 4 6 4 5 5 5 4 2 3 2 5 4 5 2 3 3 4 3 2
4 5 3 3 2 3 3 4 2 3 3 4 3 2 2 1 4 4 5 2 2 2 3 4 3
3 3 1 2 3 5 4 3 2 4 5 5 4 3 2 0 2 3 4 3 3 3 3 4 3
3 3 1 2 4 5 3 2 2 5 5 4 2 1 1 0 2 3 3 3 4 4 2 2 2
2 3 1 2 4 6 5 3 3 6 7 5 4 3 3 1 2 4 5 5 5 5 4 4 3
3 4 2 2 3 5 4 3 3 6 6 5 3 3 3 3 5 6 6 5 5 6 6 7 5
2 3 3 3 3 4 4 4 4 5 5 4 3 3 3 3 5 7 8 6 5 5 6 7 5
2 3 5 5 4 3 2 4 4 4 3 4 3 3 3 5 6 6 5 4 5 6 7 7 5
3 4 6 5 4 1 1 3 3 2 1 3 4 5 4 5 5 5 4 3 4 5 6 6 4
4 6 8 6 5 2 1 2 2 3 2 4 4 6 6 8 7 5 2 2 3 5 5 6 4
4 5 6 5 6 3 2 2 2 4 2 4 3 6 5 6 6 6 4 4 4 5 4 5 4
2 3 4 5 7 5 3 2 2 5 3 4 1 3 3 5 5 5 4 5 5 4 2 3 3
1 1 1 3 5 4 2 1 1 3 2 3 1 2 1 2 2 3 3 4 4 3 2 2 2

ソルバの解答

..#..##......#.#...#..#..
##..#.#..#.##...#..#.##.#
..#..#...#...#..#######..
..#..#.#.#.##.#.#..###.##
#.####.......#...###..#.#
#....#.#.#...#.#####.#...
#..##..#....#..##.#..##.#
#.#.....#...#..#.##.###.#
#..#..##....###.#...#.###
#.##..#.##.#.#.####....#.
#......#...##..#.###....#
#.#.#.#.####...#.#...#...
#.##..#.#..#...#.##..#...
#.....##..#.##...#.#..###
##..#.....#......#......#
....###..###.#....#.##...
#....#..#.#......#..##..#
#.#..##..###.##..###..###
.#..#.#..##.#..####.#####
...#...##........###.#..#
.####...#..##.###....####
##.#........###.##..#.##.
.###.#..#.#...####..#...#
#..###..#.#.#.#..##.##.##
....##....#.....#..##....

2022-06-11

6人の壁画の手数を短縮した話

ゲームプログラミング

こんなもん要はただのスライドパズルなので探索かければいけるやろ

趣旨
方針
IDA*実装
IDA*実行1
IDA*実行2
解についての考察1
解についての考察2
さらなる短縮へ
BFS
おわりに
おわらせない
双方向BFS
双方向BFS実行
おわり
ソースコード
- IDA*
- BFS
- 双方向BFS

趣旨

6人の壁画は61手解法(http://homentosu.blog.fc2.com/blog-entry-71.html)が知られているが、かなり人間的にわかりやすい操作をしている。これが最適でない可能性が少なからずある、というかさすがにもっと短くなるはず。実機でちまちまやるのは面倒だし間違いが起きそうなのでプログラムを書いてパソコンに頑張ってもらいましょう

方針

何種類か思いつくので検討してみる。

(1) 幅優先探索(BFS)

最初の盤面から（ありうるすべての）一手先盤面を考える。一手先の盤面（のうち初見のものすべて）から（ありうるすべての）二手先盤面を考える。目的の盤面が見つかるまで続ける。目的の盤面が見つかったらその道筋を復元する。

(2) A*

BFSするとき、何らかの評価値にもとづいて、評価値がよいものから優先して探索する。スライドパズルでいえば「このままいくと全体でX手以上はかかる」という見積もりを評価値として、Xが小さいものから探索する手法が一般的。

こういうBFS的手法は手数が多いと探索する盤面がどんどん増え、それに比例してメモリを食う。探索する盤面の総数が1000万ぐらいなら余裕だが、それより圧倒的に多くなる気がするのでちょっとまずそう。

……実はこの「まずそう」をきちんと評価しておけばよかったと後悔することになるのだが、それはまたのちほど。

次に深さ優先探索(DFS)的手法。イメージは、完成 or 手詰まりになるまで探索して、手詰まりならちょっと戻して別の方法を試していくという感じ。強みはメモリが少なくて済む点。スライドパズルに適用する場合は無限ループを回避するような工夫が必要。

(3) 反復深化深さ優先探索(IDDFS)

X手で行ける、と思ってX手先までDFS。だめならX+1手と思い直して最初からやり直す。実行時間はかさむが、基本的にX+1手のDFSで探索する盤面>X手のDFSで探索する盤面なので、最初から必要手数が分かっている場合と比較して致命的に遅くなるということはないと考えていい。

(4) 反復深化A* (IDA*)

基本はIDDFSと同じだが、X手で行ける可能性がなくなった時点で（X手まで探索せずに）早期撤退。不要な探索が減り、高速化につながる。

今回はIDA*を採用した。

IDA*実装

盤面の表現

プログラム上で盤面をどう表現すればよいだろうか。簡単に思いつくのは、「タイルに番号を振って二次元の配列上に盤面を見たまま表現する」ことだろう。

初期盤面（5人）を例にとるとこういう感じ。タイルに0から順に番号を振っていき、タイルがない場所は-1とした。

こうすると目的盤面（6人）はこのように表現できる。

というわけで、「6人の壁画を解く」という行為はプログラム上では「上の配列を変化させていって、下の配列にたどり着く」ことと言い表せる。

タイルをスライドするとは

「配列を変化させる」といっても、好き放題に変化させていいわけではなく、タイルのスライドに対応する変化しか許されない。スライドは「0から10までの好きな数字nと上下左右のうち好きな方向を1つ選んで、配列上にあるすべてのnをその方向にずらすこと」になるが、次の条件を満たす必要がある。

ずらした結果nが配列外にはみ出てはいけない（タイルは盤外に行けない）
-1（空き）でない場所をnで上書きするようなことがあってはいけない（タイルは重ならない）

よりプログラムにしやすい表現にしてみよう。簡単のためxの正の方向とする。

nが(x1, y1), (x2, y2), ... にあるとき、配列の(x1, y1), (x2, y2), ... 番目を-1で書き換えたうえで、(x1+1, y1), (x2+1, y2) ... 番目をnに書き換える。ただし

どの(x?+1, y?)も配列外になってはいけない
(x?+1, y?)をnに書き換える前に-1以外の数字が書かれていてはいけない

より効率的に配列を表現する方法

盤面サイズは5x5で25なので、配列を配列のまま保持するとint型(32bit) x 25個のデータ量になってしまう。IDA*をやるうえでそれが問題になることはほぼないのだが、もっと小さいデータ量で表現することができる。

盤面を復元するためには、各ピースの座標が分かればじゅうぶんである。

→ピースの各座標を(x, y)の配列にすると32bit変数2個 x 11個で足りる。

→(x, y)は25通りしかないので、(x, y)は32bit変数2個でなく1個で足りる。

具体的には、f (y, x) = 5y + xとすれば(0, 0)から(4, 4)までの数字の組は0から24の25個の数字に1対1で対応付けられる。またf(y, x)の値から y と x を復元することも容易で y = f(y, x) / 5, x = f(y, x) % 5

→もっというと、5bitあれば足りる（2⁵=32>25）。

→5bitの情報が11個あればいいので、55bitあれば足りる。

→64bit変数1個に圧縮できる。

評価値

IDA*をやるうえで重要になるのが評価値である。スライドパズルでは「完成までの手数がこれだけはかかる」という値を評価値とすることが多い。この値の見積もりに使う指標として最も基本的かつ簡単なものがマンハッタン距離である。

2次元空間でのマンハッタン距離は、2つの座標 (x1, y1), (x2, y2) に対して abs(x1 - x2) + abs(y1 - y2) で表される。ここで abs(x) は x の絶対値である。

スライドパズルにおける「マンハッタン距離」は現在の盤面と目的盤面を比較したときの各タイルのマンハッタン距離の総和であり、パズルを完成するにはこれ以上の手数が必要となる（証明：目的盤面とのマンハッタン距離はタイルのスライド1回で1より大きく減らない）。

よって、現在動かした手数+目的盤面とのマンハッタン距離がIDA*の深さ上限より大きくなったら、その時点で探索を打ち切ってよい。（なお、実際のところマンハッタン距離はあまりよい指標ではない）

ループ回避

タイルを右に動かして左に動かす行為は意味がない。それは何もしていないに等しいからである。IDA*で手数を決め打つとはいえ、こういう無駄な部分はなるべく避けたい。一般的な、全部のタイルが1x1のサイズで空きが1か所しかないようなパズルの場合、「1手前に戻らない」とするだけでよいのだが、6人の壁画では空きマスが多いためこれでは足りない。タイルAを右、タイルBを上、タイルAを左、タイルBを下、という一連の操作は、結局何もしていないのに、1手前に戻るような操作は行われていないため「1手前に戻らない」という処理では検出できない。ということは、今まで経由してきたすべての盤面を覚えておくのがよさそうである。

一般に「x が集合 S に含まれているか」の判定問題はC++ではsetを使って高速に行える。setで盤面を管理しながら、過去に探索した盤面にループしてきた場合は探索を打ち切る。このとき盤面の比較がsetの内部では行われるはずだが、上で書いた盤面情報を圧縮する方法を使うことで、25要素の配列の比較から整数1つの比較にでき、判定の高速化を図れる。

出力形式

1手目から順に、

タイル番号（スペース）移動方向

を必要な行数出力。（探索状況とか実行時間とかも出力するけど）

IDA*実行1

その他細かいことは実装に譲って、とにかく実行してみよう。

depth:19
0 sec.
depth:20
0 sec.
depth:21
0 sec.
depth:22
0 sec.
depth:23
0 sec.
depth:24
0 sec.
depth:25
0 sec.
depth:26
0 sec.
depth:27
12 sec.
depth:28
12 sec.
depth:29
861 sec.
depth:30
860 sec.
depth:31
^C
続行するには何かキーを押してください . . .

終わんねぇ～～～～～～

終わらな過ぎて途中で強制終了した。探索する盤面の数が爆発してしまっているようだ。

IDA*実行2

解全体の改善が無理なら、部分解だけでも改善したい。

このパズルの最も難しいところは縦に長い3番と4番を入れ替えるところであると言ってよいだろう。この2つを入れ替えるには4マスの縦幅が必要だが、それを確保するために61手解はタイルをかなり多く移動させており、とくに8番タイル（元サイトでは青1マス、と表現されている）の移動が激しい。この部分、削れる気がする。

というわけで、先に示したサイトで紹介されている解法を26手まで進めた盤面（(2)の終わりまで）に対して再試行。最適解の保証はなくなるが、ショートカットぐらいは見つかってもおかしくない。

depth:14
0 sec.
depth:15
0 sec.
depth:16
2 sec.
depth:17
2 sec.
depth:18
193 sec.
depth:19
193 sec.
depth:20
solution found
20 moves
0 L
0 D
6 U
6 U
6 L
6 L
3 U
6 L
8 U
8 L
4 L
1 D
1 D
4 D
3 R
3 R
4 U
4 U
8 R
8 D
(226 sec.)

終わった。

なんと26手を6手短縮して、20手で行けるらしい。つまり、このパズルの手数は61手から55手に短縮できる。

解についての考察1

IDA*の出力解に対して、同一タイルの連続移動をまとめて（可能なら適宜順番を入れ替えつつ）見やすくする。

0 LD
6 UULLL
3 U
8 UL
4 LD
1 DD
3 RR
4 UU
8 RD

8番のタイルを動かす回数がかなり減っている。

ところで、6 UULLLに注目したい。本来は6 UULLと6 Lに分かれていて間に3 Uがあった部分だ。位置関係から6 UULLは確かに3 Uの前に行う必要があるのだが、後ろの6 Lは3 Uの後でも問題ない操作なうえ、実際に手を動かしてみるとこの6 Lの必要性が（少なくとも20手目までは）存在しないことが分かる。この後どうしても6 Lの必要性が出てくる場合は仕方ないが、そうでないならこの操作は削ってしまえる。

というわけでこの6 Lをしなかった盤面をグッと睨んだ結果、6 Lが不要であることがわかった（ここは自分の頭で考えた）。

53手解法。

0 LD
6 UULL
3 U
8 UL
4 LD
1 DD
3 RR
4 UU
8 RD
2 RU 
0 RRDR
2 DD
4 L
3 L
1 UU
0 RU
2 RR
5 RR
6 DD
4 LLD
3 LL
1 LL
0 U
2 U
5 R
6 R
3 D
1 LL

解についての考察2

53手解中盤の 2 RU / 0 RRDR / 2 DD に注目したい。これは何をやっているかというと、0と2の順番を入れ替えているのだ。だが、ここの順番をもとから逆にしておけばこの操作は不要となり解が短縮できる。実際、最初の0 LD を 2 L / 0 Dとすればよく、51手解が発見できた。

図がちょっと誤植してる。本質変わんないから許して

51手解

2 L
0 D
6 UULL
3 U
8 UL
4 L
1 DD
4 D
3 RR
4 UU
8 RD
0 RDR
2 RRD
4 L
3 L
1 UU
0 RU
2 RR
5 RR
6 DD
4 LLD
3 LL
1 LL
0 U
2 U
5 R
6 R
3 D
1 LL

だが、IDA*の結果を見て改善できそうな部分はここまで。

6人の壁画の必要手数は51手でひとまず結論づいた。

……かに思えたが。

さらなる短縮へ

このパズルの重要な性質を見落としていた。

下半分、全然変わらない。

8番は縦長タイルを入れ替える関係で動かす。それはわかるが、そもそも頭を付け替えることで人数が増えるというギミックのこのパズルにおいて、下のほうにあるうえにデカい7番9番10番を動かす必要性はあまり感じられない。

こいつら、動かさないものと仮定してみるか。

こいつらを動かさないことに決めると、このパズルは一気に盤面数が小さくなる。見積もってみよう。

小さくなった壁画は、マス16個、2マスのタイル3個、1マスのタイル5個である。

2マスのタイルから置いていく。最初のタイルを置く方法は16通り（本当は盤外になるパターンがあるのでもっと小さい）、次は14マスから選ぶ。16x14通り（盤外とか干渉とか無視）。同様にやっていくと、

16 x 14 x 12 x 10 x 9 x 8 x 7 x 6 = 81285120

本当はもっと小さいという事実と、これを全部探索し尽くすよりはそれなりに早い段階で解が見つかるだろうという希望的観測を合わせると、幅優先探索が走れそうな問題サイズまで小さくなった。

BFS

幅優先探索のコードを書いた。枝刈りとかはしてない。

手元のパソコンのメモリを1GB以上食ったが、43手解が見つかった。もしこれが最適でない場合、下のでかい奴らのどれかを1回以上は動かす必要がある。

solution found (81 sec.)
43 moves:
0 L
1 L
1 L
1 L
4 U
6 U
6 U
6 L
8 U
8 R
3 R
2 R
1 D
1 L
2 D
3 D
4 L
6 L
4 L
3 U
3 U
8 L
8 D
2 R
2 R
2 U
5 R
5 R
5 R
4 D
4 L
3 L
3 D
6 R
0 R
1 U
4 L
3 L
6 D
0 R
0 R
0 R
6 D

きれいに書き直したものがこちら。

0 L
1 LLL
4 U
6 UUL
8 UR
3 RD
2 RD
1 DL
6 L
4 LL
3 UU
8 LD
2 RRU
5 RRR
4 DL
3 LD
6 R
0 R
1 U
4 L
3 L
6 DD
0 RRR

ちなみに、最適解で9番を動かす可能性は否定できないと思い、そのように書き直したものを実行してみたが、僕のパソコンでは無理そうだった。

時間とメモリがもっと潤沢にあればあるいは、いやどうだろう

おわりに

43手。もっと適切に枝刈りして双方向探索とかするともう少しどうにかなるかもしれないけど、今のところはその元気はないかな。IDA*を書いたのは初めてだったのでいい経験になりました。

おわらせない

43手解が見つかった。これより小さい手数の解があるとして、それは42手ではありえない（マンハッタン距離の偶奇性より）。あるとしたら41手か39手、37手……である。

動かさないタイルを設定したとはいえ、単純なBFSで43手先まで探索できたのなら、動かさないタイルの制限を取っ払っても20手ぐらいなら余裕で探索できるはず。双方向BFSが使えそう。

双方向BFS

初期盤面からBFSすると、手数が増えるにつれ盤面数は（たいてい）多くなっていくので、探索空間は下図左のような感じになる。しかし、目的盤面からもBFSして真ん中で落ち合うことにすれば、探索空間は下図右のようになり効率的である。

前から21手の盤面と後ろから20手の盤面の積集合が空でなければ、41手解が存在することがいえる。同様に前から19手の盤面と後ろから20手の盤面の積集合が空でなければ、39手解が存在することが言える。

とりあえず動かし方とかは考えないことにして、そういう盤面があるならそれを出力するプログラムを書いた。

双方向BFS実行

数字を変えて実行したが、41手、39手、37手、35手、33手、31手解は見つからなかった（それ以下の解がないことは最初のIDA*で実証済）。

おわり

いうことで6人の壁画の最短手数は43手、一例は上に示した通りです。この下にはソースコードしかないので見たい人だけ見てください。

ソースコード

IDA*

#include <iostream>
#include <vector>
#include <set>
#include <stack>
#include <algorithm>
#include <utility>
#include <chrono>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
using namespace std;

// 諸定数
const int H=5; // 縦幅
const int W=5; // 横幅
const int N=11; // ピースの数
const int space=-1; // 空白

// 方向
const string dirStr = "LURD";
const vector<int> dy = {0,-1,0,1};
const vector<int> dx = {-1,0,1,0};


// ピースの左上位置情報
const pair<int,int> undefined = {-1,-1};
vector<pair<int,int>> firstPos(N,undefined);
vector<pair<int,int>> lastPos(N,undefined);
// ピース形状
// (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
vector<vector<pair<int,int>>> pieceShape(N);

// 判定: 座標(y,x)は盤面内
inline bool inBoard(const int y, const int x){
  return 0<=y && y<H && 0<=x && x<W;
}
// 判定: ピースの左上を(y,x)としてもピース全体が盤面内に収まる
bool pieceInBoard(const int y, const int x, const vector<pair<int,int>>& dydx){
  for(const auto& p:dydx){
    if(!inBoard(y+p.first,x+p.second)) return false;
  }
  return true;
}

// 盤面のハッシュ ピースの位置から64bit整数(使うのは55bit)
const int hashDataWidth = 5;
const int mask = (1<<hashDataWidth)-1;
long long calcHash(const vector<pair<int,int>>& pos){
  long long res = 0;
  rep(i,N){
    res<<=hashDataWidth;
    int y = pos[N-1-i].first, x = pos[N-1-i].second;
    res|=(y*W+x);
  }
  return res;
}
// ハッシュ値から指定されたピースの場所を特定
inline pair<int,int> getPos(const long long hash, const int pieceID){
  int p = (hash>>(pieceID*hashDataWidth))&mask;
  return make_pair(p/W,p%W);
}
// 指定のピースを指定の方向に動かしたときのハッシュ
long long calcSlidedHash(const long long hash, const int pieceID, const int dir){
  pair<int,int> pos = getPos(hash,pieceID);
  int y = pos.first, x = pos.second;
  y += dy[dir]; x += dx[dir];
  long long yx= y*W+x;
  long long res = hash;
  int shiftSize = pieceID*hashDataWidth;
  res ^= ((res>>shiftSize)&mask)<<shiftSize;
  res ^= (yx<<shiftSize);
  return res;
}

// 盤面から各ピースの位置情報を取得
void calcPiecePosList(const vector<vector<int>>& Board, vector<pair<int,int>>& pos){
  rep(y,H){
    rep(x,W){
      int pieceID = Board[y][x];
      if(pieceID!=space && pos[pieceID]==undefined){
        pos[pieceID] = make_pair(y,x);
      }
    }
  }
}

// ピース形状の取得
// (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
void calcPieceShapeList(const vector<vector<int>>& Board, vector<vector<pair<int,int>>>& dydx){
  rep(y,H){
    rep(x,W){
      if(Board[y][x]!=space){
        int pieceID = Board[y][x];
        pair<int,int> pos = firstPos[pieceID];
        dydx[pieceID].push_back(make_pair(y-pos.first, x-pos.second));
      }
    }
  }
}

// 判定:指定の場所に指定の形状のピースを置くと他のピースと干渉する
bool conflict(const vector<vector<bool>>& NGPos, const vector<pair<int,int>>& pieceShape, const int y, const int x){
  for(const pair<int,int>& dydx:pieceShape){
    if(NGPos[y+dydx.first][x+dydx.second]) return true;
  }
  return false;
}

// 2点のマンハッタン距離
inline int MD(const pair<int,int> a, const pair<int,int> b){
  return abs(a.first-b.first)+abs(a.second-b.second);
}
// 2つの盤面の各ピースのマンハッタン距離の総和
// 基本はMDによる差分計算 乱用は控える
int MDSum(const vector<pair<int,int>>& pos1, const vector<pair<int,int>>& pos2){
  int res = 0;
  rep(i,N){
    res += MD(pos1[i],pos2[i]);
  }
  return res;
}

// 初期盤面
const vector<vector<int>> firstBoard = {
  {-1,  0, -1, -1,  1},
  {-1,  2,  3, -1,  4},
  { 5,  5,  3,  6,  4},
  { 7,  7,  7,  8,  9},
  { 7,  7, 10, 10,  9}
};

// 26手から20手に
const vector<vector<int>> lastBoard = {
  { 6, -1, -1,  4,  3},
  { 0,  2, -1,  4,  3},
  { 5,  5, -1, -1,  1},
  { 7,  7,  7,  8,  9},
  { 7,  7, 10, 10,  9}
};


// 51手を達成するための盤面 これは動く
/*
const vector<vector<int>> lastBoard = {
  {-1,  6, -1,  4,  3},
  { 2,  0, -1,  4,  3},
  { 5,  5, -1, -1,  1},
  { 7,  7,  7,  8,  9},
  { 7,  7, 10, 10,  9}
};
*/

// 完成盤面 これはさすがに終わらない
/*
const vector<vector<int>> lastBoard = {
  { 1, -1, -1, -1,  0},
  { 4,  3, -1, -1,  2},
  { 4,  3,  6,  5,  5},
  { 7,  7,  7,  8,  9},
  { 7,  7, 10, 10,  9}
};
*/


long long firstBoardHash, lastBoardHash;


int dfs(const int pred, const int depth, const int maxDepth, const long long boardHash, vector<vector<bool>>& NGPos,
  set<long long>& hashSet, stack<pair<int,char>>& operationStack){

  // 明らかにループしているなら探索不要
  if(hashSet.count(boardHash)) return 0;

  // 今の盤面をチェック済みリストにメモ
  hashSet.insert(boardHash);

  // 探索成功
  if(boardHash==lastBoardHash){
    cout << "solution found\n";
    vector<pair<int,int>> ans;
    while(!operationStack.empty()){
      ans.push_back(operationStack.top()); operationStack.pop();
    }
    reverse(ans.begin(),ans.end());
    cout << depth << " moves\n";
    for(auto a:ans){
      cout << a.first << ' ' << (char)a.second << '\n';
    }
    return 1;
  }
  rep(pieceID,N){ // 動かすピース
    pair<int,int> pos = getPos(boardHash, pieceID);
    int y = pos.first, x = pos.second;

    for(auto dydx:pieceShape[pieceID]){
      NGPos[y+dydx.first][x+dydx.second]=false;
    }
    rep(dir,4){ // 移動方向
      // ピースの移動先座標
      int nextY = y+dy[dir], nextX=x+dx[dir];
      // 盤面外には出られない
      if(!pieceInBoard(nextY,nextX,pieceShape[pieceID])) continue;
      // ほかのピースに被ってはいけない
      if(conflict(NGPos, pieceShape[pieceID], nextY, nextX)) continue;
      // 次の状態を計算する
      // 深さ = 手数
      int nextDepth = depth+1;
      // マンハッタン距離 (差分計算)
      int nowMD = pred-depth;
      int nextMD = nowMD-MD(pos,lastPos[pieceID])+MD(make_pair(nextY,nextX),lastPos[pieceID]);
      // 評価値 少なくともこれより少ない手数では不可能
      int nextPred = nextMD+nextDepth;
      // maxDepth手より多くなることが確定したら何もしない
      if(nextPred>maxDepth) continue;
      // ハッシュ
      long long nextHash = calcSlidedHash(boardHash, pieceID, dir);

      // 操作をスタックに積む
      operationStack.push(make_pair(pieceID,dirStr[dir]));
      for(auto dydx:pieceShape[pieceID]){
        NGPos[nextY+dydx.first][nextX+dydx.second]=true;
      }
      // 探索継続 成功したら探索は打ち切り
      if(dfs(nextPred,nextDepth,maxDepth,nextHash,NGPos,hashSet,operationStack)) return 1;
      // スタックに積んだ操作を削除
      operationStack.pop();
      for(auto dydx:pieceShape[pieceID]){
        NGPos[nextY+dydx.first][nextX+dydx.second]=false;
      }
    }
    for(auto dydx:pieceShape[pieceID]){
      NGPos[y+dydx.first][x+dydx.second]=true;
    }
  }
  // チェック済みリストから今の盤面を削除
  hashSet.erase(boardHash);
  return 0;
}

int main(){
  // 初期位置情報(左上)の取得
  calcPiecePosList(firstBoard, firstPos);
  // 完成位置情報の取得
  calcPiecePosList(lastBoard, lastPos);
  // ピース形状の取得
  calcPieceShapeList(firstBoard, pieceShape);
  // 初期盤面と完成盤面のハッシュ計算
  firstBoardHash = calcHash(firstPos);
  lastBoardHash = calcHash(lastPos);

  // IDA*準備
  int firstMD = MDSum(firstPos, lastPos);
  set<long long> hashSet;
  stack<pair<int,char>> operationStack;
  vector<vector<bool>> NGPos(H,vector<bool>(W,false));
  rep(pieceID,N){
    int y = firstPos[pieceID].first, x = firstPos[pieceID].second;
    for(auto dydx:pieceShape[pieceID]){
      NGPos[y+dydx.first][x+dydx.second] = true;
    }
  }
  // IDA*実行
  chrono::system_clock::time_point start, end;
  for(int maxDepth=firstMD;;maxDepth++){ // firstMD未満ではできないのでそこは省略
    cout << "depth:" << maxDepth << '\n';
    start = chrono::system_clock::now();
    if(dfs(firstMD,0,maxDepth,firstBoardHash,NGPos,hashSet,operationStack)){ // 成功
      end = chrono::system_clock::now();
      double ela = chrono::duration_cast<chrono::seconds> (end-start).count();
      cout << "(" << ela << " sec.)\n";
      break;
    }
    end = chrono::system_clock::now();
    double ela = chrono::duration_cast<chrono::seconds> (end-start).count();
    cout << ela << " sec." << '\n';
  }
  return 0;
}

BFS

#include <iostream>
#include <vector>
#include <stack>
#include <algorithm>
#include <utility>
#include <chrono>
#include <queue>
#include <map>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
using namespace std;

// 諸定数
const int H=5; // 縦幅
const int W=5; // 横幅
const int N=11; // ピースの数
const int space=-1; // 空白

// 方向
const string dirStr = "LURD";
const vector<int> dy = {0,-1,0,1};
const vector<int> dx = {-1,0,1,0};

// 判定: 座標(y,x)は盤面内
inline bool inBoard(const int y, const int x){
  return 0<=y && y<H && 0<=x && x<W;
}
// 判定: ピースの左上を(y,x)としてもピース全体が盤面内に収まる
bool pieceInBoard(const int y, const int x, const vector<pair<int,int>>& dydx){
  for(const auto& p:dydx){
    if(!inBoard(y+p.first,x+p.second)) return false;
  }
  return true;
}

// 盤面のハッシュ ピースの位置から64bit整数(使うのは55bit)
const int hashDataWidth = 5;
const int mask = (1<<hashDataWidth)-1;
long long calcHash(const vector<pair<int,int>>& pos){
  long long res = 0;
  rep(i,N){
    res<<=hashDataWidth;
    int y = pos[N-1-i].first, x = pos[N-1-i].second;
    res|=(y*W+x);
  }
  return res;
}
// ハッシュ値から指定されたピースの場所を特定
inline pair<int,int> getPos(const long long hash, const int pieceID){
  int p = (hash>>(pieceID*hashDataWidth))&mask;
  return make_pair(p/W,p%W);
}
// 指定のピースを指定の方向に動かしたときのハッシュ
long long calcSlidedHash(const long long hash, const int pieceID, const int dir){
  pair<int,int> pos = getPos(hash,pieceID);
  int y = pos.first, x = pos.second;
  y += dy[dir]; x += dx[dir];
  long long yx= y*W+x;
  long long res = hash;
  int shiftSize = pieceID*hashDataWidth;
  res ^= ((res>>shiftSize)&mask)<<shiftSize;
  res ^= (yx<<shiftSize);
  return res;
}

// 盤面から各ピースの位置情報を取得
const pair<int,int> undefined = {-1,-1};
void calcPiecePosList(const vector<vector<int>>& Board, vector<pair<int,int>>& pos){
  rep(y,H){
    rep(x,W){
      int pieceID = Board[y][x];
      if(pieceID!=space && pos[pieceID]==undefined){
        pos[pieceID] = make_pair(y,x);
      }
    }
  }
}

// ピース形状の取得
// (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
void calcPieceShapeList(const vector<vector<int>>& Board, const vector<pair<int,int>>& piecePos, vector<vector<pair<int,int>>>& dydx){
  rep(y,H){
    rep(x,W){
      if(Board[y][x]!=space){
        int pieceID = Board[y][x];
        pair<int,int> pos = piecePos[pieceID];
        dydx[pieceID].push_back(make_pair(y-pos.first, x-pos.second));
      }
    }
  }
}

// ピースがある場所の計算
void getNGPos(const long long hash, const vector<vector<pair<int,int>>>& pieceShape, vector<vector<bool>>& NGPos){
  rep(pieceID,N){
    pair<int,int> pos = getPos(hash,pieceID);
    int y = pos.first, x = pos.second;
    for(auto dydx:pieceShape[pieceID]){
      NGPos[y+dydx.first][x+dydx.second] = true;
    }
  }
}

// 判定:指定の場所に指定の形状のピースを置くと他のピースと干渉する
bool conflict(const vector<vector<bool>>& NGPos, const vector<pair<int,int>>& pieceShape, const int y, const int x){
  for(const pair<int,int>& dydx:pieceShape){
    if(NGPos[y+dydx.first][x+dydx.second]) return true;
  }
  return false;
}

struct prevHandData{
  int movedPieceID;
  int movedDir;
};

int main(){
  // 時間計測
  chrono::system_clock::time_point startTime, endTime;
  startTime = chrono::system_clock::now();
  // 初期盤面
  const vector<vector<int>> firstBoard = {
    {-1,  0, -1, -1,  1},
    {-1,  2,  3, -1,  4},
    { 5,  5,  3,  6,  4},
    { 7,  7,  7,  8,  9},
    { 7,  7, 10, 10,  9}
  };
  // 完成盤面
  const vector<vector<int>> lastBoard = {
    { 1, -1, -1, -1,  0},
    { 4,  3, -1, -1,  2},
    { 4,  3,  6,  5,  5},
    { 7,  7,  7,  8,  9},
    { 7,  7, 10, 10,  9}
  };
  // 初期/完成盤面のピースの左上位置情報
  vector<pair<int,int>> firstPos(N,undefined);
  vector<pair<int,int>> lastPos(N,undefined);
  calcPiecePosList(firstBoard, firstPos);
  calcPiecePosList(lastBoard, lastPos);
  // ピース形状の取得
  // (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
  vector<vector<pair<int,int>>> pieceShape(N);
  calcPieceShapeList(firstBoard, firstPos, pieceShape);
  // 初期盤面と完成盤面のハッシュ計算
  const long long firstBoardHash = calcHash(firstPos), lastBoardHash = calcHash(lastPos);

  queue<long long> que; que.push(firstBoardHash);
  map<long long,prevHandData> prevHand; prevHand[firstBoardHash]={-1,'X'};
  while(!que.empty()){
    const long long boardHash = que.front(); que.pop();
    // ピースがある場所の計算
    vector<vector<bool>> NGPos(H,vector<bool>(W,false));
    getNGPos(boardHash, pieceShape, NGPos);
    rep(pieceID,N){ // 動かすピース
      if(pieceID==7||8<pieceID) continue; // ここは動かさないことに決めた
      // 動かしたいピースの位置
      const pair<int,int> pos = getPos(boardHash, pieceID);
      const int y = pos.first, x = pos.second;
      // 自分を盤面からいったん外す
      for(const pair<int,int>& dydx:pieceShape[pieceID]){
        NGPos[y+dydx.first][x+dydx.second]=false;
      }
      rep(dir,4){ // 方向
        // ピースの移動先座標
        const int nextY = y+dy[dir], nextX=x+dx[dir];
        // 盤面外には出られない
        if(!pieceInBoard(nextY,nextX,pieceShape[pieceID])) continue;
        // ほかのピースに被ってはいけない
        if(conflict(NGPos, pieceShape[pieceID], nextY, nextX)) continue;
        // 次の盤面
        const long long nextBoardHash = calcSlidedHash(boardHash, pieceID, dir);
        // 探索済みならパス
        if(prevHand.count(nextBoardHash)) continue;
        // キューに突っ込み自分の一手前をメモ
        que.push(nextBoardHash);
        prevHand[nextBoardHash] = {pieceID, dir};
        // 答えが見つかったら終了
        if(nextBoardHash==lastBoardHash) goto outOfSearch;
      }
      // 自分を盤面に戻す
      for(auto dydx:pieceShape[pieceID]){
        NGPos[y+dydx.first][x+dydx.second]=true;
      }
    }
  }
  outOfSearch:
  endTime = chrono::system_clock::now();
  double ela = chrono::duration_cast<chrono::seconds> (endTime-startTime).count();
  cout << "solution found (" << ela << " sec.)\n";
  // 手の復元
  long long nowHash = lastBoardHash;
  stack<pair<int,char>> stk;
  while(nowHash!=firstBoardHash){
    const prevHandData phd = prevHand[nowHash];
    stk.push(make_pair(phd.movedPieceID,phd.movedDir));
    nowHash = calcSlidedHash(nowHash,phd.movedPieceID,phd.movedDir^2); // 反対方向に動かす
  }
  // 解の出力
  cout << (int)stk.size() << " moves:\n";
  while(!stk.empty()){
    cout << stk.top().first << ' ' << dirStr[stk.top().second] << '\n'; stk.pop();
  }
  return 0;
}

双方向BFS

#include <iostream>
#include <vector>
#include <stack>
#include <algorithm>
#include <utility>
#include <chrono>
#include <queue>
#include <set>
#include <map>
#define rep(i,n) for(int i=0;i<(int)(n);i++)
using namespace std;

// 諸定数
const int H=5; // 縦幅
const int W=5; // 横幅
const int N=11; // ピースの数
const int space=-1; // 空白

// 方向
const string dirStr = "LURD";
const vector<int> dy = {0,-1,0,1};
const vector<int> dx = {-1,0,1,0};

// 判定: 座標(y,x)は盤面内
inline bool inBoard(const int y, const int x){
  return 0<=y && y<H && 0<=x && x<W;
}
// 判定: ピースの左上を(y,x)としてもピース全体が盤面内に収まる
bool pieceInBoard(const int y, const int x, const vector<pair<int,int>>& dydx){
  for(const auto& p:dydx){
    if(!inBoard(y+p.first,x+p.second)) return false;
  }
  return true;
}

// 盤面のハッシュ ピースの位置から64bit整数(使うのは55bit)
const int hashDataWidth = 5;
const int mask = (1<<hashDataWidth)-1;
long long calcHash(const vector<pair<int,int>>& pos){
  long long res = 0;
  rep(i,N){
    res<<=hashDataWidth;
    int y = pos[N-1-i].first, x = pos[N-1-i].second;
    res|=(y*W+x);
  }
  return res;
}
// ハッシュ値から指定されたピースの場所を特定
inline pair<int,int> getPos(const long long hash, const int pieceID){
  int p = (hash>>(pieceID*hashDataWidth))&mask;
  return make_pair(p/W,p%W);
}
// 指定のピースを指定の方向に動かしたときのハッシュ
long long calcSlidedHash(const long long hash, const int pieceID, const int dir){
  pair<int,int> pos = getPos(hash,pieceID);
  int y = pos.first, x = pos.second;
  y += dy[dir]; x += dx[dir];
  long long yx= y*W+x;
  long long res = hash;
  int shiftSize = pieceID*hashDataWidth;
  res ^= ((res>>shiftSize)&mask)<<shiftSize;
  res ^= (yx<<shiftSize);
  return res;
}

// 盤面から各ピースの位置情報を取得
const pair<int,int> undefined = {-1,-1};
void calcPiecePosList(const vector<vector<int>>& Board, vector<pair<int,int>>& pos){
  rep(y,H){
    rep(x,W){
      int pieceID = Board[y][x];
      if(pieceID!=space && pos[pieceID]==undefined){
        pos[pieceID] = make_pair(y,x);
      }
    }
  }
}

// ピース形状の取得
// (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
void calcPieceShapeList(const vector<vector<int>>& Board, const vector<pair<int,int>>& piecePos, vector<vector<pair<int,int>>>& dydx){
  rep(y,H){
    rep(x,W){
      if(Board[y][x]!=space){
        int pieceID = Board[y][x];
        pair<int,int> pos = piecePos[pieceID];
        dydx[pieceID].push_back(make_pair(y-pos.first, x-pos.second));
      }
    }
  }
}

// ピースがある場所の計算
void getNGPos(const long long hash, const vector<vector<pair<int,int>>>& pieceShape, vector<vector<bool>>& NGPos){
  rep(pieceID,N){
    pair<int,int> pos = getPos(hash,pieceID);
    int y = pos.first, x = pos.second;
    for(auto dydx:pieceShape[pieceID]){
      NGPos[y+dydx.first][x+dydx.second] = true;
    }
  }
}

// 判定:指定の場所に指定の形状のピースを置くと他のピースと干渉する
bool conflict(const vector<vector<bool>>& NGPos, const vector<pair<int,int>>& pieceShape, const int y, const int x){
  for(const pair<int,int>& dydx:pieceShape){
    if(NGPos[y+dydx.first][x+dydx.second]) return true;
  }
  return false;
}

struct prevHandData{
  int movedPieceID;
  int movedDir;
};

int main(){
  // 時間計測
  chrono::system_clock::time_point startTime, endTime;
  startTime = chrono::system_clock::now();
  // 初期盤面
  const vector<vector<int>> firstBoard = {
    {-1,  0, -1, -1,  1},
    {-1,  2,  3, -1,  4},
    { 5,  5,  3,  6,  4},
    { 7,  7,  7,  8,  9},
    { 7,  7, 10, 10,  9}
  };
  // 完成盤面
  const vector<vector<int>> lastBoard = {
    { 1, -1, -1, -1,  0},
    { 4,  3, -1, -1,  2},
    { 4,  3,  6,  5,  5},
    { 7,  7,  7,  8,  9},
    { 7,  7, 10, 10,  9}
  };
  // 初期/完成盤面のピースの左上位置情報
  vector<pair<int,int>> firstPos(N,undefined);
  vector<pair<int,int>> lastPos(N,undefined);
  calcPiecePosList(firstBoard, firstPos);
  calcPiecePosList(lastBoard, lastPos);
  // ピース形状の取得
  // (a,b) in pieceShape[i] = ピース[i]の左上(y,x)とすると(y+a,x+b)もピース[i]
  vector<vector<pair<int,int>>> pieceShape(N);
  calcPieceShapeList(firstBoard, firstPos, pieceShape);
  // 初期盤面と完成盤面のハッシュ計算
  const long long firstBoardHash = calcHash(firstPos), lastBoardHash = calcHash(lastPos);

  // スタートからfrontDepth手
  const int frontDepth = 21;
  cout << "from start\n";
  queue<long long> que; que.push(firstBoardHash);
  set<long long> hashSet; hashSet.insert(firstBoardHash);
  set<long long> d21;
  rep(depth,frontDepth){
    cout << "depth: " << depth << " queue size:" << (long long)que.size() << '\n';
    queue<long long> que2;
    while(!que.empty()){
      const long long boardHash = que.front(); que.pop();
      // ピースがある場所の計算
      vector<vector<bool>> NGPos(H,vector<bool>(W,false));
      getNGPos(boardHash, pieceShape, NGPos);
      rep(pieceID,N){ // 動かすピース
        // 動かしたいピースの位置
        const pair<int,int> pos = getPos(boardHash, pieceID);
        const int y = pos.first, x = pos.second;
        // 自分を盤面からいったん外す
        for(const pair<int,int>& dydx:pieceShape[pieceID]){
          NGPos[y+dydx.first][x+dydx.second]=false;
        }
        rep(dir,4){ // 方向
          // ピースの移動先座標
          const int nextY = y+dy[dir], nextX=x+dx[dir];
          // 盤面外には出られない
          if(!pieceInBoard(nextY,nextX,pieceShape[pieceID])) continue;
          // ほかのピースに被ってはいけない
          if(conflict(NGPos, pieceShape[pieceID], nextY, nextX)) continue;
          // 次の盤面
          const long long nextBoardHash = calcSlidedHash(boardHash, pieceID, dir);
          // 探索済みならパス
          if(hashSet.count(nextBoardHash)) continue;
          hashSet.insert(nextBoardHash);
          // キューに突っ込む
          que2.push(nextBoardHash);
        }
        // 自分を盤面に戻す
        for(auto dydx:pieceShape[pieceID]){
          NGPos[y+dydx.first][x+dydx.second]=true;
        }
      }
    }
    swap(que,que2);
  }
  cout << "depth: " << frontDepth << " queue size:" << (long long)que.size() << '\n';
  while(!que.empty()){
    d21.insert(que.front()); que.pop();
  }
  // ゴールからbackDepth手
  const int backDepth = 20;
  cout << "from goal\n";
  queue<long long> reverseQue; reverseQue.push(lastBoardHash);
  set<long long> reverseHashSet; reverseHashSet.insert(lastBoardHash);
  rep(depth,backDepth){
    cout << "depth: " << depth << " queue size:" << (long long)reverseQue.size() << '\n';
    queue<long long> reverseQue2;
    while(!reverseQue.empty()){
      const long long boardHash = reverseQue.front(); reverseQue.pop();
      // ピースがある場所の計算
      vector<vector<bool>> NGPos(H,vector<bool>(W,false));
      getNGPos(boardHash, pieceShape, NGPos);
      rep(pieceID,N){ // 動かすピース
        // 動かしたいピースの位置
        const pair<int,int> pos = getPos(boardHash, pieceID);
        const int y = pos.first, x = pos.second;
        // 自分を盤面からいったん外す
        for(const pair<int,int>& dydx:pieceShape[pieceID]){
          NGPos[y+dydx.first][x+dydx.second]=false;
        }
        rep(dir,4){ // 方向
          // ピースの移動先座標
          const int nextY = y+dy[dir], nextX=x+dx[dir];
          // 盤面外には出られない
          if(!pieceInBoard(nextY,nextX,pieceShape[pieceID])) continue;
          // ほかのピースに被ってはいけない
          if(conflict(NGPos, pieceShape[pieceID], nextY, nextX)) continue;
          // 次の盤面
          const long long nextBoardHash = calcSlidedHash(boardHash, pieceID, dir);
          // 探索済みならパス
          if(reverseHashSet.count(nextBoardHash)) continue;
          reverseHashSet.insert(nextBoardHash);
          // キューに突っ込む
          reverseQue2.push(nextBoardHash);
        }
        // 自分を盤面に戻す
        for(auto dydx:pieceShape[pieceID]){
          NGPos[y+dydx.first][x+dydx.second]=true;
        }
      }
    }
    swap(reverseQue,reverseQue2);
  }
  cout << "depth: " << backDepth << " queue size:" << (long long)reverseQue.size() << '\n';
  bool found = false;
  while(!reverseQue.empty()){
    // 解が見つかったら
    if(d21.count(reverseQue.front())){
      cout << "solution found\n";
      cout << reverseQue.front() << '\n';
      found = true;
      break;
    }
    reverseQue.pop();
  }
  if(!found) cout << "solution not found\n";
  endTime = chrono::system_clock::now();
  double ela = chrono::duration_cast<chrono::seconds> (endTime-startTime).count();
  cout << ela << " sec.\n";
  return 0;
}

2	3	2	2	2	4	3	2	1	1	2	2	3	2	2	2	2	1	2	2	3	3	3	3	1
2	4	3	3	3	5	4	2	2	2	3	2	4	3	3	3	4	4	5	5	6	6	5	4	1
2	4	3	3	3	4	4	2	4	3	5	4	5	4	2	4	4	5	5	6	8	7	6	5	3
1	4	4	5	5	4	4	1	3	2	3	2	4	4	3	3	4	6	7	7	7	6	5	5	3
2	4	3	4	5	4	5	2	4	2	3	2	4	4	4	4	5	6	7	6	6	4	4	4	3
3	4	3	5	6	4	4	2	3	1	1	1	3	3	4	4	6	7	7	5	4	4	4	4	2
3	4	2	3	3	2	3	3	4	2	1	2	3	3	4	5	7	7	6	5	5	6	5	4	2
3	4	3	4	3	2	3	4	3	1	0	3	4	5	4	5	5	5	3	4	4	7	6	7	4
3	5	4	4	2	2	3	5	4	3	2	3	5	5	5	5	6	6	4	4	3	5	5	6	4
3	4	3	3	2	2	4	5	4	2	3	4	6	5	5	5	6	6	5	4	2	2	3	5	4
3	5	3	4	2	3	3	5	5	5	6	5	5	2	4	4	7	6	6	3	2	1	2	2	2
3	5	3	4	2	3	3	5	4	4	5	5	4	1	3	3	6	5	6	3	3	2	2	1	1
3	5	3	4	2	4	4	6	4	5	5	5	4	2	3	2	5	4	5	2	3	3	4	3	2
4	5	3	3	2	3	3	4	2	3	3	4	3	2	2	1	4	4	5	2	2	2	3	4	3
3	3	1	2	3	5	4	3	2	4	5	5	4	3	2	0	2	3	4	3	3	3	3	4	3
3	3	1	2	4	5	3	2	2	5	5	4	2	1	1	0	2	3	3	3	4	4	2	2	2
2	3	1	2	4	6	5	3	3	6	7	5	4	3	3	1	2	4	5	5	5	5	4	4	3
3	4	2	2	3	5	4	3	3	6	6	5	3	3	3	3	5	6	6	5	5	6	6	7	5
2	3	3	3	3	4	4	4	4	5	5	4	3	3	3	3	5	7	8	6	5	5	6	7	5
2	3	5	5	4	3	2	4	4	4	3	4	3	3	3	5	6	6	5	4	5	6	7	7	5
3	4	6	5	4	1	1	3	3	2	1	3	4	5	4	5	5	5	4	3	4	5	6	6	4
4	6	8	6	5	2	1	2	2	3	2	4	4	6	6	8	7	5	2	2	3	5	5	6	4
4	5	6	5	6	3	2	2	2	4	2	4	3	6	5	6	6	6	4	4	4	5	4	5	4
2	3	4	5	7	5	3	2	2	5	3	4	1	3	3	5	5	5	4	5	5	4	2	3	3
1	1	1	3	5	4	2	1	1	3	2	3	1	2	1	2	2	3	3	4	4	3	2	2	2

2	3	2	2	2	4	3	2	1	1	2	2	3	2	2	2	2	1	2	2	3	3	3	3	1
2	4	3	3	3	5	4	2	2	2	3	2	4	3	3	3	4	4	5	5	6	6	5	4	1
2	4	3	3	3	4	4	2	4	3	5	4	5	4	2	4	4	5	5	6	8	7	6	5	3
1	4	4	5	5	4	4	1	3	2	3	2	4	4	3	3	4	6	7	7	7	6	5	5	3
2	4	3	4	5	4	5	2	4	2	3	2	4	4	4	4	5	6	7	6	6	4	4	4	3
3	4	3	5	6	4	4	2	3	1	1	1	3	3	4	4	6	7	7	5	4	4	4	4	2
3	4	2	3	3	2	3	3	4	2	1	2	3	3	4	5	7	7	6	5	5	6	5	4	2
3	4	3	4	3	2	3	4	3	1	0	3	4	5	4	5	5	5	3	4	4	7	6	7	4
3	5	4	4	2	2	3	5	4	3	2	3	5	5	5	5	6	6	4	4	3	5	5	6	4
3	4	3	3	2	2	4	5	4	2	3	4	6	5	5	5	6	6	5	4	2	2	3	5	4
3	5	3	4	2	3	3	5	5	5	6	5	5	2	4	4	7	6	6	3	2	1	2	2	2
3	5	3	4	2	3	3	5	4	4	5	5	4	1	3	3	6	5	6	3	3	2	2	1	1
3	5	3	4	2	4	4	6	4	5	5	5	4	2	3	2	5	4	5	2	3	3	4	3	2
4	5	3	3	2	3	3	4	2	3	3	4	3	2	2	1	4	4	5	2	2	2	3	4	3
3	3	1	2	3	5	4	3	2	4	5	5	4	3	2	0	2	3	4	3	3	3	3	4	3
3	3	1	2	4	5	3	2	2	5	5	4	2	1	1	0	2	3	3	3	4	4	2	2	2
2	3	1	2	4	6	5	3	3	6	7	5	4	3	3	1	2	4	5	5	5	5	4	4	3
3	4	2	2	3	5	4	3	3	6	6	5	3	3	3	3	5	6	6	5	5	6	6	7	5
2	3	3	3	3	4	4	4	4	5	5	4	3	3	3	3	5	7	8	6	5	5	6	7	5
2	3	5	5	4	3	2	4	4	4	3	4	3	3	3	5	6	6	5	4	5	6	7	7	5
3	4	6	5	4	1	1	3	3	2	1	3	4	5	4	5	5	5	4	3	4	5	6	6	4
4	6	8	6	5	2	1	2	2	3	2	4	4	6	6	8	7	5	2	2	3	5	5	6	4
4	5	6	5	6	3	2	2	2	4	2	4	3	6	5	6	6	6	4	4	4	5	4	5	4
2	3	4	5	7	5	3	2	2	5	3	4	1	3	3	5	5	5	4	5	5	4	2	3	3
1	1	1	3	5	4	2	1	1	3	2	3	1	2	1	2	2	3	3	4	4	3	2	2	2

2	3	2	2	2	4	3	2	1	1	2	2	3	2	2	2	2	1	2	2	3	3	3	3	1
2	4	3	3	3	5	4	2	2	2	3	2	4	3	3	3	4	4	5	5	6	6	5	4	1
2	4	3	3	3	4	4	2	4	3	5	4	5	4	2	4	4	5	5	6	8	7	6	5	3
1	4	4	5	5	4	4	1	3	2	3	2	4	4	3	3	4	6	7	7	7	6	5	5	3
2	4	3	4	5	4	5	2	4	2	3	2	4	4	4	4	5	6	7	6	6	4	4	4	3
3	4	3	5	6	4	4	2	3	1	1	1	3	3	4	4	6	7	7	5	4	4	4	4	2
3	4	2	3	3	2	3	3	4	2	1	2	3	3	4	5	7	7	6	5	5	6	5	4	2
3	4	3	4	3	2	3	4	3	1	0	3	4	5	4	5	5	5	3	4	4	7	6	7	4
3	5	4	4	2	2	3	5	4	3	2	3	5	5	5	5	6	6	4	4	3	5	5	6	4
3	4	3	3	2	2	4	5	4	2	3	4	6	5	5	5	6	6	5	4	2	2	3	5	4
3	5	3	4	2	3	3	5	5	5	6	5	5	2	4	4	7	6	6	3	2	1	2	2	2
3	5	3	4	2	3	3	5	4	4	5	5	4	1	3	3	6	5	6	3	3	2	2	1	1
3	5	3	4	2	4	4	6	4	5	5	5	4	2	3	2	5	4	5	2	3	3	4	3	2
4	5	3	3	2	3	3	4	2	3	3	4	3	2	2	1	4	4	5	2	2	2	3	4	3
3	3	1	2	3	5	4	3	2	4	5	5	4	3	2	0	2	3	4	3	3	3	3	4	3
3	3	1	2	4	5	3	2	2	5	5	4	2	1	1	0	2	3	3	3	4	4	2	2	2
2	3	1	2	4	6	5	3	3	6	7	5	4	3	3	1	2	4	5	5	5	5	4	4	3
3	4	2	2	3	5	4	3	3	6	6	5	3	3	3	3	5	6	6	5	5	6	6	7	5
2	3	3	3	3	4	4	4	4	5	5	4	3	3	3	3	5	7	8	6	5	5	6	7	5
2	3	5	5	4	3	2	4	4	4	3	4	3	3	3	5	6	6	5	4	5	6	7	7	5
3	4	6	5	4	1	1	3	3	2	1	3	4	5	4	5	5	5	4	3	4	5	6	6	4
4	6	8	6	5	2	1	2	2	3	2	4	4	6	6	8	7	5	2	2	3	5	5	6	4
4	5	6	5	6	3	2	2	2	4	2	4	3	6	5	6	6	6	4	4	4	5	4	5	4
2	3	4	5	7	5	3	2	2	5	3	4	1	3	3	5	5	5	4	5	5	4	2	3	3
1	1	1	3	5	4	2	1	1	3	2	3	1	2	1	2	2	3	3	4	4	3	2	2	2