第3回 Asprova プログラミングコンテスト

第3回 Asprova プログラミングコンテストに参加しました。1位だったようです。最後にパラメータをいじって数回提出したら900万点ぐらい上がってしまった結果なので、3位から上のスコアに本質的な差を感じられず、その点は運がよかったなあという感想になりました(完全に偶然上がっただけというわけではありませんが)。

コンテスト中にやったことをまとめてみます。

問題

その前に問題の内容ですが、ここには雰囲気だけ書くことにします。もとの問題文と違うことを書いていたり情報が足りなかったりするので、正式かつ正確な内容は問題文を見てください。

あなたは工場を持っています。工場には複数の設備があり、それらを動かして複数の品目を生産しています。それぞれの品目は複数の工程を経て完成し、各品目の各工程に使用する設備は一意に決まっています。

工場ではたくさんの注文を受け、それらを受けた順に生産することになっています。注文にはそれぞれ納期が決まっており、それに遅れるとお金はもらえません。

ある日あなたは、注文を受けすぎて現状の設備では納期までに全く対応しきれないことに気付きました。そこであなたは資金を投資することで設備の生産性を上げることにしました。

適切に投資して利益を最大化する計画を立ててください。なお必ずしもすべての注文を納期に間に合わせる必要はありません。

方針が固まるまで

各工程を割り付ける設備は一意に決まっているし、同じ設備に割り付ける工程の順番も決められているので、設備の能力の上げ方を固定すればなるべく早い時間から貪欲に割り付けていくのが最適スケジュールになる(と問題文にも書いてある)。そうすると、考えるべきことはどの設備の能力をどの日に上げるかだけなる。

よくわからないけどとりあえず山登り、と思ったがこれはよくなかった。ほとんどスコアが上がらない。考えてみればそれはそうで、近傍のとり方が下手だと初期状態から遷移できる先でスコアが上がる状態がほとんど存在しない。そしてこの段階でとった近傍はだいたい下手である。

こういう状況でまず最初に思いつくのはスコア関数をいじることで、今の場合だと納期遅れを一時的に許容したい感じなので、納期遅れが発生したときは遅れの量に応じて利益に一定の割り引き係数をかけてみた(元々のスコア関数は少しでも遅れると 0 になる)。ところが、これは係数が大きすぎると遅れをあまり気にしない計画になって能力追加するメリットが少なくなるし、係数が小さすぎると能力追加コストが見返りに勝ってしまうので何もしないほうが得になってしまう。ということであまりいい結果が出ない。

ビジュアライザを見た感じだと、広範囲に対してまとめて能力追加しないとあまり効果がなさそう。追加コストの評価も一時的に小さくする、とか考えられないわけではないけど、こういうことを積み重ねるのもちょっとアドホック感が出てくるきらいがある。

そこで別の方針を考えた。能力をちょうどいいところまで上げる問題ではなく、最初に多すぎるくらい追加しておいてちょうどいいところまで下げる問題だとみなす。これがわりとよかったようだ。何しない状態だとスコアの上がる遷移があまり見つからないのに対して、上げすぎの状態は不要なところをちょっと下げるだけでスコアが上がるので、スコアが更新される遷移がたくさんあるのだ。

まず最初に、コストパフォーマンスのよさそうな(性格には、単位資金あたりの短縮稼動時間が最大の)設備を選んでちょっとずつ増やすループを回し、遅れがなくなったらそこから各設備と各日付けについて減らしても大丈夫ならちょっと減らすループを回して解を作ってみた。練習用データだといい具合になったと思ったけど、提出してみたら一部のテストケースで大幅にマイナスになって-10億点が出た。負の点数でも一応 AC は出るようだ(最高点が負だと新システムの順位表では0点になり、旧システムだと負の点数になって未提出よりも下位に表示される?)。

スコアは負になったけどこれを初期解にして山登りでもするか、という気持ちになって試してみたらすぐ局所最適解に到達したので、焼き鈍すことにした。近傍のとり方は設備と連続する何日かをランダムに選んで、その日の能力を 8〜現在値+8 の範囲でランダムに変更する、とかだったと思う。ここまでのことを実装して 957875042 が出ている。

焼き鈍しができたので、高速化する意味があるかどうかを見るために回す時間を伸ばしたところ、まだけっこう上がるようだった。

このあたりまででだいたいの方針は出ていて、あとはずっと焼き鈍しを頑張る感じになった。後半は高速化にかなりの時間を費やしていた気がする。細かいところでは 1.1 乗をメモ化するとか温度の計算を 256 回おきにするとかメモリアクセスの局所性を頑張って上げるとかいろいろよくある最適化をした。高速化はループ回数という形でわかりやすく結果が出るので楽しいのかもしれない。最終的に焼き鈍しのループは practice05 で322万回ぐらい回っている。

あと、そろそろテストケースを手元で作るべきなんだろうなと思い始めたけど、結局作らなかった。

焼き鈍し周り

近傍のとり方

日をランダムに選ぶと全く稼動してない日を選んでしまう可能性がそこそこあって無駄になってそうだと思ったので、日を直接選ぶのではなく工程を選んでそれが割り当たっている日を取得することにした(各工程の開始と終了の時刻を覚えておく)。

practice04 の局所最適解(後述)を見て、同じ品目の連続する二工程に対して変更をするのも近傍に追加してみた。効果がどれくらいあるか実ははっきりしなかったのだが、平均的にはスコアはよくなってるようにも見えるので採用することにした。

近傍は時間経過に合わせてだんだん狭める。具体的には、x = (経過時間)/(制限時間) として ±(2+log_2(x)) の幅で能力追加量を変化させる。ただし一番スコアがよかった提出では x < 0.15 のときはもっと大きく動かしている。ちなみにこの式は x = 1, 0.5, 0.25 でそれぞれ 2, 4, 6 ぐらいになってくれるとまあまあよさそうな感じが(実験的に)したのでそれっぽいのにした。別にきれいな式にする必要はないし折れ線とか階段でもよかったと思う。

最後の 0.1% の時間は、追加能力を減らす方向だけを試す山登りにする(たぶん運がいいとこれで数十万点上がる)。

差分計算

一部の追加計画を変更すると何も割り当たっていない時間に能力追加してしまって無駄が発生することがあり、それが嫌だと思ってスコア計算中にいろいろ書き変える仕様にしてしまったため(そうする必要があると思ったのかどうかよく覚えていない)、状態を変更してスコアを計算し直す際に元の状態を丸ごとコピーしなければならなくなった。ここが最初ボトルネックになっていた。結局、状態とその一時保存先へのポインタを持っておいて、保存するときは memcpy して元に戻したいときはポインタの swap をするだけで済むようにした。

利益の差分計算は、最初に工程間の依存関係を表すグラフを作り(各工程から同じ設備の次工程と同じ注文の次工程へ辺を張る)、変更の影響を直接受ける工程から到達可能な工程だけ計算し直すようにするとやや速くなる。

能力追加コストの計算は、能力追加量を変更するときに毎回更新すればよい。

スコア関数

高速化と近傍のとり方以外でスコアに効いたかなと思うのは、評価関数を少しいじったこと。最初のうちは納期遅れをやや許容する形の式にして、時間経過とともに許容度を下げて、最終的には本来のスコア関数と一致するようにする。初期に試してうまくいかなかったやつだけど、能力を上げすぎて下げる方針に変えたのでまた状況が違う。

局所最適解にはまる様子

なかなかスコアが安定しないなあと思いながらずっとやっていたのだけど、結局最後まであまり安定しなかった気がする。特に practice04 のスコアがかなりぶれるので(平均で2600万ぐらい出るが悪いときは2500万ぐらい)、解をたくさん作って様子を見てみたところ、局所解はまりポイントが(少なくとも)二つあり、それぞれ50万点ぐらいずつもっていかれる様子が観察できた。

ひとつはオーダー番号4の注文で、これは主に設備7, 8を強化すると間に合うのだが、7により比重をおくのが正解で、8のほうに比重をおきつつ10も少し強化するというのがはまる局所解のようだ。近傍をうまくとれればたぶん脱出できるのだが、複数設備にまたがった変更が必要なのが難しいところか。この状況が存在することに気付いたので二つの連続する工程の位置に同時に変更を加えるという近傍を追加してみた。

もうひとつは終盤の品目3に対する注文群で、間に合わせようとするとコストが大きくなって割に合わないので切り捨てるのが最適のようなのだが、それなりの数をまとめて切り捨ててやる必要があるためか、頑張って間に合わせるスケジュールが局所最適になって抜け出しにくいらしい。こっちは対策できていないけど、あまりこれにはまるケースは多くはなかったような気がする。

焼き鈍しのスコア変化の様子

温度を変えながらそれぞれ10回ずつ実行したときのスコア上昇の様子をグラフにしてみた。

f:id:lkozima:20190511010142p:plainf:id:lkozima:20190511010211p:plainf:id:lkozima:20190511010136p:plain
温度ごとのスコア変化の様子。右のほうが温度が高い

温度が高めのときにはあるところまで下がりぎみで、途中から上昇に転じる様子が見えたので、それなら上昇を始めるあたりの温度を初期温度にすればいいのでは、と思いそうしてみたところ具合がよさそうだった。そうすればいいと思ったのは理論的な根拠があったわけではなくてなんとなくなんだけど、これって一般的にうまくいく方法なんだろうか。

その他

終了後に、初期解をそのまま返すコードを投げてみたら 669780418 が出た。2000ms ぐらいかかっていて、わりと想定外だった。手元で真面目にテストケースを作ってみていたら事前にわかっていたのかもしれない。

でも今回は前回に比べるとある程度ちゃんとデータを見ながら正しい方向に改善を考えられた気がする。テストケース作らなかったけど。

Marathon Match 107

ラソンマッチ始めてみました。問題文

問題概要

縦 H 横 W のグリッドがあります。各マスに 9 以下の非負整数が書いてあります。マスのいくつかを塗って、塗られたマスに書いてある数の合計を大きくしてください。ただし 縦 h 横 w のどの長方形も塗られたマスを Kmin 個以上 Kmax 個以下含むようにしなければなりません。

とりあえず正の点をとる方法

h×w の subgrid の中から Kmin 以上 Kmax 以下の個数を塗る任意のパターンを作る。左上から順番に塗るとか、何でもいい。それを H×W に周期的に拡張すれば制約を満たす塗り方ができて正の点が得られる。

その後のいろいろ

初期解が作れるようになったので、とりあえずこれを動かして山登りをしてみた。何箇所か塗る/塗らないの状態を flip したものを近傍にした(最終的には二箇所だけにした)。提出してみたら 930070.87。この時点ではまだいろいろ効率のよくない部分があったが、高速化はいつでもできそうなので後回しにして、方針を改善できないか考え直すことにした。

初期解を作るときに点数が高くて端に近いやつを優先的に塗る(マスの点数をそのマスを含む subgrid の個数で割った値を優先度にする)ようにしたら、それだけで最初の山登り解と同じくらいのスコアが出ることに気付く。ただ、これだけだと問題があって、たまに制約を満たすものが見つからないことがある。中央付近があまり塗られてなくて周辺がたくさん塗られすぎている状態になってしまった場合、追加で塗るマスをどう選んで制約を満たせなくなる。例えば H = W = 4, h = w = 2, Kmin = 1, Kmax = 2 で以下のような配置になった場合とか。

 ****
 ....
 ....
 ****

そこで、うまくいかなかったらちょっとランダム性を入れて何度か作り直すことにした。それでもだめだった場合は周期的な解の中で一番スコアが高いのを作る。これは以下のようにすれば作れる。まず 0 <= i < h, 0 <= j < w の範囲にある各 (i, j) について (x, y) = (i, j) mod (h, w) になるような (x, y) のマスの値の総和を求める。和が大きいほうから Kmax 個まで (i, j) を塗って、それを周期的に全体に拡張する。

5×5 の場合は全探索で間に合うのでそうした。もう少し大きいケースも枝をちゃんと刈れば厳密解が出せたかもしれないがやらなかった。このあたりで 958728.41

高速化と焼き鈍しをやって、細かいことをいろいろやって、焼き鈍しのパラメータ調整をして 989241.86

ところで貪欲解の構成に成功した場合は焼き鈍しをやってもほとんどスコアが上がっていないケースが多いようだけどこれはどう解釈したらいいんだろうなあとずっと思っていた。

学びと反省

  • 時間を測るのに chrono を使ったら TLE が出た。調べてみると、rdtsc を使うとよいという情報が出てきたので真似してみたらまあまあいい感じになった。50ms ぐらい余裕をみておかないと TLE するときがあるみたい?
  • 盤面が制約を満たしているかどうかの判定をたくさんやってるので segment tree を使ったら速くなるだろうと思ってやってみたら遅くなったので悲しみを感じた。(落ち着いて考えてみると調べる矩形が O(HW) 個あるので全部で O(HWlogHlogW) になっていて遅いに決まっているのでは)

Asprova コンテストに参加した

年末年始は Asprova Programming Contest をやってました。なんか問題文が軽い気持ちで読むには重すぎる感じなのですが,短くまとめると

  • いくつかの設備を備えた工場があり,そこでいろんな製品を生産しています
  • 各製品はいくつかの工程を経て完成します
  • 各製品の各工程は処理できる設備とできない設備があります
  • いま,どの製品がいつまでにいくつほしいという注文が200件来ました
  • いい感じの稼動スケジュールを求めてください

という感じです。実際にはもっと条件があります。

やったこと(だいたい時系列順)

数字はその時点でのスコアです。ちなみに配布されてるサンプルコードをそのまま提出すると408307056327点です(試してないけどたぶん)。

  • 問題文が何度読んでもよくわからないので,とりあえず配布されてたサンプルコードを読んだり入力を読むコードを自分で書いたりしてみる。
  • サンプルコードを読んで出力を見てみたら,常に最初に見つかった設備に割り付けるようになっており,何も割り付けられない設備がたくさんあることに気付く。ちゃんと分散させたらそれだけでもスコアが上がるはずなので,それを実装してみることにする。
  • 割り付け可能なすべての設備について,すでに割り付けられている一番遅い工程の終了時刻が一番早いところに割り付ける貪欲を書いた(つまりサンプルをすべての設備を見るように変えた)。ついでに注文を納期の早い順にソートしておくことで早く着手したほうがよさそうなものを早い時刻に割り付けるようにした。447088233657
  • 出力を眺めていたら,一部の設備への割り付け時刻がどんどん後ろに遅延していく様子。たぶん単純に後ろに置いていくだけだと前のほうにたくさん隙間ができていそうで,そこに割り付けられる工程がありそう。そこで別の工程の間に入れられそうな隙間があったらそこに詰め込むようにした。段取り時間を計算し直す必要があったりしてこのあたりからだんだん実装がややこしくなる。448737404998
  • 基本的に納期に対して遅れがちになるみたいなので,着手遅延を狙うよりは早めに着手して納期の遅れを減らすほうがスコアはよくなる傾向がありそうだと予想し,なるべく早い時間に工程を詰め込む方針考えることにする。ただし納期遅れが発生していない場合には,最後に納期の範囲内で全体を後ろへずらして着手遅延ボーナスを稼ぎたい。工程が複数ある入力でこれをやるのが面倒だったので,input 1(マスターデータ 1 から作られたやつ)の場合だけ最後にその処理を入れることにした(たぶんスコアへの貢献はかなり小さい)。
  • 今までは開始時刻が最小のところに割り付けていたのだが,段取りのペナルティの係数をちょっといじったものを開始時刻に足した値が最小になるよう割り付けてみた。係数はいろいろ試して一番よかったのを採用する(この評価値,落ち着いて考えてみると質の異なるものを足しているのでどういう意味があるのか疑問だったのだが,最後までこれよりもよいものは見つけられなかったのでしかたなく使い続けた)。449614829656
  • 簡単なビジュアライザを書いた。input 3 では最初の1/3を除いてずっと遅れが発生しているような状態になるらしいことがわかった。この入力だと最初のほうの工程を担っている設備が1--6で他の設備と担当できる工程が disjoint なのだけど,1--6がほとんど休みなしで稼動してこれだけ遅れているので,大幅な改善は難しそうに見える。この時点で上位陣は 4498 に乗っていたので,どこに改良の余地があるのだろうかと不思議に思い始める。ただ,納期の比較的早いものがずっと後回しにされているような印象はあったので,その部分をうまく解決する方法を考えられればよかったのかもしれない。
  • 一つの注文は複数の工程からなっているのだが,注文ごとではなく複数の工程を別々に処理したほうが効率よく割り付けられるのではないか,と思い試してみる。優先度つきキューに処理すべき工程を入れて順に割り付けていくのだが,優先度のつけ方を納期の早い順にすると前よりも少しよくなった。449640963038
  • さらに納期から所要時間の見積りを引いたものを優先度にしてみたり評価関数のパラメータの動かし方を変えてみたりパラメータをひとつ追加したりした(何も割り付けれていない隙間ができることに対してペナルティを課す項を足した)。449730476138
  • 貪欲でこれ以上根本的な改善ができる気がしないので,乱数に頼るかという気持ちになり、ここまでで得られた解から始めて山登りをした。既に割り付けられている注文のいくつかをランダムに除去して割り付け直して,スコアが上がったら採用,を繰り返す。割り付けの制約が複雑なので,制約に違反しないように割り付け直すのがかなり大変で苦労する。assert をたくさん書いた。あと細かいことをいろいろ試して 449800668222
    • ちなみにこのスコアが出たとき提出したコードにはバグがあって手元で動かすとときどき segmentation fault が出るのだけど,このバグ(スコアにはさほど影響するとは思えない)を直すとスコアが1000万点ぐらい落ちるので困惑した。手元では直してもそんなに悪くなっているようには見えず,結局何が起きていたのかわからないままである。もしかしたら,たまたまジャッジサーバで実行したときにいい乱数が出るコードだったのかもしれないが,そんなことあるだろうか?

できるかもしれないと思ったけど結局やらなかったこと

  • 高速化。山登りが2秒では全然収束してなかったのでできたらよかったのだけど。割り付けられそうな隙間を探すときに,線形に探索するのではなく区間木のようなデータ構造を持っておけば効率よく見つけられたかもしれない。線形探索でも平均すると10〜20回ぐらいループを回せば割り付け箇所が見つかっているようだったし,更新のオーバーヘッドも大きくなるはずだし,本当に割り付けられるかどうかは前後の工程とか見ないと判断できないしで,全体として本当に改善するのか確信がなったのでやる気が出なかったところはある。
  • ビームサーチ。使えるタイプの問題なのかなと思ったけど,どういう探索木を考えればいいのかよくわからなかった。そもそもビームサーチを使ったことがないので勝手がわからないというのもある。
  • この問題自体はかなり複雑だけど,適当に問題を単純化したら厳密解が出て,そこから始めてうまく補正したらいい解が作れないかとちらっと思った。そこまでアルゴリズム力がないので詳しく検討せず。

感想

終わってみると,貪欲法で4497が出たあたりから先は問題の特性をうまく捉えて利用するという意味で本質的にいい改善はできなかったように思う。改善できるとしたら納期遅れペナルティを減らすことなんだろうと思ったけどどうすればそれがうまくできるのかわからなかった。

後半でいろいろ細かいことを試せたのはそれなりに勉強になったし悪くはなかったと思う一方で,根拠があるのかないのかわからない思いつきをいろいろ試していただけで意味のある洞察が新しく得られたということはなかったようにも感じる。後者の意味では後半の試行錯誤を振り返ると何をやっていたんだろうなという気もする。

二分探索

二分探索は,値の探索に使うという認識が一般的だと思いますが,実際にはもっと守備範囲が広いです。ということでなるべく一般的な形で書くとどうなるのかを整理した話を書いてみることにしました。why3 で書いて正しいことを検証してみたらけっこうすっきり理解できた気がします。

int 上の述語 f があるとします。f を真にする a と偽にする b がひとつずつわかっていて,a < b であるとしましょう(requires のところ)。そうすると,f x だが f (x+1) ではない x を求めることができる(ensures のところ,result が戻り値),というのが関数の仕様です。

中身はまあ普通の二分探索なのですが,上限 ub は必ず f を満たさないように,下限 lb は必ず f を満たすようにするとわかりやすくて不変条件が簡単に書けます。

ソートされた列 A から値 v のインデックスを求めるよくある二分探索は,f x として A[x] <= v をとった場合になります。

ちなみにコード中では f に特に条件を付けていませんが,真偽が頻繁に入れ替わるような述語に対して使うケースはあまり見た記憶がありません。ある値を境にしてそれより下で真,上で偽(あるいはその逆)になる述語を考えて,その境界の値を求めるのに使うことが多いです。また a, b は自分で与えないといけませんが,たいてい自明な上界と下界があります。

上のコードが正しいことのチェックは,例えば手元の環境では why3 prove -P cvc4 binsearch.mlw を実行すると

binsearch.mlw BinSearch WP_parameter binsearch : Valid (0.03.s)

と出て,正しいことが自動証明できました。int の長さが限られている言語だと足し算がオーバーフローした場合に正しく動かない気がしますが,そこまではチェックしてないようです。

今回は整数上でやりましたが,実数上でも適当なεに対して f x だが f (x+ε) ではない x を求める,のような形で書けるでしょう(たぶん正確にはちょっと違う)。

最短距離と自由豊穣圏

距離空間は豊穣圏であるという Lawvere の論文(http://www.tac.mta.ca/tac/reprints/articles/1/tr1abs.html)がありますが,それに似た感じで重みつきグラフから自由に生成された豊穣圏を考えると hom-object が最短経路になるらしいことに気付いて面白いと思ったので書いてみます.(よく見たら Lawvere の論文にも書いてありましたが,まあいいでしょう.)

R の monoidal structure

R は非負実数の集合に∞を追加した集合とする.これに普通と逆の順序を入れて圏とみなす.実数の普通の足し算はこの圏でのモノイド積になる(0が単位元).また,この圏における直和は inf で与えられる.特に始対象は∞である.モノイド積はすべての直和に対して分配律を満たす.

V-graph と重みつきグラフ

V を可算直和をもつモノイド圏で,分配律  A \otimes \coprod_i B_i \simeq \coprod_i (A \otimes B_i) \left(\coprod_i A_i\right) \otimes B \simeq \coprod_i (A_i \otimes B) が成り立つものとする(Lawvere の論文には詳しい条件が書かれていないけど,後でこれくらい必要になる気がするので最初から仮定しておく).

V-グラフ  \mathcal A は,集合  \mathcal A_0 と,各  u, v \in \mathcal A_0 に対する対象の割り当て  \mathcal A(u, v) \in V からなる.

特に V = R のときを考えると,V-グラフとは非負の重みつき単純有向グラフのことである.逆に非負の重みつき単純有向グラフ G が与えられたとき,辺がないところには重み ∞ の辺があると考えることで G を R-グラフとみなせる.

free V-category と最短距離

V-グラフ \mathcal Aに対して,V-category  \widetilde{\mathcal A}を対象が \mathcal A_0の要素全体で,hom-objectが  \widetilde{\mathcal A}(u, v) = \coprod_{n \ge 0} \mathcal A^n(u, v) となるようなものとする(0は始対象,1はモノイド積の単位元).ただし
 \mathcal A^0(u, v) =
  \begin{cases}
    0 & u \ne v \\ 1 & u = v,
  \end{cases}
  \quad
  \mathcal A^{n+1}(u, v) =
  \coprod_{w \in \mathcal A_0} (\mathcal A^n(w, v)\otimes \mathcal A(u, w))
とする.これは V-category になる(たぶん).Lawvere によると,この構成は V-Cat から V-graph への忘却関手の左随伴になっているそうだ(いかにもなってそうである).

重みつきグラフ G を R グラフとみなして上記の構成を行うと, \widetilde G(u, v)は u から v への最短距離になる( \otimesが足し算で \coprodがinfであることを思い出せばわかる).

Z/nZ の分解周りのメモ

 n = q_1 \dots q_k, q_i = p_i^{e_i} と書いたとき中国剰余定理として知られる同型  \mathbb Z / p^n\mathbb Z \simeq \prod_i \mathbb Z / q_i \mathbb Z の具体的な与え方が毎回思い出せないのでどこかに書いておきたいと思ってメモ。ついでに単数群の生成元も。

同型

 \mathbb Z / n \mathbb Z  \to \prod_i \mathbb Z / q_i \mathbb Z は単純にそれぞれの成分に射影すればよい。逆が少しややこしくて,
 \prod_i \mathbb Z / q_i \mathbb Z \ni (x_1, \dots, x_k) \mapsto \sum_i \left(\frac{n}{q_i} \mod q_i\right)^{-1} \cdot \frac{nx_i}{q_i} \mod n \in \mathbb Z / n \mathbb Z
となる。mod q_i での逆元は互除法を使うか  \phi(q_i) - 1 乗すると求まる。

単数群

p が奇素数のときは  (\mathbb Z / p^n \mathbb Z)^{\times}巡回群で,r を mod p での原始根とすると  (p+1)r で生成される。他にも生成元はたくさんある。

p = 2 のときは n > 2 なら巡回群ではなく, (\mathbb Z / 2^n \mathbb Z)^{\times} \simeq \mathbb Z / 2 \mathbb Z \times \mathbb Z / 2^{n-2} \mathbb Z となる。右から左への同型は  (x \mod 2, y \mod 2^{n-2}) \mapsto (-1)^x \cdot 3^y で与えられる(左辺は乗法群,右辺は加法群なので混乱しないよう注意)。他の同型もたぶんたくさんある。

Graphillion で半順序の個数を数えてみた話

Graphillion で遊んでみた。

やったこと

有限集合上の半順序の個数を数えてみた(cf. OEIS A001035)。これを効率的に求めるのは未解決問題らしい。ちなみに有限集合上の位相の個数を数える問題とだいたい同じ。

素数 n を固定して,Graphillion を使って [1..n] 上の半順序をわりと愚直に列挙した。

実行時間は n=6 で 0.3s, n=7 で 8.6s くらい。n=8 になると三時間ぐらい走らせても終わらなかったのでやめた。ちなみに OEIS を見る限り少なくとも n=18 までは既知らしいので,こういうアプローチで新しい値を見つけるのは無理がありそう。

やりかた

a <= b だったら a から b へ辺を伸ばすことにすると,半順序は有向グラフとみなせる。そこで n 点からなる完全グラフの部分グラフで,半順序になっているようなものを列挙すればよい。Graphillion は基本的な集合演算をサポートしているので,基本的には順序の条件をそのまま書くだけでよい。

ただし,Graphillion は無向グラフを扱うライブラリなのに対して,ここでは有向グラフを扱いたい。そこで,source と target を別の点にして有向グラフをむりやり無向グラフで表すことにした。

一応ちゃんと書くと以下のとおり。有向グラフ  (G, E) に対して無向グラフ  (G+G', E') (a, b') \in E' \iff (a, b) \in E となるように作る。ただし  G' G のコピーで, (-)'\colon G \to G'全単射とする。

実際のコードでは,G の点を正の整数,G' の点を負の整数で表した。あと半順序は反射的反対称的推移的関係だけど,対角線部分はあっても情報が増えないので,非対称的推移的関係(strict order; 日本語だと強順序?)にした。こっちのほうが辺の数が減って少し効率がよくなるようだ。

from graphillion import GraphSet

def posets(n):
    univ = [(i, j) for i in range(1, n+1) for j in range(-n, 0)]
    GraphSet.set_universe(univ)
    x = GraphSet({})

    # asymmetry
    for i in range(1, n+1):
        x = x.excluding((i, -i))

    # transitivity
    for i in range(1, n+1):
        for j in range(1, n+1):
            if i != j:
                for k in range(1, n+1):
                    if j != k:
                        y = x.including([(i, -j), (j, -k)]).excluding((i, -k))
                        x = x - y

    return x

import sys
print(len(posets(int(sys.argv[1]))))

実行してみた。

$ time python3 posets.py 6
130023

real	0m0.328s
user	0m0.293s
sys	0m0.028s

$ time python3 posets.py 7
6129859

real	0m8.648s
user	0m8.251s
sys	0m0.324s

出てきた数字は合ってる。

感想とか

ZDD を使うとコンパクトに表現できるかと思ったが,期待したほどコンパクトにはならなかったようだ。ZDD で効率的に表現できるのは基本的には疎な組合せ集合(正確にはちょっと違うと思う)なのだけど,考えてみれば順序全体ってそれほど疎ではないかもしれない。

強引に無向グラフにしてるところにちょっと無駄がある感じはして,Graphillion を使わないで ZDD を自分で操作すればもうちょっと効率よく列挙できそうには思う。

あとは,順序の表現を単に strict order にしたけど,推移律から推論できる辺を除くなどすればもっと辺の数を節約することはできるかもしれない。あるいは,実は BDD のほうがいいとか,(Z)SDD を使うといいとかはあるかもしれない。