mattyuuの数学ネタ集

世界は数式で溢れている

マーライオン再び〜ラグランジュの未定乗数法を用いて〜

ガロア、いいですね。私は男なのでガロアに恋心を抱くことはありませんが、ガロア理論の基本定理には恋心を抱いてしまいそうです。これまで数論を勉強してこなかったので、ガロア理論が大いに使われていそうな数論を勉強するのが今の楽しみです。

はじめに

以前下記の記事でマーライオンが一番遠くまで水を飛ばすことのできる口の角度を求めました。

mattyuu.hatenadiary.com

前回は空気抵抗を無視しましたが、今回は空気抵抗も考慮して問題を解きます。空気抵抗が入るだけで、問題がぐっと難しくなり、前回と同様の方法では解けませんでした。しかし、ラグランジュの未定乗数法という手法を用いることで、最適解\theta_{\rm max}の満たすべき方程式を得ることができましたので、問題を解きながらラグランジュの未定乗数法の強力さを紹介したいと思います。

この記事では前回必要だった高校物理の力学の知識、高校の数Ⅲで習う三角関数の微分の知識に加えて、微分方程式の基礎(変数分離法のみ)、偏微分、指数関数、対数関数の微積分の知識も前提としてしまいます(大学初年度レベルです。)。逆に言うと、数学の知識が増えると解ける問題の範囲も広がるということなので、数学を勉強するモチベーションにもなりますね。

前提知識がない方は、前回同様フジテレビのガリレオで福山雅治さんはこんな数式を書いていたんじゃないか？と想像しながら、数式を流し読みしていただけたら面白いのではないかと思います。

問題設定

では、問題設定から行いましょう。

f:id:mattyuu:20161218123227j:plain
f:id:mattyuu:20161218123235j:plain

前回と同様に、水は大きさのない質点として計算していきます。水は常に連続的に噴射され続けていますが、時刻t=0で質量mの水分子が1粒発射されたものとして読んでください。また空気抵抗は速さに比例する大きさの力が、運動の方向と逆向きに働くものとします。

では前回と同じ方法で解けるところまで解いていきましょう。口から飛び出る水の速度を水平方向(x軸右方向)の速度v_xと、垂直方向(y軸上方向)の速度v_yに分けて、ニュートンの運動方程式を立てます。

\displaystyle \begin{cases} m\frac{{\rm d}v_x(t)}{{\rm d}t} = -kv_x \rightarrow \frac{{\rm d}v_x(t)}{{\rm d}t} = -k^{\prime}v_x\\ \\ m\frac{{\rm d}v_y(t)}{{\rm d}t} = - mg -kv_y \rightarrow \frac{{\rm d}v_x(t)}{{\rm d}t} = -g -k^{\prime}v_y \end{cases}

ここで、記載の簡略化のため、k^{\prime}=k/mと置いています。出発点となる運動方程式からk、mが消えたため、今後k、mは用いず、k^{\prime}のみを用いることになります。k^{\prime}が大きいほど、空気抵抗の影響が大きくなります。

v_x(t)、v_y(t)の導出

v_xとv_yの微分方程式をそれぞれ変数分離法で解き、v_x(t)、v_y(t)を導出します。

まずは、v_x(t)です。

\displaystyle \begin{split} & \frac{{\rm d}v_x}{v_x}=-k^{\prime}{\rm d}t \\ & \\ & \rightarrow \int_{v_x(0)}^{v_x(t)} \frac{{\rm d}v_x}{v_x} = -k^{\prime}\int_0^t{\rm d}t \\ & \\ & \rightarrow \log{v_x(t)}-\log{\left(v\cos{\theta}\right)}=-k^{\prime}t \\ & \\ & \rightarrow v_x(t) = v\cos{\theta}\cdot e^{-k^{\prime}t} \end{split}

次に、v_y(t)を求めます。

\displaystyle \begin{split} & \frac{{\rm d}v_y}{g+k^{\prime}v_y}=-{\rm d}t \\ & \\ & \rightarrow \int_{v_y(0)}^{v_y(t)} \frac{{\rm d}v_y}{g+k^{\prime}v_y} = -\int_0^t{\rm d}t \\ \\ & \rightarrow \frac{1}{k^{\prime}}\log{\left(g+k^{\prime}v_y(t)\right)}-\frac{1}{k^{\prime}}\log{\left(g+k^{\prime}v\sin{\theta}\right)}=-t \\ \\ & \rightarrow g+k^{\prime}v_y(t)=\left(g+k^{\prime}v\sin{\theta}\right) \cdot e^{-k^{\prime}t} \\ \\ & \rightarrow v_y(t)=v\sin{\theta}\cdot e^{-k^{\prime}t} -\frac{g}{k^{\prime}}\left(1-e^{-k^\prime t}\right) \end{split}

無事、v_x(t)、v_y(t)が求まりました。

\displaystyle \begin{cases} v_x(t) = v\cos{\theta}\cdot e^{-k^{\prime}t}\\ \\ v_y(t) = v\sin{\theta}\cdot e^{-k^{\prime}t} -\frac{g}{k^{\prime}}\left(1-e^{-k^\prime t}\right) \end{cases}

x(t)、y(t)の導出

次に、v_x(t)、v_y(t)をそれぞれtで積分して、時刻tでの水の位置(x(t),y(t))を求めましょう。

まずは、x(t)から。

\displaystyle \begin{split} & x(t) = \int_0^t v_x(t){\rm d}t = \int_0^t v\cos{\theta}\cdot e^{-k^{\prime}t} {\rm d}t \\ & \\ & \qquad = \frac{v\cos{\theta}}{k^{\prime}}\left(1-e^{-k^{\prime}t}\right) \end{split}

次に、y(t)です。

\displaystyle \begin{split} & y(t) = h + \int_0^t v_y(t){\rm d}t = \int_0^t \left(v\sin{\theta}\cdot e^{-k^{\prime}t} -\frac{g}{k^{\prime}}\left(1-e^{-k^\prime t}\right)\right) {\rm d}t \\ & \\ & \qquad = h + \frac{v\sin{\theta}}{k^{\prime}}\left(1-e^{-k^{\prime}t}\right)-\frac{g}{k^{\prime}}t+\frac{g}{{k^{\prime}}^2}\left(1-e^{-k^{\prime}t}\right) \end{split}

前回よりはるかに複雑ですが、x(t)、y(t)も求まりました。

\displaystyle \begin{cases} x(t) = \frac{v\cos{\theta}}{k^{\prime}}\left(1-e^{-k^{\prime}t}\right) \\ \\ y(t) = h + \frac{v\sin{\theta}}{k^{\prime}}\left(1-e^{-k^{\prime}t}\right)-\frac{g}{k^{\prime}}t+\frac{g}{{k^{\prime}}^2}\left(1-e^{-k^{\prime}t}\right) \end{cases}

空気抵抗がない場合は、このx(t)、y(t)を連立させtを消去することで、xとyの関係をy=ax^2+bx+cという二次関数(=放物線)の形で表すことができましたが、式の形を見ればわかる通り空気抵抗がある場合では二次関数で表すことができません。マーライオンの口から発射された水は、空気抵抗が入ることで放物線とは異なる軌道を描くことなります。

せっかくx(t)、y(t)が求まりましたので、ここで一休みしてその軌道を見てみましょう。空気抵抗を大きくする程その軌道は歪曲され、遠くまで飛ばなくなることがわかりますね。

f:id:mattyuu:20161218112702g:plain

地面に落ちる時間が求められない

前回はy(t)=0を解いて、水が地面に落ちる時刻Tを求め、そこから水平到達距離Lを\thetaの関数で表し、それを微分することで最適解\theta_{\rm max}を求めました。しかし、今回は最初の関門であるy(t)=0が解ける気がしません。y(t)の中にはtの一次関数と指数関数が共に現れています。例えば私は、次のようなシンプルな方程式も解析的に解くことができません。*1

\displaystyle x+e^{x}=0

そして、解析的に解くことを諦めて、y(t)=0をコンピュータを使って数値的に解いたとしても、数値解は離散的にしか求められないため、その先の水平到達距離Lを\thetaで微分をするところで行き詰まるでしょう。

私はここでこの問題を解くことを諦めかけました。

しかし、諦めきれず式をいじくり回しているときに、だんだんと頭の片隅から蘇ってくるものがありました。それがラグランジュの未定乗数法だったんです。

ラグランジュの未定乗数法

ここでは、一瞬マーライオンを忘れて、ラグランジュの未定乗数法の紹介をします。(変数の記号は今回の問題に合わせてtと\thetaで記載しています。)

ラグランジュの未定乗数法
g(t,\theta)=0という条件の下、f(t,\theta)を最大化したいとき、
F(t,\theta,\lambda)=f(t,\theta)-\lambda g(t,\theta)
とおくと、下記の連立方程式の解が、その最適解の候補を与える。

\displaystyle \begin{cases} \frac{\partial F}{\partial t}(t,\theta,\lambda) = 0 \\ \\ \frac{\partial F}{\partial \theta}(t,\theta,\lambda) = 0 \\ \\ \frac{\partial F}{\partial \lambda}(t,\theta,\lambda) = 0 \\ \end{cases}

2変数の問題なのに、わざわざ変数\lambdaを付け加えて3変数にしているのがミソです。大学初年度で習う陰関数定理を用いれば証明が出来ますが、ここでは証明は行いません*2。

さて、このラグランジュの未定乗数法がなぜ今回の問題に関係あるのでしょうか。

先ほどx(t)、y(t)と、あたかもtだけの関数のようにx、yを書き下しましたが、この問題では一番遠くまで飛ばす\thetaを求めるという意味で、\thetaも変数と考えた方が都合がよく、x(t,\theta)、y(t,\theta)と書いた方がより適切です(h、v、g、k^\primeも値を変更することができますが、ここでは定数という位置付けです。パラメータと呼んだ方がいいかもしれません。)。

先ほど解くことを諦めたy(t,\theta)=0という方程式に注目しましょう。これはある角度\thetaで発射された水が地面に落ちる時刻tを求める方程式と考えてきましたが、逆に時刻tで地面に落ちるように口の角度\thetaを求める方程式と考えることもできます。その意味でtと\thetaは全く対等です。

y(t,\theta)=0を満たす(t^{\ast},\theta^{\ast})は、口の角度が\theta^{\ast}のとき、時刻t^{\ast}で地面に水が地面に落ちているということを意味します。このような(t^{\ast},\theta^{\ast})の組の中でx(t^{\ast},\theta^{\ast})を最大にするものを求めると、その\theta^{\ast}が求めたい最適の角度\theta_{\rm max}になります。つまり、y(t,\theta)=0という条件の下、x(t,\theta)を最大化したいのです。これはラグランジュの未定乗数法の前提そのものですね。

簡単な例

今回の問題をラグランジュの未定乗数法で解く前に、とても簡単な例題を解くことで前回の空気抵抗なしの場合の解法と、この後に行う今回の解法の比較をしてみましょう。

問題
g(t,\theta)=t-\theta^2=0という条件の下、f(t,\theta)=2\theta-t+4を最大化しましょう。

ではまず、前回の解法から。

直接的な解法
g(t,\theta)=t-\theta^2=0よりt=\theta^2となる。これをf(t,\theta)=2\theta-t+4に代入して、
f(t,\theta)=-\theta^2+2\theta+4=-(\theta-1)^2+5
となるから、\theta=1、t=1のとき、最大値5をとる。

この解法はg(t,\theta)=0の条件からt=\tilde{g}(\theta)のようにtを\thetaで表して、\tilde{f}(\theta)=f\left(\tilde{g}(\theta),\theta\right)という風にfを\thetaだけの関数で表すことで、通常の1変数関数の最大値問題に帰着させています。この解法の肝はt=\tilde{g}(\theta)とtを\thetaの関数で表せたことになります。しかし空気抵抗がある場合のマーライオンの問題のようにそれはいつも可能とは限りません。

では、ラグランジュの未定乗数法を用いた解法を見てみましょう。

ラグランジュの未定乗数法による解法
F(t,\theta,\lambda)\equiv f(t,\theta)-\lambda g(t,\theta)=2\theta-t+4 - \lambda(t-\theta^2)とする。

ラグランジュの未定乗数法より、最適解を得るには下記の連立方程式を解けばよい。

\displaystyle \begin{cases} \frac{\partial F}{\partial t}(t,\theta,\lambda) = -1-\lambda = 0 \\ \\ \frac{\partial F}{\partial \theta}(t,\theta,\lambda) = 2+2\lambda\theta = 0 \\ \\ \frac{\partial F}{\partial \lambda}(t,\theta,\lambda) = -t+\theta^2 = 0 \\ \end{cases}

第1式より\lambda=-1、これを第2式に代入して\theta=1、そして第3式よりt=1を得る。
最大値はf(1,1)=5となる*3。

先の解法とは全く違うのに、これで同じ解が得られることに驚きます。

ラグランジュの未定乗数法を用いて解く

さて、いよいよ本来の問題を解くときがきました。まず、条件g(t,\theta)=0のg(t,\theta)を次で定義します。

\displaystyle g(t,\theta)\equiv {k^{\prime}}^2y(t,\theta)={k^{\prime}}^2h + k^{\prime}v\sin{\theta}\left(1-e^{-k^{\prime}t}\right)-{g}k^{\prime}t+g\left(1-e^{-k^{\prime}t}\right)

g(t,\theta)=y(t,\theta)としてもよいのですが、定数{k^{\prime}}^2をかけることで記載が簡易になるようにしています。

最大化するf(t,\theta)もx(t,\theta)にk^{\prime}をかけて、次式で、

\displaystyle f(t,\theta)\equiv k^{\prime}x(t) = v\cos{\theta}\left(1-e^{-k^{\prime}t}\right)

また、関数F(t,\theta,\lambda)を、次式で、

\displaystyle \begin{split} & F(t,\theta,\lambda)\equiv f(t,\theta)-\lambda g(t,\theta) \\ & \\ & \qquad = v\cos{\theta}\left(1-e^{-k^{\prime}t}\right)-\lambda\left({k^{\prime}}^2h + k^{\prime}v\sin{\theta}\left(1-e^{-k^{\prime}t}\right)-{g}k^{\prime}t+g\left(1-e^{-k^{\prime}t}\right)\right) \end{split}

それぞれ定義します。

関数の準備が整いました。ラグランジュの未定乗数法を用いると、解くべき連立方程式は下記になります。

\displaystyle \begin{cases} \frac{\partial F}{\partial t}(t,\theta,\lambda) = k^{\prime}v\cos{\theta}\cdot e^{-k^{\prime}t} -\lambda\left({k^{\prime}}^2v\sin{\theta}\cdot e^{-k^{\prime}t}-{g}k^{\prime}+gk^{\prime}e^{-k^{\prime}t}\right) = 0 \\ \\ \frac{\partial F}{\partial \theta}(t,\theta,\lambda) = -v\sin{\theta}\left(1-e^{-k^{\prime}t}\right) -\lambda k^{\prime}v\cos{\theta}\left(1-e^{-k^{\prime}t}\right)= 0 \\ \\ \frac{\partial F}{\partial \lambda}(t,\theta,\lambda) = -\left({k^{\prime}}^2h + k^{\prime}v\sin{\theta}\left(1-e^{-k^{\prime}t}\right)-{g}k^{\prime}t+g\left(1-e^{-k^{\prime}t}\right)\right) = 0 \\ \end{cases}


それでは、解いていきましょう。まずはt\neq 0に注意して、第2式より、

\displaystyle \lambda = -\frac{\tan{\theta}}{k^{\prime}}

を得ます。いきなり\lambdaがわかりました。幸先のよいスタートです。

次に、得られた\lambdaを第1式に代入して、

\displaystyle k^{\prime}v\cos{\theta}\cdot e^{-k^{\prime}t} +\frac{\tan{\theta}}{k^{\prime}}\left({k^{\prime}}^2v\sin{\theta}\cdot e^{-k^{\prime}t}-{g}k^{\prime}+gk^{\prime}e^{-k^{\prime}t}\right) = 0

となり、この両辺に\cos{\theta}をかけて整理すると、

\displaystyle \begin{split} & k^{\prime}v\cos^2{\theta}\cdot e^{-k^{\prime}t} +k^{\prime}v\sin^2{\theta}\cdot e^{-k^{\prime}t}-g\sin{\theta}+g\sin{\theta}\cdot e^{-k^{\prime}t} = 0 \\ & \\ & \rightarrow k^{\prime}ve^{-k^{\prime}t} +g\sin{\theta}\cdot e^{-k^{\prime}t}=g\sin{\theta} \\ & \\ & \rightarrow e^{-k^{\prime}t} = \frac{g\sin{\theta}}{k^{\prime}v+g\sin{\theta}} \\ & \\ & \rightarrow t=\frac{1}{k^{\prime}}\log{\frac{k^{\prime}v+g\sin{\theta}}{g\sin{\theta}}} \end{split}

形は複雑ですが、tを\thetaの式で求めることができました。では、最後にこれを第3式に代入して整理しましょう。

\displaystyle \begin{split} & {k^{\prime}}^2h+k^{\prime}v\sin{\theta} \frac{k^{\prime}v}{k^{\prime}v+g\sin{\theta}} - g\log{\frac{k^{\prime}v+g\sin{\theta}}{g\sin{\theta}}} + \frac{gk^{\prime}v}{k^{\prime}v+g\sin{\theta}} = 0 \\ & \\ & \rightarrow {k^{\prime}}^2h\left(k^{\prime}v+g\sin{\theta}\right)+{k^{\prime}}^2v^2\sin{\theta}-g\left(k^{\prime}v+g\sin{\theta}\right)\log{\frac{k^{\prime}v+g\sin{\theta}}{g\sin{\theta}}}+gk^{\prime}v = 0 \end{split}

あとは、これを\thetaについて解くだけなのですが、見るからに無理そうです。この方程式を満たす(0,\pi/2)の範囲の解が、求めたい最適解\theta_{\rm max}になります*4。

最後の最後で数値解になってしまうことは口惜しいですが、仕方ありません。

最後の式の左辺を\thetaではなく、一般の変数xの関数f(x)とみなして、demosでy=f(x)のグラフを描いてみましょう。

f:id:mattyuu:20161218110018p:plain

k^{\prime}、h、v、gのパラメータをいじることで、関数のグラフは形を変えますが、確かに(0,\pi/2)の範囲に1つだけ解がありそうです。

そして、 k^{\prime}を大きくしていくと、最適解\theta_{\rm max}は小さくなっていきます。これはh=0の場合によく知られている事実と合致します。
f:id:mattyuu:20161218110457g:plain

demosで確かめる

残念ながらdesmosには方程式をsolveする機能はありませんので、パラメータを連続的に動かした時の最適解\theta_{\rm max}をdesmos上だけで変動させて表示させることはできません。パラメータを固定して\theta_{\rm max}を数値計算で求め、desmosを使って実際に最適解が得られていそうということを確認します。前回同様、マーライオンにいろいろな角度で水を発射してもらい、飛距離が最大になる最適解の角度の方向にはオレンジ色の点線を引き、最大到達地点x(\frac{1}{k^{\prime}}\log{\frac{k^{\prime}v+g\sin{\theta_{\rm max}}}{g\sin{\theta_{\rm max}}}},\theta_{\rm max})=\frac{v^2\cos{\theta_{\rm max}}}{k^{\prime}v+g\sin{\theta_{\rm max}}}には赤い旗を置いておきましょう。そして、比較のため空気抵抗がない場合の軌道を紫色の点線で表示します。

下記の例では、 k^{\prime}=0.2、h=1、v=2、g=1としました。その時の\theta_{\rm max}はおよそ0.56となります。

f:id:mattyuu:20161218115515g:plain

うまくいってますね!

まとめ

最後の最後は数値解になってしまいましたが、ラグランジュの未定乗数法を用いて曲がりなりにも最適解\theta_{\rm max}を求めることができました。最後に得られた式を使って、空気抵抗がある場合は最適な角度は常に下方に修正されるのか、といったことを調べたら面白いと思うのですが、今回はこの辺りにしておきます。

空気抵抗なし版の問題もラグランジュの未定乗数法を用いて、前回と同じ解析解を得ることができます。やってみるとラグランジュの未定乗数法を用いる方が遥かに簡単でした。一度計算してみたら面白いと思います。

ラグランジュの未定乗数法は大学で習っており、「あぁ、そんなものもあるんだなぁ、便利そうだなぁ」くらいの印象でした。しかし自分が問題設定してラグランジュの未定乗数法を使わないと解けない問題に出会うとその凄さに驚嘆せずにはおれません。「あぁ、楽しい計算だったなぁ」

f:id:mattyuu:20161218123247j:plain

*1:ここで、解析的に解くとは、「誰もが知っている関数で解を書き下す」という意味で使っています。

*2:https://www.amazon.co.jp/これならわかる工学部で学ぶ数学-千葉-逸人/dp/490381419X参照

*3:ラグランジュの未定乗数法は最大値の候補(極値)を与えるだけなので、最大値となっていることは別途議論が必要です。

*4:実際は極値になるだけですので、最大値になるかの議論は必要です。

麻雀×プログラミング×数学

この記事は、伝道師になろう！ Advent Calendar 2016 - Adventar 11日目の記事です。10日目は狡猾な狐さんの「ARとしての知識を実装して出かけよう！」(予定)です。

はじめに

東京では今年から「伝道師になろう」というイベントが定期的に開催されています。伝道師になろうってなんだろう?」という方は、狐さんの下記の記事を参照してください。


wreck1214.hatenablog.com


簡単に言えば「自分の好きなことを、知らない人にもわかりやすく、なぜ好きかを伝えて、面白さを伝道しよう」という企画です。僕も参加したことがあります*1。

コンセプトが本当に素晴らしいので、来年はこのイベントが全国に広がればいいですね。狐さんによると、ビブリオバトルのように誰が始めたかわからないくらいとのことです。

数学以外に好きなもの

僕が好きなものは絶対的に数学です。伝道できる(と自分で思っている)ものは数学しかありません。しかし、伝道師になろうは日曜数学会から派生したイベントということもあって、このアドベントカレンダーに数学ネタだけを書くのは少し気が引けます*2。

僕は数学以外で好きなものを聞かれると「麻雀」、「プログラミング」、「京都」、「熊野寮」、「料理(←最近やっていない)」、「小籠包」、「ヒカルの碁」から選出して答えていますが、数学とこれらには、好きレベルでいうと、


数学>>>>>>|超えられない壁|>>>>>>麻雀、プログラミング、京都、熊野寮、料理、小籠包、ヒカルの碁(順不同)


くらいの差があります。数学に関しては永遠に熱く語る自信がありますが、その他の好きなものは最高でも3時間程度しか語れない気がします。ただその他の好きなものも数学とコラボさせることでそれなりのボリューム、ネタになるのではないかと思い、今回「麻雀」、「プログラミング」を選抜して、投稿しました。

まず最初に、伝道師になろうのルールに則って、「麻雀」、「プログラミング」、「数学」がそれぞれなぜ好きなのかを語っていくことにしましょう。

麻雀が好きな理由 (半分は熊野寮の話)

麻雀が好きな理由を2点挙げます。

まず1点目は運の流れを感じることができる、そして運がいい時は最高に気持ちがいいということです。

普段生活の中で運がいいとか、悪いとかを考えることは少なくなってきましたが、麻雀をやっていると何をやってもうまくいく無双状態というものを体感できます。とにかく配牌(最初に配られる牌)がいい、とにかくいい牌をツモる(牌を引くこと)、そしてとにかく点数が高い状態です。まさに下記の女子プロの豊後さんのような状況です。


【麻雀】 ブンゴ無双 【役満】

無双状態は滅多に発生するものではありませんが(感覚的には5半荘に一回)、無双状態になった時の快感を味わうために麻雀をやっていると言っても過言ではない気がするくらい楽しく気持ちいいです。

2点目は仲間と語り合いながらやるとやっぱり楽しい、というところです。

(ここからは熊野寮の話を交えつつ説明していきます。)

僕が麻雀を覚えたのは大学生時代の時です。先に挙げた数学以外の好きなものの中に「熊野寮」というものがありましたが、僕は学生時代に熊野寮に6年間住んでいて、そこで先輩に麻雀を教わりました。熊野寮は古いですが最高の寮です。素数大富豪の大会が開催された最初の場所でもあり、最近のトレンドをいち早く吸収するあたりはさすがと思います。Googleで"熊野寮"を画像検索するとその凄さをわかっていただけると思いますので、是非調べてみてください。なお、家賃は700円です。ただ電気代や食堂のお姉さんたちの保険料等が追加され、実際は4,100円でした。それでも安い!

熊野寮にはA棟、B棟、C棟があり、それぞれ各階に談話室があります。僕はA棟3FのA3談話室に住んでいました。これは談話室に入り浸っていたという比喩的な表現ではなく、談話室の片隅に本棚で人一人分寝れるスペースを作って、そこに住んでいました。談話室には大学生の娯楽が詰まっています。ゲーム、漫画、お酒、麻雀。。。みんな好きな時に談話室に来て、好きなことをして語り合っていました。今思っても楽しい空間です。そして下記理由から麻雀こそ仲間との語り合いにピッタリな遊び、ゲームだと思うのです。

  • 手を動かすだけで疲れないから長時間できる
  • ご飯食べながらでもできるから長時間できる
  • 時々「ポン、チー、ロン、ツモ」という以外は発声の必要がないいからお話しができる

実際、談話室でやる麻雀は娯楽でやる麻雀ですので、いつも楽しく語り合っていました。朝まで仲間*3と語りながら麻雀をやることも珍しくなかった僕としては、麻雀は語り合いのゲームで楽しいんだ!ということを伝道したいと思います*4。

女子用の雀荘(タバコは絶対禁止、フードメニューにスイーツ、カクテル多め)とか作れば、女子会にピッタリな気がしますね。

では麻雀が好きな理由はこれくらいにして、先に進みましょう。

プログラミングが好きな理由

これに関しては、数学の計算ができるからという一言に尽きます。コンピュータというのは和訳すると「計算機」ですし、プログラミングはコンピュータに命令を与えるものですので、「計算機に命令を与えるもの」がなぜ好きかと聞かれて、「計算ができるから」と答えてるだけになります。何の回答にもなっていないように思えますが、その言葉以上に計算できることは素晴らしいと思うのです。

例えば、あなたがお望みであれば僕は4,294,967,297や147,573,952,589,676,412,927が素数でないことをプログラミングで簡単に示すことができますが、一昔前はそれは大変難しいものでした。

4,294,967,297に関してはtsujimotterさんの下記の記事を、
tsujimotter.hatenablog.com

147,573,952,589,676,412,927に関しては、下記のフランク・ネルソン・コールの逸話をご参照ください。
メルセンヌ数 - Wikipedia

上記の記事にあるように、昔の人は20桁程度の素因数分解ですら難しかったのです*5。円周率の計算に一生を費やした人がいます。そして死後に誤りを指摘されてしまった悲しい人もいます。その人達が一生を費やして計算した円周率の値は、おそらく僕は30分程度プログラミングをすると手に入れることができます。

17世紀のケプラーはブラーエが残した大量の観測データから、地球が楕円軌道で太陽の周りを回っていることを突き止めました。また、18世紀のオイラーが\zeta(2)が\frac{\pi^2}{6}と予想したのも*6、19世紀にガウスが小惑星ケレスを発見できたのも全部大量の計算の賜物でしょう。

僕自身も前回のブログで書いた素数大富豪の会心の一撃のシミュレーションは、プログラミングをしなければ一生かかっても終わらなかったでしょうし(というか絶対途中でやめる)、プログラミングをせずにゼータ関数の零点を計算するのがいかに大変であるかを経験したこともあります。

ネイピアが対数を発見するために引きこもった強い動機となったように、昔の人の計算との戦いは想像を絶するものだったと思います。もしオイラーやガウスにプログラミングを教えてあげれば、二人は泣いて喜ぶのではないでしょうか*7。プログラミングをすることでそのような計算をいとも簡単に行うことができる。直接的な理由ですが、これがプログラミングを好きな理由です。

数学が好きな理由

これは語ると切りがありませんし、語り切れるとも思えません。そしてまだまだ自分でもなぜ好きがかわからない状況です。ちょうど半年くらい前に自分がなぜ数学が好きなのかをマインドマップで殴り書きしていました。狂気じみていますが、こちらを貼り付けることで好きな理由に代えさせていただきます。
f:id:mattyuu:20161211024737j:plain

ついでに同じタイミングで数学ってどんな分野があるかも殴り書きしていました。
f:id:mattyuu:20161211024754j:plain

どちらも何も見ずに5分くらいで書き上げたものです。汚い字で読めないですが、それだけ好きなんです、

ただ、これじゃ全然伝道していないので、1点だけ言わせてください。(ここから先はあくまで持論です。)

「数学は紙とペンがあればどこでもできる*8」ということを好きな理由として強調したいです。これは「紙とペンがあれば計算できるしね」みたいな単純な話ではありません。数学的実体なんてこの世のどこにも存在しません。 虚数はその名の通り現実に存在しないような名前が不本意に与えられていますが、実数の1という数はどこに存在していますか？紀元前にピタゴラスによってピタゴラスの定理が発見されましたが、ピタゴラス以前にもピタゴラスの定理は成り立っていたはずです。ピタゴラスの定理はどこにあったのでしょうか?

1もピタゴラスの定理もこの世界のどこにも存在しません。でもこの世界のどこにでも存在しているのです。そして紙とペンを使うことでその数学的実体の化身が目の前に現れ、時代、人種、性別、年齢を超えて他の人とも寸分違わず同じものを認識しあえるのです。それは本当に奇跡のようなことだと思います*9。こんな奇跡を体感できるなんて。。。はぁはぁ。。数学ってすごいですよね。

この話はエドワード・フレンケルさんの「数学の大統一に挑む」から影響を受けていると思います。

数学の大統一に挑む

数学の大統一に挑む

また、日曜数学者の始祖であるtsujimotterさんも過去にtwitterで「旅行に行っている暇がない。今は数学的世界を探索するのに忙しい。」のようなことをつぶやいておられました。過去のツイートなので探しきれませんでした。まさに同感で、僕らは紙とペンを使うことで、この世界よりもずっと広大な数学的世界を探検できるのです。。。はぁはぁ。。数学ってすごいですよね。

いよいよコラボ まずは麻雀の基本説明

なぜ好きかの説明が終わりました。「麻雀」×「プログラミング」×「数学」のコラボをしたいと思います。

そのためには最低限麻雀のことを知る必要がありますので、説明します。

まずは麻雀牌にどのようなものがあるかを説明します。麻雀牌は大きく数牌と字牌に分かれるのですが、この後の考察に字牌は不要ですので、数牌だけを考えます。数牌には萬子(マンズ)、筒子(ピンズ)、索子(ソーズ)の3種類があり(これを色と言います)、それぞれ1から9の数字が書かれています。色と数字が同じである牌はそれぞれ4枚あります。つまり数牌は3(色)×9(数字)×4=108枚存在します。

f:id:mattyuu:20161211024426j:plain

次にアガリの形を説明します。麻雀は14枚の牌を「3枚1組のグループ4つ」と「2枚1組の雀頭(ジャントウ)１つ」の形にすることでアガることができます。「3枚1組」は全く同じ牌3枚でもいいですし、同じ種類の牌で2,3,4のように3つの連続する整数の組でも良いです。雀頭は全く同じ牌の2枚組になります。

下記がアガリの形の例になります。

f:id:mattyuu:20161211024500j:plain

麻雀は基本的には手持ちの牌は13枚になっています。13枚の状態で1枚牌を引いて、14枚の状態になった時に、アガリの形になっていればアガることができます。アガれない場合は1枚引いたので、手持ちの14枚の牌どれでもいいので1枚捨てて再び13枚の状態にします。このようにアガリの形になるまで牌を交換し続けるゲームになります。アガリの形が特定の条件を満たせば役というものが付きます。役が珍しいものであればあるほどアガった時にもらえる点数が多くなります。

最後に、待ち牌について説明します。次に良い牌を引けばアガれる状態のことをテンパイしているといいます。テンパイしている時に引くことでアガれる牌、つまりアガるために待っている牌を待ち牌といいます。下記の例のように待ち牌は1種類とは限りません。

f:id:mattyuu:20161211024510j:plain
f:id:mattyuu:20161211024532j:plain

いろんなパターンの待ち牌がありますね。

以上で麻雀の基本説明を終わります。

第二の九蓮宝燈探索

唐突ですが、13枚の手持ちの牌が次の形の時、1〜9の牌全てが待ち牌になります。つまり1〜9のどの数牌を引いてもアガることができます。

f:id:mattyuu:20161211024518j:plain

実際に確かめてみると、確かにアガリの形になっています。
f:id:mattyuu:20161211024543j:plain
f:id:mattyuu:20161211024556j:plain
f:id:mattyuu:20161211024550j:plain

なお、この形でアガると九蓮宝燈(チューレンポートウ)という役の中でも特別な役満という役になります。九蓮宝燈はアガったら死んでしまうという迷信があるくらい珍しい役になっております。

「1〜9のどの牌がきてもアガれるなんて九蓮宝燈はすごいなあ。」と思うところですが、知られていないだけで、他にも1〜9のどれがきてもアガれる13枚の手牌が存在するも知れません。そんな未だ日の目を見ない第二の九蓮宝燈を見つけてあげたい、見つかったら僕だけでもそいつのことを役満と呼んであげたい。しかし、麻雀の手牌の組み合わせはそれなりに大きい数になりそうだし、それぞれの手牌でちゃんと待ち牌を漏れなく見つけられるか不安。。。。。

そんな悩みを解決するのがプログラミングなんです！

コンピュータはプログラミングした通りに高速に、嫌とも言わず、間違うことなく動いてくれます。正しくプログラミングすれば正しい結果が得られます。今回はどのようにプログラミングしたかの紹介はしませんが*10、実際にプログラミングして実験した結果、九蓮宝燈以外には1〜9の全てが待ち牌になるような手牌は存在しませんでした。

第二の九蓮宝燈が見つからなかったのは悲しいですが、それだけ九蓮宝燈はすごいものなんだ!ということがわかりました!

「プログラミングって強力だなあ。家に麻雀牌なくてもこういうことがわかるんだからなあ。」と思います。

一定の成果を得ましたが、まだ「麻雀」×「プログラミング」止まりです。

数学的に意味がある待ち牌

さて、プログラミングによって色々な手牌の待ち牌を調べられるようになりました。では最後に数学的エッセンスを加えて本記事を締めくくりましょう。

スライドを見ていただければ、説明は不要です。なお、不思議なことに下記の4種類すべてちょうど2パターンずつでした。

f:id:mattyuu:20161211024611j:plain
f:id:mattyuu:20161211024617j:plain
f:id:mattyuu:20161211024623j:plain

最後はみんな大好き素数です!

f:id:mattyuu:20161211024637j:plain

明日はHaru Negamiさんの「好きすぎて食べちゃいたい◯◯」です(^_^)

ありがとうございました。

f:id:mattyuu:20161211024649j:plain

*1:今回は学生時代の話をするため、一人称は「僕」を使いたい気分です。

*2:数学に限らず好きなものをどんどん発表していこう!というのが元々の誕生理由な気がするため。

*3:熊野寮のメンバーは家族でも友達でもない仲間という表現がよく当てはまる気がします。

*4:しかし熊野寮のような特殊な環境はほとんどありません。誰かのアパートで麻雀しても多分隣の部屋の人からうるさいと怒られるでしょう。雀荘に行ったらお金がかかります。仲間が集い会えるような場所が欲しいところです。

*5:本質的には今も難しいです。RSA暗号の基本原理です。でもそれはまた別のお話。

*6:オイラーは後にちゃんと証明しています。

*7:しかし、この二人は数学界の偉人中の偉人で常任にはその底は計り知れないので、コンピュータどころでは泣いて喜ばない気もします。オイラーは人が呼吸をするように、鳥が空を飛ぶように計算していたと言われています。

*8:実際は教科書も欲しいところ。

*9:かなりスピリチュアルなことを言っていますが、哲学的な話は大の苦手です

*10:汚いソースを見られるのは超恥ずかしいため。

素数大富豪〜会心の一撃の成功確率〜

はじめに

今素数大富豪が数学好きの間で静かなブームになっています。「素数大富豪って何ですか?」という方は、本アドベントカレンダー1日目のにせいさんの記事を読んでください。

さて、今年もアドベントカレンダーの季節がやってきました！ この記事は、素数大富豪アドベントカレンダー 1日目の記事です。

素数大富豪アドベントカレンダー！ - にせいの日記


私は素数大富豪を1回しかやったことがなく、戦略やゲームの中で出会った素数達について語れるレベルではありませんが、せっかくの機会ですので何か投稿したいと思い、素数大富豪 AdventCalendar 2016に2回分登録させていただきました。本記事はその1つ目の投稿です。2つ目は12月24日です。

では、よろしくお願いします。

問題

突然ですが、素数大富豪には一発大逆転の秘策があり、例えば対戦相手の手札が1枚、自分の手札が10枚で圧倒的に不利な状況でも、うまくヒットすれば逆転勝利が可能です。これは通常の大富豪にはない素数大富豪の特徴です。

その秘策とは"秘"密にするほど大したものではなく、自分からカードを出せる場面で、手持ちの全カードを一列に並べて場に出すことです。場に出したカードから決まる数が素数であれば、手札が全てなくなるため勝利となります。しかし合成数だった場合は、ペナルティとして山札から出したカードと同じ枚数だけ引かないといけないためハイリスク、ハイリターンの策といえるでしょう。

ここでは、自分の手札全てを一気に場に出すこの策を「会心の一撃」と呼ぶことにします(上述の通り会心の一撃に失敗すると大ダメージをくらうため、メガンテと呼んだ方が妥当な気もします。)。

先日のMathPowerで(第一回)素数大富豪大会が開かれましたが、MathPowerの司会も務められたタカタ先生が会心の一撃を狙って場を盛り上げていた記憶があります。今回の記事では、会心の一撃について考察します。タイトル通りですが、問題は下記です。

f:id:mattyuu:20161206014435j:plain

答えの確率をNを使った簡単な数式で表すことができれば大変よいのですが、ただでさえ不規則な素数が考察対象になる上に、トランプという自然界から見たらなんの必然性もないよく分からない意味不明で超人工的な存在(数は1～13の13種類で、数字は人間の指が10本あるためか10進法で、スートというよく分からない記号があるため各数のカードの枚数は4)も考察対象になることから、それは不可能に思えます。今回はコンピュータを用いて頑張ることにします。

シミュレーションによる解法

N枚のトランプで会心の一撃を行うときの成功確率は下記で求めることができます。


\displaystyle \frac{\mbox{N枚のカードを一列に並べた時に素数となるパターン数}}{\mbox{N枚のカードを一列に並べる全パターン数}}


さて、まずは分母の全パターンを計算してみましょう。
それは手札にN枚のカードが配られていて、それを一列に並べるパターンとなります。
まず52枚のカードからN枚を選出するパターンが{}_{52}{\rm C}_N通り、その各々に対してN枚のカードを一列に並べるパターンはN!通りになるので、

\displaystyle {}_{52}{\rm C}_{N}*N!=\frac{52\cdot 51 \cdots (52-N+1)}{N!}\cdot N!=52\cdot 51\cdots (52-N+1) \bigl(={}_{52}{\rm P}_N\bigr)

となります。(最初から{}_{52}{\rm P}_Nで計算すればよかったですね。)

{}_{52}{\rm P}_Nのパターンの数のうち、何パターンが素数になるのかをコンピュータに計算させることで確率の分子を得ることができます。そこから手札がN枚の時の確率を求めますることができるのですが、この分母の{}_{52}{\rm P}_N、つまり素数判定しないといけないパターンの数が、Nと共にどのくらいの数になるのかというと、

N=1:52\mbox{通り}
N=5:311,875,200\mbox{通り}
N=10:57,407,703,889,536,000\mbox{通り}
N=20:306,532,583,223,109,543,994,300,006,400,000\mbox{通り}
N=30:71,759,895,859,162,404,648,138,410,532,134,505,676,800,000,000\mbox{通り}
N=40:168,388,112,212,869,181,588,664,081,406,834,062,715,634,990,448,640,000,000,000\mbox{通り}
N=50:40,329,087,585,471,939,285,830,318,428,201,883,487,644,752,720,441,638,912,000,000,000,000\mbox{通り}
N=52:80,658,175,170,943,878,571,660,636,856,403,766,975,289,505,440,883,277,824,000,000,000,000\mbox{通り}

となります。N=52では8,065不可思議通りもの素数判定をしなければいけません。そして、その判定する桁数は68桁です(1桁*9(1～9)*4(スート)+2桁*4(10～K)*4(スート))。そんな巨大な数を素数判定するだけでも時間的なコストがかかりそうです。。。

実際N枚のカードで全通りのパターンに対して素数判定を行う方法では、N=8(←うろ覚え)で私のコンピュータと私の我慢の限界に達しました。
1日かかっても終わりません。N=20くらいまでいくと、現在のスーパーコンピュータでも宇宙誕生から今日までの日数をかけても計算できないかもしれません。

そこで、今回は各Nに対してそれぞれ10万回ずつ会心の一撃を繰り出し、その成功確率を算出しました。(実際は10万回のシミュレーションを10セット行っていますので、各Nに対して100万回シミュレーションしています。)

コンピュータを使って問題を解かないといけない時点で悔しいですが、その上正確な値ではなく数値実験の結果に止まることは甚だ遺憾です。しかし、世の中は解ける問題よりも解けない問題の方がはるかに多いのです。あきらめるしかないでしょう。

一つ気がかりがあります。N=52では8,065不可思議通りものパターンがあります。そんなにパターンが多いのにたった10万回の施行でそれなりに正しい確率は出せるのでしょうか? 10万回1セットのシミュレーションの結果の確率が上下すれば、それは真の確率に近いものとは言えないでしょう。

この辺りは専門的な確率統計の知識があれば定量的に評価できるものと思われますが、実際に10万回やれば結果の確率にぶれは無く、確率はある値に近づいていくことがわかりました。関東地方の視聴率は関東地方1600万世帯のうち600世帯だけでしか調査していないとのことなので、同じ感覚で10万回やれば十分なのでしょう!

それにしても回数が全然違いますが、、、この辺りは今後統計検定の資格を狙って勉強していきたいと思います。
80,658,175,170,943,878,571,660,636,856,403,766,975,289,505,440,883,277,824,000,000,000,000
100,000

結果

ではシミュレーション結果をお見せします。各Nに対して、10万回の会心の一撃を繰り出しその成功回数を表にしております。上述の通り10万回1セットだけでは本当に有効なシミュレーションができているかわからないため10セット繰り返しています。一番右の列が各Nに対してトータル100万回の施行での確率です。

f:id:mattyuu:20161205220325p:plain

各Nで会心の一撃が成功した回数はセットによらずおおよそ一定となっており、それなりに有効なシミュレーションになっていそうということはわかっていただけると思います。Nに対して確率がどのように変動するのかがわかるグラフも書いておきます。

f:id:mattyuu:20161206014449j:plain

見ていただければわかるかと思いますが、Nが大きくなるにつれて会心の一撃の成功確率は急激に下がっていきます。
N=1では成功確率46%(13種類の数のうち、素数は2,3,5,7,11,13の6種類なので6/13≒46%)ですが、N=10で成功確率は4%に落ちます。
素数定理によると、整数x以下の素数の個数\pi(x)は、

\displaystyle \pi(x)\simeq\frac{x}{\log x}

で表せるので、整数xが素数である確率p(x)は

\displaystyle p(x)\simeq\frac{\pi(x)}{x}=\frac{1}{\log x}

のように考えることができます。xが大きくなればなるほどp(x)はどんどん小さくなります。つまり整数が大きくなっていくにつれて素数はまばらになっていくのです。Nを大きくするということはより大きな整数で会心の一撃を行うわけで、出した数が素数になる確率はどんどん下がっていくのです。このp(x)と今回得られた確率のグラフと比べるのは乱暴ですが、似ているような気がしませんか?

f:id:mattyuu:20161206014508j:plain

さて、得られた結果をよく見ると、N=52で会心の一撃の成功確率が上がっています。これは一体どうしたことでしょうか? 答え(と予想するもの)は注釈に書きますので、みなさん考えてみてください。*1

曲がりなりにも答えの確率を求めることができましたが、不満に思われている方もいるでしょう。素数大富豪に限らず、2の倍数、3の倍数、5の倍数は簡単に判定できます。会心の一撃を繰り出す時に、その判定を行わない人はいないでしょう。

ということで、会心の一撃で出したトランプから決まる数が2、3、5のいずれかの倍数の場合は、ノーカウントととしてシミュレーションしたグラフも載せておきましょう。

f:id:mattyuu:20161206014524j:plain

グラフだけでは分かりにくいですが、N=52でも確率は2.5%程度までしか落ちません。40回に1回は成功する計算になります。

まとめ、感想

冒頭に記載した、素数大富豪アドベントカレンダー5日目のにせいさん記事を見ると、上級者は会心の一撃を狙ってないことがわかり、少し意外でした。確かに今回の結果を見ても会心の一撃の成功率は高くありません。大会となると人の目も気になるので、会心の一撃を出しにくいのかもしれませんね。

ただ、素数大富豪は究極の素数出会い系カードゲームですので、今後街中で会心の一撃に失敗した人を見かけたら、笑ったり、哀れんだりせずに、「全力で未知の素数に出会おうとしている肉食系の人だ」と心の中で褒めてあげましょう。

おまけ

最後に今回コンピュータを使ってどのようにシミュレーションしたのかを、超簡単にお伝えします。

私は普段Javaを使っているのですが、Javaで素数大富豪で出てくるような巨大な数を扱うのは難しいと思っています。BigIntegerというクラスを使わないといけないですし、BigIntegerに素数判定のために用意されているメソッドも大きい整数に対しては確率的にしか判定できないようです。(当初BigInteger用の素数判定メソッドを作りかけていたのですが、絶対遅いだろうと思って途中でやめました。)

そのため、以前から少し触っていたSageというアプリを使うことにしました。インストールしたのが昔ですので、どのサイトを参照してどうインストールしたかは覚えていません。インストールしたい方は、「Sage 数学」で検索してみてください。

Sageを使えば\zeta(2)や1000!も簡単に求めることがきます。

f:id:mattyuu:20161206015221p:plain

is_prime()というメソッドが用意されており、次のような巨大な数でも素数判定を一瞬で行ってくれます。

f:id:mattyuu:20161206015205p:plain

Sageでプログラミングをする時はPythonというプログラミング言語を使うようです。今回初めてPythonでプログラミングしました。今回のシミュレーションのソースは下記の通りです。

import random

# n枚のトランプでkaisuu回会心の一撃を繰り出して、その成功回数を出力する
def sosuDaifugo(n,kaisuu):
	# 成功した回数を格納する変数
	counter = 0

	# kaisuu回繰り返す
	for i in range(0,kaisuu):
		# トランプの定義
		trump = [ "1", "1", "1", "1", "2", "2", "2", "2", "3", "3", "3", "3", "4", "4", "4", "4", "5",
				"5", "5", "5", "6", "6", "6", "6", "7", "7", "7", "7", "8", "8", "8", "8", "9", "9", "9", "9", "10", "10",
				"10", "10", "11", "11", "11", "11", "12", "12", "12", "12", "13", "13", "13", "13"]

		num = ""
		for j in range(0,n):
			# ランダムにトランプを配るイメージ
			card = trump[random.randint(0,len(trump)-1)]
			# 配られたトランプを右に並べていくイメージ。文字列結合であることに注意
			num = num  + card
			# 配られたトランプは山から無くすイメージ
			trump.remove(card)

		# 素数判定
		if is_prime(int(num)):
			# ここのコメントを外すと出会えた素数も出力される
			#　print int(num)
			counter = counter + 1

	print counter


# メイン部分 1枚から52枚までのループ
for n in range(1,53):
	print n
	# 100000回繰り返す
	sosuDaifugo(n,100000)

このソースファイルをsosu.py (.pyはPythonのソースコードの拡張子)として保存し、Sage上でloadすれば実行されます。

f:id:mattyuu:20161206015118p:plain

素数大富豪のシミュレーションと聞くと大変そうなイメージを持つ方が多いかもしれませんが、ソースは意外に短いでしょう?こんな簡単なソースでシミュレーションできるわけですから、素数大富豪に限らず身の回りで気になることがある方はどんどんシミュレーションしていきましょう!

f:id:mattyuu:20161206015636j:plain

*1:「52枚のカードを全て出すと、各桁の和は256となり、3の倍数にならない。つまり52枚全て使ったときの数も3の倍数にならない。3の倍数が出なくなる分、心持ち素数が出やすくなる」と予想します。

有限体Fpの乗法の可視化

twitterを始めて数学好きの方をフォローするとp進数体(\mathbb{Q}_p)という言葉をよく聞くようになりました。p進数体は数論をやる中では空気のように出てくる数学的対象のようですが、数論を勉強して来なかった私は今年5月の第13回数学カフェ(素数回)でその存在を初めて知りました。

ずっと数学は好きだったので、読み物とかでp進数体は何度か目にしていたとは思うのですが、おそらくその度に頭の中で有限体{\mathbb F}_pに変換されていたものと思います。({\mathbb F}_pは体(たい)だし、1,2,\cdots,p-1,p(=0)という風にp進むと0になって、なんか10進法に似てるので、{\mathbb F}_pをp進数体と勘違いしても仕方ないような気がします。。。)

p進数体は今後本格的に勉強して何かブログを書いてみたいと思いますが、今回は{\mathbb F}_pに関して、その乗法(掛け算)、除法(割り算)の構造を可視化してみたので紹介したいと思います。(今回は可視化した絵を紹介したい!ということが目的ですので{\mathbb F}_pが何なのかといった数学的背景の説明は必要最低限に留めます。というかほとんどないです。)

では可視化の対象となる{\mathbb F}_pとその乗法、除法が何者であるかを説明していきましょう。

以下pを素数として、次のp個の元からなる集合を{\mathbb F}_pと定義します。(厳密には加法(足し算)、乗法をセットで定義したものを{\mathbb F}_pとします。)

\displaystyle \mathbb {F}_p = \{0,1,2,\cdots,p-1\}

さて、{\mathbb F}_pの元a、bに対してその積a*bを、

\displaystyle a*b=(ab\mbox{を}p\mbox{で割った余り})

で定義します。通常の積とは異なるので*を付けています。

いかめしく「積を定義する」と書いてしまいましたが、何のことはございません。{\mathbb F}_pを整数の部分集合と考えて通常の積を行ってしまうと、その結果が{\mathbb F}_pからはみ出してしまうので、pで割った余りを積としているだけです。
例えば 、p=5の時、{\mathbb F}_5の元3と4を通常の積で計算すると12となりますが、{\mathbb F}_5=\{0,1,2,3,4\}であり、12は含まれません。それでは困るので、5で割った余りを考えて3*4=2としているのです。

実際にp=5の場合の九九表を書いてみましょう。(正確には五五表です。)

f:id:mattyuu:20161107003050j:plain

次に除法を定義しますが、実際にどうやって計算するのかを示した方が早そうです。積と同様に{\mathbb F}_5にて3\div 4を計算してみましょう。{\mathbb F}_5を有理数の部分集合と考えて通常の割り算を行うと、3\div 4=3/4=0.75となりますが、積の時と同様こちらも{\mathbb F}_5に含まれません。有理数での考えは捨てて、3\div 4の答えが満たす条件を考えましょう。

\displaystyle 3\div 4=X

とした時、両辺に4をかけて、

\displaystyle 3=4*X

となりますので、3\div 4の答えは、4を掛けると3になる{\mathbb F}_5の元ということになります。

そこでもう一度{\mathbb F}_5の九九表を見てみると、

f:id:mattyuu:20161107003108j:plain

2に4を掛けると3になることが分かりますので、{\mathbb F}_5の世界では3\div 4=2となります*1。日常生活の計算を実数の世界で行っている我々には不思議に思えますが、これは{\mathbb F}_5の世界ではこれは立派な除法になっているのです。

では除法も九九表で表しておきましょう。

f:id:mattyuu:20161107003121j:plain

簡単ではありましたが、可視化の対象である{\mathbb F}_pとその乗法、除法の説明を終わります。

さて、任意の素数pに対して、{\mathbb F}_pの乗法と除法の九九表を書けるようになったのですが、数字だけの九九表だけを見てもあまり面白みが湧きません。

f:id:mattyuu:20161107003728j:plain

幸い私たち人間は色を認識することができますので、この九九表を決まった規則で色分けして可視化してみましょう(九九表になっている時点で乗法、除法は可視化できているとも言えるのですが、この記事では色分けを以って可視化するいうことにしています)。

色分けすると言っても、適当な色を塗ってしまっては台無しです。数が大きくなるにつれてだんだん薄くなるように色を塗りましょう(厳密には{\mathbb F}_pは順序集合ではないため、数が大きい、小さいという概念はありません)。

f:id:mattyuu:20161107003205j:plain

{\mathbb F}_pの九九表は人間が勝手に決めた表ではなく、自然の規則だけで作成した表であり、塗る色もできるだけ人為的でない規則で選んだつもりですので、きっと美しい絵図が現れるはずです！

なお、乗法の0が入る行、列は全部0で面白みがありません。また除法も0の列は禁断の零除算で定義されませんし、0の行は0が並ぶだけでこちらも面白みがありません。そのため乗法、除法共に0の行と列は可視化対象外とします。

では、まずは乗法から。最初にp=19の乗法の九九表をお見せしましょう！

f:id:mattyuu:20161106224054j:plain

！！w(*ﾟoﾟ*)w！！
最初に見たとき、びっくりしました！何か綺麗な図が現れるだろうと予想してたら、本当に綺麗な図が現れたんです！

p=97、p=397とpを大きくしてくと、次のようにパターンが複雑になっていきます！

p=97
f:id:mattyuu:20161106224104j:plain

p=397
f:id:mattyuu:20161106224116j:plain

p=397まで大きくすると、何か星のようなものがはっきり見えます。私の予想では星の境界線はおおよそ[xy=1]の双曲線を伸縮したり、並行移動したりしたものになっているものと思います。今回は可視化した絵を紹介したいだけということで考察は行っていません。

次に除法です。同じルールでp=19を可視化してみましょう！ワクワク！

f:id:mattyuu:20161106224051j:plain

Σ（￣Д￣;）がーんっ!
き、綺麗じゃない(←あくまで主観)、、、切なかったです。。。でもこれも自然が決めた図柄です。「除法は乗法と違って構造が複雑なんだよ」と自然が教えてくれているのではないでしょうか。

乗法のときと同じようにp=97、p=397とpを大きくしてみましょう。

p=97
f:id:mattyuu:20161106224058j:plain

p=397
f:id:mattyuu:20161106224109j:plain

何やら、左上から右下にかけてピンク色の線、左下から右上にかけて白色の線が現れています。やはり除法に関してもまだ考察は行っていませんので、この線が何なのかはまだわかりません。

ここで乗法、除法の可視化を通して特筆したいことがあります。実はどの横一列、どの縦一列を選んでも同じ色が重複しているものはありません*2。どの一列を選んでも、その列の色を並び替えることで白からピンク色へ単調なグラデーションを作ることができます。

そんな厳しい条件がある中で目の前に現れてくれた乗法、除法の図に「ありがとうございましたm(_ _)m」と言いたいところですが、実はこの条件はそんなに厳しいものではありませんでした。

9は素数ではありませんが、パズルで流行っている数独(←ナンプレともいう)もこの条件を満たしますね。この条件がとても厳しいものであれば、世の中の数独の答えのパターンは多く存在することができないため、「この問題どこかで解いたことあるぞ、つまんね」となって、数独はここまで流行らなかったと思われます。
しかも数独は3\times 3の小さい区画に関しても1から9が重複しない(=色が重複しない)というより厳しい条件を満たしていますので、なおさらです。

f:id:mattyuu:20161107003228j:plain

では、n\times nの表で、どの横一列、どの縦一列にも1からnの数(もしくはn色の色)を重複しないように並べる方法は何通りあるのでしょうか。私は小一時間この問題を考えて諦めました。ネットで調べたところこのような条件を満たす表をラテン方陣、またはラテン方格というようで、今現在nを使った解析解は求まっていないようです。あのオイラーも挑戦したことがあるようで、私が少し考えるだけで答えに辿り着けるわけありませんでした。

ちなみにですが、数独のn=9ではラテン方陣は15,224,734,061,278,915,461,120通りあるようです。組み合わせ爆発の凄さを考えると、おそらくn=p=19ではこの宇宙の全星の数、というか全原子の数を優に超えているのではないでしょうか。乗法、除法の図に対しては「想像もできないような多くのパターンの中から、私の目の前に現れてくれてありがとうございましたm( )m」と言うべきなのでしょう。

最後にp=2から397まで素数pを大きくするに従って、可視化した乗法、除法の絵図がどのように変化するかのアニメーションをお見せして終わります。
乗法の方は何やら星がきらめいているような感じがします。

乗法
f:id:mattyuu:20161106224144g:plain

除法
f:id:mattyuu:20161106224205g:plain

ありがとうございました。

f:id:mattyuu:20161107003243j:plain

          • 以下、数学的補足

今回数学的な背景の説明はほとんどしておりません。念のために記載しておいた方が良いと思われる、2点に関して補足コメントします。

【 {\mathbb F}_pにおいて割り算ができること】
{\mathbb F}_5において、4をかけて3になるものを見つけることによって、3\div 4の答えとしましたが、そのような答えはいつでも見つかるのでしょうか。
実はpが素数であることが大いに影響しています。

素数とは限らない一般の正の整数nに対して、今回考えたような乗法を同様に定義することができます。加法も適当に定義することで、{\mathbb Z}_nという足し算、掛け算のできる数学的対象が得られます(数学的には環と言います)。実は{\mathbb Z}_nはnが素数のときに限り0以外の任意の元で割り算を行うことが可能になります(数学的には体(たい)と言います)。nが素数pのとき、割り算ができるという理由から{\mathbb Z}_pを特別に{\mathbb F}_pと書いていたのです。

{\mathbb Z}_nに関して勉強したい方はtsujimotterさんの下記の動画がオススメです！


2-2. 時計の世界の整数論 - 2015/3/27

最後の「ありがとうございました」の画像にはこちらの動画で紹介されている表を可視化し、オイラー基準とフェルマーの小定理が成り立っていそうということを確認しております。

【 乗法、除法の九九表のどの一列を選んでも数字、もしくは色が重複しないこと】
実は割り算ができることと、任意の列(0の列行は除く)の中で数字、色が重複しないことは同値になります。
ここはただの補足コメントということで、割り算ができれば、数字、色が重複しないことだけ簡単に示します。
例えば、乗法のある横一行で数が被ってしまったとします。つまり、b_1\neq b_2に対して、a*b_1=a*b_2が成り立ったとすると、両辺をaで割ることでb_1=b_2となりこれは仮定に矛盾します。

先のコメントで{\mathbb F}_pは割り算ができることを(証明していませんが)伝えていますので、どの一列でも数字、色が重複しないことがわかります。

                                              • -

*1:数学的補足参照

*2:数学的補足参照

マーライオンと初めてのdesmos

10月1日の第7回日曜数学会で数式お絵描きのパイオニアであるもっちょさんとコロちゃんぬさんが初めて出会う場面、そして数式お絵描きの手法について語り合う場面に立ち会うことができました。私はTwitterアカウント(@mattyuu123)を持っていてもつぶやいたことがなかったのですが、この時初めてつぶやきました。私の中では記念の瞬間でした。


私が思うもっちょさんとコロちゃんぬさんの代表作をそれぞれ載せておきます。


この数式お絵描きはdesmosというアプリを用いておられます。私もdesmosを使って何かしてみたいと思い、今回、高校物理の力学、高校数学の範囲でdesmos遊んびしてみたので発表いたします。題材としてはマーライオンを選びました。(数式お絵描きはしておりません)

なお、この記事では高校物理の力学の知識、高校の数Ⅲで習う三角関数の微分の知識を前提としています。前提知識がない方でも、フジテレビのガリレオで福山雅治さんはこんな数式を書いていたんじゃないか？と想像しながら、数式を流し読みしていただけたらと思います。

今回の問題は下記です。

f:id:mattyuu:20161023014914j:plain

マーライオンが水を噴射するときに、一番遠くまで水を飛ばすには口をどの角度にして発射したらよいかという問題になります。

この問題ではマーライオンの身長(正確には口の高さ)をhというパラメータで含めております。実はh=0、つまり地面から直接水を発射する場合であれば、この問題はとても簡単です。高校物理の力学の最初の方に必ず習うと言ってもいいかもしれません。答えは皆さんのご想像通り45°になります。しかし、身長が加わるだけでこの問題は難しくなります。立式は容易ですが、解けるかがわからないのです。数値解であれば容易に出せますが、解析解を出せるかは、冷や汗ものでした。

解く前に一点注意として、この問題では水は大きさのない質点として計算していきます。また空気抵抗も考えません。水は常に連続的に噴射され続けていますが、時刻t=0で水分子が1粒発射されたものとして読んでください。

では解いていきましょう。
(しばらく単調な計算が続きますので、流れだけ簡単に説明すると、
「時刻tでの水の位置を求める」→「水が地面に落ちる時刻Tを求める」→「水が地面に落ちた時の位置を求める」→「水が地面に落ちた時の位置を発射角\thetaの関数と思って微分する」→「微分を0とする\thetaが答えとなる」)

解くためには、座標軸の設定が必要ですので、次の通りとします。

f:id:mattyuu:20161023015005j:plain

まずは、口から飛び出る水の速度を水平方向(x軸右方向)の速度v_xと、垂直方向(y軸上方向)の速度v_yに分けて考えましょう。水平方向には力が加わりませんので、口から飛び出た瞬間の速度のまま変化しません。しかし垂直方向の速度に関しては、重力の影響を受け、時刻tが1増えるごとに重力加速度gずつ小さくなっていきます。
よって、時刻tでの速度はそれぞれ下記になります。

\displaystyle \begin{cases} v_x(t) = v\cos{\theta}\\ \\ v_y(t) = v\sin{\theta} - gt \end{cases}

次に時刻tでの水の位置(x(t),y(t))を出しましょう。高校力学では公式として学んでおり式を導くことができますが、実は、速度を時間で積分することで位置を出すことができます。逆に位置を時間で微分することで速度を出すこともできます。ここでは積分で位置を導きます。公式で覚えている方は結果が一致していることを確認してもいいかもしれません。

\displaystyle \begin{cases} \int_0^t{v_x(t)}{\rm d}t = \int_0^t {v\cos{\theta}}{\rm d}t \rightarrow x(t) = v\cos{\theta}\cdot t + x(0) = v\cos{\theta}\cdot t \\ \\ \\ \int_0^t{v_y(t)}{\rm d}t = \int_0^t {(v\sin{\theta} - gt)}{\rm d}t \rightarrow y(t) = v\sin{\theta}\cdot t -\frac{1}{2}gt^2 + y(0) = v\sin{\theta}\cdot t -\frac{1}{2}gt^2 + h \end{cases}

x(0)、y(0)はそれぞれ左辺の積分を行った結果出てくる、水の時刻t=0での初期位置になります。水は(0,h)より発射されているため、x(0)=0、y(0)=hとなっております。このx(t)、y(t)を連立させtを消去することで、xとyの関係をy=ax^2+bx+cという二次関数の形で表すことができます。これは、二次関数が放物線と言われる所以(ゆえん)になります。今回の計算では二次関数の形にする必要がないのでこれ以上は立ち入りません。

次に、水が地面に落ちる時刻Tを求めます。これはy(t)=0となるtを求めればよいということになります。y(t)=0とすると、

\displaystyle \begin{split} & v\sin{\theta}\cdot t - \frac{1}{2}gt^2 + h = 0 \\ & \rightarrow gt^2-2v\sin{\theta}\cdot t - 2h = 0 \\ & \rightarrow t = \frac{v\sin{\theta} \pm \sqrt{v^2\sin^2{\theta}+2gh}}{g} \end{split}

中学校で習った二次関数の解の公式を、それなりに現実っぽい問題を解く時に使えるのは嬉しいですね。t\geq 0より、水が地面に着く時間Tは、

\displaystyle T = \frac{v\sin{\theta} + \sqrt{v^2\sin^2{\theta}+2gh}}{g}

となります。

水が地面に落ちた位置を(L,0)としましょう。つまり、(x(T),y(T))=(L,0)です。Lは次の式で表すことができます。

\displaystyle \begin{split} L = x(T) = v\cos{\theta}\cdot T & = \frac{v\cos{\theta}}{g}\left(v\sin{\theta} + \sqrt{v^2\sin^2{\theta}+2gh}\right) \\ & = \frac{v^2}{g}\left(\cos{\theta}\sin{\theta} + \cos{\theta}\sqrt{\sin^2{\theta}+\alpha}\right) \end{split}

ここで、今後計算を進めるにあたって記載を容易にするために、\alpha\equiv 2gh/v^2としています。今考えている変数は\thetaだけであり、g、h、vはそれぞれ定数ですので、このように置いてしまって中身が分からなくなってもこの後の計算には影響しません。

このLを\theta変数の関数L(\theta)と思って、その最大値を求めていきます。最大、最小を求める方法は微分が定石です。微分して得られる導関数を0とする\thetaが、最大、または最小を与えます。


\displaystyle \frac{{\rm d}L}{{\rm d}\theta} = \frac{v^2}{g}\left(-\sin^2{\theta}+\cos^2{\theta} - \sin{\theta}\sqrt{\sin^2{\theta}+\alpha}+\cos{\theta}\frac{\sin{\theta}\cos{\theta}}{\sqrt{\sin^2{\theta}+\alpha}}\right)

微分を0にする\thetaを求めたいため、{\rm d}L/{\rm d}\theta=0とおき、両辺をv^2/gで割ることで、下記のようになります。

\displaystyle -\sin^2{\theta}+\cos^2{\theta} - \sin{\theta}\sqrt{\sin^2{\theta}+\alpha}+\frac{\sin{\theta}\cos^2{\theta}}{\sqrt{\sin^2{\theta}+\alpha}} = 0

\sqrt{}が出てきおります、、、これを満たす\thetaを求めることはできるのでしょうか。。。

まずは\sqrt{}を含むものと、含まないもので分けましょう。

\displaystyle -\sin^2{\theta}+\cos^2{\theta} = \sin{\theta}\sqrt{\sin^2{\theta}+\alpha}-\frac{\sin{\theta}\cos^2{\theta}}{\sqrt{\sin^2{\theta}+\alpha}}

両辺に\sqrt{\sin^2{\theta}+\alpha}をかけましょう。

\displaystyle \left(-\sin^2{\theta}+\cos^2{\theta}\right)\sqrt{\sin^2{\theta}+\alpha} = \sin{\theta}\left(\sin^2{\theta}+\alpha-\cos^2{\theta}\right)

両辺を二乗しましょう。

\displaystyle \begin{split} & \left(-\sin^2{\theta}+\cos^2{\theta}\right)^2\left(\sin^2{\theta}+\alpha\right) = \sin^2{\theta}\left(\sin^2{\theta}+\alpha-\cos^2{\theta}\right)^2 \\ & \rightarrow \left(\sin^4{\theta}-2\sin^2{\theta}\cos^2{\theta}+\cos^4{\theta}\right)\left(\sin^2{\theta}+\alpha\right) \\ & \qquad = \sin^2{\theta}\left(\sin^4{\theta}+\alpha^2+\cos^4{\theta}+2\alpha\sin^2{\theta}-2\alpha\cos^2{\theta}-2\sin^2{\theta}\cos^2{\theta}\right) \\ \\ & \rightarrow \sin^6{\theta}-2\sin^4{\theta}\cos^2{\theta}+\cos^4{\theta}\sin^2{\theta}+\alpha\sin^4{\theta}-2\alpha\sin^2{\theta}\cos^2{\theta}+\alpha\cos^4{\theta} \\ & \qquad = \sin^6{\theta}+\alpha^2\sin^2{\theta}+\sin^2{\theta}\cos^4{\theta}+2\alpha\sin^4{\theta}-2\alpha\sin^2{\theta}\cos^2{\theta}-2\sin^4{\theta}\cos^2{\theta} \\ \\ & \rightarrow -\sin^4{\theta}+\cos^4{\theta}-\alpha\sin^2{\theta}=0 \end{split}

途中計算間違いしそうな、ややこしい式変形を経て、なんとも綺麗な形になりました。さて、これ以上先に進めますでしょうか？



進めるんです！

まず、最初の2項は次のように変形できます。

\displaystyle \begin{split} & -\sin^4{\theta}+\cos^4{\theta}=\left(\sin^2{\theta}+\cos^2{\theta}\right)\left(-\sin^2{\theta}+\cos^2{\theta}\right) \\ &\qquad =-\sin^2{\theta}+\cos^2{\theta}=\cos{2\theta} \end{split}

ここで、「和と差の積は二乗の差」と中学時代に呪文のように習った展開(因数分解)の公式と、高校で習う\sin^2{\theta}+\cos^2{\theta}=1、及び\cosの加法定理である\cos{(\alpha+\beta)}=\cos{\alpha}\cos{\beta}-\sin{\alpha}\sin{\beta}が役に立っております。

なんと基礎学力が大事かがわかる場面でしょうか！

私が中学校の先生であったなら、「この因数分解はマーライオンが吐き出す水を解析するために必要なんだよ」と生徒たちに教えたいところです。

後ろの1項はどうでしょうか?

こちらも高校で習う、\cosの加法定理を応用することで、

\displaystyle -\alpha\sin^2{\theta}=-\alpha\frac{1-\cos{2\theta}}{2}

と変形できます。(\cos{2\theta}=\cos^2{\theta}-\sin^2{\theta}=1-2\sin^2{\theta}を用いてます。)

以上から、{\rm d}L/{\rm d}\theta=0は、

\displaystyle \begin{split} & \cos{2\theta}-\alpha\frac{1-\cos{2\theta}}{2}=0 \\ & \rightarrow \cos{2\theta}=\frac{\alpha}{2+\alpha} \end{split}

となります。ここで、\alpha=2gh/v^2であったことを思い出すと、最終的に、

\displaystyle \cos{2\theta}=\frac{\frac{2gh}{v^2}}{2+\frac{2gh}{v^2}}=\frac{gh}{v^2+gh}

なんとも、綺麗にまとまった式でしょうか！

高校生の皆さんは、「求めたい答えは、\cosをとると\frac{gh}{v^2+gh}となる角度を半分にしたものです。」と答えれば正解になります。

大学では\arccosという関数が登場し、\cosをとるとxになる角度を\arccos{x}と表すことができます。

よって、
\displaystyle 2\theta=\arccos{\frac{gh}{v^2+gh}}

とできるので、最終的な答え(求めたい角度を\theta_{\rm max}とでも名前をつけましょう)をまとめると、次のようになります。

f:id:mattyuu:20161023022325j:plain

前回記事に書いたトーナメントの数式同様、今回も無事解析解を得ることができたことはなんとも不思議であります。うまく式変形ができなければ、コンピュータによる数値解しか得られなかった可能性は大いにあるのです。

さて、忘れかけていたdesmosを使う時がついにきました。desmosを使って、マーライオンにいろいろな角度で水を発射してもらいましょう。そして今回求めた飛距離が最大になる角度の方向にはオレンジ色の点線を引き、最大到達地点(L=x(T)=L(\theta_{\rm max})には赤い旗を置いておきましょう。

そうすると下記のようになり、たしかに今回求めた解析解が正しそうということがわかっていただけるかと思います。

f:id:mattyuu:20161023025014g:plain

初速vをもっと大きくしても、下記のように問題なさそうですね。

f:id:mattyuu:20161023025102g:plain

今回私が作成したdesmosは下記リンクから遊ぶことができます。

https://www.desmos.com/calculator/vhvcloibx9

もっちょさんやコロちゃんぬさんに比べたらとても簡単なものですが、このように解析した結果をアニメーションで見ることができることには感動を覚えます。私が学生時代であれば、Javaを使って水の軌道を数値計算し、gnuplotを使って静的な画像を何枚も作っていたことでしょう。

また、この数式は小便小僧やサクランボの種飛ばし世界選手権(あるんですね、、)でもそのまま使えます(さくらんぼの種は質点ではなく剛体なので、別の力学を用いる必要がある点には目をつむります)。数学の力を感じますね。

読んでいただき、ありがとうございました。

f:id:mattyuu:20161023021715j:plain

トーナメントは運か実力か(解析解編)

先日の第7回日曜数学会で福原和朗さんが「トーナメントは運か実力か」という大変面白いLTをしておられました。今回は福原さんのLTを受けて、私が見つけた数式について発表させていただきます。

福原さんのLTに関しては下記のリンクのブログにまとめられております。

発表

第７回日曜数学会で発表してきました。 - kazurof weblog

福原さんは後述するある確率の数値解(←コンピュータ等を用いた計算による数値で与えられる解)を発表されたのですが、LTの後に、「必ず確率を文字式で解を表せるよね(つまり解析解(←問題を特徴付ける変数で表された代数的な解)が得られるよね)」と2人で会話させていただきました。先日のMathPowerで行われた素数大富豪のトーナメント表を見てそのことを思い出し、無事に解析解を見つけることができましたので、ここに発表したく思った次第です。

本記事の数学面の説明は、高校レベルの確率(特にコンビネーション\mathrm Cの使い方、赤い玉と白い玉が入っている袋から玉を取り出す系の確率の問題の解き方、何人かを一列に並べて特殊な並べ方の場合の数を求める方法)を習得している人を想定して記載しております。ただ高校確率の知識がない方でも問題の設定は理解していただけると思いますので、こういった現実的な問題を数学を使って解くことができるということ、つまり数学の力を感じていただけたら幸いです。

始まり、始まり
f:id:mattyuu:20161016085420j:plain
↑
「トーナメントは運か実力か(解析解編)」というタイトルで発表させていただきます。

夏の高校野球やリオオリンピックの柔道は負けたら終わりのトーナメント方式ですね(オリンピックの柔道に関しては敗者復活戦があるため完全に負けたら終わりというわけではありませんが)。トーナメントでは第1回戦で実力が1位と2位の人が戦ったら、2位の人は初戦敗退となってしまいます。また、実力がない人でもたまたま対戦相手が自分よりも下の実力であれば、第二回戦、第三回戦と勝ち進むことができます。理不尽だな、と思った経験はないでしょうか。

(余談ですが、駒大苫小牧高校は2004年、2005年と夏の甲子園を連覇しております。そしてなんと2006年の夏も甲子園の決勝まで行っているんです。この2006年の決勝こそ歴史に残る引き分け再試合となった、「ハンカチ王子こと斎藤佑樹」VS「マー君こと田中将大」なんです。夏2連覇して、その翌年の夏も決勝まで進むとは、駒大苫小牧の実力やばいですね。)

数学の確率を使えば理不尽なことがどれほど起こり得るのか、ということを定量的に表すことができます。では問題設定を行います。

f:id:mattyuu:20161016085425j:plain
↑
nを正の整数として、 2^n人が参加するトーナメントを考えます。ただしトーナメントは上図のような綺麗な形のトーナメントを考え、シードがあったり、参加者によって優勝するまでに必要な勝利数が異なったりしないものとします。

参加者を 2^n人(つまり2人、4人、8人、16人、、、)という人数に縛るのは綺麗な形のトーナメントにするために必要な条件となります。 2^nという制約を課すことで第1回戦、第2回戦、第3回戦、、、という風に第 k回戦終わるごとに勝ち残っている参加者の数が半分になり続け、第 n回戦で1人になる(つまり優勝者が決まる)綺麗な形のトーナメント表を作成することができます。綺麗な形のトーナメントという数学的ではない表現が出てしまっておりますが、その意味するところはわかっていただけると思います(ちなみに数学的には完全2分木という言葉で表します)。

f:id:mattyuu:20161016085429j:plain
↑
参加者には「強さ」を示す数値があらかじめ与えられており、各試合では「強さ」が大きい参加者が必ず勝つものとします。
「強さ」は1から 2^nまでの整数値をとり、参加者の「強さ」はそれぞれ異なるものとします。つまり、1から 2^nまでの各整数値が 2^n人の参加者一人一人に対して重複なく割り振られているものとします。

「強さ」 2^nを持つ参加者はランキング1位、「強さ」 (2^n-1)を持つ参加者はランキング2位、「強さ」1を持つ参加者はランキング 2^n位(つまりビリ)となります。毎回『「強さ」 xを持つ参加者』と書くのは大変ですので、今後は『「強さ」 xを持つ参加者』のことを xさんと呼ぶことにしましょう。

f:id:mattyuu:20161016085449j:plain
↑
上述した設定の元、本記事の問題はこのようになります。

トーナメント表をランダムに組むという意味は、 2^n人の参加者をランダムに一列に並べ、トーナメント表の一番下に据えるということになります。

各試合では「強さ」が大きい参加者が必ず勝つため、 2^nさん(つまり「強さ」 2^nを持つ参加者)が必ず優勝しますし、1さん(つまり「強さ」1を持つビリの参加者)は必ず初戦敗退します。2さんは初戦が1さんでなければ第2回戦に進みますが、第2回戦では必ず負けるでしょう。 P_{a,b}で aさんがちょうど b回勝つ確率を表しますから、上で話した簡単な考察ですぐに下記のことがわかります。問題が解けた暁には、下記が成り立っているかを確認することでその妥当性を判定できそうです。

2^nさんは絶対優勝する。必ずちょうど n回勝つ。逆にちょうど k回( 0\leq k \leq n-1)勝つことはない。
P_{2^n,n}=1
P_{2^n,k}=0

1さんは必ず初戦敗退。逆にちょうど k回( 1 \leq k\leq n)勝つことはない。
P_{1,0}=1
P_{1,k}=0

2さんは初戦が1さんでなければ第1回戦勝てるが、第2回戦は必ず負ける。初戦の相手が1さんでなければ初戦敗退する。
(トーナメント表はランダムのため、例えば初戦の相手が1さんになる確率は、参加者から2さんを除いた (2^n-1)人のうち1さんだけということで、 1/(2^n-1)となります。)
P_{2,0}=(2^n-2)/(2^n-1)
P_{2,1}=1/(2^n-1)

f:id:mattyuu:20161016085452j:plain
↑
具体的な例を見つつ、何が問われているかを確認しましょう。上図は n=3の例になります。

このトーナメントではランキング2位の7さんが初戦敗退している上、ランキング5位の4さんが決勝戦に進んでいるなど、かなり理不尽なことが発生しております。このような理不尽なことがどれほど発生するのかということを、確率という道具を使って評価しよう、ということになります。

以上で問題設定を終わります。イメージは掴めましたでしょうか。

f:id:mattyuu:20161016085456j:plain
↑
「 aさんがちょうど b回勝つ確率 P_{a,b}を求めよ」と言われてもどこから手をつけてよいかわかりません。 2^n人が参加するトーナメント表は単純に考えると 2^n!(=2^n \cdot (2^n-1)\cdots3\cdot2\cdot1)通りのパターンがあり、考えようとするだけでも頭がこんがらがりそうです。
今、 aさんがちょうど b回勝つ確率を考えたいので、まずは aさんだけに注目することで考察を始めるための入り口を模索していきたいと思います。

トーナメント表をランダムに組むと aさんの最初の位置は様々です。しかし上図のように aさんを隣の zさんと入れ替えた2つのトーナメント表は本質的に同じものであるといえます。同様に aさんと zさんの2人分の組で構成されたトーナメントのブロックを隣のブロックと入れ替えてもトーナメント表は同じものであると言えるでしょう。トーナメント表は参加者が勝ち進んでいくと、どこを勝ち進んだ勝者と戦うことになるのかを示すものであり、ブロックの左右の入れ替えには影響を受けないため、このような操作をしてもトーナメント表の本質は変わりません。 aさんができるだけ左に行くように各階層でブロックを入れ替えていくと、元のトーナメント表と本質的に同じトーナメント表で aさんが一番左下に配置されているものを得ることが可能になります。つまり aさんは始めからトーナメント表の一番左下に固定してよいということがわかります(数学的にはこのことを「 aさんを左下に固定しても一般性を失わない」と表現します。)。

aさんを一番左下に置くトーナメント表を考えることで何かいいことがあるでしょうか。まず最初に感じられるメリットは問題のサイズを落とすことができたということになります。元々 2^n!通りのトーナメント表の可能性がありましたが、 aさんを左下に固定することで残りの (2^n-1)人の並べ替えである、 (2^n-1)!通りに問題のサイズを落とすことができました。つまり問題のサイズが 1/2^nになっています。例えば n=5の場合を考えると 1/32の時間で結果を出すことができるわけで、これは数値シミュレーションを行う場合には大変都合の良いメリットになります。

一般性を失わないようにもっと場合の数を減らすことは可能でしょう(例えば第一回戦の対戦では左側に強さが小さい人が来るように配置する等)。ただこれ以上の一般化は数値シミュレーションをする上でメリットになっても、今回のように解析的に解を求めたい場合は必ずしもメリットになりません。場合分けが多く発生し、論理的に考える思考の妨げになるからです。実は今回の一般化は問題のサイズを落とす以上のメリットがあります。 aさんを左下に持っていくことで、 aさんが b回勝つための必要条件を簡単に考察することができるのです。

f:id:mattyuu:20161016085500j:plain
↑
本格的な考察に入る前に、今後の考察、記載を容易にするため、上図のようにトーナメント参加者にラベルをつけておきます。

f:id:mattyuu:20161016085508j:plain
↑
では aさんが b回勝つための必要条件を求めましょう。 aさんをトーナメントの左下に固定した威力が発揮されるときがきました。

必要条件というのは aさんが b回勝つために最低限満たさなければいけない条件となります。つまり必要条件を満たしても aさんがちょうど b回勝つとは限りません。 (b+1)回以上勝つ可能性もあるわけです。「 aさんが b回勝つための必要条件」というのは「 aさんが最低 b回勝つための必要十分条件」ということになります。なおしばらくの間、「 aは1ではない」、かつ「 bは0でない」として話を進めます。 a=1または b=0は、 aさんが初戦敗退するという少し趣が異なる事象になるため、考察がしにくくなります。

まずは aさんが1回勝つ必要条件(つまり初戦突破する条件)を求めましょう。数学では文字式を使った一般的な議論の前に具体的な数値でイメージをつかんで一般化するという手法が有効です。 aさんが1回勝つためには初戦の相手である X_1が aさんよりも弱い(=「強さ」の小さい)1さんから (a-1)さんであることが明らかに必要条件になります。つまり aさんが1回勝つための必要条件は「 X_1が1さんから (a-1)さんのいずれかである。」ということになります。

次に aさんが2回勝つ必要条件を求めましょう。 X_1が aさんよりも弱い参加者であることは明らかに必要ですが、 aさんは第二回戦ですぐ隣のブロックの X_2と X_3の勝者と戦うわけですから、 X_2、 X_3のいずれも aさんよりも弱い1さんから (a-1)さんのいずれかであることが必要条件になります。つまり aさんが2回勝つための必要条件は「 X_1、 X_2、 X_3が1さんから (a-1)さんのいずれかである」ということになります。

同様に考察を行うと aさんが b回勝つための必要条件は「 X_1〜 X_{2^b-1}が1さんから (a-1)さんのいずれかである」ということが分かるかと思います。

めでたしめでたしとしたいところですが、ここでもう一つ考えないといけない条件があります。例えば5さんは4回勝つことができるでしょうか。5さんが4回勝つためには X_1〜 X_{2^5-1}(=X_{31})が1さんから4さんである必要があります。 X_1から X_{31}に重複なく1さんから4さんを割り振ることは不可能です。つまり aさんが b回勝つためには a>2^b-1も必要になるわけです。

例えば5さんがb回勝つためには
5>2^b-1
が満たされる必要がありますから、
2^b<6
となり、これを満たす bは1、2のみとなります。つまり5さんはどれだけ頑張っても3回以上勝つことができないのです。

以上まとめるとaさんがb回勝つための必要条件は、
「 a>2^b-1」 かつ 「 X_1〜 X_{2^b-1}が1さんから (a-1)さんのいずれかである」
となります。

f:id:mattyuu:20161016085514j:plain
↑
最終的に求めたいものは aさんがちょうど b回勝つ確率の P_{a,b}ですが、その一歩手前として、 aさんが最低 b回勝つ確率 \tau_{a,b}を求めましょう。( \tauという文字に特別の意味はなく、最初に私がノートで計算した際にたまたま用いた文字で、ここでもそのまま使わさせていただきます。)

前章で aさんが b回勝つための必要条件、つまり aさんが最低 b回勝つための必要十分条件が得られております。これを丁寧に確率の計算に当てはめてあげれば良いです。上図に記載した通り、確率の定義を思い出しておきましょう。

まず、求めたい確率の分母である全体の場合の数を求めます。今左下に aさんを固定しているため、トーナメント表の組み合わせとして起こりうる全体の場合の数は、 aさんを除く (2^n-1)人を一列に並べる並べ方の総数である (2^n-1)!に等しくなります。一列に並べて左の人から順に X_1、 X_2、 X_3、、、 X_{2^n-1}と割り振っていけばよいのです。

f:id:mattyuu:20161016085519j:plain
↑
次に確率の分子になる aさんが最低 b回勝つトーナメントの総数を求めましょう。

まず、 X_1から X_{2^b-1}に aさんより弱い1さんから (a-1)さんまでを割り振る場合の数を求めます。これは (a-1)人の中から (2^b-1)人を選出し、選出された (2^b-1)人を一列に並べる場合の数の総数になります。 (a-1)人から (2^b-1)人を選出する組み合わせは {}_{a-1}{\mathrm C}_{2^b-1}通り、その各々に対して選出された (2^b-1)人を一列に並べる総数は (2^b-1)!ですから、この場合の数は {}_{a-1}{\mathrm C}_{2^b-1}\cdot(2^b-1)!となります。
そして、その各々に対し参加者の割り振られていない X_{2^b}から X_{2^n-1}へ (2^b-1)人に選出されなかった残りの \{(2^n-1)-(2^b-1)\}人、つまり (2^n-2^b)人を割り振ります。この並べ方は aさんが最低 b回勝つという事柄に影響しませんので、単純に一列の並べ方の総数を出せばよく (2^n-2^b)!となります。以上まとめると、 aさんが最低 b回勝つトーナメントの総数は {}_{a-1}{\mathrm C}_{2^b-1}\cdot (2^b-1)!\cdot(2^n-2^b)!となります。

f:id:mattyuu:20161016085524j:plain
↑
確率の分母、分子が求まりましたので、 aさんが最低 b回勝つ確率 \tau_{a,b}は上図の通りとなります。 {}_n{\mathrm C}_r=n!/(n-r)!r!を使うことで、なんともシンプルで綺麗な表記になりました。

なお、最後の \prodは \piの大文字であり、高校で習ったシグマ( \sum)の掛け算版の記号として使われます。
つまり、
\displaystyle \prod_{k=1}^na_k=a_1a_2\cdots a_n
です。

a<=2^b-1のときは、 aさんは絶対に b回勝てないので、 \tau_{a,b}=0となります。

今、 a\neq 1で考察していました。 a=1つまり、1さんは最弱のため1回も勝つことができません。つまり b>0となる bに対して、明らかに \tau_{1,b}=0です。
そして b>0なる bに対して、 a=1<=2^b-1ですから、今回求めた \tau_{a,b}は a=1に対しても成り立つことがわかります。

f:id:mattyuu:20161016085528j:plain
↑
さて、本題であった aさんがちょうど b回勝つ確率 P_{a,b}を求めましょう。ただし、ここでも b\neq 0とします。 b=0に対しては次章で求めたいと思います。

以前にも注意した通り、 \tau_{a,b}は aさんが最低 b回勝つ確率であり、ちょうど b回勝つ確率ではありません。しかし下記の考え方から、 \tau_{a,b}から、簡単に P_{a,b}を求めることができます。

トーナメントが一つ与えられたとき、 aさんが初戦敗退するか、ちょうど1回勝つか、ちょうど2回勝つか、ちょうど3回勝つか、、、ちょうど n回勝つか(つまり、優勝するか)はどの2つの事象も同時に起こることはありません。「ちょうど3回勝つ」と「ちょうど4回勝つ」という事象は同時に起こらないのです。このことを確率の言葉で排反事象といいます。

「 aさんが最低 b回勝つトーナメントの総数」は、「 aさんがちょうど b回勝つトーナメントの総数」と「 aさんがちょうど (b+1)回勝つトーナメントの総数」と「 aさんがちょうど (b+2)回勝つトーナメントの総数」、、、「 aさんがちょうど n回勝つトーナメントの総数」の和になります。

このことから(両辺をトーナメントの総数で割ることで)、「 aさんが最低 b回勝つ確率」=「 aさんがちょうど b回勝つ確率」+「 aさんがちょうど (b+1)回勝つ確率」+「 aさんがちょうど (b+2)回勝つトーナメントの確率」+・・・+「 aさんがちょうど n回勝つ確率」となるのです。

今求めたいものは「 aさんがちょうど b回勝つ確率」 P_{a,b}ですが、今具体的に求まっている \tau_{a,b}には、上図の通り邪魔な確率「 aさんが (b+1)回以上勝つ確率」が含まれてしまっております。しかし、「 aさんが (b+1)回以上勝つ確率」= \tau_{a,b+1}ですので、 \tau_{a,b}から邪魔な確率 \tau_{a,b+1}を引いたものが、 P_{a,b}になるのです。

f:id:mattyuu:20161016122736j:plain
↑
これまでは b\neq 0を考えておりました。最後に b=0つまり、 aさんが初戦敗退する確率を求めましょう。これまでの議論に比べればとても簡単なものになります。

上図では少し丁寧に確率を求めましたが、 aさんが初戦敗退する必要十分条件は、「初戦の対戦相手の X_1が aさんより強い」になります。 aさんの初戦の相手となり得る参加者の総数は、 2^n人の参加者から aさんを除いた 2^n-1です。そのうち aさんより強い人は (2^n-a)人おりますし、初戦の相手が誰になるかは aさんを除く参加者全てに対して同様に確からしいため、 aさんが初戦敗退する確率は、
\displaystyle P_{a,0}=\frac{2^n-a}{2^n-1}
になります。

f:id:mattyuu:20161016085537j:plain
↑
以上の議論をまとめると、最終的な解析解は上図の通りになります。誰かが注文したわけではないのに、このような綺麗な形になるところに数学の美しさを感じてしまします。

a、 bに具体的な値を入れて計算すると、とてもよくできた数式だと感心します。例えば、1さんが必ず初戦敗退することや、 2^nさんが必ず優勝するということを上記の式から計算して確認してみてください。一体この数式は誰が考え、生み出したのでしょうか。とても神秘的な気持ちになります。

今回の所感として、問題が解けたことに驚きを感じています。この世の中には未だ解けない難問が数多く存在します。一見簡単そうな問題でも、解けていないものが多くあります。今回の問題もたまたま解くことができましたが、解けるという保証はありませんでした。簡単に解ける問題、解けない問題の違いは何か、なぜ今回解けたのか、考えるのも怖いくらい難しい問題な気がします。

今回得られた解析解からどんな面白いことがわかるかという考察はしておりません。それはまた機会があれば考察したいと思います。例えば「一回戦敗けの人の言い訳はどれだけ信用してもよいか？」や「 n\rightarrow \inftyにして、 x=b/nを [0,1]区間の連続変数と考えた時の P_a(x)の分布はどうなるか?」等々。後者の分布を出せば、福原さんが日曜数学会で発表されていた、「2位の人は運ゲーである」ということについても目に見える形でわかるのではないかと思います。

f:id:mattyuu:20161016132135j:plain

リーマンゼータ関数のゼロ点を手計算してみた

10/1に第7回日曜数学会で掲題の件でLTをさせていただきました。
拙い発表でしたが、自分の好きな数学に関して発表する機会を作っていただき、日曜数学会の幹事のキグロさん、tsujimotterさん、塩鯖さんに厚く御礼を申し上げます。

また、休憩時間や懇親会にてその他多くの参加者の皆様とも話させていただき、とても有意義な時間が過ごせました。
数学の話題しか出ない懇親会はとてもいいものですね。
皆様ありがとうございました。

私の発表に対してtwitterでつぶやいてくださった方がいらっしゃり大変嬉しかったです。
ありがとうございます。
私はつぶやいたことがなかったのですが、つぶやいてもらえると嬉しいということがわかりましたので、今後日曜数学会に参加した際は全力でつぶやかさせていただきますのでよろしくお願いします。

さて、昨日の発表ですがせっかくなので公開したくブログを作成しました。
日曜数学会のLTは5分という短い時間で伝えたいことを伝えなければいけません。
5分に全力をかける、5分で伝えたいことを凝縮するという点でとても勉強になりましたが、本当は説明したかったところを泣く泣く削ることがありました。
この記事ではそういったところも含めて発表資料のスライドを紙芝居のように使いながら説明させていただこうと思います。
今後リーマンゼータ関数のゼロ点を手計算したい方がいらっしゃいましたら、ぜひ読んでいただけたらと思います。

(本来こういった発表は前提とする数学的知識を明確にすべきなのですが、今回はあまり考えておりませんでした。。。世間一般では複素数とはどれだけ浸透しているのでしょうか。。複素数とは何かといった説明は全くありません。。)
(実際に紙芝居のようにまとめてみると文字数が1万字を超えてしまい、相当のモチベーションがある人じゃないと最後まで読んでもらえないだろうと思います。ただこういう風にまとめておくのは今後の自分のためにもなると思い、公開します。)

始まり、始まり
f:id:mattyuu:20161002134126j:plain
↑
「リーマンゼータ関数のゼロ点を手計算してみた」というタイトルで発表させていただきます。

f:id:mattyuu:20161002134143j:plain
↑
リーマンゼータ関数とは何かというと、18世紀にオイラーが研究し、1859年にリーマンの記念碑的論文「与えられた数より小さい素数の個数について」で{\zeta(s)}という記号で名付けられた関数で、自然数の s乗の逆数の無限和で定義されます。 sが変数で、 \zeta(←ゼータと読みます)が関数名となります。注意したいことはこの関数は通常中学校で学ぶ関数とは異なり、複素変数で複素数値を返す関数になるところです。つまり、変数 sと\zeta(s)は複素数になります。

例えばs=2とすると、\zeta(2)は平方数の逆数の無限和になりますが、なんとこの無限和は\frac{\pi^2}{6}と等しくなります! これは誠に驚くべきことで、小学生でも理解できて、円とは全く関係のないように見える分数の足し算の先に、\piが現れるのです。この級数が無限大にならず収束することは容易に分かるのですが、誰もその値に辿り着けずにいました。1735年に人類で初めてオイラーがこの級数が\frac{\pi^2}{6}と等しくなることを発見、証明しました。とても美しいですね。私は学生時代に学んだ数式の中で2番目にこの数式が好きです(1番はコーシーの積分公式です。)。

関数の記号は\zetaなのに、この関数の名前には"リーマン"もしくは"リーマンの"という風に接頭語がつきます(接頭語という表記が正しいかはわかりませんが。)。これは後の時代に様々な\zeta関数が見つかっているからのようです。私はまだ勉強中の身であるためこれ以上は語ることができません。

f:id:mattyuu:20161002134148j:plain
↑
今回なぜリーマンゼータ関数のゼロ点をテーマに選んだかをお伝えします。
それは、リーマンゼータ関数のゼロ点は、人類が抱える問題の中で間違いなく最難関の部類に入る未解決問題リーマン予想と密接に関係しているためです。リーマンゼータ関数のゼロ点とは\zeta(s)=0となるようなsのことであり、マイナスの偶数がゼロ点になることは容易に証明できます。これは自明なゼロ点と呼ばれます(スライドの水色の点となります。)。

\zeta(-2)=0 ?と不思議になられた方もいらっしゃるかもしれません。それは\zeta(s)を自然数の逆数和で表記した時に、明らかにゼロにならないように見えるからです。
{ \displaystyle \zeta(-2) = \sum_{n=1}^{\infty} \frac{1}{n^{-2}}=1+2^2+3^2+4^2+\cdots = 0 (?) }
実は\zeta(s)を自然数の逆数和で表記できるのは\rm{Re}(s)>1の時だけであり、\rm{Re}(s)\leq1に対しては解析接続という技術を使って関数を表記しなおさなければいけません。ここでは「解析接続が必要なんだよ」というところで説明は省かせていただきます。

さて話を戻すと、マイナスの偶数という自明なゼロ点以外にもリーマンゼータ関数のゼロ点が存在します。それらは非自明なゼロ点と呼ばれます(スライドのピンク色の点となります。)。リーマン予想とは「リーマンゼータ関数の非自明なゼロ点は一直線上に並んでいるはずだ」という予想になり、この一直線とは実部が\frac{1}{2}の直線になります(スライドのピンク色の点線になります)。

リーマン予想が生まれて160年弱、世界中の大天才達が挑み、未だ証明できていません。このことを考えるだけでもゼロ点が凄いものということがわかっていただけると思います。

「天才達が挑んだ歴史なんてどうでもいい、そもそもなぜゼロ点が重要か? ゼロ点がどこかわかって何かいいことある?」と思われる方もいらっしゃるかもしれません。すみません、私の力不足で詳しく色々語ることができないのですが、世の中にはリーマン予想が正しければ正しいという定理が何千もあるようです。それだけリーマン予想は数学の根底に関わる問題なのでしょう。リーマン予想を証明することは多くの数学者が最後のピースを埋められずに証明できていない、数々の定理を証明することになるのです。

f:id:mattyuu:20161002134153j:plain
↑
次に、なぜゼロ点を手計算したくなかったを説明します。
2009年のNHKスペシャル「魔性の難問　リーマン予想・天才たちの闘い」を見たところ、「実際のゼロ点の位置を、まず4つほど求めたリーマンは予想外の事実に気づきます。」というナレーションと共に、リーマンが驚いている描写がされておりました。また当時のリーマンの計算用紙も紹介されておりました。リーマンゼータ関数は無限和や無限区間の積分で定義されます。1859年というコンピュータのない時代にどうやってゼロ点を手計算したんだろう? と疑問に思ったことがきっかけです(今回、実際にはリーマンの計算方法ではなくグラームの計算方法で計算しております。)。


f:id:mattyuu:20161002134159j:plain
↑
リーマンゼータ関数のゼロ点を手計算する方法は私が探した限りネット上では見つけることができませんでした。

なお、「リーマン予想 零点 手計算」でgoogle検索を行うと、tsujimotterさんのブログが一番に出てきました。

ジーゲルのZ関数を数値計算する

ジーゲルのZ関数を数値計算する - tsujimotterのノートブック

私がご紹介するのもおこがましいのですが、tsujimotterさんは日曜数学者の生みの親であり、私の憧れであります。数学会では3Dプリンタでリーマンゼータ関数を可視化したことで有名です。ゼロ点がなぜ重要か? ということに関しても、下記の記事にて、リーマンの素数公式の右辺に出てくるためであるという新たな知見を私に与えてくれました。(ここまで不思議と素数について記載していなかったのですが、リーマンゼータ関数のゼロ点は数学者たちが愛してやまない素数の情報をすべて持っているとも言えるのです。)

リーマンの素数公式を可視化する

リーマンの素数公式を可視化する - tsujimotterのノートブック

ただtsujimotterさんのブログではゼロ点をruby(←プログラミング言語)を用いて計算しております。私はコンピュータのない時代にどうやって計算したかを知りたいため、tsujimotterさんのブログだけでは疑問を解消することができませんでした。

f:id:mattyuu:20161002134205j:plain
↑
リーマン予想とゼロ点計算の歴史を簡単に紹介します。リーマン予想が1859年に世に生み出されたことは既にお伝えした通りですが、リーマン予想は1900年のヒルベルトの23の問題、2000年のミレニアム懸賞金問題に選出されております。
ヒルベルトの23の問題とは、1900年にドイツのヒルベルトという有名数学者が20世紀に解くべき数学の23個の未解決問題をまとめたものです。同様にミレニアム懸賞金問題とは、2000年にクレイ研究所が定めた数学上の7つの重要未解決問題のことであり、解決した数学者には100万ドルの懸賞金が与えられるとのことです。1900年、2000年と100年経っても数学会の最重要未解決問題であり続けるリーマン予想のラスボス感を思い知らされます。

驚くことに1900年のヒルベルトの23の問題が発表された当時、人類は非自明なゼロ点に関する数値情報を1つも持っていなかったようです(リーマンただ一人を除いて)。初めて(リーマンを除く)人類が非自明なゼロ点を目にしたのは、1903年グラームの研究結果となります。今回は鹿野健先生編著の「リーマン予想」(日本評論社)を参照して、グラームの方法でゼロ点を手計算しました。

なお、リーマンは研究結果を論文として公表しますが、その結果に至った証明、計算過程のメモを公表することはなかったようです(←リーマン予想に関してだけかもしれません)。そのため1900年代にリーマン予想を証明しようとしてできなかった数学者達に、「リーマン予想は成り立たないのであろう、リーマンは計算していないのだから、リーマン予想は夢想だろう」という風に思われたことがあるようです。1932年ジーゲルというドイツの数学者がその驚くべき洞察力で、リーマンが残した断片的なメモを解読し、リーマンがゼロ点を実際に手計算していたという事実、またゼロ点を計算するのにグラームの方法より有用な公式(←今日ではリーマン・ジーゲルの公式と呼ばれています)を導いていたことを発見します(NHKスペシャルのナレーションには嘘がなかったことになります。)。リーマンは頭脳だけで定理を見つけるような天才ではなく、実際には膨大な手計算を元に定理を導くエレファントな一面もあったわけです。

今回、リーマンが実際に用いたリーマン・ジーゲルの公式は用いず、グラームの方法でゼロ点を計算しています。これはリーマン・ジーゲルの公式がリーマンの天才的頭脳を持ってして思いつく積分路の変更と、エレファントで複雑な計算を経て得られる異形な式であり、ぱっと見なぜこれを計算することでリーマンゼータ関数のゼロ点が計算できるかわからないためです。それに比べグラームの方法は初等的と言え、リーマンゼータ関数を計算しているという実感を得ることができる方法になっており、日曜数学会の発表ではグラームの方法の方が適切と考え、そのようにいたしました。

f:id:mattyuu:20161002134213j:plain
↑
リーマンゼータ関数のゼロ点を計算する中で困難が2つあります。まず1つ目はリーマンゼータ関数は無限個の足し算で定義されていることです。上記のように足し算を\frac{1}{999999999999^s}まで計算したとしても、まだまだ足し算を止めるわけにはいきません。

f:id:mattyuu:20161002134226j:plain
↑
そこでグラームが採用したのがオイラー・マクローリンの和公式です。この公式自体はオイラーが導いたものでリーマンゼータ関数とは直接関係ありません。ここではこういう便利な公式があるというくらいの認識で問題ありません。

f:id:mattyuu:20161002134233j:plain
↑
オイラー・マクローリンの和公式にて、N=\infty、f(z)=1/{z^s}として、
\displaystyle \zeta(s)=\sum_{n=1}^{M-1}\frac{1}{n^s}+\sum_{n=M}^{\infty}\frac{1}{n^s}
とすることで、上記の式を導くことができます。

まず主要項と呼ばれる部分ではリーマンゼータ関数を最初のM-1まで真面目に計算します。その次の補正項と呼ばれる部分で計算していない残りの部分のk次の近似が与えられます(kを大きくすることで誤差が小さくなります)。近似の精度はM、kを大きくすることで上げることができます。スライドでは説明していませんが、誤差項R_{2k}は下記のように上から評価できます(\sigma={\rm Re}(s)です)。

\displaystyle \left|R_{2k}\right| \leq \left|\frac{s+2k+1}{\sigma+2k+1}\right|\left|\frac{s(s+1)(s+2)\cdots(s+2k)}{(2k+2)!}\right| \left|B_{2k+2}\right|M^{-\sigma-2k-1}

なお、計算したい非自明なゼロ点(s=1/2+it)のtが大きくなるにつれて、十分な精度でゼロ点を計算するためにはM、kをどんどん大きくしていく必要があります。例えばt=10であれば、M=10、k=10とすることで、誤差項の絶対値を10^{-15}よりも小さくできますが、t=50となると、M=10、k=10としても、誤差項の絶対値が10^{-2}程度まで大きくなってしまいます。実際に手計算をやってみればわかりますが、M、kが1大きくなるだけで、計算をする気力がどっと落ちてしまいます。これはグラームの方法の欠点となっています。

とはいえ、オイラー・マクローリンの和公式を使うことで無限個の和という掴みどころのなかったリーマンゼータ関数を、誤差の大きさを評価できる有限回の計算に落とし込むことができました。これは大きな進歩と言えます。

f:id:mattyuu:20161002134240j:plain
↑
ここまで数式のみの続きイメージがしにくかったですが、実際にs=1/2+5i(つまりt=5)、M=5、k=3とした場合の展開式を紹介します。式は長いですが、これを実直に計算することで\zeta(1/2+5i)を手計算できるのです。

f:id:mattyuu:20161002134248j:plain
↑
ところで、
\displaystyle \frac{1}{2^{\frac{1}{2}+5i}}
はどうやって計算するのでしょうか。
実はまたまたオイラーが見つけた公式である、
\displaystyle e^{i\theta}=\cos{\theta}+i\sin{\theta}
を用いて、上記のように計算できるのです。
なお、オイラーの公式にて\theta=\piとして右辺を左辺に移行した、
\displaystyle e^{i\pi}+1=0
は、ネイピア数e、円周率\pi、虚数単位i、和の単位元0、積の単位元1が1つの式で結ばれており、その美しさから人類の至宝と呼ばれています。

思いがけないところで\sin、\cos、\logが出てきました。これは簡単には手計算できないぞ、手計算諦めよう、とも考えましたが、コンピュータのなかったリーマンの時代でも数表と呼ばれるものがあり、代表的な関数の値を数表から得ることができます。過去に数表を作ってくださった方々に多大な敬意を払いつつ、発表を続けます。

f:id:mattyuu:20161002134253j:plain
↑
次の困難は複素数値関数であるリーマンゼータ関数では中間値の定理が使えないことです。私はあまり数値計算法に明るくはないのですが、中間値の定理は関数のゼロ点を求めるために、とても重要なものであろうと容易に想像がつきます。例えば上の例のようにf(1)=-1、f(2)=1であれば、xが1から2まで増加する間にその符号が変わっているため1<\alpha<2となる\alphaが存在して、f(\alpha)=0になるでしょう。f(1)、f(2)の符号が異なるだけで、ゼロ点が存在することを知ることができるのです。

f:id:mattyuu:20161002134302j:plain
↑
しかしリーマンゼータ関数は実数値関数ではなく複素数値関数です。実数と違い複素数は2次元の平面と同一視できます。2次元の複素平面上では-1と1を0を通らずに結べるため、上記の例のようにf(1)=-1、f(2)=1であっても、f(\alpha)\neq 0 (1<\alpha<2)ということも考えられるのです。

f:id:mattyuu:20161002134308j:plain
↑
そこで登場するのが、ジーゲルのZ関数です。複素数値関数で中間値の定理が使えないのであれば、実数値関数にしてしまおうという考えです。ここで細かく説明はできないのですが、{\rm Re}(s)=1/2上で、Z(t)と\zeta(1/2+it)はその絶対値が一致します。つまり、
\displaystyle Z(t)=0 \Longleftrightarrow \zeta(1/2+it)=0
が成り立つのです。複素数値関数で中間値の定理が使えない困った\zeta(1/2+it)ちゃんのゼロ点を計算する代わりに、実数値関数のZ(t)のゼロ点を計算すれば良いということです。これも大きな進歩です。

なお、このZ関数にはジーゲルのZ関数と名付けられておりますが、その存在はリーマンが論文に書いていた関数等式から初等的に導出することができます。また、ジーゲルより前の時代のグラームがゼロ点の計算に用いていることから、ジーゲルはこの関数の発見者ではないと思われます。ジーゲルがリーマンのメモを解読する中で何かこの関数に重大な意味づけした等の功績で名前が付いているものと私は考えております。

f:id:mattyuu:20161002134315j:plain
↑
Z(t)の計算のために必要な\theta(t)もガンマ関数のGaussの公式、
\displaystyle \Gamma(s)=\lim_{N\rightarrow\infty}\frac{N!N^s}{s(s+1)\cdots(s+N)}
やオイラー・マクローリンの和公式を用いることで上記のように近似計算ができます。実際にはその後e^{\theta(t)}=\cos{\theta(t)}+i\sin{\theta(t)}を計算しないといけないのですが、数表を使えばなんとかできそうです。

f:id:mattyuu:20161002134320j:plain
↑
ついに武器は揃いました。tに対して、オイラー・マクローリンの公式から\zeta(1/2+it)を計算し、その後Z(t)を計算します。いくつかのtに対してZ(t)を計算して、例えばZ(t_1)<0、Z(t_2)>0であれば、複素平面上でs=1/2+it_1とs=1/2+it_2を結ぶ線分上に\zeta(s)=0となるゼロ点が存在するのです。

f:id:mattyuu:20161002134326j:plain
↑
今回はコンピュータのなかった時代にどうやってゼロ点を計算したのかを知りたいため、実際に計算する前に上記のルールを策定しておきます。

f:id:mattyuu:20161002134331j:plain
↑
こちらが私が実際に計算したtになります。最初はt=1,2,3,4,\cdots,20と計算していこう、その後符号が変わったあたりをt=X.1,X.2,\cdots,X.9と計算していき
とゼロ点をそれなりの精度で求めたかったのですが、実際にやってみるとその大変さからたった3つのtに留まってしまいました。。。

なお、それぞれのtに対して、M、kを上記のように選べば、その誤差を10^{-5}程度まで小さくすることができます。

f:id:mattyuu:20161002134337j:plain
↑
計算結果をお伝えする前に白状しないといけないことがあります。計算を進めている中で突然のルール改定が行われました。計算を続け大量の筆算で嫌気がさしてきた頃に、急遽45.9727÷101が現れました。「÷100にしたいよ、÷101とほとんど結果一緒でしょ、、もう嫌だ電卓使っちゃおう」という流れでルール改定を決心しました。電卓を使っていいと言っても1ステップの計算に限ります。そこは譲れません。

「昔の人は偉いな、すごいな。自分たちは便利なものに囲まれて何か大事なものを見失っている気がする。」そういう風に思いました。

f:id:mattyuu:20161002134344j:plain
↑
こちらが最終結果になります。\zeta(1/2+it)、Z(t)それぞれに対して、私の手計算の結果と、Wolframというアプリで計算した正確な値を比較しております。自分でも驚くほど計算結果が近いものになっておりました。面白いと感じたのはZ(t)は実数値関数になるため、私の手計算の結果の虚部もかなり小さい値になっていたことです。

計算時間に関してはt=5の途中で突然のルール改定が行われたため、t=10、t=15に対しては計算時間が短くなっております。

Z(10)<0、Z(15)>0となっていることから、10<\alpha<15となる\alphaが存在して、\zeta(1/2+i\alpha)=0となる非自明なゼロ点が存在することがわかりました！！！！

今回はこれ以上精度を上げて計算はできませんでした。。。
今回はと書きましたが、今後もやりません。こんな計算はコンピュータを使えば一瞬です。便利なものは使います。

f:id:mattyuu:20161002134355j:plain
↑
計算をしていて気づいたことは私の計算用紙とリーマンの計算用紙の見た目がなんとなく似ているということです。リーマンという知の巨人に少し近づけた気がしました。

f:id:mattyuu:20161002134401j:plain
↑
こちらは今回大変参考になった参考文献になります。私が使った式は全てこの参考書からの引用です。

f:id:mattyuu:20161002134409j:plain
↑
お読みいただきありがとうございました。


最後にちょっとした感想です。

実は虚部が一番小さい非自明なゼロ点はs=1/2+14.1\cdots i程度になるのですが、最初のゼロ点が私たち人類が日常的に手の届く範囲に存在していることは奇跡ではないかと思いました。最初のゼロ点の虚部が14程度である必然性はありません(少なくとも私は知りません)。ということは数学の神様がいたずらをして例えば最初のゼロ点が1/2+10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000iであれば、コンピュータを使って計算を続けても非自明なゼロ点を1つも見つけられないということがあり得たのではないかと思います。非自明なゼロ点を具体的に1つも挙げられない中、全ての非自明なゼロ点が一直線上に並んでいることを証明するのは相当心が参るのではないかと。

もう一つ最後に、実は下のようなスライドも用意していました。わかる人にはわかる数学ネタになっております。
f:id:mattyuu:20161002184153j:plain

ありがとうございました。