Holograph1c Attract0r

日々の数学やプログラミングに関係する話。

今から15回の反復計算で円周率を20億桁求めます

みなさん、お手元に無限の精度*1の電卓のご用意はできていますでしょうか。

2つ変数を用意します。

f:id:mikan_alpha:20200920220608p:plain

初期条件

以下の漸化式に従って、1 / z_{15} を求めます。15回繰り返しましょう。

f:id:mikan_alpha:20200920221241p:plain

漸化式

そうしたらば、きっとあなたの手の中には20億個の見覚えのある数たちが列を成していることでしょう。

*1:もしなければ、10^10桁程度の精度があるものが望ましいです

項別微分っていつできる?

自分用の備忘録も兼ねて、項別微分について書こうと思います。

ある関数列 \{ f_n(x) \}_{n \geq 1} について、\sum_{n=1}^{\infty} f_n(x) という関数を考えます。

この関数を微分したいと思ったとき、普通は(微分の線形性的に)項別に微分しようとしてしまうと思うのですが、実はこれができないケースが存在します。
(これは微分の定義にも極限が絡んでいるために、極限の交換が発生するためです)

具体的には、以下のような定理が知られています。

定理 1.
区間 I 上で微分可能でかつその導関数が連続である(C^1 級)関数列 \{ f_n(x) \}_{n \geq 1} に対して、

(i) \sum_{n=1}^{\infty} f_n(x)IS(x)各点収束

(ii) \sum_{n=1}^{\infty} f_n'(x) が IT(x)一様収束

の2つを満たすとき、S(x) もまた I 上の各点で連続で微分可能になり S'(x) = T(x) が成立する。

すなわち、\frac{d}{dx} \sum_{n=1}^{\infty} f_n(x) = \sum_{n=1}^{\infty} \frac{d}{dx} f_n(x) となる。(微分と総和が交換可能)

今日はこの証明のために、一様収束などの話をしたいと思います。

一様収束とは?

まず、関数の一様normというものを定義します。

定義 1.
ある区間 I 上で定義された関数 f について、

\begin{align}
\| f \| := \sup_{x \in I} | f(x) |
\end{align}

f一様normという。

このノルムは、一般的な絶対値と同じ性質を持ちます。(三角不等式など)

この一様normを用いて、一様収束という概念を定義します。

定義 2.
区間 I 上で定義された関数列 \{ f_n(x) \}_{n \geq 1} を考える。
このとき、\{ f_n(x) \}_{n \geq 1}I 上の関数 f(x) に一様収束することを、\lim_{n \to \infty} \| f_n - f \| = 0 により定義する。

ある区間 I 上の関数列 \{ f_n(x) \}_{n \geq 1} が 関数 f(x) に一様収束するとき、I 上の各点 x_0 において f_n(x_0)f(x_0) に収束します。
(これは直感的にも明らかだと思うので証明は省きます)

この I 上の各点 x_0 において f_n(x_0)f(x_0) に収束、が成り立つとき、\{ f_n(x) \}_{n \geq 1}f(x)各点収束するといいます。一様収束するならば各点収束しますが、逆は一般に成り立ちません。

また、一様収束する連続な関数列の収束先は、また連続になります。

一様収束の性質

色々な性質がありますが、ここでは冒頭の定理の証明に必要なものだけピックアップしました。

定理 2.
区間 I:[ a,b ] 上で一様収束する連続関数列 \{ f_n(x) \}_{n \geq 1} に対して、以下が成り立つ。

\begin{align}
\lim_{n \to \infty} \int_a^b f_n(x) \, dx = \int_a^b \lim_{n \to \infty} f_n(x) \, dx
\end{align}

定理 3.
区間 I 上で一様収束する連続な導関数列 \{ f_n'(x) \}_{n \geq 1} に対して、f_n(x) もまた各点収束するとき以下が成り立つ。

\begin{align}
\lim_{n \to \infty} \frac{d}{dx} f_n(x) = \frac{d}{dx} \lim_{n \to \infty} f_n(x)
\end{align}

定理2の証明

f(x) := \lim_{n \to \infty} f_n(x) とする。仮定より、f(x) もまた連続関数である。

\begin{align}
0 &\leq | \int_a^b f_n(x) \, dx - \int_a^b f(x) \, dx | \\
&= | \int_a^b \left( f_n(x) - f(x) \right) \, dx | \\
&\leq \int_a^b | f_n(x) - f(x) | \, dx \\
&\leq \int_a^b \| f_n - f \| \, dx \\
&= \| f_n - f \| (b - a)
\end{align}

従って、以下を得る。

\begin{align}
0 \leq | \int_a^b f_n(x) \, dx - \int_a^b f(x) \, dx | \leq \| f_n - f \| (b - a)
\end{align}

ここで各辺において n \to \infty とするとはさみうちの原理から以下が従う。

\begin{align}
\lim_{n \to \infty} | \int_a^b f_n(x) \, dx - \int_a^b f(x) \, dx | = 0
\end{align}

よって、定理が示された。

\begin{align}
\lim_{n \to \infty} \int_a^b f_n(x) \, dx = \int_a^b f(x) \, dx
\end{align}

定理3の証明

g(x) := \lim_{n \to \infty} f_n'(x) とする。仮定より、g(x) もまた連続関数である。

x, x_0 \in I に対して、

\begin{align}
\int_{x_0}^x g(t) \, dt &= \int_{x_0}^x \lim_{n \to \infty} f_n'(t) \, dt \\
&= \lim_{n \to \infty} \int_{x_0}^x f_n'(t) \, dt \\
&= \lim_{n \to \infty} \left[ f_n(t) \right]_{x_0}^x \\
&= \lim_{n \to \infty} \left( f_n(x) - f_n(x_0) \right) \\
&= f(x) - f(x_0)
\end{align}

従って、以下を得る。

\begin{align}
f(x) = \int_{x_0}^x g(t) \, dt + f(x_0)
\end{align}

ここで、両辺を x を変数として微分する。

\begin{align}
\frac{d}{dx} f(x) &= \frac{d}{dx} \left( \int_{x_0}^x g(t) \, dt + f(x_0) \right) \\
&= g(x) \\
&= \lim_{n \to \infty} f_n'(x)
\end{align}

f(x) = \lim_{n \to \infty} f_n(x) から、定理が示された。

\begin{align}
\frac{d}{dx} \lim_{n \to \infty} f_n(x) = \lim_{n \to \infty} \frac{d}{dx} f_n(x)
\end{align}

冒頭の定理の証明

上で証明した定理3を使うと、すぐに示すことができます。

定理1の証明

\begin{align}
\frac{d}{dx} S(x) &= \frac{d}{dx} \sum_{n=1}^{\infty} f_n(x) \\
&= \frac{d}{dx} \left( \lim_{n \to \infty} \sum_{k=1}^{n} f_k(x) \right) \\
&= \lim_{n \to \infty} \frac{d}{dx} \sum_{k=1}^{n} f_k(x) \\
&= \lim_{n \to \infty} \sum_{k=1}^{n} f_k'(x) \\
&= \sum_{n=1}^{\infty} f_n'(x)
\end{align}

よって、S'(x) = T(x) が成立。
ただし、2行目から3行目の微分と極限の交換に定理3を用いた。

具体例とか

最後に、項別微分不可能な例を紹介したいと思います。

以下のような関数を考えてみます。

\begin{align}
f(x) = \sum_{n=1}^{\infty} a^n \cos (b^n x)
\end{align}

a,b はそれぞれ正の定数で、a < 1b > 1ab > 1 とします。

級数判定法というのを用いると、この(関数項)級数は一様絶対収束することが分かります。これで各点収束の条件は満たされています。
また、各項は明らかに連続な関数なので、その収束先である f(x) は連続であることも分かります。

あとは、f_n(x) := a^n \cos (b^n x) としたとき、\sum_{n=1}^{\infty} f_n'(x) が一様収束するかが重要になります。

実際に計算してみると f_n'(x) = - a^n b^n \sin (b^n x) となります。ここで、ab > 1 の仮定から、これの総和は明らかに発散してしまいます。

従って、導関数列の和は一様収束しないので項別微分不可能ということになります。

次回予告とか

夏も真っ盛りで、暑苦しい日々が続いていますね。私も毎日溶けそうになっています。

まあ実際はほとんど家の中でゲームしてるんですけど……。

次回は、ちょっと計算機科学よりの内容か、自己同形数の話を書いてみるか、どっちかかなあという気分です。全然予定は立ててないので確定ではないです。

一様収束の話とか、自分も最近初めて学んだのですが、他の数学を学ぶ人の役に立ったらいいなと思います。というわけで、お疲れさまでした。

N以下の素数の逆数和の発散速度がlog log Nくらいなことのお気持ち

なんか適当に検索かけても出なかったのでお気持ちだけ書きます。(厳密な証明は英語Wikipediaあたりにあります)

以下、特に断りが無ければ p素数とします。

さて、以下のような事実が知られています。

\begin{align}
\sum_{p} \frac{1}{p} = \infty
\end{align}

素数の逆数和は、調和級数(全ての自然数の逆数和)と同様に発散します。感覚的には、平方数の逆数和は \frac{\pi^2}{6} に収束するので、素数は全ての自然数よりは少なくて、でも平方数よりはずっと多いといった感じです。

一方で、上の無限和は発散することには発散するのですが、その速度がめちゃくちゃ遅いことでも有名です。具体的には下のようになります。

\begin{align}
\sum_{p \leq n} \frac{1}{p} = O(\log \log n)
\end{align}

対数のさらに対数ということでかなり遅いことが見て取れます。気になる方は、お手元の電卓か、Wolfram|Alphaで大きい数で計算してみましょう。
ざっくり計算したければ、N という整数の対数 \log N はだいたいその桁数くらいなので、それを2回すると良いです。かなり小さい数字になることが分かると思います。

お役立ち情報

この発散速度が役に立つ場面として、計算機科学における素数判定アルゴリズムの一つである、エラトステネスの篩の(時間)計算量があります。

エラトステネスの篩とは、ある与えられた整数 N 以下の素数を見つけるためのアルゴリズムです。詳しい話は端折りますが、計算量解析をしてあげると、このアルゴリズムの計算量は O(N \log \log N) 程度になります。愚直に整数を1つずつ N まで素数判定をすると O(N \sqrt{N}) になってしまうので、かなりの計算量の改善になっています。

お気持ち

それでは本題です。使う道具は、ゼータ関数オイラー積表示と関数 \log (1+x)テイラー展開、そして調和級数の発散速度が対数オーダーであることです。

\begin{align}
\zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s} = \prod_{p} \left( 1 - p^{-s} \right)^{-1}
\end{align}

\begin{align}
\log(1+x) = x - \frac{1}{2} x^2 + \frac{1}{3} x^3 - \frac{1}{4} x^4 + \cdots
\end{align}

\zeta(1) の対数をとります。

\begin{align}
\log \zeta(1) &= \log \left( \sum_{n=1}^{\infty} \frac{1}{n} \right) \\
&= \log \left( \prod_{p} \left( 1 - \frac{1}{p} \right)^{-1} \right) \\
&= - \sum_{p} \log \left( 1 - \frac{1}{p} \right)
\end{align}

総和の中身の \log ( 1 - 1/p ) をテイラー級数を使って展開します。

\begin{align}
- \sum_{p} \log \left( 1 - \frac{1}{p} \right) &= - \sum_{p}  \left( - \frac{1}{p} - \frac{1}{2} \cdot \left( \frac{1}{p} \right)^2 - \frac{1}{3} \cdot \left( \frac{1}{p} \right)^3 - \cdots \right) \\
&= \sum_{p}  \left( \frac{1}{p} + \frac{1}{2} \cdot \left( \frac{1}{p} \right)^2 + \frac{1}{3} \cdot \left( \frac{1}{p} \right)^3 + \cdots \right) \\
&= \sum_{p} \frac{1}{p} + \frac{1}{2} \sum_{p} \frac{1}{p^2} + \frac{1}{3} \sum_{p} \frac{1}{p^3} + \cdots
\end{align}

したがって、

\begin{align}
\log \left( \sum_{n=1}^{\infty} \frac{1}{n} \right) = \sum_{p} \frac{1}{p} + \frac{1}{2} \sum_{p} \frac{1}{p^2} + \frac{1}{3} \sum_{p} \frac{1}{p^3} + \cdots
\end{align}

が得られました。

左辺はもう既に完成されてますので、右辺に注目します。

1つ目の項以外は、それぞれ以下のように上から評価できるので、有限値に収束することが分かります。(\zeta(s)s > 1 で絶対収束するため)

\begin{align}
\sum_{p} \frac{1}{p} + \frac{1}{2} \sum_{p} \frac{1}{p^2} + \frac{1}{3} \sum_{p} \frac{1}{p^3} + \cdots \quad \leq \quad \sum_{p} \frac{1}{p} + \frac{1}{2} \zeta(2) + \frac{1}{3} \zeta(3) + \cdots
\end{align}

有限の値は収束発散には影響を及ぼさないので、有限値の部分を全てまとめて S とおいてあげると、

\begin{align}
\log \left( \sum_{n=1}^{\infty} \frac{1}{n} \right) = \sum_{p} \frac{1}{p} + S
\end{align}

となります。ざっくり、有限項の和でもこの関係が成り立っていると考えてあげると、以下が得られます。

\begin{align}
\sum_{p \leq N} \frac{1}{p} \quad \approx \quad \log \left( \sum_{n=1}^{N} \frac{1}{n} \right) \quad \approx \quad \log \left( \log N \right)
\end{align}

Cauchyの平均値の定理とTaylor展開

Cauchyの平均値の定理を用いてTaylorの定理を証明します。

Cauchyの平均値の定理

関数 g,h は閉区間 [a,b] 上で連続かつ開区間 (a, b) 上で微分可能であって、h(b) \not= h(a) かつ区間内の各点 x において (g'(x), h'(x)) \not= (0, 0) とする。このとき、

\begin{align}
\frac{g(b) - g(a)}{h(b) - h(a)} = \frac{g'(c)}{h'(c)}
\end{align}

を満たす c \in (a,b) が存在する。

証明

f(x) := (g(x) - g(a))(h(b) - h(a)) - (g(b) - g(a))(h(x) - h(a)) とおく。

このとき、f(a) = f(b) = 0 である。これは定義式に代入することで直ちに確かめられる。

ここで(Lagrangeの)平均値の定理から、

\begin{align}
\frac{f(b) - f(a)}{b-a} = f'(c)
\end{align}

を満たす c \in (a,b) が存在する。ここで、左辺は 0 となることに注意する。

f導関数

\begin{align}
f'(x) = g'(x) \cdot (h(b) - h(a)) - h'(x) \cdot (g(b) - g(a))
\end{align}

であるから、ここで x=c とすることで定理の式を得る。

Taylorの定理

区間 a, b 上で連続かつ開区間 (a,b) 上で n微分可能な関数 f に対して、

\begin{align}
f(b) = f(a) + f'(a)(b-a) + \frac{f''(a)}{2!} (b-a)^2 + \cdots + \frac{f^{(n-1)}(a)}{(n-1)!} (b-a)^{n-1} + \frac{f^{(n)}(c)}{n!} (b-a)^{n}
\end{align}

を満たす c \in (a,b) が存在する。

証明(n=2の場合)

平均値の定理より、n=1 のケースは既に示されている。

\begin{align}
f(b) = f(a) + f'(c) (b - a)
\end{align}

g(x) := f(x) + f'(x)(b-x), h(x) := (b-x)^2 とおく。

ここで、Cauchyの平均値の定理を適応すれば、

\begin{align}
\frac{f(b) - (f(a) + f'(a)(b-a))}{0 - (b-a)^2} = \frac{f''(c)(b-c)}{-2 (b-c)}
\end{align}

となる c \in (a,b) が存在する。

この式を整理することで n=2 のケースが示される。

また一般の n についても、同様に証明ができる。

Taylor展開

Taylorの定理を書いたのでわざわざ別に書くまででもないですが、Taylorの定理において関数が無限回微分可能で、剰余項が n \to \infty としたとき0に収束するならばTaylor展開可能です。

微分方程式の美味しい炊き方、そして黄金比を食べることによるその効果。

お久しぶりです。今日は数学の話をします。

今回は、この微分方程式を解いていきたいと思います。

\begin{align}
f'(x) = f^{-1}(x)
\end{align}

ある関数の導関数逆関数が等しい、と言っています。

どうでしょう、ぱっとすぐに解が思いつきますか?
一体どんな関数がこの方程式の解なのでしょうか。

タイトルでネタバレしてないか?

観察 1. f'(x) = f(x)

まずは簡単なものから考えてみましょう。

最初にこれを解いてみます。

\begin{align}
f'(x) = f(x)
\end{align}

これはすぐに解ける人も多いことでしょう。

何回微分しても変わらない、これは、指数関数の特徴ですね。
上の微分方程式の解は、f(x) = e^x です。

また、x に依存しない定数をかけておいてあげても微分には影響しないので、もっと一般的に、f(x) = Ce^xC は定数)が解となります。

ちょっと天の邪鬼な人は、f(x) = 0 も解であることに気がついているでしょう。ただ、自明な解はあまり面白くないので、今日はこれは無しにしておきます。以下も同様です。

ここで、f(x) = Ce^x としたときの逆関数 f^{-1}(x) についても見ておきます。

式変形してあげると、f^{-1}(x) = \log \frac{x}{C} です。
これはどうも、最初の微分方程式を満たしそうにありません。

観察 2. f''(x) = -f(x)

さて次はこれを見てみます。

\begin{align}
f''(x) = - f(x)
\end{align}

2回微分すると、自身の符号を反転したものに等しい関数を探せ、ということです。
これも、分かる人にはすぐ解かれてしまうでしょう。

例えば、f(x) = \sin x としてあげれば解の1つになります。また、別の人は f(x) = \cos x としたかもしれません。

これもまた一般化してあげると、f(x) = A \sin x + B \cos xA,B は定数)という形(三角関数の線形結合)の関数は全て解になります。

さて、また逆関数についても考えておきたいですが……線形結合した形を変形していくのは少し面倒です。

r = \sqrt{A^2 + B^2} とおいてあげれば、f(x) = r \left( \frac{A}{r} \sin x + \frac{B}{r} \cos x \right) ですから、三角関数を合成してあげて、f(x) = r \sin (x + \alpha) のような形になるでしょう。

これを x について解いてあげれば、f^{-1}(x) = \arcsin(\frac{x}{r}) - \alpha でしょうか。

あまり綺麗じゃない気がしますが、とりあえず逆三角関数であることだけ分かれば十分です。三角関数と逆三角関数は、名前こそ似ていても違う性質の関数であることに注意しておきます。

少なくとも、これも最初の微分方程式の解にはならなそうです。

観察から分かること

2つほど例を上げてみましたが、何か気づくことはありましたでしょうか。私は気が付きませんでした。

ある関数 f(x)f'(x) = f^{-1}(x) を満たす時、どんな条件が必要でしょうか。

それは、f'(x)f^{-1}(x)共に同じ性質種類であることです。
数学らしからぬふわっとした言葉ですが……。

つまりは、例えば f'(x)指数関数ならば同時に f^{-1}(x)指数関数でなければならなくて、f'(x)三角関数ならば同時に f^{-1}(x)三角関数でなければならないわけです。
しかしながら、この2つの関数についてはこれは有り得ないということを上で述べました。指数関数の逆関数は対数関数になってしまいますし、三角関数逆関数は逆三角関数です。

では、これが成り立つような関数にはどんなものがあるでしょうか?

 

私たちは、どうやら最も簡単な例を忘れているようです。

多項式関数

ここで、多項式関数について考えてみます。多項式と言いつつも、ここでは1項のみからなる関数を考えましょう。

f(x) = A x^rA は定数、r \not= 0)とおきます。指数は整数に限りません。

この関数の導関数は、明らかに f'(x) = rA x^{r-1} です。

逆関数はどうでしょうか?

少し計算すれば、f^{-1}(x) = (\frac{x}{A})^{1/r} となることが分かります。

そしてこの2つは、実は上であげた解となる条件を満たしています
どちらも多項式関数の形を保っているのがお分かりでしょうか。

計算のために、少し式を整理します。

\begin{align}
f'(x) &= (rA) x^{r-1} \\
f^{-1}(x) &= A^{-1/r} x^{1/r}
\end{align}

これでもう明らかです!

あとは係数部分、指数部分について r,A連立方程式を解いてあげましょう。

指数に着目してあげれば、

\begin{align}
r-1 = \frac{1}{r}
\end{align}

という式が出てきます。
これを整理してあげれば、

\begin{align}
r^2 - r - 1 = 0
\end{align}

です。おや、いやに見覚えがありませんか?
解の公式などでこれを解けば、

\begin{align}
r = \frac{1 \pm \sqrt{5}}{2}
\end{align}

となりました。

これは、黄金比 \phi とその共役無理数 \bar{\phi} です!

\begin{align}
\phi = \frac{1 + \sqrt{5}}{2} \quad , \quad \bar{\phi} = \frac{1 - \sqrt{5}}{2}
\end{align}

あとはこれを使って係数部分についても方程式を解いてあげましょう。

\begin{align}
\phi A = A^{-1 / \phi}
\end{align}

式を整理します。

\begin{align}
A^{1 + 1 / \phi} = \frac{1}{\phi}
\end{align}

ここで黄金比の性質 \phi^{-1} = \phi - 1 を使うと、

\begin{align}
A^{\phi} = \frac{1}{\phi}
\end{align}

こうです!後はもう簡単です。

\begin{align}
A = \sqrt[\phi]{\phi^{-1}}
\end{align}

これで微分方程式を解くことができました!

まとめ

今日はこんな微分方程式を扱いました。

\begin{align}
f'(x) = f^{-1}(x)
\end{align}

ここで解を多項式関数と仮定してあげると、解(の1つ?)である以下の関数が出てきました。

\begin{align}
f(x) = \sqrt[\phi]{\phi^{-1}} x^{\phi}
\end{align}

自分もこの微分方程式を解いていて「まさか黄金比が出てくるとは……」という感じでした。

微分方程式の解き方的には、他の解の存在とかが分からないので良くないのかもしれませんが、ともかく解の一例を与えることができました。

この微分方程式の解って一意なのでしょうか?微分方程式の解の一意性等について詳しい方がいたら、情報をお待ちしています。

以上、お疲れさまでした。

参考


A very interesting differential equation.

 

AtCoderで水色になりました

緑の色変記事を書いてから早くも3ヶ月ほどが経ちましたが……ABC168にて、無事水色Coderになりました!!

f:id:mikan_alpha:20200518105542p:plain

うれしい

順位は557位で、パフォーマンスは1781でした!これが初めての青perfです。

f:id:mikan_alpha:20200518110108p:plain

ABC168の戦績

さて色変記事なので「水色になるまでやったこと」的なのを書きたいのですが……やったことと言えばコンテストに出て、復習して、たまに水diffを解くくらいのものなのであまり参考にならない気がしています、申し訳のなさ(AtCoderではまだ440問くらいしかAC出してません)
最近はApexかMinecraftしかしていません、FAKEだろ

ただ言えることは、やっぱりコンテストに出るのって最大の精進だと思うのでCodeforcesとかでコンテストに出まくるといいと思います。AtCoderと必ずしも問題傾向が等しいわけではないので、AtCoderのレーティングに直結することは保証しかねますが、逆に言うとAtCoderでレートが上の人に勝つチャンスがあります
何もAtCoderだけが競プロでは無いので、競プロが好きな人こそCodeforcesはおすすめです。ある程度英語が読めれば十分なので。どうしても無理ならDeepLを使いましょう。

あとE8さんのこちらの精選100問はとてもおすすめです。既に読んで取り組んでいる方も多いと思います。 

qiita.com

まだ半分も終わっていないので自分も青になるために解いていきたいと思います……。

 

最後になりますが、vtuberにハマると競プロとの相性が悪すぎる(精進時間の確保的な意味で)ので気をつけたほうがいいです、はい。趣味が基本競プロと相性悪くて泣ける。(最近は追加で学業が競プロを侵食しています)

ちなみに私はホロライブ推しです。


【ホロライブ】だいたい清楚なデビュー初期と現在の配信を比較してみた【ホロろぐ001】

 

高校卒業までに青目指して頑張りたいと思います。

競プロで悩んだ計算量の話

競プロの精進をしていて、個人的に悩んだことがあったので書いておきます。
誰かの助けになれば、と思います。

※各問題の解法ネタバレ有り

ABC097 C - K-th Substring

問題概要:
文字列 s が与えられる。s の部分文字列のうち、辞書順で K 番目に小さいものを出力せよ。
|s| \leq 5000

N = |s| とします
基本方針としては、制約が小さいのでC++ならsetにでも突っ込んで順番に取り出せばよいです。setはデフォルトでソートされた状態で保持してくれますからね。

というわけで一番最初に提出したコードがこちら。

Submission #11947075 - AtCoder Beginner Contest 097

うん、ループは O(N ^ 2) だしsetへのinsertは対数オーダーだから通るはず……って言ってると、なぜかTLEに。さて、なぜでしょうか。

強い人に質問してみると、これ実は O(N ^ 3 \log N) になってるらしいとのこと。

どこで O(N) が発生しているかというと、文字列の比較でした。
文字列の比較なんてどこでした?という感じなのですが、よく考えてみるとsetへinsertするときに毎回発生しています。ソートされた状態で要素を保持するためには比較をしなければなりませんからね……。これが盲点でした。
ソートをする際、比較そのものの回数は O(\log N) なわけですが、その比較1回あたりに O(N) かかっていたわけです。(正確には比較する文字列の長さのmin?)

そんなわけで、もうちょい計算量を削らなければなりませんが、この問題の場合には K の制約が小さいおかげで、長さ6以上の部分文字列はsetへinsertしないことによりTLEを回避できました。

ACコード(1行追加しただけ):
Submission #11947828 - AtCoder Beginner Contest 097

ABC162 E - Sum of gcd of Tuples (Hard)

問題概要:
1 以上 K 以下の整数からなる長さ N の数列の全てについてgcdを求め、その和を計算せよ。
2 \leq N \leq 10 ^ 5, 1 \leq K \leq 10 ^ 5

基本方針は公式解説の通りで、1 \leq X \leq K に対して最大公約数が X となるような数列の数を求めてその和を求めます。

こちらが書いたコード。
Submission #11975132 - AtCoder Beginner Contest 162

注目してほしいのは、もちろんループの部分。

for (int x = k; x > 0; x--) {
    lint t = k / x;
    t = modpow(t, n);
    lint x2 = x + x;
    while (x2 <= k) {
        t -= c[x2];
        x2 += x;
    }
    t = ((t % MOD) + MOD) % MOD;
    c[x] = t;
    res += t * x;
    res %= MOD;
}

内側のwhileで、倍数の分の数を引いています。
ここで気になるのが、このwhileループ、計算量はどうなるの?ということです。 ACを出したときはよく分かっていないまま通ってしまったのですが、別にTLEとなってもおかしくなかったと思います。
ここで今一度、ループの実行回数を考えてみます。

whileループの回数は、コード中の変数 x の値に依存しています。具体的には、x の倍数の個数の分だけ、ループが実行されるはずです。

ということは、x = k のときは \frac{k}{k} 回、x = k-1 のときは \frac{k}{k-1} 回、…という風になるので(本当は小数点以下切り捨てですがあえてこう書きます)、全体でのループの回数は以下の和になるはずです。

\frac{k}{k} + \frac{k}{k-1} + \frac{k}{k-2} + \cdots + \frac{k}{2} + \frac{k}{1}

もうこう書いてしまったら分かりやすいのですが、これを k でくくります。
そうすると、こうです。

k \cdot \left( \frac{1}{k} + \frac{1}{k-1} + \frac{1}{k-2} + \cdots + \frac{1}{2} + \frac{1}{1} \right)

右の分数の和は第 K 調和数 H_K なので、これは \log K に近似することが出来て、最終的な計算量は O(K \log K) ということになります。

コンテスト中にここまで思考を詰めれる人、すごいな……。