読者です 読者をやめる 読者になる 読者になる

mattyuuの数学ネタ集

世界は数式で溢れている

有限体Fp上の楕円曲線'のパズル

はじめに

先日職場の勉強会でRSA暗号楕円曲線暗号について発表をしました。面白いことに話の全体を通してフェルマー(17世紀のフランスのアマチュア数学者)が登場しました。

フェルマーパスカルと共に確率論を創始するなど、上記の暗号関連の話以外にも重要な仕事を行なっております。フェルマーは17世紀の人ですが、現代社会の根っこの部分に彼が与えた、与えている影響は大きそうです。ただ、今回はこれ以上フェルマーの話はしません。勉強会の資料作成を通して、有限体上の楕円曲線を可視化したのですが(←グラフ書いただけですが)、可視化した楕円曲線を眺めたり、いじったりしているうちに、「パズルとして売り出せるんじゃない?」というような面白い性質を見つけたのでその紹介をしたいと思います。

まず、楕円曲線とは下記で定義されます。

楕円曲線の定義

abを定数とし、方程式
\displaystyle
y^2=x^3+ax+b
を満たす(x,y)全体の集合(←グラフと呼ぶ)を楕円曲線という。

注1:この記事では上記の一般形を得るために、標数が2、3以外の体K上で考えています。abKの元です。
注2:(右辺のxの3次式)=0は重解を持たないものとします。

中学校や高校で習う関数のグラフは、実数体(\mathbb R)という足し算、引き算、掛け算、割り算ができる世界*1の上で考えています。実数上で楕円曲線のグラフを描くと、例えば下記の図のようになります。

a=1b=1の場合

f:id:mattyuu:20170129155536p:plain

a=-1b=0の場合

f:id:mattyuu:20170129155547p:plain

2つ目のy^2=x^3-x (a=-1b=0)の例でグラフの描き方を簡単に説明します。

x=1のときは、y^2=0となりますからy=0です。(x,y)=(1,0)に赤い点を描きます。同様にx=0-1のときもy=0となるので、(x,y)=(0,0)(-1,0)にも赤い点を描きます。

x=2のときは、y^2=6となるので、y=\pm\sqrt{6}です。(x,y)=(2,\sqrt{6})(2,-\sqrt{6})に赤い点を描きます。先程と違い、1つのxに対して2つのyが対応しました。

x=1/2のときは、y^2=-3/8となるので、これを満たす実数yは存在しません。そのため実数上のグラフではx=1/2のときは描く点はありません。

x-\inftyから\inftyまで動かしながら、各xに対して方程式を満たすyを同様に求めて(x,y)に点を描いていくと、先程の2つ目のグラフになるのです。

有限体{\mathbb F}_p上の楕円曲線

楕円曲線暗号では、楕円曲線実数体(\mathbb R)上ではなく、有限体{\mathbb F}_p上で考えます。有限体{\mathbb F}_pが何かわからない方は以前書いた下記の記事を読んでみてください。

mattyuu.hatenadiary.com

{\mathbb F}_pは簡単に言えば集合\{0,1,2,\cdots,p-1\}に適切に足し算、引き算、掛け算、割り算を定義した世界になります。今回の記事では足し算と、掛け算しか使わないため、例として{\mathbb F}_5の足し算表、掛け算表を載せておきます。

足し算表

+ 0 1 2 3 4
0 0 1 2 3 4
1 1 2 3 4 0
2 2 3 4 0 1
3 3 4 0 1 2
4 4 0 1 2 3

掛け算表

× 0 1 2 3 4
0 0 0 0 0 0
1 0 1 2 3 4
2 0 2 4 1 3
3 0 3 1 4 2
4 0 4 3 2 1

有限体{\mathbb F}_pの復習が終わりましたので、実際に{\mathbb F}_p上の楕円曲線を描いてみましょう。ここでは例としてp=5a=1b=1とします。

f:id:mattyuu:20170129155625j:plain

楕円曲線というのに、グラフは繋がってないし、曲がった感じもしませんが、楕円曲線の定義に反してませんし、描き方も実数体上の時と全く一緒です。実数体では-\inftyから\inftyまでの無限個のxに対して、yを求めないとグラフは描けませんでしたが、有限体{\mathbb F}_p上ではp個のxに対してyを求めればよいことが本質的に異なります。

では実際に上のグラフがy^2=x^3+x+1のグラフになっていることを確認してみましょう。

x=0のとき、y^2=1となるので、2乗して1になる数求めればよいことになります。{\mathbb F}_5上では0^2=01^2=12^2=43^2=44^2=1より、x=0のときはy=14となります。よって(x,y)=(0,1)(0,4)のマスにピンク色を塗っております。

同様に例えばx=4のとき、y^2=4となります。{\mathbb F}_5上では2^2=43^2=4より、(x,y)=(4,2)(4,3)のマスにピンク色を塗ります。

様子が違うのはx=1のときです。このときy^2=3となりますが、{\mathbb F}_5上では2乗して3になる数は存在しません。そのためx=1の縦列のセルにはピンク色に塗るセルはありません。これは先述した実数体上の楕円曲線の例で、y^2=-3/8の解が存在せず描く点がなかったことと同じです。

今回有限体上の楕円曲線javascriptでプログラムを作って描画しています。プログラムの良いところは一回プログラムを作ってしまえばパラメータ(ここではpab)を変えたグラフも容易に描画できることです。手計算、手描きでグラフを作成すると、パラメータを変えるごとに改めて計算を行う必要があります。またpを大きくすると計算自体が大変になります。

さあプログラムの利点を活かして{\mathbb F}_5上の楕円曲線'を全て描画してみましょう。abはそれぞれ\{0,1,2,3,4\}の5通りから選ぶことができるので、全部で25通りの楕円曲線'が得られます。

{\mathbb F}_5上の楕円曲線'たち
f:id:mattyuu:20170129155643p:plain

p=23にしても描けます!
{\mathbb F}_{23}上の楕円曲線'たち(出力結果を一部抜粋、実際には23^2=529通りの楕円曲線'が現れます。)
f:id:mattyuu:20170129155655p:plain

はぁはぁ、、楽しい。。

ここで一点白状します。上記には楕円曲線以外のグラフも含まれています。楕円曲線の定義の注2で記載しましたが、y^2=x^3+ax+bの右辺のxの3次式について、(右辺のxの3次式)=0という方程式を作ると、この方程式は重解を持たないという前提があります。例えば{\mathbb F}_5上ではx^3+2x+2=(x-1)^2(x-3)となり、重解を持ってしまうので、{\mathbb F}_5上でy^2=x^3+2x+2楕円曲線ではありません。

しかし、今回考えたパズルは全てのabを対象にしないと成立しません。誠に情けないですが楕円曲線の集合に、楕円曲線とは言ってはいけない{\mathbb F}_5上のy^2=x^3+2x+2といった曲線も追加して楕円曲線'という名前で広い範囲の集合を考えます。

気づいた性質

さて、楕円曲線'のグラフを並べてみると気づくことがあります。それはpaを固定して、b0からp-1まで動かしたp個のグラフを重ね合わせると全マスがちょうど1回ずつピンク色に塗られるのです!

言葉だけではわからないと思いますので、つまりこういうことです。
f:id:mattyuu:20170129155713j:plain

つまり、こういうことです。
f:id:mattyuu:20170129155805g:plain

つまり、こういうことです。
f:id:mattyuu:20170129155831g:plain

最後の例ではピンク色の部分はセロハンのような素材と思ってください。セロハンを重ねるごとにピンク色が濃くなります。aが同じ値の時は各セルが足並みをそろえて濃くなっていく様子が確認できます。

p=23でも、綺麗に成り立ちます。(ファイルサイズが大きくなってしまうので、アニメは途中で止めています。)
f:id:mattyuu:20170129155851g:plain

はぁはぁ、、面白い。。

p=97とか、p=691とかもっと大きい素数では楕円曲線'のグラフもより複雑になります。複雑なグラフがp個あり、それを重ね合わせると綺麗なピンク一色!何かに役立つわけではないですが、数学の力を感じさせられます!

この性質が成り立つ証明は超簡単です。残念ながらこの性質は楕円曲線ではない方程式でも成り立ちますし、もっと言えば体{\mathbb F}_p上じゃなく剰余類環{\mathbb Z}_n上でも成り立ちます。

証明

{\mathbb F}_p上の楕円曲線':y^2=x^3+ax+bに対してa=a_0\in {\mathbb F}_pを1つ固定する。この時、下記2点を示せばよい。

  1. \forall (x_0,y_0) \in {\mathbb F}_p \times {\mathbb F}_pに対して、\exists b_0 \in {\mathbb F}_p{y_0}^2={x_0}^3+a_0x_0+b_0となる。
  2. b_0\ne b_1 \in {\mathbb F}_pのとき、y^2=x^3+a_0x+b_0y^2=x^3+a_0x+b_1は共通の解を持たない。
  • 1の証明:b_0={y_0}^2-{x_0}^3-a_0x_0とすればよい。
  • 2の証明:(x_0,y_0)が共通の解とすると、{y_0}^2={x_0}^3+a_0x_0+b_0{y_0}^2={x_0}^3+a_0x_0+b_1がそれぞれ成り立つ。辺々を引くと、b_0=b_1となり、これは前提に矛盾。

証明はいかにも当たり前ですが、可視化してみないと気づかなかった性質なので、満足感はあります。

応用例

楕円曲線'に関して面白い性質を見つけたので、何か応用したくあります。

まずは、大人向けのパズルとして活躍の余地はないでしょうか。p=5のときは、25種類の楕円曲線'が描かれたカードを用意します。色が塗られていない部分は透過するようにしたいためアクリル板で作成しましょう。25枚のカードから適切に5枚ずつ選んで、それぞれの組で各カードを適切に回転、反転させて重ね合わせることで、ピンク一色のカードの組が5セットできたら成功です。*2

f:id:mattyuu:20170129155916j:plain

p=5が余裕でできるようになったら、p=7p=11と難易度を上げていけば素数の数だけ楽しめます!

次に思いつくのは子供向けの絵合せパズルです。2〜3歳の時はaを固定して、5枚のカードでくまさんの絵合せをしましょう。

f:id:mattyuu:20170129155959g:plain

幼稚園に行きだしたら、ゴレンジャーを題材にしましょう。この頃になると25枚あっても大丈夫なはずです。

f:id:mattyuu:20170129160032g:plain

このように子供の時から素数楕円曲線'に慣れ親しみ純粋に育った男の子は、プロポーズにもメッセージカードと称して楕円曲線'カードを使用してしまう恐れがあります。

f:id:mattyuu:20170129160100g:plain

おまけ

今回は有限体上の楕円曲線を可視化しました(←グラフ書いただけ)。可視化すると楕円曲線の方程式を満たす解の個数も一目瞭然ですが、この解の個数がフェルマーの最終定理の証明に深く関係しています。また、この解の個数はハッセ・ヴェイユ境界と呼ばれる、より高度な数学にも関連しています。私はステートメントはなんとなくわかるレベルです。

今回楕円曲線とは何かということに関しては何も触れていません。楕円曲線google検索すると、やはりtsujimotterさんの記事がwikipediaの次に出てきましたので、ここにも貼り付けておきます。楕円曲線の面白いネタがたくさん書かれています。興味がある方は一緒に楕円曲線の勉強をしましょう。
tsujimotter-sub.hatenablog.com

ついでながら、今回ゴレンジャーのアニメはjavascriptcanvasで作成したのですが、画像を回転させるのにこれまたtsujimotterさんの下記の記事が大変役に立ちました。原点を中心としてcanvas全体を回転させるのは簡単にできるのですが、画像を中心として、その画像だけを回転させる方法はネットで調べてもよくわからず、ネットサーフィンを続けているうちに下記の記事に行き着きました。参考にこちらも貼り付けておきます。

tsujimotter.hatenablog.com

f:id:mattyuu:20170129160143j:plain

*1:数学的には四則ができる世界を体(たい)と呼ぶ

*2:このパズルを売り出すためには解の一意性を証明しておく必要があります。

蔵本予想とは何か? 〜予想紹介編〜

1/7に第8回日曜数学会にて「蔵本予想とは何か?」というタイトルで力学系理論に関するLTをさせていただきましたのでご紹介します。私の不勉強のため今回は蔵本予想の紹介に留まりましたが、近いうちに証明も体得して概説したいと思い、本記事では〜予想紹介編〜というサブタイトルを付けております。

f:id:mattyuu:20170108202050j:plain

蔵本予想とは同期現象に関する数学的な予想になります。"予想"と書いていますが、九州大学の千葉逸人准教授によって証明されています。今回は蔵本予想の概説を行いたいと思います(不勉強ゆえ間違った記載があるかもしれません。あくまでざっくり蔵本予想がどういうものなのかを知っていただければと思います。)。


f:id:mattyuu:20170108202111j:plain

まず同期現象とは何かを説明します。多数の蛍が集まっている場面、または、多数のメトロノームを集めた状況(←台車の上に置く等、相互に作用するようにする。固い台の上においてはダメ。)を想像してください。最初はそれぞれの蛍が各々のタイミングで発光し、メトロノームが各々のペースで針が振れているのに、十分な時間が経つとまるで一匹の蛍のように、まるで1つのメトロノームのように発光のタイミング、針の触れ方、ペースが揃ってくることが観測されます。これを同期現象といいます。

youtubeに蛍の美しい同期現象がありましたので、リンクを貼っておきます(開始35秒くらいから同期現象の動画が始まります。)。

www.youtube.com


f:id:mattyuu:20170108202123j:plain

こういった同期現象を語る数学的なモデルに蔵本モデルがあります。これは京都大学蔵本由紀教授により提案されたものです。

蔵本モデルは上記のような微分方程式で表されます。3ステップでこの式が意味するところを解説したいと思います。


f:id:mattyuu:20170108202137j:plain

まずステップ1です。蛍の点滅やメトロノームの針が左右に振れる運動のように、同期現象の元となる単体の運動は周期的な運動です。周期的な運動ということで、質点が一定速度で円周上を回る運動をモデルの出発点としましょう。

このモデルの運動を記述する変数として質点が円周の中心となす角\thetaをとります。\thetaは時間tの経過と共に変化していきますので、変数\thetaはまた時間tの関数とも思えます。その意味で\theta\theta(t)とも書きます。

次に回るスピード(自然振動数と呼びます。単位時間当たりどれだけの角度回るかという値です。)を\omegaとすると、このモデルを記述する微分方程式は、

\displaystyle
\frac{{\rm d}\theta}{{\rm d}t}=\omega

となります。

自然振動数\omegaが大きければ大きいほど質点が円周を回るスピードが速くなっていきます。\omegaは固体によって異なる値を取り得ます。例えば短い間隔で点滅を繰り返すせっかちな蛍の\omegaは大きく、長時間かけて点滅を繰り返すのんびり屋の蛍の\omegaは(絶対値が)小さくなります。

この微分方程式は簡単に解けて下記のように表せますが、力学系理論ではシステム(系)の状態を微分方程式で表すことが通常ですので、次のスライド以降も微分方程式の形で進めていきます。

\displaystyle
\theta(t)=\omega t + \theta_o

ステップ1では1つの質点の円周上の運動を考えましたが、質点が1つでは同期現象を考えようがありません。では、次のステップに進みます。


f:id:mattyuu:20170108202151j:plain

ステップ2では同期現象を発生させるために質点をたくさん、ここではN個集めます。

質点をN個に変更することに伴い、変数\thetaも増やす必要がありますので、添え字iを付けて\theta_i(i=1,2,\cdots,N)とします。そして自然振動数も固体によって異なってもよいものでしたので、変数\theta_iに対応する自然振動数も\omega_iのように添え字を付けて表します。

ステップ2で質点の数をN個に増やしましたが、\theta_iの運動を定める微分方程式には\theta_iの運動の自然振動数\omega_iしか出てこず、\theta_j(j\neq i)は登場しません。つまりそれぞれの質点は自分のペース(自然振動数\omega_i)で回り続けるだけであり、まだ同期現象は発生しません。では、次のステップに進みます。

f:id:mattyuu:20170108202205j:plain

ステップ3では質点の運動を定める微分方程式に上図の項を加えます。\theta_iの運動を決める微分方程式\theta_1から\theta_Nが登場することで、初めて変数間の相互作用が発生します。これで蔵本モデルの完成です。

この相互作用項の主要な部分は\sin{(\theta_j-\theta_i)}になります。

\theta_i微分方程式の右辺の相互作用項を考察していることに注意して、\theta_j(j\neq i)\theta_iより少しだけ大きい場合*1\sin{(\theta_j-\theta_i)}は正になりますので、\theta_i\theta_jから回るスピードが速くなるように相互作用を受けます。

逆に\theta_j\theta_iより少しだけ小さい場合は、\theta_i\theta_jから回るスピードが小さくなるように相互作用を受けます。

\sin{(\theta_j-\theta_i)}という相互作用により\theta_iさんと\theta_jさんがお互いにお近づきになりたい気持ちがあるということがわかっていただけると思います。実際は\theta_1から\theta_N全てから相互作用を受けるので、\theta_jだけの相互作用を考えるだけでは、回るスピードが速くなるのか遅くなるのかはわかりません。

\sin{(\theta_j-\theta_i)}\sum_jで和を取っています。質点の個数Nが大きくなっても相互作用項が大きくなり過ぎないように相互作用項の分母にNが入っています。

また、Kは相互作用項の大きさを決めるパラメータになります。このKが大きければ大きいほど同期現象が起こりやすくなります。つまりKは、(次スライドに続く)


f:id:mattyuu:20170108202228j:plain

つまりKは、空気(Kuki)を読む力です。

図のように\theta_jさんに比べて少しだけ遅れている\theta_iさんは\theta_jさんにおいて行かれないように頑張って速く回ろうとします。逆に\theta_jさんは\theta_iさんが遅れていることに気づいて、ペースを落としてあげようとします。

このように空気を読んで相手を思いやる力の大きさ、それがKなのです。

Kを大きくするに従って、同期が発生しやすくなる。これが蔵本予想の要となります。


f:id:mattyuu:20170108211232j:plain

ここで、同期の発生具合を定量的に把握するための指標を定義しておきます。

N個の質点が半径1の円周上を運動している時に、N個の質点の重心*2と円の中心との距離rを同期具合と呼ぶことにします*3

図では青い点が質点の重心となっています。全く同期が発生していない場合の同期具合r0ですが、同期現象が現れると同期具合rはだんだん1に近づいていきます。

数式ばかり見ていてもイメージが湧かないので、空気を読む力Kを大きくしていくと、同期具合rも大きくなっていくことを実際にシミュレーションをして確認したいと思います。そのためのモデルとしては、(次スライドに続く)


f:id:mattyuu:20170108202451j:plain

シミュレーションのモデルとしてはクマさんの着ぐるみのバイトの子を使いましょう*4。図のバイトの子はクマさんの被り物を一定のペースで付けたり外したりする趣味を持っています。そんなバイトの子を30人集めて、各々自分のペースで趣味を楽しんでもらいましょう。ただいつもと違うのは周りに同じ趣味を持つ子が29人いることです。空気を読む力Kを大きくしていくと、どうなるでしょうか?

【空気を読む力K=0
みんな自分のペースで楽しんでいます。
f:id:mattyuu:20170108194310g:plain

【空気を読む力K=1
変わりはないように見えます。
f:id:mattyuu:20170108194359g:plain

【空気を読む力K=2
みんなのペース、タイミングが同じになってきています。
f:id:mattyuu:20170108194433g:plain

【空気を読む力K=10
何かのお遊戯会のように足並みが揃っています。
f:id:mattyuu:20170108194458g:plain

確かにKを大きくしていくと同期具合が大きくなっていっていることがわかりますね。では最後に蔵本予想とは何かをざっくり説明します。


f:id:mattyuu:20170108202508j:plain

蔵本予想
N\rightarrow \inftyの極限で、空気を読む力Kと同期具合rの関係は図のグラフのようになる。

ポイント1:
空気を読む力Kが小さい間は同期具合rは0のまま。空気を読む力がある値K_c(←自然振動数の確率密度関数だけに依る値)を超えると、同期具合rが大きくなり始める。

ポイント2:
空気を読む力KK_cに近い時、同期具合r\sqrt{K-K_c}に比例する。

ポイント1のように状態の変化具合がK_cを境に急に変わること(微分が不連続になること)を相転移といいます。固体が液体に、液体が気体になる現象と同じ定義です。みんなが多少空気を読んだところで全く同期が起きません。空気を読む力を上げていくと、急に同期が起こり始め「え?さっきまでみんな同期してなかったのに急にどうしたんだよ?」という状態になります。これが蔵本予想なのです。(先ほどのシミュレーションではK=1では同期が起きてなさそうでしたが、K=2では明らかに同期が発生していました。)

今回の記事はインターネットに公開されている千葉逸人准教授の下記論説の最初の2ページくらいを日曜数学会発表用に自分なりに解釈してまとめたものです。

http://www2.math.kyushu-u.ac.jp/~chiba/paper/sugaku2016.pdf

私の理解が至らず説明不足の点*5もありますし、解釈、説明が間違っている箇所もあるかもしれません。今後、何かの文献で証明を体得できたあかつきには、正しいステートメントで記事を書きたいです。まずは、「同期現象に関する蔵本予想というものがあって、それは空気を読む力と同期具合に関するものなんだ。そして、千葉逸人准教授が証明されたんだ。」ということを知っていただければ幸いです。

おまけ

今回行ったバイトの子のシミュレーションをするためのソースコードは下記になります。
ソースコードと言ってもただのhtml(javascriptも含みますが)なので、Internet ExplorerSafariなどのブラウザがあれば自分のパソコンに閉じた世界で遊べます。(Safariでしか動作確認していません。ご了承ください。)

ソースの内容をkumamoto.html等の"html"拡張子でローカルに保存し、ダブルクリックで実行することで動かせます。(ただし、バイトの子と、クマさんの画像ファイルに代わる画像ファイルを何か2つ用意し、htmlファイルと同階層においてください。その後、ソース内の画像ファイル名を指定する箇所をそのファイル名で置換してください。)

<html lang="ja">
  <head>
    <meta charset="utf-8" />
    <title>Kuramoto Simulator</title>

    <!-- javascriptで実装する。 -->
    <script type="text/javascript">
      const NUM = 30; // 質点の数(変更してもいいが、バイトの子の描画は対応できていない。)
      const WIDTH = 400;  // 描画領域の横幅
      const HEIGHT = 400;  // 描画領域の縦幅
      var omega = new Array(NUM); // 自然振動数
      var theta = new Array(NUM); // 変数θ
      var theta_next = new Array(NUM);  // 次の時刻でのθを計算するための一時変数
      var r =  new Array(NUM);  // 質点の色を保持する変数 red
      var g =  new Array(NUM);  // 質点の色を保持する変数 green
      var b =  new Array(NUM);  // 質点の色を保持する変数 blue
      var ctx_circle;  // Canvasの2Dコンテキスト(使い方しかわかっていない)
      var ctx_kuma;

      var RAD = 150;  // 円周の半径
      var delta_t = 0.01; // 計算するときの1ステップの時間幅
      var K = 0.0;  // 空気を読む力 K

      var timer;  // startボタン、stopボタンで制御するためのタイマー変数(使い方しかわかっていない)

      var img_woman = new Image(); // バイトの子の画像
      var img_kuma = new Image(); // クマさんの画像


      // startボタンを押した時の処理
      // 一度止めて、Kの値を更新して再び動かす。
      function start(){
        clearInterval(timer);
        K = document.getElementById('K_VALUE').value;
        // 33msec間隔でdraw(描画)メソッドを呼び出す
        timer = setInterval(draw, 33);
      }

      // stopボタンを押した時の処理
      function stop(){
        clearInterval(timer);
      }

      // ネットで見つけた正規分布に従う乱数を発生させる関数
      /**
        * 正規分布乱数関数 参考:http://d.hatena.ne.jp/iroiro123/20111210/1323515616
        * @param number m 平均μ
        * @param number s 分散σ^2
        * @return number ランダムに生成された値
       */
      function seiki(m, s) {
        var a = 1 - Math.random();
        var b = 1 - Math.random();
        var c = Math.sqrt(-2 * Math.log(a));
        if(0.5 - Math.random() > 0) {
          return c * Math.sin(Math.PI * 2 * b) * s + m;
        }else{
          return c * Math.cos(Math.PI * 2 * b) * s + m;
        }
      }

      // 初期処理 htmlが読み込まれたタイミングで処理が走る
      function init(){
        // Canvasの2Dコンテキストの取得
        ctx_circle = document.getElementById('circle').getContext('2d');
        ctx_kuma = document.getElementById('kuma').getContext('2d');

        // 各質点の色、自然振動数、初期値の設定
        for(var i = 0; i < NUM; i++){
          
          // r,g,bいずれも乱数にすると見た目が良くないので、rは255に固定する。
          r[i] = 255;//Math.floor(Math.random() * 256);
          g[i] = Math.floor(Math.random() * 256);
          b[i] = Math.floor(Math.random() * 256);

          // 自然振動数は平均1、分散1になるように設定
          omega[i] = seiki(1.0,1.0);
          
          theta[i] = Math.random()*2*Math.PI;
        }

        // 画像ファイルの設定
        img_woman.src = "pose_kiri_woman.png"; // TODO なんでもいいので画像ファイルに差し替えてください
        img_kuma.src = "kuma.gif"; // TODO なんでもいいので自分が持っている画像ファイルに差し替えてください

        // 描画をスタートする。
        start();
      }

    // 描画メソッド(この中で微分方程式の差分の計算もしている。。。メイン部分) 
    function draw(){

      // 微分方程式の計算
      keisan();
      // 円周上の運動の描画
      draw_circle();
      // バイトの子の描画
      draw_kuma();
    }

    // 微分方程式の計算
    function keisan(){
      // 10ステップを計算する。1ステップごと描画していた遅々とアニメが進まなかったため。
      for(var l=0;l<10;l++){
        // 微分方程式の差分の計算部分
        for(var i = 0; i < NUM; i++){
          // 相互作用のsin部を足すための一時変数
          var sum_sin = 0.0;
          for(var j=0;j<NUM;j++){
            sum_sin += Math.sin(theta[j]-theta[i]);
          }
          theta_next[i] = theta[i] + (omega[i]+K*sum_sin/NUM)*delta_t;
        }
        // 次のθの値で現状のθを更新
        for(var i = 0; i < NUM; i++){
          theta[i] = theta_next[i];
        }
      }
    }

    // 円周上の運動の描画
    function draw_circle(){
      // 画面描画の更新仕様に関するおまじない
      ctx_circle.globalCompositeOperation = "source-over";
      ctx_circle.fillStyle = '#FFFFFF';
      ctx_circle.fillRect(0, 0, WIDTH, HEIGHT);

      // 円周の描画
      ctx_circle.beginPath();
      ctx_circle.fillStyle = '#000000';
      ctx_circle.arc(convertX(0), convertY(0), RAD, 0, Math.PI*2.0, true);
      ctx_circle.stroke();
 
      // 円周上の質点の描画部分
      for(var i = 0; i < NUM; i++){     
        //更新した座標で円を描く
        ctx_circle.beginPath();
        ctx_circle.fillStyle = 'rgb(' + r[i] + ',' + g[i] + ',' + b[i] + ')';
        ctx_circle.arc(convertX(Math.floor(RAD*Math.cos(theta[i]))), convertY(Math.floor(RAD*Math.sin(theta[i]))), 10, 0, Math.PI*2.0, true);
        ctx_circle.fill();
      }

      // 重心の描画部分 まず平均を出して、描画する。
      var heikinX = 0.0;
      var heikinY = 0.0;
      for(var i = 0; i < NUM; i++){
        heikinX += RAD*Math.cos(theta[i]);
        heikinY += RAD*Math.sin(theta[i]);
      }
      heikinX/=NUM;
      heikinY/=NUM;
      ctx_circle.beginPath();
      ctx_circle.fillStyle = '#0000FF'; // 色は青固定
      ctx_circle.arc(convertX(Math.floor(heikinX)), convertY(Math.floor(heikinY)), 10, 0, Math.PI*2.0, true);
      ctx_circle.fill();
    }

    // バイトの子の描画
    function draw_kuma(){

      // 画面描画の更新仕様に関するおまじない
      ctx_kuma.globalCompositeOperation = "source-over";
      ctx_kuma.fillStyle = '#FFFFFF';
      ctx_kuma.fillRect(0, 0, WIDTH, HEIGHT);

      for(var i=0;i<6;i++){
        for(var j=0;j<5;j++){
          ctx_kuma.drawImage(img_woman, convertX(-175+60*i), convertY(150-70*j),40,40);
          ctx_kuma.drawImage(img_kuma, convertX(-185+60*i), convertY(170-70*j+15*Math.sin(theta[i+6*j])),50,40);
        }
      }
    }

      // htmlの座標の考え方が画面左上を原点、画面右に向かってx軸、下に向かってy軸になっており、わかりにくい。
      // 画面描画の中心を原点、画面右にx軸、画面上にy軸を取る数学で通常使う座標系から
      // htmlの座標系へ変換する関数
      function convertX(x){
        return x+WIDTH/2;
      }
      function convertY(y){
        return -y+HEIGHT/2;
      }

    </script>

    <style type="text/css">
      canvas { background-color:#FFF; border: 1px solid #999; }
    </style>
  </head>

  <!-- 画面表示部分 -->
  <body onload="init();">
    <!-- 各種ボタン、Kの値を設定するテキストボックス -->
    <input type="button" value ="start" onclick="start()">
    <input type="button" value ="stop" onclick="stop()">
    <input type="button" value ="reset" onclick="init()">
    K=<input type="text" id="K_VALUE" value ="0">

    <!-- javascriptのCanvasという仕様で描画する。
        円周上の運動と、バイトの子の運動の2つ描画領域を作成する。
        ブラウザのウィンドウの大きさによらず横に並べたいのでtable内に作成する。-->
   <table>
      <tr>
        <td>
          <!-- 制御を行うため"circle"というidを付けたCanvas部品 -->
          <canvas id="circle" width="400" height="400">
          </canvas>
        </td>
        <td>
          <!-- 制御を行うため"kuma"というidを付けたCanvas部品 -->
          <canvas id="kuma" width="400" height="400">
          </canvas>
        </td>
      </tr>
    </table>
</html>

なお、作成にあたって下記サイトを参考にしました。
http://yoppa.org/taumedia10/2065.html

いろいろいじって遊んでみてください。

f:id:mattyuu:20170108202522j:plain

*1:実際には2\piの整数倍の違いを無視して考えているので、「絵的に少しだけ前にいる場合」という表現の方が良い気がします。

*2:N個の質点の質量は等しいものとします。

*3:同期具合という言葉はこの記事で勝手に作った言葉です。

*4:日曜数学会ではポケモンのナッシーをモデルにしましたが、ブログでナッシーを使うと面倒なことになりそうなので、フリー素材を使います。

*5:自然振動数の確率分布に関すること。漸近安定の概念。そもそもシミュレーションのNは有限であったが、蔵本予想はN\rightarrow \infty

素数大富豪〜素数大富豪数判定機を作ってみた〜

この記事は、素数大富豪 Advent Calendar 2016 - Adventar 24日目の記事です。23日目はicqk3さんの素数大富豪ver2です。

icqk3さんの素数大富豪 ver.2にて「ラマヌジャン革命」が追加実装されていますが、「ラマヌジャン革命」と共に新たに公式ルールに加わることが決まった「合成数出しにおける指数表記出し」も実装されているようです。仕事が早い!

はじめに

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

nisei.hatenablog.com

私は素数大富豪を1回しかやったことがなく、戦略やゲームの中で出会った素数達について語れるレベルではありません。その観点では本アドベントカレンダー素数大富豪の玄人の方たちが存分に語っていますので、ぜひ読んでください。

私はコンピュータを使って面白いことわからないかなぁという観点で素数大富豪 AdventCalendar 2016に2回分登録させていただきました。本記事はその2つ目の投稿で、1つ目は12月6日に投稿しています。こちらも読んでいただけたら幸いです。

mattyuu.hatenadiary.com

素数大富豪数とは

素数の観点から考えると、1より大きい任意の整数は素数合成数かのいずれかになります。そして素数大富豪の観点から考えると、1より大きい任意の整数は素数大富豪数素数大富豪数のいずれかになります。

この記事では素数大富豪数を下記で定義しています。

素数大富豪数
素数大富豪にて、ペナルティを受けることなく手として出しうる数を「素数大富豪数」と呼ぶ。また、1より大きい整数のうち「素数大富豪数」でない数を「非素数大富豪数」と呼ぶ。

まず、素数大富豪素数*1はいずれも素数大富豪数になります。最大素数大富豪素数より大きい素数は、素数大富豪にて手として出すことができないため、全て非素数大富豪数になります。

素数大富豪には「合成数出し」というルールがあるため、素数大富豪数は合成数も含まれることに注意が必要です。例えば公式ルールで解説されている46,793という合成数73\times 641となるため、73641のカードを捨てることで場に出すことができます。そのため46,793素数大富豪数になります。

しかし、128128=2\times 2\times 2\times 2\times 2\times 2\times 2ですが、ジョーカー2枚を2として使っても7枚の2を場に捨てることができません。そのため128は非素数大富豪数になります(実は最小の非素数大富豪数です。)。

なお、本アドベントカレンダー22日目のせきゅーんさんの下記の記事にて、次回のルール改版時に「合成数出しにおける指数表記出し」ルールが公式ルールとして追加されることが予告されました。

integers.hatenablog.com

せきゅーんさんの記事でも言及されていますが、指数表記出しを行うことで128素数大富豪で出せるようになり、晴れて素数大富豪数に昇格する見通しです。なお、指数表記出しルール追加後の最小の非素数大富豪数は2000になります。突然の指名で2000もさぞかし驚いていることでしょう。今度はどうか2000も救ってあげてください。

素数大富豪数判定機

1より大きい整数を素数の観点、素数大富豪数の観点で分類すると、下記の4分類になります。1より大きい任意の整数は下記の分類のどこかただ1つだけに所属しています。

素数であるか否かは、(簡単かどうかは別にして)素因数分解を行うことで判定できます。しかし、素数大富豪数かどうかの判定は簡単でしょうか。

例えば、
\displaystyle
45,838,577,369,109,312

という数を見てこれが素数大富豪数かどうかを判定するのは難しいと思います。

数字列の最後の12の箇所で1のカードと2のカードを出すべきか、12のカードを出すべきかパッとわかりません。素因数に22,22,219が出てくるのであれば2のカードは温存しておくべきで12を出す必要がありそうです。しかし、素因数に121,212,127が出てくるのであれば逆に12を温存しておかないとカードが足りなくなります。

今回、1より大きい整数を入力すると、

  • その数が4つの分類のうちのどの分類になるのか
  • 素数大富豪数の場合に、どのようにカードを出せばいいのか

を出力する素数大富豪数判定プログラムを作成しました*2

本当はどのようなアルゴリズム素数大富豪数を判定しているのかを紹介したかったのですが、「合成数出しにおける指数表記出し」ルール追加後のプログラム修正時等の別の機会にさせていただきます。今回新ルールでの判定ロジックも作りたかったのですが、12月22日に追加されたルールということで時間が足りませんでした。新ルールでは分割数とか考える必要がありちょっと難しそうです。

この素数大富豪数判定機があれば、身の回りに溢れるあらゆる数たちが素数大富豪で活躍できる数であるかどうかがわかります。今回は素数大富豪数判定機を使って、いろいろな数の判定をしましたので、その紹介を行い締めさせていただきます。以下乱文が続きます。。

明日は本アドベントカレンダーの作成者にせいさんによるまとめの記事です。ありがとうございました。

最大素数大富豪メルセンヌ素数

2^n-1の形をした素数メルセンヌ素数といいます。8番目のメルセンヌ素数(2,147,483,647)までは素数大富豪素数であることがわかりました。そして9番目以降のメルセンヌ素数素数大富豪素数でした。そもそもメルセンヌ素数nに対して爆発的に大きくなっていくので、すぐに最大素数大富豪素数より大きくなってしまいます。

最大の素数大富豪数かつメルセンヌ数(最大素数大富豪メルセンヌ素数)は、8番目のメルセンヌ素数である

\displaystyle
2,147,483,647

で、最小の素数大富豪数メルセンヌ素数は(最小のエデンの園メルセンヌ素数)、は9番目のメルセンヌ素数である

\displaystyle
2,305,843,009,213,693,951

でした。なお、エデンの園素数とは素数大富豪で出すことができない素数です。

integers.hatenablog.com

最小のエデンの園双子素数

双子素数というのは差が2である素数の組のことです(せきゅーんさんによると、「数の組は素数とは言えないため、素数双子と呼ぶべき」とのことです。)。

最小の素数大富豪数の双子素数(最小のエデンの園双子素数)は

\displaystyle
70001


\displaystyle
70003

でした。

最小のエデンの園非正則素数

非正則素数とはなんですか?という方は下記をご参照ください。

integers.hatenablog.com

最小の素数大富豪数の非正則素数(最小のエデンの園非正則素数)は

\displaystyle
70009

でした。

フィボナッチ数列素数大富豪数性

フィボナッチ数列11235813というふうに直前の2項の和を次の項とする数列になります。

最小の素数大富豪フィボナッチ数(場に出すことができないフィボナッチ数)は
\displaystyle
75,025=5\times 5\times 3001
でした。

タクシー数

2番目のタクシー数までは素数大富豪数でしたが、それ以降は素数大富豪数合成数でした。

素数大富豪タクシー数

\displaystyle
\begin{split}
{\rm Ta}(1)&=2\\
\\ \qquad
{\rm Ta}(2)&=1729 = 7\times 13 \times 19\\
\end{split}

(1729合成数出ししなくても「ラマヌジャン革命」の数として使用可能です。)

素数大富豪数かつタクシー数

\displaystyle
\begin{split}
{\rm Ta}(3)&=87,539,319 = 3\times 3\times 3\times 7\times 31\times 67\times 223\\
\\ \qquad
{\rm Ta}(4)&=6,963,472,309,248 = 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 3\times 3\times 3\times 7\times 13\times 19\times 31\times 37\times 127\\
\\ \qquad
{\rm Ta}(5)&=48,988,659,276,962,496=2\times 2\times 2\times 2\times 2\times 2\times 3\times 3\times 3\times 7\times 7\times 7\times 7\times 13\times 19\times 43\times 73\times 97\times 157\\
{\rm Ta}(6)&=24,153,319,581,254,312,065,344 = 2\times 2\times 2\times 2\times 2\times 2 \times  3\times 3\times 3 \times 7\times 7\times 7\times 7 \times  13 \times  19 \times  43 \times  73 \times  79\times 79\times 79 \times  97 \times  157
\end{split}

タクシー数とはなんですか?という方は下記をご参照ください。

integers.hatenablog.com

ちなみにせきゅーんさんのタクシー数は素数大富豪数合成数でした。素数大富豪で使用可能です。さすがですね。

\displaystyle
6,058,655,748=2\times 2\times 3\times 3\times 31\times 151\times 157\times 229


integers.hatenablog.com

最大素数大富豪フェルマー

フェルマー数は2^{2^n}+1の形をした整数です。素数とは限りません。

5番目のフェルマー数までは素数大富豪素数でした。

\displaystyle
\begin{split}
{\rm F}(0) &= 2^1 + 1 = 3\\
\\ \qquad
{\rm F}(1) &= 2^2 + 1 = 5\\
\\ \qquad
{\rm F}(2) &= 2^4 + 1 = 17\\
\\ \qquad
{\rm F}(3) &= 2^8 + 1 = 257\\
\\ \qquad
{\rm F}(4) &= 2^{16} + 1 = 65,537\\
\end{split}

そして6番目のフェルマー数は素数大富豪数
\displaystyle
\begin{split}
{\rm F}(5) &= 2^{32} + 1 = 4,294,967,297= 641\times 6,700,417
\end{split}

しかし、7番目のフェルマー数は素数大富豪数合成数になります。

\displaystyle
\begin{split}
{\rm F}(6) &= 2^{32} + 1 = 18,446,744,073,709,551,617=274,177 \times 67,280,421,310,721
\end{split}

それ以降のフェルマー数で素数大富豪素数になるものはありません。よって最大の素数大富豪フェルマー数は、

\displaystyle
4,294,967,297

で、最小素数大富豪フェルマー数(最小の、素数大富豪数かつフェルマー数)は、

\displaystyle
18,446,744,073,709,551,617

になります。

最小の素数大富豪数かつ完全数


\displaystyle
33,550,336=2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 2\times 8,191

でした。

完全数って何?という方は下記をご参照ください。

integers.hatenablog.com

山の高さ

世界的に有名な山の標高はいずれも素数大富豪数合成数でした。意外に実戦で使えそう?

エベレスト(8,848m)
\displaystyle
8,848= 2\times 2\times 2\times 2\times 7\times 79

K2(8,611m)
\displaystyle
8,611=79\times 109

キリマンジャロ(5,895m)
\displaystyle
5,895=3 \times 3\times 5\times 131

モンブラン(4,809m)
\displaystyle
4,809=3\times 7\times 229

富士山標高(3,776m)
\displaystyle
3,776=2\times 2\times 2\times 2\times 2\times 2\times 59

【その他】
下記のような特徴的な合成数素数大富豪数でした。

\displaystyle
\begin{split}
&1,234,567,89=2\times 3\times 3\times 5\times 3,607\times 3,803 \\
\\ \qquad
&123,456,789=3\times 3\times 3607\times 3,803 \\
\\ \qquad
&12,345,678,910,111,213=113\times 125,693\times 869,211,457 \\
\end{split}

私の誕生日も素数大富豪数でした。。。やはり0がネック。

\displaystyle
19,860,207=3\times 67\times 98,807

*1:素数大富豪」に手として出しうる素数を「素数大富豪素数」と呼ぶ。by tsujimotterさん

*2:いつかjavascript等で誰もが使えるように公開できれば良いのですが、javascriptだけでは巨大な整数の素数判定が難しい気がします

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

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

はじめに

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

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と置いています。出発点となる運動方程式からkmが消えたため、今後kmは用いず、k^{\prime}のみを用いることになります。k^{\prime}が大きいほど、空気抵抗の影響が大きくなります。

v_x(t)v_y(t)の導出

v_xv_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を消去することで、xyの関係を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だけの関数のようにxyを書き下しましたが、この問題では一番遠くまで飛ばす\thetaを求めるという意味で、\thetaも変数と考えた方が都合がよく、x(t,\theta)y(t,\theta)と書いた方がより適切です(hvgk^\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=1t=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)=0g(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}hvgのパラメータをいじることで、関数のグラフは形を変えますが、確かに(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.2h=1v=2g=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組の雀頭(ジャントウ)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}_pp進数体と勘違いしても仕方ないような気がします。。。)

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の元abに対してその積a*bを、

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

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

いかめしく「積を定義する」と書いてしまいましたが、何のことはございません。{\mathbb F}_pを整数の部分集合と考えて通常の積を行ってしまうと、その結果が{\mathbb F}_pからはみ出してしまうので、pで割った余りを積としているだけです。
例えば 、p=5の時、{\mathbb F}_5の元34を通常の積で計算すると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

24を掛けると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=97p=397pを大きくしてくと、次のようにパターンが複雑になっていきます!

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=97p=397pを大きくしてみましょう。

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}_nn素数のときに限り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:数学的補足参照