2024-03-18

拡張GCDについて

スナネコです。拡張GCDについての質問に答えます。

問題：
以下のように実装された extgcdがある。
整数 $a,b$ が $(a,b)\neq(0,0)$ を満たし、extgcd(a, b) の返り値を $(x,y)$ としたとき、 $\max(|x|,|y|)\leq \max(|a|,|b|)$ が成り立つことを示せ。

def extgcd(a, b):
  # ax+by=gcd(x,y) を満たす (x,y) を返す
  if b == 0: return (1, 0)
  X, Y = extgcd(b, a % b)
  x = Y
  y = X - a // b * Y
  return (x, y)

extgcd(a, b) が「 $|ax+by|=\gcd(a,b)$ を満たす整数 $x,y$ を返すこと」の証明は省略します。

最初に補題を示しておきます。
補題★：
0 でない整数 $a,b$ と整数 $x,y,g$ について、 $|ax+by|=g$ かつ $|g|< |b|$ かつ $|x|\leq|b|$ が成り立つならば、 $|y|\leq |a|$ である。
証明：
$ax+by=\pm g$ を変形すると、 $|by|=|\pm g-ax|\leq|g|+|ax|$ となり、この両辺を $|b|$ で割ると $|y|\leq|g/b|+|ax/b|$ となります。
ここで $|x|\leq |b|$ の仮定から $|x/b|\leq 1$ なので $|ax/b|\leq |a|$ であり $|y|\leq |g/b|+|ax/b|\leq |g/b|+|a|$ です。
仮定から $|a|,|y|$ は整数で $|g/b|<1$ なので、小数部分を無視してよく $|y|\leq |a|$ が成り立つことがわかりました。

もとの問題の証明をします。
証明：
$a=0$ または $b=0$ のとき、extgcd(a,b) の返り値は $(1,0)$ か $(0,1)$ になることが実際にアルゴリズムの挙動を追うことでわかるのでOKです。
$a\neq 0$ かつ $b\neq 0$ のときを考えます。
$|x|\leq |b|$ かつ $|y|\leq |a|$ が成り立つことを示します。このことが示せれば $|x|\leq |b|\leq \max(|a|,|b|)$ かつ $|y|\leq |a|\leq \max(|a|,|b|)$ となって $\max(|x|,|y|)\leq \max(|a|,|b|)$ が成り立つことが言えます。
再帰の深さに関する帰納法で証明します。
まず再帰の深さが1回で終わるケース、つまり $a\%b=0$ のときを考えます。
このときは実際にアルゴリズムの挙動を追うことで $(0,1)$ が返されることがわかるのでOKです。
次に再帰の深さが2回以上になるケース、つまり $a\%b\neq 0$ のときを考えます。
このとき、帰納法の仮定を使いつつ示すべきことをまとめると
「 $|bX+(a\%b)Y|=\gcd(a,b)$ かつ $|X|\leq |a\%b|$ かつ $|Y|\leq |b|$ を満たす整数 $X,Y$ が与えられる。 $x=Y$ 、 $y= X - \left\lfloor\frac{a}{b}\right\rfloor * Y$ としたとき、 $|x|\leq |b|$ かつ $|y|\leq |a|$ となることを示せ」です。
$x=Y$ なので $|x|=|Y|\leq |b|$ は明らかです。
$a\%b\neq 0$ なので $|\gcd(a,b)|<|b|$ になることに注意すると、「 $a,b\neq 0$ かつ $|ax+by|=\gcd(a,b)$ かつ $|\gcd(a,b)|< |b|$ かつ $|x|\leq |b|$ 」が成り立っているので補題★を使うことができて、 $|y|\leq |a|$ が成り立つことがわかります。
よって帰納法により証明できました。

2023-08-25

Difficlutyについて

パークガイド「長くなったのでブログの方でお答えしました」https://t.co/5iCqTXMgkJ #マシュマロを投げ合おう https://t.co/2YaT6kUXWU pic.twitter.com/EDONHtZPHV
— 競技プログラミングをするフレンズ (@kyopro_friends) August 25, 2023

Difficultyに関するご質問を頂いたので私たち競プロフレンズとしての見解を述べておきます。
本来であればアライさんを筆頭とした、精力的に競プロをされているアニマルガールの方々に説明をしてもらった方がよいことかもしれませんが、kyopro_friendsとしての表明ということで、私、ミライが代表してお話しさせていただきます。
これから説明する内容はあくまでkyopro_friendsの見解であり、他のwriterの方々や、AtCoder社の見解ではないことに注意してください。

・Difficultyの逆転は良いことだと思うか
良いことではありませんが、避けられるものでもありませんし、必要以上に避けようとするものでもありません。
避けられるものではない、ということについてはまたあとでご説明します。
避けようとするものではない、というのは例えば、310ではD,Eの難易度が逆転していますが、この結果を知った上で再度同じ問題を出すとしても、私たちはこの2問の順序を入れ替えようとは思わないでしょうという意味です。

・Difficultyと配点のどちらを見るべきか
基本的にこの2つには強い相関があるので、どちらを見ても結果はほとんど同じになるはずです。ほとんどの問題では解く方の得意不得意次第で変わる程度の差しかないと思います。
難易度とdifficultyが大きく離れている問題に関しても基本的には配点の方を信じてほしいところですが、本当にまれにある313Fのような難易度想定ミスの問題や、ABC288のように難易度が高くなっているセットなどもあるので、difficultyを目安に解いた方が精進はスムーズだと思います。

さて、difficultyの逆転は避けられるものではない、という話をしましょう。

まずはコンテストでの立ち回りについて考えてみます。
ABCが始まって開始70分時点でE問題までをACしました。実力的に残り30分でFGを両方解くのは無理です。配点はF問題が500点でG問題が600点です。
G問題の方が配点が大きいので、G問題が解けるならG問題を解いた方が得です。つまり
・FGどちらも解けなければレート-20
・Fが解ければレート±0
・Gが解ければレート+20
となります。FとGどちらに手をつけますか？

もちろんコンテスト本番中には「Gが解ければレートがいくら上がるか」を正確に予測することは困難ですし、自分がF,G問題を解ける可能性を正しく見積もることも困難なので、どちらを選ぶのが正しいとも一概には言えません。

ではもしF問題が500点、G問題が525点で他は先程と同じ状況だとどうでしょうか？

一見すると全く同じ状況に見えますが、「writerは、F問題とG問題の難易度の差は小さいと思っている」という情報が増えています。
ということは、先程の状況よりもG問題を解こうとする人が増えるでしょう。

このように、配点を細かくしたことで問題間の(想定)難易度の差が示されることは、「1つ飛ばして次の問題を解く」ことにインセンティブを与えます。
したがって、仮に問題が正確に難易度順に並んでいたとしても、difficultyが逆転する可能性は以前に比べ高くなります。
また、中盤までにFGの逆転が起こってしまえば、E問題までを解いた人にとって「F問題を解いても順位が対して上がらない」という状況になるため、G問題に手を付ける理由が増し逆転は加速します。
よって、特に配点の差が小さい問題が並んでいる箇所では、今後もdifficlutyの逆転は頻繁に起こると考えられます。

最後に、ABC301以降で逆転が起きた回について確認しておきましょう

301～315(15回)
315 D > E
313 F > H (FGは配点が同じ)
312 E > G (EFは配点が同じ)
310 D > E
307 C > DE
(304 D>E) 差100未満

このうち明確に難易度想定ミスと言えるのは313Fくらいでしょう。
315Dも簡単ではありませんが、E問題とここまで大きな難易度差がある問題には思えませんし、307もD問題が比較的簡単だったので、CDが逆転することまでは理解の範疇ですが、E問題と逆転するような問題にはとても思えませんね。

2023-05-30

パークガイドが教える！　けもフレ入門編

パークガイド「気合を入れて書きました。是非読んでください」https://t.co/oPugj95QiL #マシュマロを投げ合おう https://t.co/BLH5Yfv3zp pic.twitter.com/axPvtf1Hor
— 競技プログラミングをするフレンズ (@kyopro_friends) 2023年5月30日

パークガイドのミライです。
けものフレンズ入門記事を書いて欲しいとのことでしたので、ご紹介をしていきます。現実世界での出来事・役割よりも作品そのものの説明に重きをおいて紹介します。
けものフレンズという作品自体の紹介ということで、物語の登場人物としての私ではどうにも難しいところがあるので、今日はメタな視点からお話しします。

要約

現行のコンテンツを追いたい→けものフレンズ３
私たち (@kyopro_friends) が準拠する世界観を知りたい→ネクソン版
そもそもけものフレンズってなんで流行ったんだっけ？→アニメ1期

けものフレンズとは

公式サイトの説明が一番わかり易いですね

この世界のどこかに造られた超巨大総合動物園「ジャパリパーク」。そこには世界中の動物が集められ、研究・飼育が行われていました。ところがある日、神秘物質「サンドスター」の影響で、動物たちが次々とヒトの姿をした「アニマルガール」へと変身！いつしか彼女たちは“フレンズ”と呼ばれるようになり、ジャパリパークの新たな主として、にぎやかに暮らすようになりました。しかしジャパリパークにはフレンズの登場とほぼ同時期に「セルリアン」と名付けられた謎の生物も出現するようになり……。

ジャパリパークを舞台に、アニマルガールのみなさんの冒険や日常を描いた物語がさまざまなメディアで展開されています。

主要コンテンツ

発表された順に並べています。

種別	タイトル	年代	通称	備考
ソシャゲ	けものフレンズ	2015～2016	ネクソン版
漫画	けものフレンズ ‐ようこそジャパリパークへ!‐	2015～2017	フライ版
アニメ	けものフレンズ	2017	アニメ1期
舞台	舞台「けものフレンズ」	2017	舞台1
ソシャゲ	けものフレンズぱびりおん	2018～2021	ぱびりおん	現在閲覧不可
舞台	あにてれ×=LOVE ステージプロジェクト「けものフレンズ」	2018	イコラブ版
アニメ	ようこそジャパリパーク	2018～2020	ようジャパ	現在閲覧不可
舞台	舞台「けものフレンズ」2 〜ゆきふるよるのけものたち〜	2018	舞台2
アニメ	けものフレンズ２	2019	アニメ2期
漫画	けものフレンズ２	2019～2020	漫画版２
舞台	舞台けものフレンズ「JAPARI STAGE!」～おおきなみみとちいさなきせき～	2019	JS
ソシャゲ	けものフレンズ３	2019～	３

この他にもソシャゲを中心にいくつかのコンテンツがありますが、主要なものは以上の作品群です。
重要な注意点として、けものフレンズの作品は、「私たちヒトが普通にいるパーク」を舞台にしたものと、「ヒトがいないパーク」を舞台にしたものの2種類のに大別され、この2つでは作品のテイストが大きく異なります。
ここで挙げた作品のうち、『ネクソン版』『フライ版』『３』、及び『ぱびりおん』のごく一部、そしてもう1作品(ネタバレのため詳細省略)が「ヒトがいるパーク」、それ以外が「ヒトがいないパーク」での物語になります。
「ヒトがいるパーク」は現代を舞台に、「ヒトがいないパーク」はそれらよりいくらか未来を舞台にしたものであると考えられています。
近年では、全ての物語が同一世界にあるとする考え方が一般的です。ただし、ゼルダの伝説シリーズのように、「大きな出来事の結果に基づいて分岐した世界での出来事である」とする解釈や、パワポケシリーズのように「"正史"は存在するが、それが作中で語られているとは限らない」とする解釈もあります。

各コンテンツについて

一部の作品は前の作品の直接の続編などであり、履修の順序に注意する必要があります。以下の個別解説で具体的に述べます。

けものフレンズ（ソシャゲ、『ネクソン版』、2015～2016）

私たち @kyopro_friends が準拠する物語です。
新米パークガイドである私ミライと主人公である"あなた"が、サーバルさんたちとともにパークを冒険します。
サーバルさんそっくりの姿をしたセルリアンとして生まれたセーバルさんと友達になるまでの物語であり、セルリアンの手からパークを取り戻す戦いの物語でもあります。
ゲームはすでにサービス終了していますが、有志による動画により、現在でもストーリーを知ることができます。

けものフレンズ ‐ようこそジャパリパークへ!‐（漫画、『フライ版』、2015～2017）

新米飼育員であるナナさんが、キタキツネさんはじめとするアニマルガールのみなさんに振り回される日常ものです。
「ヒトがいるパーク」の日常を描いた数少ない作品の1つです。

けものフレンズ（アニメ、『アニメ1期』、2017）

正体不明の人物かばんさんが、自分を知るため、パークを知るために、サーバルさんとともにパークを冒険します。
『ネクソン版』と繋がりがあるため、『ネクソン版』を知っていた方が物語をより深く理解できますが、知らなくとも最低限の理解はできる作りになっています。

舞台「けものフレンズ」（舞台、『舞台1』、2017）

オカピさんとサーバルさんが、ジャパリパークのアイドルグループ"PPP"のライバルとしてアイドルグループを結成してライブをします。
全体的に『アニメ1期』の色が濃いですが、直接の繋がりはありません。

けものフレンズぱびりおん（ソシャゲ、『ぱびりおん』、2018～2021）

現在閲覧不可
『フライ版』とは逆に「ヒトがいないパーク」の日常ものです。ショートストーリー(数百文字～数千文字)が2000本以上あります。特定の条件を満たしこれらのストーリーを解禁することがゲームの目的の1つでした。
サービス終了時に全てのストーリーを閲覧可能としたオフライン版がしばらく公開されていましたが、現在では入手不可能となっています。

あにてれ×=LOVE ステージプロジェクト「けものフレンズ」（舞台、『イコラブ版』、2018）

ツチノコさんとニホンオオカミさんを主人公に、アニマルガールのみなさんがお芝居に挑戦するお話です。お二人の終盤のやりとりには胸を打たれます。
リョコウバトさん、オーストラリアデビルさんなど、出演メンバーは意味深なチョイスとなっています。

ようこそジャパリパーク（アニメ、『ようジャパ』、2018～2020）

現在閲覧不可
『ネクソン版』を原作にしたショートアニメ(5分程度×36話)です。
おそらく尺の都合とは思いますが、ストーリーが大幅に簡略化されているので、時間が許すのであれば原作である『ネクソン版』の方を見ていただきたいところです。
もっとも、2021年のあにてれサービス終了に伴い、現在この作品を見る手段はないようですが……。

舞台「けものフレンズ」2 〜ゆきふるよるのけものたち〜（舞台、『舞台2』、2018）

ゆきやまに住むギンギツネさんキタキツネさんを主人公にした物語です。姉妹愛がテーマの1つでしょうか。
『舞台1』と同一世界観のようですが、特に前作を知らなくとも楽しめます。

けものフレンズ２（アニメ、『アニメ2期』、2019）

正体不明の人物キュルルさんが、おうちに帰るために、サーバルさんカラカルさんとともにパークを冒険します。
『アニメ1期』の直接の続編ですが、『アニメ1期』とは異なる雰囲気と強いメッセージ性を持っています。

けものフレンズ２（漫画、『漫画版２』、2019～2020）

『アニメ2期』の漫画化です。
後半ではアニメと異なる展開を見せ、アニメとは一部の重要な点が異なる結末を迎えます。

舞台けものフレンズ「JAPARI STAGE!」～おおきなみみとちいさなきせき～（舞台、『JS』、2019）

オオミミギツネさんとシベリアンハスキーさんを主人公に、これまでのけものフレンズ作品ではあまり扱われなかった展開を巧みなシナリオで扱います。
『舞台1』『舞台2』と同一世界観のようですが、特に前作を知らなくとも楽しめます。

けものフレンズ３（ソシャゲ、『３』、2019～）

『ネクソン版』より少し後のパークを舞台に、主人公である"あなた"が探検隊の隊長として、ドールさんとともにパークを冒険します。
ジャパリパークを巡って人々の思惑が交錯するシーズン1、新種のセルリアンとその秘密に迫るシーズン2、そして、新たな来園者を迎えて現在展開中のシーズン3と、物語は今なお続いています。
この作品は『ネクソン版』と『フライ版』の直接の続編です。最低限の要素は作中で説明されますが、『ネクソン版』の事前履修を強くおすすめします。またシーズン2に触れる際には『フライ版』の事前履修も強くおすすめします。

おわりに

ここでは紹介しませんでしたが、ゲームとして「ぱずるごっこ」「ピクロス」、ソシャゲとして「FESTIVAL」「KINGDOM」、ツールとして「あらーむ」、『３』に関連するショートアニメとして「ちょこけも」「カレンダレコード」、書籍として『アニメ2期』のノベライズ、Vtuberとして「けものフレンズ Vぷろじぇくと」など、けものフレンズにはこの他にも様々なものがあります。
この記事が皆さんがけものフレンズに触れるきっかけとなれば幸いです。

2023-04-24

除原理について

スナネコです。
競プロ界隈で"除原理"と呼ばれているものについて説明します。
この言葉は厳密な定義があるものではないので、これから説明することは単に「私はこうだと思っている」というお話として聞いてください。

定義

「 $g:=\zeta f$ と $x$ が与えられるので、 $f(x)$ を求めてください」という問題を、
メビウス変換の定義に基づいて $f(x)=\sum_{y\leq x}\mu(x,y)g(y)$ と求めるのが包除原理、
小さい方から順に $f(x)=g(x)-\sum_{y< x}f(y)$ と求めるのが除原理です。
ただし、シグマの中身を定義通りに計算するものだけでなく、和をとる順序を工夫したり状態をまとめたりするものもそれぞれ包除原理、除原理に含めることにします。

例

EDPC Y問題を考えてみます。
(1,1) から (H,W) への経路それぞれに対して、i 番目の壁を通らない、通る、わからないを表す o, x, ? の長さ N の文字列を対応させます。
また、「その文字列になるような経路の個数」をその文字列そのもので表すことにします。
N=3で考えます。求めるものは ooo です。

包除原理は定義通りに書くと ooo=???-(x??+?x?+??x)+(xx?+x?x+?xx)-xxx で、
除原理は ooo=???-oox-oxo-oxx-xoo-xox-xxo-xxx ですね。
このままだとどちらも項の数が $O(2^N)$ になってしまうので計算をまとめる必要があります。
ここでは詳しい説明はしませんが、
包除原理は ooo=(xが偶数個)-(xが奇数個) に
除原理は ooo=???-x??-ox?-oox にまとめるとうまく計算できます。
kiwamachanさんの解説やけんちょんさんの解説は包除原理で、snukeさんの解説やフェネックさんの解説が除原理ですね。読み比べて見てください。

例2

「 $1\leq x,y \leq N$ であって、 $\mathop{gcd}(x,y)=1$ を満たす組 $(x,y)$ の個数 $f(N)$ を求めよ」という問題を考えてみます。
包除原理で考えると、gcdが1の倍数になるケース、2の倍数になるケース、3の倍数になるケース、……を係数の辻褄をあわせながら足し引きすればいいので $\sum_{1\leq d \leq n}\mu(d)\left(\frac{N}{d}\right)^2$ になります。エラトステネスの篩のように計算すれば $O(N\log\log N)$ 、ディリクレ級数の積を考えると $\tilde{O}(N^{2/3})$ で計算できます。
除原理で考えると、全てのケースから、gcdが2になるケース、3になるケース、……を引けばいいので $N^2 - \sum_{2\leq d \leq N}f\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$ になります。これは定義通りにメモ化再帰で計算すると $O(N\log N)$ 、商列挙を組み合わせると $O(N^{3/4})$ 、小さい $x$ での $f(x)$ を別の方法で先に計算しておくことにすれば $\tilde{O}(N^{2/3})$ で計算できます。

計算量

例2で見たとおり、特別な工夫をしない限り除原理の方が包除原理よりも計算量が悪くなることが多いです。
一方で、除原理は排反な事象に分割して計算できればあとは簡単なのに対し、包除原理だと足すのか引くのか何倍するのか係数を考えるところが難しいです。
包除原理がよくわからないと思うときは、いったん除原理で式を立ててゼータ関数(順序集合のなす束の構造)を見極めてからメビウス関数を書き下して包除原理に書き直す、ということもできます。私は包除原理を考えるのは苦手なので、最初は除原理で考えて、計算量が間に合わないときだけ包除原理で考え直すことにしています。

高速メビウス変換

最初の例での ooo=???-x??-ox?-oox という計算方法は、集合に関する高速メビウス変換と同じです。
（もしわからないなら、3次元の立方体の8つの頂点を000,……,111に対応させて、高速ゼータ変換/高速メビウス変換の外側のループを1回進めるごとに、どこの値がどう足し算/引き算されているか図に描いてみてください）
さらに、集合に関する高速メビウス変換は多次元累積和の逆操作でイメージできるものです。（さっきの図は描いてみましたか？）
つまり、多次元累積和の逆変換と高速メビウス変換とこのタイプの除原理はだいたい同じものということになりますね。
約数ゼータ変換/約数メビウス変換も、「素数が2と3しかない世界」を考えてみると、2次元累積和とその逆変換になっていることがわかります。

除原理？

ここまでで説明した私のイメージでは、ABC162E「Sum of gcd of Tuples (Hard)」の公式解説のO(NlogN)解法は除原理だと思うのですが、実際には「約数包除原理」と呼ばれることが多いようです。
"除原理"という言葉を私とは違う意味で使う子もたくさんいるみたいですね。

2023-04-12

Grundy数がXORになることの証明

スナネコです。
Grundy数がxorで計算できることの説明をしてほしいとサーバルさんに頼まれたのでします。

ざっと検索してみましたが、「Grundy数が0でないとき0にできる」くらいしか証明が書いてあるものは見つけられないですね。これではXORになることの説明には不十分です。
（追記）この記事にかかれていました　Grundy数（Nim数, Nimber）の理論（追記終わり）

私が以前 Nimber の解説を書いたときも、方針だけ書いて細かい証明は書いていませんでした。
せっかくなのでちゃんと証明しましょう。
ここから先では XOR を $\oplus$ で表すことにします。

"山"の個数に関する帰納法から、"山"が2個のときを示せば十分です。
grundy数が a と b の2つの"山"を合わせた局面のgrundy数を $g(a,b)$ と表すことにします。
定義から分かっていることは $g(a,b) = {\mathrm mex}(\{g(a,b') \mid b'< b\}\cup \{g(a',b) \mid a'< a \})$ で、示したいことは $g(a,b)=a\oplus b$ です。
$a,b$ に関する帰納法で示します。 $a+b$ の小さい順だと思ってもいいですし、 $(a,b)$ の辞書順だと思ってもいいです。
定義から $g(0,0)={\mathrm mex}(\emptyset)=0$ であり、 $(a,b)=(0,0)$ のとき成り立っています。
$(a,b)\neq(0,0)$ のとき帰納法の仮定から
$g(a,b) = {\mathrm mex}(\{a\oplus b' \mid b'< b\}\cup \{a'\oplus b \mid a'< a \})$ です。なので示すべきことは
・ $a\oplus b$ は $\{a\oplus b' \mid b'< b\}\cup \{a'\oplus b \mid a'< a \}$ に含まれない
・0 以上 $a\oplus b$ 未満の任意の数は $\{a\oplus b' \mid b'< b\}\cup \{a'\oplus b \mid a'< a \}$ に含まれる
の2つです。
1つ目は $x\neq b \Rightarrow a\oplus x \neq a \oplus b$ と $x\neq a \Rightarrow x\oplus b \neq a \oplus b$ が成り立つので言えます。
2つ目を示すのはちょっと大変です。
示したいことは「 $n < a \oplus b$ ならば『 $b \oplus n < a$ または $a \oplus n < b$ 』」です。
$b \oplus n < a$ なら、 $a'=b \oplus n$ とすれば $a'\oplus b=n$ にできて、 $a \oplus n < b$ なら、 $b'=a \oplus n$ とすれば $a\oplus b'=n$ にできますからね。
証明には、「 $x < y$ ⇔ $x\oplus y$ の最上位bitの位置では、 $x$ のbitが0で、 $y$ のbitが1」という事実を使います。
$n < a \oplus b$ ⇔ $a\oplus b\oplus n$ の最上位bitの位置では、 $n$ のbitは0で、 $a \oplus b$ のbitは1。つまり(n,a,b)は(0,0,1)か(0,1,0)。
$b\oplus n < a$ ⇔ $a\oplus b\oplus n$ の最上位bitの位置では、 $b\oplus n$ のbitは0で、 $a$ のbitは1。つまり(n,a,b)は(0,1,0)か(1,1,1)。
$a\oplus n < b$ ⇔ $a\oplus b\oplus n$ の最上位bitの位置では、 $a\oplus n$ のbitは0で、 $b$ のbitは1。つまり(n,a,b)は(0,0,1)か(1,1,1)。
となることから、たしかに「 $n < a \oplus b$ ならば『 $b \oplus n < a$ または $a \oplus n < b$ 』」が成り立っていることがわかりました。

これで $g(a,b)=a\oplus b$ が証明できました。

2022-11-15

実装例について

パークガイドのミライです。

ABCの解説の実装例について議論があるようなので、私たち競プロフレンズとしての見解を述べておきます。
本来であればアライさんを筆頭とした、精力的に競プロをされているアニマルガールの方々に説明をしてもらった方がよいことかもしれませんが、kyopro_friendsとしての表明ということで、私、ミライが代表してお話しさせていただきます。

まず最初にお伝えしておくべきこととして、解説のクオリティーに対してAtCoder社からの要求は一切ありません。どのような解説を書くかは完全にwriterの裁量に任せられています。
なので、これから説明する内容はあくまでkyopro_friendsの方針であり、他のwriterの方々や、AtCoder社の見解ではないことに注意してください。

本題に入ります。
前提として、C問題以降の解説に実装例は不要だと私たちは考えています。
writerは、わざわざ実装を読まなくとも理解できるような解説を書くべきですし、読者は、そのような解説を読んで理解できないのであれば理解できない部分は自分で考えたり調べたりするべきです。
解説がPDFで書かれていた時代には実装例を載せることはほとんどありませんでしたが、「作問に使ったwriter解を公開して欲しい」という声は何度も聞こえていました。現在私たちが解説に実装例を載せているのはこの声を受けてのものです。
動機がこのようなものであるため、基本的に実装例にはwriter解がそのまま貼ってあります。writer解というのはつまり「私たちにとって読みやすいコード」です。
なので、
「人に見せるコードは整えて欲しい」に対する答えは「人に見せるコードではありません」
「人に見せるコードを用意して欲しい」に対する答えは「その必要はないと考えています」
「人に見せるコードでないなら解説に載せないで欲しい」に対する答えは「その通りだと思います」
になります。

現在のスタイルが不評なようであれば、解説にwriter解へのリンクを載せるのはやめ、実装を読みたい人が各自でkyopro_friendsの提出一覧を漁るようにするのが妥当な解決策だと思います。
みなさんはどうお考えですか？

2022-03-08

等比数列に単項式が掛かっているものの和を計算するyosupo judgeの問題

スナネコです。yosupo judgeにある $\sum_{i=0}^{n-1}r^i i^d$ を解きました。
https://judge.yosupo.jp/problem/sum_of_exponential_times_polynomial

yosupo judge の github を見ると参考文献が載ってます。
https://github.com/yosupo06/library-checker-problems/issues/98

そういうわけで、この記事は Min_25 さんによる以下の記事を私の理解で書き直したものです。
ひとによってはオリジナルの方がわかりやすいかもしれません。
https://web.archive.org/web/20170622095007/http://min-25.hatenablog.com/entry/2015/04/24/031413

話の流れはこうです。
・数列 $a_n=\sum_{i=0}^{n-1}r^i i^d$ を漸化式で表す
・係数行列の固有多項式を求める
・固有多項式で割ったあまりを explicit な式で表す

ステップ1. 漸化式を出す

本題に入る前に、単項式 $n^d$ の値を $n$ に関する漸化式で求めるにはどうすればよいか考えてみましょう。
全ての $i\leq d$ についての $n^i$ を持てばできますね。具体的には、 $(n+1)^d$ が $n^i$ と二項係数を使って書けることから言えます。例えば
$\pmatrix{ 1 & 0 & 0 & 0\\ 1 & 1 & 0 & 0\\ 1 & 2 & 1 & 0\\ 1 & 3 & 3 & 1\\ } \pmatrix{ 1\\ n\\ n^2\\ n^3\\ }= \pmatrix{ 1\\ (n+1)\\ (n+1)^2\\ (n+1)^3\\ }$
となります。二項係数を並べたこの下三角行列のことを下三角パスカル行列というそうです。 $n\times n$ の下三角パスカル行列をwikipediaでは $L_n$ と表しているので、ここでもそう書くことにします。

本題に戻ります。ここまでわかれば元の問題を漸化式で書くのは簡単です。
$\pmatrix{ 1r & 0r & 0r & 0r & 0\\ 1r & 1r & 0r & 0r & 0\\ 1r & 2r & 1r & 0r & 0\\ 1r & 3r & 3r & 1r & 0\\ 0 & 0 & 0 & 1 & 1 } \pmatrix{ r^n\\ r^n n\\ r^n n^2\\ r^n n^3\\ \sum_{i=0}^{n-1} r^i i^3 }= \pmatrix{ r^{n+1}\\ r^{n+1} (n+1)\\ r^{n+1} (n+1)^2\\ r^{n+1} (n+1)^3\\ \sum_{i=0}^{(n+1)-1} r^i i^3 }$
こうですね。左上の $(d+1)\times(d+1)$ の小行列を $L_{d+1}$ とまとめて書くことにして
$M:=\pmatrix{ & 0\\ rL_{d+1} & \vdots\\ & 0\\ 0\ 0\ \ldots \ 0\ 1 & 1 }$
$A_n:=\pmatrix{ r^n\\ r^n n\\ \vdots\\ r^n n^d\\ \sum_{i=0}^{n-1} r^i i^d }$
とします。求めたいものは $A_n$ の最後の成分であり、 $A_n=M^nA_0$ 、 $A_0=(1,0,0,\ldots,0)$ です。
ここまでくれば、線形漸化式を持つ数列の第 $n$ 項を求めるのと同じようにできそうな気がしてきましたね。

ステップ2. 係数行列の固有多項式を求める

なぜ固有多項式を求めたいのか一応説明しておきましょうか。一般に、行列 $A$ の固有多項式 $p_A$ は $p_A(A)=O$ を満たします。このことから、 $f(x)=(x^n \bmod p(x))$ とすると、 $M^n=f(M)$ となります。*1
そうすると、小さい次数の項しか残らないので、いろいろ話が簡単になるわけです。

それでは本題に戻ります。さっき求めた係数行列 $M$ の固有多項式 $p$ を求めましょう。
$M$ は下三角行列で、対角成分は $r$ が $d+1$ 個、 $1$ が $1$ 個なので、 $p(x)=(x-r)^{d+1}(x-1)$ です。……はい、これだけです。

詳しくは省略しますが、ここまでわかれば、線形漸化式のときと同じように高速きたまさ法で $O(d\log d\log n)$ で元の問題を解くことができます。

ステップ3. 固有多項式で割ったあまりを求める

きたまさ法を使った解法では $f(x)$ を求めるために $O(d\log d\log n)$ の時間がかかりますが、実はこの問題の場合は $f(x)$ の係数を直接計算することができるので、よりよい計算量で解くことができます。
ここから先はひたすら計算が続きます。ちゃんとついてきてくださいね。

$f$ の形から、定数 $c$ と、 $d$ 次以下の多項式 $g$ を用いて、 $f(x)=c(x-r)^{d+1}+g(x)$ と書けます。この $c$ と $g$ を頑張って求めます。

$g$ について

この $g$ は要するに $(f\bmod (x-r)^{d+1})=(x^n\bmod (x-r)^{d+1})$ です。 $n< d+1$ なら自明です。そうでないときは $g(x+r)=( (x+r)^n\bmod x^{d+1})=\sum_{k=0}^{d}\binom{n}{k}r^{n-k}x^k$ になるので、あとは $x$ をずらすだけですね。
$\begin{eqnarray} g(x)&=&\sum_{k=0}^{d}\binom{n}{k}r^{n-k}(x-r)^k\\ &=&\sum_{k=0}^{d}\binom{n}{k}r^{n-k}\sum_{i=0}^{k}\binom{k}{i}(-r)^{k-i}x^i\\ &=&\sum_{i=0}^{d}r^{n-i}x^i(-1)^i\sum_{k=i}^{d}\binom{n}{k}\binom{k}{i}(-1)^k \end{eqnarray}$
$k$ に関する和が残ってしまいました。二項係数なので階乗に書き下して変形してもよいですが、組合せ的な言い換えをした方がわかりやすいことも多いです。今回、二項係数の積の部分は「 $n$ 個から $k$ 個選び、さらに選んだ $k$ 個から $i$ 個選ぶ方法の数」だと思えます。これは「 $n$ 個から $i$ 個選び、さらに選ばなかった $n-i$ 個から $k-i$ 個選ぶ方法の数」と等しいので、こういうふうに変形できますね。
$\begin{eqnarray} & &\sum_{k=i}^{d}\binom{n}{k}\binom{k}{i}(-1)^k\\ &=&\sum_{k=i}^{d}\binom{n}{i}\binom{n-i}{k-i}(-1)^k\\ &=&\binom{n}{i}(-1)^i\sum_{k=0}^{d-i}\binom{n-i}{k}(-1)^k\\ &=&\binom{n}{i}(-1)^i\binom{n-i-1}{d-i}(-1)^{d-i} \end{eqnarray}$
最後の等号は、 $\binom{x}{y}=\binom{x-1}{y-1}+\binom{x-1}{y}$ を使って $\binom{n-1}{*}$ を全部 $\binom{n-1-i}{*}$ に置き換えて計算するとわかります。そういうわけで、 $g$ は具体的に次のような形で書けました。
$g(x)=\sum_{i=0}^{d}r^{n-i}(-1)^{d-i}\binom{n}{i}\binom{n-i-1}{d-i}x^i$

$c$ について

すでに $g$ が求まっているので、適当な値を $x$ に代入して、恒等式 $f(x)=c(x-r)^{d+1}+g(x)$ を $c$ に関する $1$ 次方程式だと思って解くだけですが、適当な $x$ に対して $f(x)$ を求めるのはそんなに簡単ではないので、簡単に計算できる $x$ を探す必要があります。どんな $x$ なら簡単に求められるでしょうか。実は、 $(x-1)|p(x)$ と $f(x)=(x^n\bmod p(x))$ をあわせると $f(1)=(f(x)\bmod (x-1))=(x^n \bmod (x-1))=1$ になるので、 $x=1$ を選べばよいです。このとき、 $1=c(1-r)^{d+1}+g(1)$ を得ます。 $r\neq 1$ ならこれで話は終わりです。 $r=1$ のときは $c$ の係数が $0$ になってしまうので、これではうまくいきません。 $r=1$ のとき、元の問題は単なる単項式の和なので、最初から場合分けをして取り除いておいてもいいんですが、せっかくなのでいまやったのと同じような方法で求められないか考えてみましょう。

$r=1$ のときに $c$ を求める

$f(x)=c(x-r)^{d+1}+g(x)$ という恒等式を思い出すと、 $g$ は $d$ 次以下だったことから、 $c$ は $f$ の $x^{d+1}$ 次の係数と等しくなります。 $r=1$ のとき $p(x)=(x-1)^{d+2}$ なので、「 $f(x)=(x^n \mod (x-1)^{d+2})$ の $d+1$ 次の係数を求めよ」という問題になりますが、上で $g$ を求めるときにやったのと全く同じ議論で $\binom{n}{d+1}$ になることがすぐわかります。

ステップ4. $a_n$ を求める

これも線形漸化式の第 $n$ 項を求めるときと全く同じですが、一応説明しておきましょう。
求めたいものは $A_n=f(M)A_0$ の最後の成分でした。 $f(x)=\sum_{i=0}^{d+1}c_ix^i$ とおくと、 $A_n=\sum_{i=0}^{d+1}c_iM^iA_0=\sum_{i=0}^{d+1}c_iA_i$ となります。つまり $a_n=\sum_{i=0}^{d+1}c_ia_i$ なので、 $a_0,\ldots,a_{d+1}$ を予め愚直に求めておけば計算できます。

計算量解析

計算する必要があるものは $a_0,\ldots,a_{d+1}$ と $f(x)=c(x-r)^{d+1}+g(x)=\sum_{i=0}^{d+1}c_ix^i$ です。
計算量は $O(d \log {\rm MOD}+\log n)$ に見えますが、線形篩・線形逆元列挙・N乗数の高速列挙・二項係数の差分計算など、細かい工夫をたくさんすると $O(d+d \log {\rm MOD}/\log d+\log n)$ になります。

私の実装例はこれです。線形篩ではなく普通の篩を使っているので、計算量は $O(d\log\log d+d \log {\rm MOD}/\log d+\log n)$ になっています。
https://old.yosupo.jp/submission/81322

ラグランジュ補間でも解けるという話を聞いたのですが、私はまだわかっていません。誰か教えてくれますか？

*1:二項演算 mod は、剰余環に送るのではなく、割り算して余りを求める演算子とします

要約