DivineJK’s diary

Z cjikf c evf'

AtCoder青になったので水色から青までの一年半とついでに今まで書きそびれた黒から水までも全部書く

自己紹介

DivineJKです。地元の友達からは苗字で呼ばれています。

10月2日に、AtCoder青になりました。数学が強いタイプです。

音楽は主に洋楽メタルが好きです。最近はK-POPも聴き始めました。

あとふわふわツインテールとか好きです。なんというか、ツインテールに限らずハーフアップとかポニーテールとかとりあえずなんでも好きです。嘘かもしれません。

ゆめロリとか甘ロリとかも好きです。というか、原宿系の趣味があります。その影響でメイクを勉強中です。

あと、たまにイラストを描いたりしてます。最近は女子の制服を描くのとかマイブームです。正直イラストの方の精進もしたいんですよね。

このご時世で最近はなかなかできないのですが、地下アイドルのヲタをやったりもしてます。推しのアイドルは、Twitter（@realDivineJK）できいてくれたら教えます。

趣味の話はさておき、軽く自分の経歴を紹介します。

大学では物理学を専攻し、もう一回遊んでから卒業しました。ママ、パパ、ごめんな。まじで朝起きれないし電車乗ると吐き気する。

で、そのもう一回遊ぶときに、単位数確保のためにフラッと別の学部のプログラミングの授業でPythonに触れたら、見事にハマってしまいました。ハマりすぎてBCH公式の各項の係数計算や、スピン合成の固有状態の計算なんかも自主的に組んでました。そんなわけでPythonで遊びながらTwitterサーフィンしていると、顔の前に手を置いて何かを考えているアイコンのとある社長のツイートで、競技プログラミングやAtCoderの存在を知りました。

僕は今まで数学オリンピックとか情報オリンピックとかに無縁の人生だったので、競技学問ってちょっとハードルが高いなーとか、凡人の自分には無関係かなーと思っていました。僕は科学オリンピックや開発や研究ですごい成果をあげたとか、中学高校がとんでもない進学校だったとかそういうわけでもなかったので。

しかし、実際にAtCoderの問題を開いてみると、意外に解ける！となって、それからどんどん問題を解く楽しみを味わっていきました。多分、僕の探していたものはこれだったんじゃないか、僕の適正はここなんじゃないかとか思えてきて、もう少し早く出会いたかったなって思いました。

とりあえず自己紹介はこれで以上です。他何か聞きたいことがあったら僕のTwitter（@realDivineJK）のリプとかDMにでも質問してきてください。答えたくない質問はスルーしますけど。

びふぉーあふたー

これが、

f:id:DivineJK:20211003105920p:plain

こう。

f:id:DivineJK:20211003110001j:plain

ジャスト+400！

およそ1年半かかりました。Ratedコンテスト初参加からなんと2年1ヶ月（！）

水になるまでのハイライト

1997/08/27 誕生

生まれました。

やったこと：

  • 生まれる

2019/09/01 黒から灰になる（ABC139）

成長して、その過程で世の中の色々な理不尽や矛盾と闘いながら、そんな中でも家族や友達と笑い合える束の間の幸せを噛み締める日々を過ごしているときにAtCoderに出会い、黒から灰になりました。

やったこと：

  • AtCoderを知ってコンテストに出る

2019/11/04 灰から茶になる（AGC040）

AGCに参加したら茶色になりました。この当時は灰色もAGC ratedだったのだ。

やったこと：

  • コンテストに出る

2020/01/11 茶から緑になる（第6回ドワンゴからの挑戦状 予選）

スマホコーディングで青Perfとって緑！w 8級から7級すっ飛ばして6級！w

バイト帰りにスマホから参加。A問題でRE出たときは初めての冷えを覚悟しました。スマホの充電残り7%でB問題が通ったときは街中でスキップしてました。AtCoderはじめた当初は緑ぐらいを目指して終わろうかなと思ったのですが、競プロの知識は結構面白かったので続けることにしました。

やったこと：

  • mod逆元

2020/03/28 緑から水になる（ABC160）

3ヶ月で緑から水になりました。あとワニが死にました。

2020年になってから人が変わったように水Perfがとれるようになりました。ちなみに初水Perfより初青Perfのほうが先という謎の現象が起こっています。

やったこと：

  • Union-Find
  • BFS
  • FFT
  • 形式的冪級数を少し

水から青のふりかえり

ABC161〜ABC168（2020/04/04〜2020/05/17）

Rating 1203 -> 1311

水色になってもだいたい緑のときと同じようなPerf（1300くらい）をとっていて、これ1300になったらどうしようとか思いながらコンテストに参加していたら、いつのまにか1300台にのってました。

なんというか、このときは自分が青になる世界が想像できてなかったです。きっとレートは上がるけど青にはなれないだろうな、みたいな。でも青にはなりたかったです。

AGC044〜ABC171（2020/05/23〜2020/06/21）

Rating 1311 -> 1245

AGC044で0完太陽をとってしまい、およそ半年ぶりの茶Perfで激冷えしました。それ以来、NoSub戦略というものを選択肢にいれつつコンテストに参加していました。

また、ABC169で全完のチャンスだったのですが、自分の実装力の無さというか、実力不足で結局冷えになってしまったのが悔しかったです。

ABC172〜AGC048（2020/06/27〜2020/10/18）

Rating 1245 -> 1349

ABC172で初めてtourist出しというものをやってみました。これが大成功（Perf 1988）して以来、先のNoSub戦略と合わせて、企業コン以外のABCは「Eが解けたら出す、解けなかったら出さない」と決めてコンテストに参加するようになりました。

ARC106〜ABC186（2020/10/24〜2020/12/19）

Rating 1349 -> 1494

ARCで2回黄Perfとって一気にHighest1500台になりました（は？）。ただ、ABC181で緑diffのEが通せなくて、ここで自分の実装力の無さやあいまいさと直面することになりました。

ABC186以降はNoSubを繰り返すなどして空白の3ヶ月になりました。

ARC112〜ABC199（2021/02/13〜2021/04/24）

Rating 1494 -> 1477

復活します。もうこの時期からNoSub戦略に飽き出したというか、レートを落とす覚悟で行くか、と考え出しました。

ここから3回青以上のPerfが取れて、Highestが1576になったときは青見えた！やばい！青見えた！と思いました。

しかし、ABC198とABC199で激冷えして戦意喪失しました。このときは、もう自分は競技者としての適性がないんじゃないかとも思えて、引退も考えました。また3ヶ月近くコンテスト不参加の日々が続きました。

ちなみに、3月に大学を卒業して、4月頭からプログラミング系のアルバイトを始めました。

ABC212〜ABC221（2021/07/31〜2021/10/02）

Rating 1477 -> 1603

やっぱり問題が解けない、解けない問題が多いのは嫌だったので、人々の影響もあってABCバチャをやることにしました。

8問ABCになるので出てみるか、という軽い気持ちで出たらABCDGの5完で激温まりました。嬉しかったので、またコンテスト参加をはじめることにしました。

レートが1579になったとき、体調不良で参加できなかったABC219のバチャに参加したら、

本番なら黄Perfの成績でした。せつない。

このときから僕は「実質青」を名乗っていました。

そして、レート1563で迎えたABC221。

僕はこの回で青になれるとは思ってなくて、まああったまったらいいかな、青Perf取れたらいいな、ぐらいの気持ちでした。

そしたら、ABCDまで順調に進み、E問題の解法もすぐに浮かびました。BITに数を入れる処理でmodとり忘れてTLEでペナしたけど、ペナ含め47分でEが解けたのはわりと順調だなと思い、順位表を見たら...。

...ええ、250位とかになってました。

これ、このまま行けば黄Perf出るじゃないですか。AtCoder参加して初めてペナルティの5分に対する後悔を覚えました。

というか、ペナ消化した後の順位の下がり方がえらい鈍くて、え、これワンあるのでは？青コーダー誕生ちゃいますか？そう思い、ac-predictorを見に行ったら、1563->の矢印の先に青い数字が表示されていて、僕は落ち着きを失いました。

もうそのあとはPerfが1888以下にならないことを祈るのみ。でもやっぱり落ち着かないので、FとかHを解いて順位を解いて順位を少しでも上げようと頑張りました。結局解けずに終わったんですが。

コンテスト終了後も順位表は見れませんでした。Perfが1888になる順位が大体わかっていたので。必死にAnthrax聴いて気を紛らわせていました。

で、Twitterの人々がコンテスト結果をツイートした瞬間に自分のプロフィール画面を観にいきました。自分の名前が濃い青になってて思わず手を叩いてコロンビアになってました。

レートが1579で元Highestのときですら青になれる実感が湧かなかったので、まさか青が実現できるなんて驚きと嬉しさが隠し切れないですね。

あと心の余裕がなかったので画像が準備できなくてすみません。

水から青になるまでにやったこと

過去問解くやつ

AGC044からABC199までは停滞期と呼んでいいかもしれません。そして停滞期の時期は、コンテストをレート上げのためのものだと思っていたかもしれません。レートが上がれば嬉しいし、レートが下がったらカス以外の感想がなかったです。

しかしあるとき、「僕のヒーローアカデミア」1巻を読みました。デクくんが雄英合格に向けて積み上げを繰り返してるシーンは親指を立てながら涙なしにはみられませんでした。

そして、このシーンで僕はあることを思いました。

「レートは一時的に上下するかもしれないけど、解いた問題数は下がることはない。そして、解いた問題数が多いほどレートが上がる確率が高くなるから、過去問を積み上げない手はない！」

というわけで、一旦レートはおいといて、解いた問題数を積み上げていくことにしました。

特に、ABCではバンバン解説を開けて知らない知識を勉強していこうという気持ちで臨みました。僕は今まで特別な天才エピソードもなく、本当にただの凡人なので、解法が思いつくのを待つよりは解説を読んで理解するほうが効率がいいかなーとも思ったのです。

なんというか、パズルを解く楽しさより知識を得る喜びを求めたというか。知らない知識に触れたり知識を積み上げていく日々は刺激に溢れて充実してました。

そうして知識積み上げの日々を過ごすと、レートが落ちてもまた上げればいいしな、みたいに思えて冷えに対する抵抗が減ったと思います。

そういう充実した日々を過ごすと、自然と人は笑顔が増えていくものです。

笑顔で明るく、穏やかな気持ちで過ごしている人とは関わりたいじゃないですか。そんなわけで、ゆったりとしたそんな気持ちで過ごしていたらおよそ3ヶ月で性別、国籍問わず友達がめっちゃ増えました。

緊急事態宣言が出てたのでリモートだったのですが、飲み会なんかもやったりしました。

そして、その中の一人の女性とお付き合いすることになりました。

出身は宮城で、去年東京にきたんだそうです。まだ東京慣れしてなくて、休みの日は原宿デートがメインですが東京探訪をしています。

彼女はとても感情豊かで、話しているととっても癒されます。この間も、コンテストで冷えたことを報告したら泣きながら僕の胸元に抱きついて慰めてくれました。

仕方ないので、僕が頭を撫でていると、気づいたら僕の胸元で寝ていました。まったく、どちらが慰めてるのか、って感じですね。

そんな感じで、青になった今も楽しく問題を解いています。今では彼女と結婚に向けた話を進めているところです。

データ構造・アルゴリズム・知識

水の時期にやったと思われる内容をほぼ全て紹介します。

  • Segment Tree
  • Lazy Segment Tree
  • Binary Indexed Tree
  • Slope Trick
  • Li Chao Tree
  • 重み付きUnion-Find
  • Euler Tour
  • 木の半径
  • 非再帰DFS
  • 形式的冪級数
  • グレイコード
  • Suffix Array（SA-IS）
  • Manacherのアルゴリズム
  • Run-Length Encoding
  • Z-algorithm
  • 偏角ソート
  • 円の共有点、包含関係
  • 全方位木DP
  • 凸包（O(N2)）
  • 2次元アフィン空間
  • SCC
  • LCA
  • 木の重心
  • ベルマンフォード法
  • ダイクストラ法
  • クラスカル法
  • TSP
  • bitDP
  • 牛ゲー
  • ワーシャルフロイド法
  • NTT
  • 中国剰余定理
  • Garnerのアルゴリズム
  • 高速ゼータ変換
  • 高速メビウス変換
  • 高速アダマール変換
  • 離散対数
  • レピュニット数
  • 約数系包除
  • FloorSum
  • and畳み込み
  • or畳み込み
  • xor畳み込み
  • スターリング数
  • カタラン数
  • 分割数
  • 多点評価
  • 多項式補間
  • LIS
  • 半開区間は便利
  • 桁DPをやるときは母国語でコメントを書いた方がいい
  • 中央値を求めたいときは、ある数以上のものを1、そうでないものを0として二分探索して調べる
  • 二次元配列は一次元に直した方が気持ち軽くなることがある
  • Union-Findではサイズが大きい集合に繋ぎ直すと計算量が落ちることがある
  • インドゾウの足の裏の面積からインドゾウの体積を求める方法（Sreekumar Function）
  • 日焼け止めは将来のためにどんな理由にせよ塗っておいた方がいい
  • コンドームを付けて致すと手は洗ってもベタベタだけどゴミが少なくなるしにおいも汚れも少ない
  • 半蔵門線は水天宮前駅で日比谷線の人形町駅と接続するというが、それは真っ赤な嘘
  • アヴ様の誕生日は8/27ではなく9/27
  • セグ木の二項演算は抽象度を上げるために外部から入れた方がいい
  • 睡眠は大事

イベント

  • 卒論執筆
  • 大学卒業
  • 大学院入試
  • プログラミング系アルバイト
  • 正社員登用

ちなみに、過去問のとこで述べた彼女の話は全て作り話です。

まあこんな下手くそな嘘で騙されるような純粋な方はさすがにこの世に存在しないとは思いますが、万が一騙されたという人がいれば、今度インチキ量子論セミナーに送り込んでやるので覚悟しておいてください。

ファッション

スカートパンツを買いました。結構良いんですけど、自転車に引っかかりそうで怖くて着れてないです。

あと化粧水と乳液でのケアを始めました。

アイドルヲタ

推しが卒業とグループ解散を経験し、現在活動が第3章です。

イラスト

カイリューさんにファンアートをプレゼントしたり、きりみんちゃんのイラストを描きました。

蟹江もなみさんときりみんちゃんの4コママンガも描きました

これからやる予定のこと

これからの目標というか、これできたら強そうとか強くなれそうだなとかさまざま。

  • AtCoder Rating1800 → 達成できたらARC出る
  • 任意modFPS
  • フロー
  • 2-SAT
  • Segtree beats
  • Codeforces復帰、Rating1900
  • MojaCoder作問
  • 問題解説記事執筆

さいごに

青になるのは難しいですが、青diff以上の問題は面白いです。解説ACでもやってみる価値ありです。

ごめん言い忘れた

少なくとも水色安定したいのであれば、BITとセグ木のライブラリは忘れずに整備しておきましょう。

それでは、今日はこの辺で。

黄色になったときの色変記事に続きます。

ABC180F - Unbranchedを形式的冪級数を使って解く

形式的冪級数（FPS）でやる解法で解いたので、まとめてみます。

問題リンク↓

atcoder.jp

考察

グラフに課される条件のうち、最初の2つの条件

  • 自己ループを持たない
  • すべての頂点の次数が 2以下である

から、考えるべきグラフは直線グラフか頂点数2以上のサイクル（を集めたもの）です。そして、3つめの条件

  • 各連結成分のサイズを並べたとき、その最大値がちょうど Lである

についてですが、典型知識として、最大値がちょうど Lのものの数え上げは、（最大値が L以下であるものの個数）-（最大値が L-1以下であるものの個数）を求めれば良いです。ここまでは公式解説の通りです。

条件を満たすグラフで、連結成分のサイズの最大値が L以下であるものについて考えます。まずは、各連結成分のサイズの振り分け方を決めます。サイズ i(1\le i \le L)の直線グラフの個数を x _ {i}、サイズ i(1\le i \le L)のサイクルの個数を y _ {i}とおくと、頂点数と辺の数の条件から、

 \sum ^ {L} _ {i = 1}ix _ {i} + \sum ^ {L} _ {i = 2}iy _ {i} = N, \quad \sum ^ {L} _ {i = 1}(i - 1)x _ {i} + \sum ^ {L} _ {i = 2}iy _ {i} = M

という関係があります。この2式を整理して、

 \sum ^ {L} _ {i = 1}ix _ {i} + \sum ^ {L} _ {i = 2}iy _ {i} = N, \quad \sum ^ {L} _ {i = 1}x _ {i} = N - M

としておきます。

次に、各連結成分の頂点に番号を振ることを考えましょう。サイズ iの直線グラフの各頂点に与える番号を i個とってきます。これらを各頂点に並べる方法は、反転して同じになるものを考慮して i\ge 2のとき i!/2通り、 i = 1のとき 1通り、 i\le 0のとき 0通りです。これを一般に f _ {1}(i)とおきます。

サイクルについても同様に、 i\ge 3のとき (i-1)!/2通り、 i = 2のとき 1通り、 i\le 1のとき 0通りです。これを一般に f _ {2}(i)とおきます。

さらに、辺にラベルがついてないことから、同じ頂点数、同じタイプの連結成分は区別しません。これをふまえて、求める個数は

 N!\sum \prod ^ {L} _ {i=0}\prod ^ {L} _ {j=0}\dfrac{\{f _ 1(i)\}^{x _ i}}{x _ {i}!(i!)^{x _ {i}}}\dfrac{\{f _ 2(j)\}^{y _ j}}{y _ {j}!(j!)^{y _ {j}}}

です。和の範囲は x _ {i}, y _ {j}が先ほどの条件を満たすようにとります。

形式的冪級数を使う

さて、この和の中にある積の部分ですが、なにか見覚えがないでしょうか？そうですね、expの級数展開の項に出てくるやつですね。以下のものを考えます。

 \sum \prod ^ {L} _ {i=0}\prod ^ {L} _ {j=0}\dfrac{\{f _ 1(i)X ^ {i}Y\}^{x _ i}}{x _ {i}!(i!)^{x _ {i}}}\dfrac{\{f _ 2(j)X ^ {j}\}^{y _ j}}{y _ {j}!(j!)^{y _ {j}}}

これの X, Yの指数は、それぞれ条件より N, N - Mとなります。そしてこれは、

 \prod ^ {L} _ {i=0}\prod ^ {L} _ {j=0}\left(\sum ^ {\infty} _ {k=0}\dfrac{1}{k!}\left(\dfrac{f _ 1(i)X^{i}Y}{i!}\right)^k\right)\left(\sum ^ {\infty} _ {k=0}\dfrac{1}{k!}\left(\dfrac{f _ 2(j)X^{j}}{j!}\right)^k\right)

の X^{N}Y^{N-M}の項にあたります。それぞれの級数はまさにexpの形そのものなので、それを用いて整理すると、

 \exp\left(Y\sum ^ {L} _ {i=0}\dfrac{f _ 1(i)X^{i}}{i!} + \sum ^ {L} _ {i=0}\dfrac{f _ 2(i)X^{i}}{i!}\right)

というとてもすっきりした形になります。これの X^{N}Y^{N-M}の項を求めたいのですが、変数が2つあると複雑なので、偏微分的な要領で Xを定数とみなして Y^{N-M}の項を求めましょう。exp型なので、この場合は Yで N-M回偏微分すれば良いです。結局、

 \dfrac{1}{(N-M)!}\left(\sum ^ {L} _ {i=0}\dfrac{f _ 1(i)X^{i}}{i!}\right)^{N-M}\exp\left(\sum ^ {L} _ {i=0}\dfrac{f _ 2(i)X^{i}}{i!}\right)

の X^{N}の項を求めればよいことになります。これを L-1でも同様に行い、その結果を Lのときの結果から引けば答えが求められます。

\bmod 10^{9} + 7ですが、 N, M, Lが小さいことと、形式的冪級数の逆数などが登場しないことから、FPSのライブラリを持ち出さずにexpの冪級数展開を用いて計算したりしても間に合います。

ABC193E - Oversleepingの解説を書いてみます

コンテスト中に通せたので、自分の解法を書いてみます。

問題

概要

電車が X秒で一方の街からもう一方の街に移動し、街に着いたときに Y秒停車する。高橋君は P秒寝て Q秒起きるを繰り返す。高橋君は街 Bで降車することができるか？もし可能ならば、最短何秒で降車できるか？ Tケースについて求めよ。

制約

  • 1\le T\le 10
  • 1\le X\le 10^{9}
  • 1\le Y\le 500
  • 1\le P\le 10^{9}
  • 1\le Q\le 500

厳密な問題の定義はもとの問題文を読んでください。

E - Oversleeping

解法

状況整理

電車が街 Aを出発してから、戻ってきて次に街 Aを出発するまで 2X+2Y秒かかります。つまり、電車の運動は周期 2X+2Y秒の周期的な運動です。同様に、高橋君の睡眠に関する行動の周期は P+Q秒です。ところで、高橋君が街 Bで電車を降りるためには、

  1. 電車が街 Bに停車している間に高橋君が目を覚ます。
  2. 高橋君が起きている間に電車が街 Bに到着する。

のいずれかが成立していれば良いです。そして、このいずれかが成り立つ最小の時刻が求める答えとなります。

周期的な運動であることから、電車の運動の状態は、時刻を 2X+2Yで割った余りによって決まります。つまり、時刻を tとすると、 t = (2X+2Y)n+k\quad(n, k\in \mathbb{Z}, n\le 0, 0\le k\lt 2X+2Y)で表したときの kによって電車の状態が決まります。一方、高橋君についても同じように t = (P+Q)m+l\quad(m, l\in \mathbb{Z}, m\le 0, 0\le l\lt P+Q)の lによってのみ決まります。 tを消去して、 (2X+2Y)n+k = (P+Q)m+lの形にすることで、 n, mについての二元一次方程式が出来上がります。このそれぞれの時刻の表示を用いて、高橋君が街 Bで電車を降りる条件を定式化してみると、

  1. l = Pかつ X\le k\lt X+Y、つまり、方程式 (2X+2Y)n - (P+Q)m = P - k
  2. k = Xかつ P\le l\lt P+Q、つまり、方程式 (2X+2Y)n - (P+Q)m = l - X

となります。それぞれの条件において、固定されていない方の kおよび lのとりうる範囲は高々 500なので、全ての場合について方程式を解いてみれば良いです。

商とあまりの定義

整数 a, b(b \neq 0)について、整数 q, rを用いて、 a = bq + rと表します。 aを bで割ったあまりとは、 aをこのような形で表したときの rであって、 0\le r\lt |b|となるような唯一の整数のことをいうものとします。また、 aを bで割ったときの商とは、 rに aを bで割ったあまりをとったときの qのことを指すものとします。

以下に例を示します。

a b あまり 商
24 5 4 4
-24 5 1 -5
24 -5 4 -4
-24 -5 1 5
40 -8 0 -5
-13 26 13 -1
-17 17 0 -1
-11 -5 4 3

拡張ユークリッドの互除法

いきなり一般的な問題を解くのは難しいと思うので、具体例を考えます。 8n - 27m = 1という方程式を解いてみましょう。後で示しますが、 n, mの係数(ここでは 8, -27)の最大公約数が右辺の約数であるとき、またその時に限り、この形の方程式は整数解を持ちます。よって、この具体例では整数解が存在します。

2つの変数に対して方程式が1つなので、整数解は無数に存在します。たとえば、 n = 17, m = 5は解のひとつです。他にも n = -10, m = -3も解のひとつです。このような解のひとつを機械的に求めるために、ユークリッドの互除法を使います。以下のように進んでいきます。

  1. 8 = -27\times 0 + 8
  2. -27 = 8\times (-4) + 5
  3. 8 = 5\times 1 + 3
  4. 5 = 3\times 1 + 2
  5. 3 = 2\times 1 + 1
  6. 2 = 1\times 2 + 0
  7. あまりが 0になったので終了。

よって、最大公約数が 1であることが求まりますが、これらの各式を a = bq + rの形でみたときに、 a, bがどのように変化するかを行列表示で見てみたいと思います。

行列表示を知らない方もいると思うので必要なことだけ簡単に説明すると、行列とは要素を矩形に並べたもので、行列同士の積として

 \begin{pmatrix}a & b \\\ c & d\end{pmatrix} \begin{pmatrix}u \\\ v\end{pmatrix} = \begin{pmatrix}au + bv \\\ cu + dv\end{pmatrix}
 \begin{pmatrix}a & b \\\ c & d\end{pmatrix} \begin{pmatrix}p & q \\\ r & s\end{pmatrix} = \begin{pmatrix}ap + br & aq + bs \\\ cp + dr & cq + ds\end{pmatrix}

がなりたちます。また、2つの行列が等しいとは、両方の行列のサイズが等しく、かつ同じ場所にある要素が等しいことを意味します。

a _ 0 = 8, b _ 0 = -27として、 n + 1番目の式を a _ {n} = b _ {n}q _ {n} + r _ {n}と表します。ユークリッドの互除法のアルゴリズムのもっとも重要な事実として、 a _ {n}と b _ {n}の最大公約数は b _ {n}と r _ {n}と等しい、というものがあります。よって、 a _ {n+1} = b _ {n}, b _ {n+1} = r _ {n} = a _ {n} - b _ {n}q _ {n}です。このことを行列表示すると

 \begin{pmatrix}a _ {n+1} \\\ b _ {n+1}\end{pmatrix} = \begin{pmatrix}0 & 1 \\\ 1 & -q _ {n}\end{pmatrix} \begin{pmatrix}a _ {n} \\\ b _ {n}\end{pmatrix}

となります。これは

 \begin{pmatrix}a _ {n} \\\ b _ {n}\end{pmatrix} = \begin{pmatrix}0 & 1 \\\ 1 & -q _ {n-1}\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & -q _ {n-2}\end{pmatrix} \cdots \begin{pmatrix}0 & 1 \\\ 1 & -q _ {0}\end{pmatrix}\begin{pmatrix}8 \\\ -27\end{pmatrix}

のように書けるので、

 \begin{pmatrix} x _ {n} & y _ {n} \\\ u _ {n} & v _ {n}\end{pmatrix} = \begin{pmatrix}0 & 1 \\\ 1 & -q _ {n-1}\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & -q _ {n-2}\end{pmatrix} \cdots \begin{pmatrix}0 & 1 \\\ 1 & -q _ {0}\end{pmatrix}

とおきます。すると、

 \begin{pmatrix} x _ {n+1} & y _ {n+1} \\\ u _ {n+1} & v _ {n+1}\end{pmatrix} = \begin{pmatrix}0 & 1 \\\ 1 & -q _ {n}\end{pmatrix}\begin{pmatrix} x _ {n} & y _ {n} \\\ u _ {n} & v _ {n}\end{pmatrix} = \begin{pmatrix} u _ {n} & v _ {n} \\\ x _ {n} - q _ {n}u _ {n} & y _ {n} - q _ {n}v _ {n}\end{pmatrix}

となります。この行列を求めることや計算していくことは、漸化式にしたがってその都度一つずつ順番に行えます。そして、初期値によらず最終的には b _ {n}が 0、 a _ {n}が a _ 0, b _ 0の最大公約数になって終了するので、

 \begin{pmatrix}1 \\\ 0\end{pmatrix} = \begin{pmatrix}x _ {f} & y _ {f} \\\ u _ {f} & v _ {f}\end{pmatrix}\begin{pmatrix}8 \\\ -27\end{pmatrix}

となります。これの上側の要素を比較すると 8x _ {f} - 27 y _ {f} = 1という式を得ます。この x _ {f}, y _ {f}は元の方程式 8n - 27m = 1の解になっています。具体的には、

 \begin{pmatrix} x _ {f} & y _ {f} \\\ u _ {f} & v _ {f}\end{pmatrix} = \begin{pmatrix}0 & 1 \\\ 1 & -2\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & -1\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & -1\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & -1\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & 4\end{pmatrix}\begin{pmatrix}0 & 1 \\\ 1 & 0\end{pmatrix} = \begin{pmatrix} -10 & -3 \\\ 27 & 8 \end{pmatrix}

であるので、解 n = -10, m = -3が見つかりました。この解を利用することで、 8n - 27m = 2の解はこれを2倍した n = -20, m = -6がみつかることもわかります。

また、 a _ 0, b _ 0を a _ {f}, 0に変えるまでの最悪計算量は、 a _ 0, b _ 0がフィボナッチ数列の隣り合う項のときなので \mathrm{O}(\log \max(a _ 0, b _ 0))です。

さらに、 n = -10 + 27u, m = -3 + 8u(u \in \mathbb{Z})も方程式 8n - 27m = 1を満たすことがわかると思います。このことは一般に、 an + bm = cに整数解が存在するなら、 nが 0以上 |b|未満になるように解を決めることが可能であることが言えます。そのように実装しておくと、解が存在しないときに -1, -1を返すことで判定できるので便利です。

改めて、方程式 ax + by = cの整数解を求める手順をまとめてみます。ここで、 k/lを kを lで割った商、 k \% lを kを lで割ったあまりとします。

  1. x = 1, y = 0, u = 0, v = 1, k = a, l = bとする。
  2. lが 0ならば手順4. へ。そうでなければ手順3. へ。
  3. x = u, y = v, u = x - u \times (k / l), v = y - v \times (k / l), k = l, l = k\% lとする。手順2. にもどる。
  4. この時点で kは aと bの最大公約数になっている。もし cが kで割り切れなければ、 (x, y) = (-1, -1)を返し終了。
  5. a = a / k, b = b / k, c = c / kとする。
  6. x, yを c倍する。
  7. (x + ub, y - ua)も方程式の解としてとれる。 u = \lfloor \dfrac{b - x - 1}{b}\rfloorとすると、 0\le x + ub\lt |b|となる。
  8. (x + ub, y - ua)を返して終了。

本題

拡張ユークリッドの互除法で条件1. と条件2. を満たす n, mを求めます。ただし、解が存在するならばどちらも 0以上でなくてはならないので、その条件も確かめる必要があります。 nは自動的にこの条件を満たします。 mが負になってしまう可能性があるように思えますが、それは起こらないことが示せます。

0\le n \lt P+Qであったとすると、式変形により、条件1. では \lfloor \dfrac{Q-1 + k}{P+Q}\rfloor \le m \le a + \lfloor \dfrac{Q-1 + k}{P+Q}\rfloor、条件2. では \lfloor \dfrac{P + Q + X - l -1}{P+Q}\rfloor \le m \le a + \lfloor \dfrac{P + Q + X - l -1}{P+Q}\rfloorですが、 1\le Q + X - 1\le Q - 1 + k \lt Q + X + Y - 1、 0 \le X - 1\lt P + Q + X - l - 1\le Q + X - 1であるので、自動的に mが 0以上になることがわかります。

さらに、以上の方法で求めた nは最小の非負整数解であるので、これで欲しいものは求まりました。あとは、解が存在するような (2X+2Y)n + kのもとでの最小値を求めれば良いです。

全体的な計算量は \mathrm{O}((P+Q)\log (\max(P+Q, 2(X+Y))))です。

コンテスト中のACコード

atcoder.jp

余談

解が存在しないときの出力が'infinity'になっていることに気をつけましょう。僕はそれで2ペナ食らいました。

ABC183E - Queen on Grid：自分用メモ

問題

E - Queen on Grid

H行 W列のマス目がある。左上のマス目を (1, 1)として、上から i番目、左から j番目のマス目をマス (i, j)とする。1手ごとに、そのマスから下、右、右下方向の好きなマスに移動できる。ただし、障害物に移ったり、障害物をを飛び越える移動や、マス目の外に出る移動は許されない。マス (H, W)に移動する方法は何通りか？

初手の考察

素朴な遷移式

マス目 (x, y)にいく方法を dp[x][y]とする。マス目 (x, y)について、左、上、左上方向のそれぞれで最も (x, y)に近い障害物のある座標を (x, l _ {x, y}), (u _ {x, y}, y), (d _ {x, y}, y - x + d _ {x, y})とする。遷移式は

dp[x][y] = \sum^{y-1} _ {j = l _ {x, y}}dp[x][j] + \sum^{x-1} _ {i = u _ {x, y}}dp[i][y] + \sum^{x-1} _ {i = d _ {x, y}}dp[i][y-x+i]

となる。

実装

メモ化再帰で実装します（思考停止）。

結果

落ちます。

考察をやり直す

L[x][y] = \sum^{y-1} _ {j = l _ {x, y}}dp[x][j]\\ U[x][y] = \sum^{x-1} _ {i = u _ {x, y}}dp[i][y]\\ D[x][y] = \sum^{x-1} _ {i = d _ {x, y}}dp[i][y-x+i]

という値を定義します。こうすると、遷移式は、

dp[x][y] = L[x][y] + U[x][y] + D[x][y]

です。それぞれの値について、特徴を整理してみます。

dp[1][1] = 1

L[x][y] = \begin{cases} L[x][y-1] + dp[x][y-1] \quad (\mathrm{Grid}[x][y] = ".") \\ 0 \quad (\mathrm{otherwise})\end{cases}

U[x][y] = \begin{cases} U[x-1][y] + dp[x-1][y] \quad (\mathrm{Grid}[x][y] = ".") \\ 0 \quad (\mathrm{otherwise})\end{cases}

D[x][y] = \begin{cases} D[x-1][y-1] + dp[x-1][y-1] \quad (\mathrm{Grid}[x][y] = ".") \\ 0 \quad (\mathrm{otherwise})\end{cases}

トポロジカル順序を考えると、どの式も W+3\le i \le H(W+2)+Wについて、 iの昇順で (i//(W+2), i\% (W+2))を更新していく方針で問題なさそうです。

正しい実装

atcoder.jp

群論と離散対数のお時間です - 3

さて、ここからが本題です。 X^{K}\equiv Y\bmod Mで、 Xと Mが互いに素とは限らない一般の場合の離散対数問題を解きます。

実験

Aから AX\bmod Mの向きに辺を結んでできた有向グラフを考えます。鳩の巣原理より、どのような Xについても、少なくとも M乗するとより小さい指数乗の値に戻ることがわかるので、いつかは同じ頂点に戻ってきます。この巡回する部分を「サイクル」と呼ぶことにします。また、サイクルに含まれない部分を「テール」と呼ぶことにします。いくつかの Mについて、 X^{k}がどのようになるか見てみたいと思います。

mod 10

Mの素因数分解 \varphi(M) \varphi(M)の約数
2\times 5 4 \{1, 2, 4\}
X サイクルの長さ テールの長さ
0 1 1
1 1 0
2 4 1
3 4 0
4 2 1
5 1 1
6 1 1
7 4 0
8 4 1
9 2 0

f:id:DivineJK:20210204135707p:plain
mod 10でのふるまい

mod 8

Mの素因数分解 \varphi(M) \varphi(M)の約数
2^{3} 4 \{1, 2, 4\}
X サイクルの長さ テールの長さ
0 1 1
1 1 0
2 1 3
3 2 0
4 1 2
5 2 0
6 1 3
7 2 0

f:id:DivineJK:20210204141617p:plain
mod 8でのふるまい

mod 30

Mの素因数分解 \varphi(M) \varphi(M)の約数
2\times 3\times 5 8 \{1, 2, 4, 8\}
X サイクルの長さ テールの長さ
0 1 1
1 1 0
2 4 1
3 4 1
4 2 1
5 2 1
6 1 1
7 4 0
8 4 1
9 2 0
10 1 1
11 2 0
12 4 1
13 4 0
14 2 1
15 1 1
16 1 1
17 4 0
18 4 1
19 2 0
20 1 1
21 1 1
22 4 1
23 4 0
24 2 1
25 1 1
26 2 1
27 4 1
28 4 1
29 2 0

f:id:DivineJK:20210204142358p:plain
mod 30でのふるまい

mod 12

Mの素因数分解 \varphi(M) \varphi(M)の約数
2^{2}\times 3 4 \{1, 2, 4\}
X サイクルの長さ テールの長さ
0 1 1
1 1 0
2 2 2
3 2 1
4 1 1
5 2 0
6 1 2
7 2 0
8 1 1
9 1 1
10 1 2
11 2 0

f:id:DivineJK:20210205180500p:plain
mod 12でのふるまい

予想

以上の結果から、いくつかの予想ができます。

  1. Xの素因数分解が X = p _ 1^{e _ 1}\times p _ 2^{e _ 2}\times\cdots\times p _ m^{e _ m}であったとする。すべての i(1\le i \le m)について、 X^{k}の素因数分解における p _ iの指数が Mのそれ以上であれば、 X^{k} \bmod Mはサイクルに含まれる。そうでなければ、テールに含まれる。
  2. サイクルの長さは \varphi(M)の約数である。
  3. サイクルに含まれる数と Mとの最大公約数は全て等しい。

証明

前提

  • \mathrm{gcd}(A, B) = \mathrm{gcd}(A\%B, B)
  • \mathrm{gcd}(A, B) = 1 \Rightarrow \varphi(AB) = \varphi(A)\varphi(B)

1. と2. の証明

Xを素因数分解して現れる全ての素数について、 Mに含まれるものをすべて掛け合わせたものを Pとし、 M' = M / Pとします。累乗の計算によって素数の種類数は変化しないので、 Pと M'は Xの何乗かに依存しません。

予想1. での仮定から、 X^{k}\equiv 0 \bmod Pです。よって、任意の非負整数 lに対して X^{k+l}\equiv X^{k} \equiv 0 \bmod Pです。

また、 M'の定義より、 Xと M'は互いに素です。よって、フェルマーの小定理から X^{\varphi(M')}\equiv 1\bmod M'が成り立ちます。両辺に X^{k}を掛け合わせて、 X^{k+\varphi(M')}\equiv X^{k}\bmod M'を得ます。

以上より、中国剰余定理を用いることで X^{k+\varphi(M')}\equiv X^{k}\bmod Mであることが示されました。さらに、予想2. のより強い主張として、「サイクルの長さは \varphi(M')の約数である。」が得られました。

また、ある i(1\le i \le m)について、 X^{k}の素因数分解における p _ iの指数が Mのそれ未満であるときを考えてみます。このとき、 X^{k}は Pで割り切れないためサイクル内のどの要素とも一致しません。よって、サイクルには含まれません。そして、 Xを掛けていくことで一度サイクルに入ってしまうとそれ以降はサイクルから抜けることはありません。よってこのとき、 X^{k}はテールに含まれます。

3. の証明

先程の証明の M', Pを引き継ぎます。サイクルに含まれる任意の要素 Aについて、 A\equiv 0 \bmod Pであり、 Aと M'は互いに素です。よって、 \mathrm{gcd}(A, M) = Pです。

実装への考察

テールの部分を別に調べて、そこで解が見つかれば終了、そうでなければサイクルをBSGSの要領で調べれば良さそうです。その際、BSGSで見つかった解にはテールの長さを足す必要があります。

テールの長さ

Xの素因数分解を X = p _ 1^{e _ 1}\times p _ 2^{e _ 2}\times\cdots\times p _ m^{e _ m}とします。また、 P = p _ 1^{f _ 1}\times p _ 2^{f _ 2}\times\cdots\times p _ m^{f _ m}とします。予想1. の帰結より、 X^{k}がサイクルに含まれる条件は、全ての i(1\le i \le m)が存在して ke _ i \ge f _ iとなることです。よって、テールの長さは \min _ {1\le i\le m}(\lceil f _ i / e _ i\rceil)になります。

サイクルの性質

サイクルの長さを Cとします。サイクルの中の任意の数 Aについて、 AX^{C}\equiv A\bmod Mです。よって、このサイクルの単位元「のようなもの」は X^{C}です。「のようなもの」という表現は、 X^{C}がこのサイクルの中に含まれていないことがあるためです（例： \bmod 40, X = 14のとき、 C=2, X^{C}=36となりますが、これはサイクルに含まれません）。しかし、 Xをかけていくといつかはサイクルの中に入ることと、 X^{tC}もまた AX^{tC}\equiv A\bmod Mを満たすので、サイクルの中から単位元を見つけることができます。実装上はどちらでも問題ないです。 X^{C}の方が計算量は少ないですが、サイクルの中の要素で演算を閉じさせたいという場合はそのような要素を探すこともできます。

そして、「 Xを掛ける」という操作の逆演算、つまり、サイクルの進む向きとは逆に進む操作を考えることができます。これが X'を掛けることによって実現するとしましょう。言い換えると、サイクル内の任意の数 Aについて AXX'\equiv A\bmod Mが成り立ちます。先程の単位元「のようなもの」を考えたときの式から推察できることとして、たとえば XX'\equiv X^{C}\bmod M、すなわち X'\equiv X^{C-1}とするといつでもこの等式が成り立ちます。この X'を逆元「のようなもの」と呼び、 X^{-1}と表記します。

別の方法でこの X^{-1}を求める方法があります。単位元「のようなもの」を eとします。満たすべき等式は XX'\equiv e\bmod Mです。ところで、競技プログラミングでおなじみのmod逆元というものがありますが、大きく2通りの求め方がありました。ひとつはフェルマーの小定理を使って求めるもの、もうひとつは拡張ユークリッドの互除法を使って求めるもの、です。 X^{C-1}で求めた方法はフェルマーの小定理的です。実は、拡張ユークリッドの互除法によって X^{-1}を求めることができます。

等式 ax \equiv b\bmod mを満たす整数 xが存在する条件として、 b\% \mathrm{gcd}(a, m) = 0です。今回の場合、 e\% \mathrm{gcd}(X, M) = 0です。予想1. の証明で出てきた Pを用いると、 \mathrm{gcd}(e, m)は Pまたは X^{C}に等しく、どちらも \mathrm{gcd}(X, M)で割り切れます。よって、 X^{-1}の存在が保証されたと同時に、拡張ユークリッドの互除法で X^{-1}を求めることができることがわかりました。

実装

実装方針

以上のことをまとめると、実装方針は次のようになります。

入力: M, X, Y: 整数 (M\ge 1, 0\le X, Y\lt M)

出力: X^{K}\equiv Y\bmod Mを満たす最小の非負整数 K（存在しなければ -1）

  1. X = 0のとき、別で処理する。
  2. Xを素因数分解する。 X = p _ 1^{e _ 1}\times p _ 2^{e _ 2}\times\cdots\times p _ m^{e _ m}とする。
  3. 全ての 1\le i\le mについて、 Mが p _ iで最大 f _ i回割り切れるとして、 ke _ i \le f _ iとなる最小の kを求める。この kがテールの長さで、 tailとおく。
  4. 0\le i \lt tailについて、 X^{i}\equiv Y\bmod Mとなる iがあれば、 iを返して終了。
  5. bas = X^{tail}とする。もし Yが \mathrm{gcd}(bas, M)で割り切れなければ、 -1を返して終了。
  6. サイクルの長さが \varphi(M)の約数であることから、 bas \times X^{C}\equiv bas\bmod Mである最小の約数 Cをさがす。この Cがサイクルの長さである。
  7. 以下、サイクル内の要素でBSGS。

計算量を見積もってみると、BSGSと素因数分解がネックとなり、 \mathrm{O}(\sqrt{M} + \sqrt{X})です。

コード

github.com

verify

Discrete Logarithm：Library Checker

群論と離散対数のお時間です - 2

前回の続きです。一般化されたBSGSはまとまりをよくするために次回に回し、今回は、

  • Pohlig-Hellman algorithm
  • BSGS, Pohlig-Hellman algorithmの実装

について書いていきたいと思います。

Pohlig-Hellman algorithm

\varphi(M)を素因数分解したときの指数が大きい場合に有効なアルゴリズムであるPohlig-Hellman algorithmを紹介します。ようやくここで群論が登場します。

ここでも Mと Xが互いに素である必要があります。ところで、これは X^{\varphi(M)}\equiv 1\bmod Mが成り立つことを意味しますが、 \varphi(M)は X^{k}\equiv 1\bmod Mを満たすものの中で最小ではない可能性があります。しかし、この kは \varphi(M)の約数であることが示せます。そこで、 \varphi(M)の約数のうちで等式を満たすような最小なものを探しておきます。そして、これが巡回群の位数となります。

以下、この巡回群の位数を nとします。そして、 n = p_1^{e_1}\times p_2^{e_2}\times\cdots\times p_m^{e_m}のように素因数分解できるとします。ここで注意するべきなのは、今考えているのはあくまで \bmod Mでの話題なので、 nではなく Mでmodをとらなくてはなりません（私が実装でハマった原因はこれでした）。

まず、Pohlig-Hellman algorithmの特別な場合の部分問題として、 n = p^{e}のように、素因数分解してひとつの素数のみのべき乗の形になっている場合を考えます。解く問題としては、 X^{K}\equiv Y\bmod Mなる 0\le K\lt p^{e}を求める、という問題です。位数が p^{e}なので、解があるとすればこの範囲に必ずおさまることが保証されます。この Kは、 K=a_0 + a_1 p +\cdots +a_ {e-1} p^{e-1} \, (0\le a_0 ,a_1 ,\cdots , a _ {e-1} \lt p)と表せるので、この a_0, a_1, \cdots, a _ {e-1}を求めたいです。

a_0, a_1, \cdots, a _ {k-1}までが求まっているとして、 R = a_0 +a_1p+ \cdots + a _ {k-1}p^{k-1}とします（何も求まっていないときは R = 0とします）。このとき、方程式を変形して X^{a _ k p^{k} + \cdots + a _ {e-1} p^{e-1}} \equiv Y(X^{-1})^{R}\bmod Mを得ることができます。この両辺を p^{e-k-1}乗します。今考えている巡回群の位数は p^{e}ですから、 X^{cp^{e}} \equiv 1\bmod M\quad (c=0,1,2,\cdots)となります。よって、左辺は X^{a _ k p^{e-1}} = (X^{p^{e-1}})^{a _ k}のみ残ります。

ここで、 \gamma = X^{p^{e-1}}とおくことで、 \gamma^{a _ k} \equiv (YX^{-R})^{p^{e-k-1}}\bmod Mです。 \gammaの生成する巡回群の位数は pであるので、これを満たす a _ kをBSGSで求めることができます。もし解が存在しなければその旨を報告し、終了です。そうでないとき、 Rに a _ kp^{k}を足し、 a _ {k+1}を求める計算に進みます。

この操作を k=0から順に行い、 K=Rが求める答えになります。 0\le k\lt eで、各 kあたりの計算量は \mathrm{O}(\sqrt{p})なので、全体的な計算量は \mathrm{O}(e\sqrt{p})です。

もとの問題に戻ります。位数が n = p_1^{e_1}\times p_2^{e_2}\times\cdots\times p_m^{e_m}と表され、 X _ i = X^{n/p _ i^{e _ i}}は p _ i^{e _ i}の巡回群を生成します。同時に、 Y _ i = Y^{n/p _ i^{e _ i}}とします。問題の方程式の両辺を n/p _ i^{e _ i}乗することで、方程式 X _ i^{K} \equiv Y _ i\bmod Mを得ます。この方程式は先程の部分問題の解き方で解くことができます。また、この方程式の解を K _ iとすると、これは Kの \bmod p _ i^{e _ i}での値になります。

そうして得られた位数 p _ i^{e _ i}での方程式の解 K _ iをもとに、中国剰余定理から Kを求めることができます。

全体の計算量は \mathrm{O}(\sum_i e _ i(\sqrt{p _ i} + \log n))です。

実装例

BSGSの実装例

def divisors(n): # nの約数
    L = []
    p = 1
    while p * p <= n:
        if n % p == 0:
            L.append(p)
        p += 1
    for i in range(len(L), 0, -1):
        if L[i-1] * L[i-1] != n:
            L.append(n // L[i-1])
    return L # nの約数が昇順に並んでいる。
def totient(n): # nと互いに素なn以下の正整数の個数
    a = n
    p = 2
    res = n
    if a % p == 0:
        res >>= 1
        while 1 - a % 2:
            a >>= 1
    p += 1
    while a != 1:
        if a % p == 0:
            res //= p
            res *= p - 1
            while a % p == 0:
                a //= p
        p += 2
        if p * p > n and a != 1:
            res //= a
            res *= a - 1
            break
    return res
def babystep_giantstep(n, a, b): # a^k = b (mod n)なる最小の非負整数kを求める（存在しなければ-1）。
    baby_step = {}
    candidate = divisors(totient(n)) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], n) != 1: # a^k = 1 (mod n)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない（一般化した時に扱う）。
    p = candidate[pnt] # 巡回群の位数
    l, r = 0, p
    m = (l + r) // 2 # m = ceil(sqrt(p))
    while r - l > 1:
        if m * m <= p:
            l = m
        else:
            r = m
        m = (l + r) // 2
    if m * m < p:
        m += 1
    x, y, u, v, k, l = 1, 0, 0, 1, a, n
    while l:
        x, y, u, v, k, l = u, v, x - u * (k // l), y - v * (k // l), l, k % l
    bas = pow(x, m, n) # bas = a^(-m) (mod n)
    e = 1
    for i in range(m): # baby_stepに(a^i, i)を記録していく。
        baby_step[e] = i
        e *= a
        e %= n
    y = b
    for i in range(m):
        if y in baby_step:
            return i * m + baby_step[y] # これが最小の解
        y *= bas # giantstep
        y %= n
    return -1 # 見つからなければ-1を返す。

Pohlig-Hellman algorithmの実装例

def divisors(n): # nの約数
    L = []
    p = 1
    while p * p <= n:
        if n % p == 0:
            L.append(p)
        p += 1
    for i in range(len(L), 0, -1):
        if L[i-1] * L[i-1] != n:
            L.append(n // L[i-1])
    return L # nの約数が昇順に並んでいる。
def totient(n): # nと互いに素なn以下の正整数の個数
    a = n
    p = 2
    res = n
    if a % p == 0:
        res >>= 1
        while 1 - a % 2:
            a >>= 1
    p += 1
    while a != 1:
        if a % p == 0:
            res //= p
            res *= p - 1
            while a % p == 0:
                a //= p
        p += 2
        if p * p > n and a != 1:
            res //= a
            res *= a - 1
            break
    return res
def babystep_giantstep(n, a, b): # a^k = b (mod n)なる最小の非負整数kを求める（存在しなければ-1）。
    baby_step = {}
    candidate = divisors(n-1) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], n) != 1: # a^k = 1 (mod n)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない（一般化した時に扱う）。
    p = candidate[pnt] # 巡回群の位数
    l, r = 0, p
    m = (l + r) // 2 # m = ceil(sqrt(p))
    while r - l > 1:
        if m * m <= p:
            l = m
        else:
            r = m
        m = (l + r) // 2
    if m * m < p:
        m += 1
    x, y, u, v, k, l = 1, 0, 0, 1, a, n
    while l:
        x, y, u, v, k, l = u, v, x - u * (k // l), y - v * (k // l), l, k % l
    bas = pow(x, m, n) # bas = a^(-m) (mod n)
    e = 1
    for i in range(m): # baby_stepに(a^i, i)を記録していく。同じkeyは小さい方をとる。
        if e not in baby_step:
            baby_step[e] = i
        e *= a
        e %= n
    y = b
    for i in range(m):
        if y in baby_step:
            return i * m + baby_step[y] # これが最小の解
        y *= bas # giantstep
        y %= n
    return -1 # 見つからなければ-1を返す。
def extgcd(a, b, c): # ax + by = cを満たす整数の組(x, y)を探す。
    if b == 0:
        if a == 0:
            if c == 0:
                return 0, 0
            return None
        if c % a == 0:
            return c // a, 0
        return None
    if b < 0:
        a, b, c = -a, -b, -c
    tk, tl = a, b
    while tl:
        tk, tl = tl, tk % tl
    if c % tk:
        return None
    a //= tk
    b //= tk
    c //= tk
    x, y, u, v, k, l = 1, 0, 0, 1, a, b
    while l:
        x, y, u, v = u, v, x - u * (k // l), y - v * (k // l)
        k, l = l, k % l
    x *= c
    k = x // b
    x -= k * b
    y *= c
    y += k * a
    return x, y
def CRT(num, a_list, m_list): # 複数modの条件を満たす整数を見つける。
    r = 0
    bas = 1
    x, y = 0, 0
    for i in range(num):
        x, y = extgcd(bas, -m_list[i], a_list[i]-r)
        r += bas * x
        bas *= m_list[i]
    return r % bas
def pohlig_hellman(mod, a, b):
    candidate = divisors(totient(mod)) # 巡回群の位数の候補
    pnt = 0
    while pow(a, candidate[pnt], mod) != 1: # a^k = 1 (mod)となる最小のkを探す。それが巡回群の位数となる。
        pnt += 1
        if pnt == len(candidate):
            return -1 # nとaが互いに素でない（一般化した時に扱う）。
    n = candidate[pnt] # 巡回群の位数
    factors = [] # nの素因数分解
    res_seed = [] # 素数のべき乗での答えを格納
    mod_list = []
    r = 0
    cnt = 0
    tmp = n
    # 素因数分解パート
    while tmp % 2 == 0:
        cnt += 1
        tmp >>= 1
    if cnt:
        r += 1
        res_seed.append(0)
        mod_list.append(0)
        factors.append((2, cnt))
    pr = 3
    while tmp != 1:
        cnt = 0
        while tmp % pr == 0:
            tmp //= pr
            cnt += 1
        if cnt:
            r += 1
            res_seed.append(0)
            mod_list.append(0)
            factors.append((pr, cnt))
        pr += 2
        if pr * pr > n and tmp != 1:
            res_seed.append(0)
            mod_list.append(0)
            r += 1
            factors.append((tmp, 1))
            break
    # 素数のべき乗についての部分問題を解く。
    for i in range(r):
        p = factors[i][0]
        e = factors[i][1]
        p_i = pow(p, e)
        mod_list[i] = p_i
        g = pow(a, n//p_i, mod) # gは位数p^eの群をなす。
        h = pow(b, n//p_i, mod)
        # g^res_seed[i] = h (mod) となるres_seed[i]を求める。
        # res_seed[i] = x_0 + x_1 * p + ... + x_{e-1} * p^{e-1} (0 <= x_j < p)の形で表せる。
        bas = 1
        inv_g = extgcd(g, mod, 1)[0] # gの逆元
        gamma = pow(g, pow(p, e-1), mod) # これは位数pの巡回群をなす。
        for k in range(e):
            # (x_0, x_1, .... x_{k-1})までわかっているとする。
            # 式変形により、g^(x_k*p^k+...+x_{e-1}*p^{e-1}) = h * g^(-(x_0+x_1*p+...+x_{k-1}*p^{k-1}))
            # 両辺p^{e-k-1}することで、gamma^(x_k) = ...の形になるので、x_kをbsgsで求める。
            hk = pow(h*pow(inv_g, res_seed[i], mod)%mod, pow(p, e-k-1), mod)
            xk = babystep_giantstep(mod, gamma, hk)
            if xk < 0:
                return -1 # x_kが見つからなければ-1
            res_seed[i] += xk * bas
            bas *= p
    # CRTでそれぞれの部分問題の答えをまとめる。
    res = CRT(r, res_seed, mod_list)
    return res

verify

atcoder.jp Xと Pが固定なので、辞書や位数を一度計算してメモしておくように実装すると、計算量を減らすことができます。

群論と離散対数のお時間です - 1

BSGSを知っていますか？私は知りません。

ということでBSGSのWikipediaを読み漁っていたら、何やらさらにつよいPohlig-Hellman's algorithmなるアルゴリズムがあるじゃないですか！ちょうどその数日前にSA-ISを実装して調子に乗っていたため、アルゴリズムがわりとシンプルじゃけん実装できるじゃろ！と思い実装したところ、世の中はそんなに甘くなく、めちゃくちゃバグり散らしました。世の中ってどうしたら甘くなるんでしょうね。この世界に真っ赤なジャムでも塗れば良いんでしょうか。ところで私ははちみつ派です。おすすめのはちみつがあればご一報ください。

さて、このアルゴリズムですが、なにやら群論の知識をたくさん使うらしいので、メモ化再帰的にわかるところまで進んでから一つずつ回収していく、という手法で理解しようと思いました。また、一応私は大学での専攻が物理学で、世の中の99.7%くらいの人には理解されないので詳細は省きますが、群論など抽象的な代数学を多用する分野に興味があるので、まとめてやっちゃおう、という勢いでやってます。個人用のメモでも良かったんですが、世の中には私と同じようにBSGSを一般化したい人がいるかもしれないと思い、公開します。

材料としてはおもに昔受けた群論に関する講義での手書きの講義ノートです。このこともデジタルな文書として残しておきたい理由で、私のページをケチる癖でよくわからないことになっている部分があるのと、後で見返したときにページ内検索ができるためです。この記事を読んでいる若い学生さんはちゃんと整理整頓されたわかりやすいノートを書きましょうね。

目標

  • ラグランジュの定理についてやる。
  • Pohlig-Hellman's algorithm を実装する。

基本

定義がいっぱいです。

群の公理

群の公理は以下のようなものです。 Gを群とし、 Gの元 x, yの積を単に xyと書くことにします。

  1. 演算が閉じている。 \forall x, y \in G; xy \in G
  2. 結合則が成り立つ。 \forall x, y, z \in G; x(yz) = (xy)z
  3. 単位元が存在する。 \exists e\in G \,\, \mathrm{s.t.}\,\,\forall x\in G; ex = xe = x
  4. 逆元が存在する。 \forall a\in G\, \exists b\in G\,\,\mathrm{s.t.}\,\, ab = ba = e

群の一例として、整数の集合に群の演算として足し算の構造を入れたものがあります（どの隣り合う二つの数の足し算から実行しても結果は同じで、単位元0と整数 aに対する逆元 -aがあります）。また、1.と2.のみを満たすものを半群、1.と2.と3.のみを満たすものをモノイドといいます。モノイドはセグ木に乗ることでよく知られていますね。 n\times n行列の集合に積の演算を入れたものはモノイドの一例です。

部分群

群 Gの部分集合 Hであって、 Gの演算で群のとき、 Hは Gの部分群であるといいます。

群の生成系

群 Gの部分集合 Sであって、 Gからいくつか（無限個でもよい）元をとってきたものを考えます。つまり、 S = \{ g_1, g_2, \cdots \} (g_i\in G)です。 Sの生成する Gの部分群とは、この集合の（重複を許して）有限個の元およびその逆元の積全体からなる集合、すなわち、 nを選んだ元の個数として、 g_1^{\pm 1}g_2^{\pm 1}\cdots g_n^{\pm 1}がなす集合です。この集合は、 \langle S\rangleと書かれます。

すこし説明が難しいので例を考えます。 S = \{2, 3\}とすると、 \mathrm{gcd}(2, 3) = 1であるので任意の整数 nについて 2x + 3y = nを満たす整数 x, yが存在します。よって、 Sは \mathbb{Z}を生成します。ここで、 4 = 2 + 2のように、同じ元を複数回選んでも良いです。

巡回群

一つの元 g\in Gから生成される群 \langle g\rangleは巡回群とよばれます。そして、もし g^{m} = e となるような最小の自然数 mが存在するなら、この mを位数といい、 \langle g\rangleは位数 mの有限巡回群といいます。存在しなければ無限巡回群と呼ばれます。

有限巡回群の例として、 \mathrm{mod}\, m上での足し算があります。 1を m個足していくと単位元 0になります。

群準同型

f: G\rightarrow G'な写像を考えます。例によって G, G'は群です。そして、 \forall g_1, g_2\in Gに対し f(g_1g_2) = f(g_1)f(g_2)となるとき、 fは準同型写像であるといい、このような写像が存在すれば Gと G'は群準同型であるとします。

準同型写像というのはどういうものかを説明します。もともと群には二項演算が入っています。これは、 g_1, g_2\in Gに対し g_1g_2\in Gという群の要素を定めるので、この二項演算は G\times G \rightarrow Gです。さらに、写像 fによってこの要素を写すと、 f: G\rightarrow G'なので f(g_1g_2)\in G'です。つまり、 Gの二つの要素に対して二項演算をしてから写像 fで G'に写す、という操作です。この操作はまとめると G\times G\rightarrow G'です。

一方、 Gの要素を写像で写してからそれらの二項演算を行っても G\times G\rightarrow G'となります。このとき、 g_1, g_2\in Gに対し f(g_1)f(g_2)\in G'となります。つまり、 fが準同型写像であれば、この2通りの結果が等しいということなのです。

さらに、この写像が全単射であれば、 fは同型写像、または、 Gと G'が同型であるといい、 G\simeq G'と表されます。

同値類

集合 Aについて R: A\times A\rightarrow \{0, 1\}というような関係を二項関係といい、 a, b, c\in Aに対して以下の3つを満たす二項関係 Rを同値関係といいます。

  1. a\,R\,a （反射律）
  2. a\,R\,b\, \Rightarrow b\,R\,a （対称律）
  3. a\,R\,b \wedge b\,R\,c\Rightarrow a\,R\,c （推移律）

例えば等号の =はこの条件を満たします。意味的にはもっと広く、「 aと bがある基準を以って仲間ですか」という質問にYesであれば同値、Noであれば同値ではない、という感じです。以降、同値関係は \simで表します。

そして、集合 Aに対して、 [x] = \{y\in A | x\sim y \}と定義します。この [x]は、 xと同値な Aの要素の集合で、同値類と呼ばれます。そして、 [x]からなる集合 \{[x]|x\in A\}を Aに対する商集合といい、 A/\simと表します。

剰余類

Gを群として、 Hを Gの部分群とします。そして、 Gの元 g, g'に対し、 g^{-1}g'\in H\Rightarrow g\sim g'であるとします。この二項関係が同値なことは非自明ですが、がんばって説明してください。そしてあなたが大嫌いな人に高らかに自慢しましょう。

ここで、 [ g]がどういう集合かを考えてみたいと思います。

まず、 g^{-1}g'\in H \Rightarrow \exists h\in H\, \mathrm{s.t.}\, g^{-1}g' = h \Rightarrow g' = ghとなります。同値類 [g]の定義は、 gと同値なものの集合という定義でした。 g'は [g]の要素であるので、 [g] = \{g, gh_1, \cdots\}であることがわかります。これは、 Hの各要素の左から gをかけたものを集めたもの、と解釈できます。そして、この [g]を gHと表し、 gHを要素とした商集合は G/H = \{gH|g\in G\}と書かれ、左剰余類と呼びます。右剰余類も似たように定義できます。 G/Hの要素数は、 Gにおける Hの指数、といいます。

例として、 [0, 28)の整数の集合に \bmod 28の足し算の構造を入れたものを考えます。これは群であることがわかると思います。この部分群として、集合 \{0, 7, 14, 28\}に \bmod 28の足し算の構造を入れたものを考えられます。すると、 [i] = \{i+h|h\in \{0, 7, 14, 28\}\}(0\le i\le 6)として、左剰余類は \{[0], [1], [2], [3], [4], [5], [6]\}です。

ところで、 g_i, g_j\in Gという二つの元が同値でなかったとき、 [g_i], [g_j]の関係性はどうなるでしょうか。1300234241で割ったあまりを求めてください。具体的には、この二つの集合の共通部分はどうなるでしょう。そこで、 [g_i], [g_j]が共通部分をもつ、つまり、 h_k, h_l\in Hについて g_ih_k = g_jh_lなる要素があると仮定します。 h_lは部分群の要素なので、 Hに逆元を持ちます。よって g_ih_kh_l^{-1} = g_jです。そして、 h_kh_l^{-1}という積はハンターハンターなので H\times Hなので、 Hのある元です。この元を h_kh_l^{-1} = h\in Hとおくと、 g_j = g_ihですので、 g_i^{-1}g_j\in Hです。 g_i^{-1}g_j\in H\Rightarrow g_i\sim g_jですが、これは g_i, g_jが同値でないという仮定に反します。よって、 [g_i], [g_j]は共通部分をもちません。

ラグランジュの定理

一応、準備は整ったので、本題のひとつ「ラグランジュの定理」の証明と応用に移ります。

定理の内容

Gを有限群とし、 Hを Gの部分群としたときに、 Hの位数は Gの位数の約数である。より具体的には、 Gの位数は G/Hの要素数（ Gにおける Hの指数）と Hの位数を掛け合わせたものである。

証明

g\in Gとします。 Hの位数を mとして、 H = \{h_1, h_2, \cdots, h_m\}とします。すると、左剰余類 gH = \{gh_1, gh_2, \cdots, gh_m\}です。また、 f: H\rightarrow gHを f(h) = gh\, (h\in H)という写像とします。

まず、任意の1\le i, j\le mの h_i, h_j\in Hについて f(h_i) = f(h_j) \Leftrightarrow gh_i = gh_jであれば、両辺に g^{-1}を作用させることで h_i = h_jが言えるため、 fは単射です。

加えて、 gHは Hの要素に左から gをかけて作った集合なので、 f(h)=gh\in gH\Rightarrow \exists h\in Hであるはずです（そうでなければ私が病み散らかしてしまいます）。よって、 fは全射です。

以上より、fが全単射であることがわかりました。よって Hと gHが同型であることから、 gHの要素数は mです。

さて、 Gの互いに同値でないものを集めた集合を、要素数を kとして \{a_1, a_2, \cdots a_k\}とします。 Hとはそもそも Gの部分集合ですから、 a_iHのどの元も Gに含まれています。そして、類別を行うときに全ての元の対で同値関係を決めたため、 Gの元であってどの a_iHにも含まれていないものはありません。つまり、 Gのどの元もa_iHのどこかに属しています。したがって、 G = a_1H \cup a_2H\cup\cdots \cup a_kHです。

また、先に示したように、任意の iについて a_iHの要素数は mです。さらに、剰余類のところで見たように、任意の 1\le i, j\le k, i\neq jについて a_iH\cap a_jH = \emptysetです。 Gの位数は包除原理を用いるまでもなくただそれぞれの剰余類の要素数を足し合わせれば良いです。

以上より、 Gの位数は kmであるので、定理は示されました。

応用

では、ラグランジュの定理を利用して何ができるかをみていきたいと思います。

巡回群について

g\in Gが生成する巡回群 \langle g\rangleは Gの部分群です。これにラグランジュの定理を適用すると、 \langle g\rangleの位数は Gの位数の約数であることが言えます。このことから、 g^{| \langle g\rangle|}=eより g^{|G|} = eです。

フェルマーの小定理

フェルマーの小定理とは、素数 pと、 pと互いに素な整数 aについて、 a^{p-1}\equiv 1\bmod pであるというものです。二項定理を用いた証明が簡単ですが、ラグランジュの定理を用いても示せます。

pが素数であることから 1, 2, \cdots, p-1のすべての数は pと互いに素です。このことは、 \bmod pでの逆元の存在を保証します。 [i]を、 iを含む剰余類とすると、 pの倍数を除いた集合は \bmod p上の積で群となり、 \{[1], [2], [3], \cdots, [p-1]\}と類別されます。この位数は p-1なので、 a^{p-1}\equiv 1\bmod p\,(1\le a\le p-1, m\in \mathbb{Z})です。

フェルマーの小定理はよく知られているようにさらに一般化できます。一般の2以上の整数 nと、 nと互いに素な整数 aについて、 n以下の正の整数で nと互いに素なものの個数 \varphi(n)を用いて a^{\varphi(n)}\equiv 1\bmod nです。これも同じようにラグランジュの定理を用いて示せます。ただし、今度は n未満の正の整数で nと互いに素ではないものがあるかもしれないので、それを除いた剰余類で群を作らなくてはなりません。これが指数の \varphi(n)の由来です。

小休止

あー疲れた。群論つら。ちなみにこの記事をここまで書くのに3~4日くらいかかりました。

ところでおすすめの原宿系ファッションのお店の情報は随時お待ちしております。どんな方法でも情報送ってくれたら僕が喜びます。

離散対数

問題

正の整数 Mと、 0以上 M-1以下の整数 X, Yが与えられる。 X^{K}\equiv Y\bmod Mを満たす最小の非負整数 Kが存在すればそれを求め、存在しなければその旨を報告せよ。 (1\le M\le 10^{9})

Baby-step giant-step（BSGS）

\bmod Mで Xに逆元があるとき、すなわち、 Xと Mが互いに素のときは、baby-step giant-step（略してBSGS）というアルゴリズムが使用できます。

Mと Xが互いに素であれば、剰余類 [X]はフェルマーの小定理の証明に出てきた群の要素です。よって、問題の等式を満たす Kが存在するなら、それは \varphi(M)-1以下の非負整数のどれかです。しかし、 1\le M\le 10^{9}なので、まともに全探索したら間に合いません（無理矢理高速化したりするとコンピュータが泣いてショートしてしまいます）。

そこで、 m = \lceil \sqrt{\varphi(M)}\rceilとして、 K = im + j\, (0\le i, j \le m-1)とおいてみます。すると、問題の等式はX^{im+j}\equiv Y \Leftrightarrow X^{j}\equiv Y(X^{-m})^{i}と変形できます。ここで、 Xと Mが互いに素であることから、 Xに逆元が存在することを使いました。

ここで、あらかじめ 0\le j\le m-1について v \equiv X^{j}となるような最小の jを求めておくことで、もし v \equiv Y(X^{-m})^{i}なる jが存在すれば im+jが答えになります。ちょっとわかりにくいですが、以下に示す実装の手順を見ればわかると思います。

【実装の手順】

入力： X, Y, M

出力： K（等式 X^{K}\equiv Y\bmod Mを満たす最小の非負整数、存在しなければ -1）

  1. m = \lceil \sqrt{M}\rceilを求める。
  2. b = X^{-m}とする。
  3. 辞書などを用いて 0\le j\le m-1について組 (X^{j}, j)を持っておく。
  4. 0\le i\le m-1について、 X^{j}=(Yb)^{i}なる jが存在すれば im+jを返して終了。
  5. 終了しなければ、 -1を返して終了。

実装上の注意点として、辞書に値を格納するときに同じkeyに対して異なる jが複数対応することがあり、そのとき jの昇順のループでは jの最大値で更新されてしまいます。 jの降順でループするか、すでにkeyが辞書に入っているときは値の更新をしない、などの方法で対応できます。

時間計算量、空間計算量ともに \mathrm{O}(\sqrt{M})です。

お知らせ

すみません、ここにPohlig-Hellman algorithmの解説を書く予定でしたが、ちょっと想定より長くなったので次回の自分に投げます。次回は

  • Pohlig-Hellman algorithm
  • 一般modでの離散対数の考察

について書く予定です。余裕があれば具体的な実装を書く予定ですが、分量が多くなればまた次回以降に回します。