流れる空の中で数学を。

とある数学好きの「手作りすうがく」と「気ままな雑記」。

【クリスマス】ウラムの螺旋描いてみた

ウラムの螺旋

 ウラムの螺旋とは、自然数の螺旋上に並べ素数だけを塗りつぶしたものである。詳しくは、ウラムの螺旋 - Wikipediaを参照。今回描くウラムの螺旋は、数学ガールの秘密ノート整数で遊ぼうを参考にして、0から始まるものにしてある。

 

クリスマスカラーのウラムの螺旋

 今回は、自分でウラムの螺旋を描いてみた。使った言語は、勉強中のPython。せっかくなので、季節に合わせてクリスマスカラーにしてみた。他の配色のウラムの螺旋が欲しい人がいたら、いってくれたら作ります(rgbを指定してくれると作れます)。

 まずは、100までの素数のウラムの螺旋。素数は赤色、合成数は緑で表示してある。

f:id:FoxQ:20181210125828p:plain

次は、1000までのもの。直線上に素数が並んでるように見える。

f:id:FoxQ:20181210125837p:plain

10000までのもの。いくつかの直線に並んだ模様が見える。

f:id:FoxQ:20181210125943p:plain

最後に、100000までのもの。点が細かすぎてよくわからなくなってきたので、ここまでにする。

f:id:FoxQ:20181210130106p:plain

再配布・ご利用について

 完全にフリー素材としますので、好きなものがあったら、ご自由に使用してくださって構いません。

【近似解の見直し】ぺるせんたげさんの数学コンテスト【問Ⅱ】

ぺるせんたげさんの数学コンテストの問Ⅱの解と近似解

 ぺるせんたげさん(@percentage011)が前回取り上げた問題の答えを公開してくれた。

  そこで、前回作った近似解と比較してみることにした。

sky-time-math.hatenablog.jp

  まず、P(n)で比較してみる。オレンジが前回作った近似解で、青がぺるせんたげさんの厳密解だ。

f:id:FoxQ:20181104123654p:plain

次に、nP(n)で比較すると、値のずれが目立ってくる。

f:id:FoxQ:20181104123753p:plain

 特に、nが大きくなると、近似解の方が明らかに大きな値を持つようになってきてしまっている。
 このずれがどこからきたのかを考えていたのだが、おそらく前回の記事の(4)式の規格化定数を(6)式で計算したのがまずかったようだ。そこで、今回は規格化定数を本来の正規分布のものに戻すと、どれくらい近似が改良されるのかを見てみることにする。

近似の改良

 前回と同様にして計算するが、サイコロの目の平均値Xが従う分布関数f(X;/mu,\sigma^2)

f(X;/mu,\sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp( - \frac{(X-\mu)^2}{2\sigma^2})……(1)

と近似する。これ以降の計算は、前回の記事と同様である。
 まず、規格化定数の値を見てみると、

 \frac{1}{\sqrt{2\pi \sigma^2}} \simeq 0.234……(2)

と前回の1.36より5~6倍ほど小さくなる。*1これを見て、近似は(1)式の正規分布でした方が良いと反省した。改めて、X=1,2,3,4,5,6での値を見積もると、

f(1;\mu,\sigma^2)=f(6;\mu,\sigma^2)\simeq 3.43 \times 10^{-17}……(3)f(1;\mu,\sigma^2)=f(6;\mu,\sigma^2)\simeq 4.66 \times 10^{-7}……(4)f(1;\mu,\sigma^2)=f(6;\mu,\sigma^2)\simeq 0.0543……(5)

  前回の記事と同様に、これらの値から求める確率P(n)は、

P(n)\simeq \frac{0.543}{5n+1}……(6)

となる。これにパラメータa,\kappa,\gammaを用いた因子をかけて補正することを考える。つまり、

P(n) \simeq \frac{0.543 [1+a \exp (-\kappa(n-1)^\gamma ]}{5n+1}……(7)

として、nの小さい値からパラメータを決定して、最終的な近似式とする。

 P(1)=1より、

a \simeq 10.0……(8)

 P(2)=\frac{1}{2}より、

\kappa \simeq 0.0911……(9)

 P(3) =\frac{1}{3}より、

\gamma \simeq 0.460……(10)

 以上より、

P(n) \simeq \frac{0.543 [1+10 \exp (-0.0911(n-1)^{0.46} ]}{5n+1}……(11)

と近似式が求まった。 

改良版近似式と厳密解の比較

 改めて作った式と厳密解とでP(n)を比較してみる。

f:id:FoxQ:20181104132700p:plain

 前回と同様ほとんど違いが見られない。

 次に、nP(n)を見てみる。オレンジが前回の近似式で、緑が今回の近似式、青が厳密解である。

f:id:FoxQ:20181104134155p:plain

 nが大きいところでの値が若干良くなっているくらいで、思ったほど改善されなかった。さらに大きいn=100くらいまでを見てみる。

f:id:FoxQ:20181104135023p:plain

 厳密解との差がさらに大きくなった。nP(n)n\rightarrow \inftyでの収束が近似解の場合かなり遅く、厳密解だとかなり早いことがわかった。これは(7)式の形で近似していることからくる限界だろう。
 そして、厳密解の場合の\lim_{n \rightarrow \infty} nP(n)の値はどうやら0.1よりも小さいようなので、前回の近似で得られる

\lim_{n \rightarrow \infty} nP(n) =0.632……(12)

は明らかに不適だとわかる。今回の近似では、十分大きなnをとる必要があり、収束も遅いものの、

\lim_{n \rightarrow \infty} nP(n) = 0.1086……(13)

なので、まだマシだと言えるだろう。厳密解の値を見ていても、私には本当の\lim_{n \rightarrow \infty} nP(n)の値がわからないので真実は謎のままだ。

*1:この記事では、有効数字3桁で計算していく。

【問Ⅱ】ぺるせんたげさんの数学コンテスト【近似解】

ぺるせんたげさんの数学コンテストの問Ⅱ

 今回はぺるせんたげさん(@percentage011)の数学コンテストの問Ⅱの「近似解」を作ったので、それを紹介したいと思います。問題文は、以下のツイートを参照してください。

  私の実力では近似解*1を出すだけで精一杯だったので、厳密な解き方がわかった方がいたらぜひ教えてください。

n個のサイコロと平均値

 n個のサイコロをふって、出た目の総和がnの倍数になるのは、総和がn,2n,3n,4n,5n,6nの6通りの場合である。よって、総和がそれぞれの値になる確率を見積もって、それらを全て足し合わせれば、それが求める確率P(n)になる。
 サイコロの個数が多い場合を考えて、総和の値を取り扱う代わりに平均値に着目することにする。つまり、サイコロの目の平均値がちょうど1,2,3,4,5,6になる確率をp(1),p(2),p(3),p(4),p(5),p(6)とすると、それらの和が求める確率になる。

P(n)=p(1)+p(2)+\cdots+p(6)……(1)

 1個のサイコロを多数回ふった時の期待値を\mu、分散を\sigma^2として、定義より計算すると、

\mu=3.5……(2)
\sigma^2=\frac{35}{12}……(3)

となる。各々のサイコロは独立同分布に従うとすると、中心極限定理より、サイコロの個数nが十分多いとき、サイコロの目の平均値Xの分布は正規分布で近似できる。その分布関数をf(X;\mu,\sigma^2)とすると、

f(X;\mu,\sigma^2) \simeq A \exp(-\frac{(X-\mu)^2}{2\sigma^2})……(4)

と近似できる。ここで、Aは確率分布関数の規格化定数である。ただし、nが有限である限り、平均値は

1 \le X \le 6……(5)
の範囲になければいけないので、正規分布の定義域を本来のX \in (-\infty,\infty)からX \in (1,6)に制限しなければいけない。この意味で、いくらかの任意性の残る近似をする必要があり、ここでは規格化定数の値を変更することにした。*2すなわち、規格化定数A

\int_1^6 f(X;\mu,\sigma^2) \mathrm{d} X = 1……(6)

によって定める。数値計算でこの値を求めると*3

A \simeq 1.36……(7)

 となる。

平均値がちょうど1,2,3,4,5,6に一致するとは

 n個のサイコロの目の平均値がとりうる値の場合の数は、目の総和の取りうる場合の数に等しい。*4すなわち、

6n-n+1=5n+1……(8)

通りの場合がある。これより、1 \le X \le 6区間5n+1等分して取り扱う。すなわち、ある1つの区間は、平均値がある1つの値になる場合に相当し、各区間の幅の値を\deltaとおくと、

\delta=\frac{6-1}{5n+1}=\frac{5}{5n+1}……(9)

となる。そして、1,2,3,4,5,6の値を含む区間におけるf(X;\mu,\sigma^2)積分値を用いて、p(1),p(2),\cdots,p(6)と見積もることができる。すなわち、k=1,2,3,4,5,6として、その確率は、

p(k) \simeq \int_{k-\delta/2}^{k+\delta/2} f(X;\mu,\sigma^2) \mathrm{d} X……(10)

 nが十分大きいとき、分布関数が連続なため、ある1つの着目した区間内におけるf(X;\mu,\sigma^2)の値は一定値と見なして近似できる。つまり、

f(X;\mu,\sigma^2) \simeq f(k;\mu,\sigma^2) \quad (k-\frac{\delta}{2} \le X \lt k+\frac{\delta}{2} )……(11)

と近似できる。したがって、(10)(11)式より、

p(k) \simeq \delta f(k;\mu,\sigma^2)……(12)

が得られる。これらの値を数値計算で見積もると、

f(1;\mu,\sigma^2)=f(6;\mu,\sigma^2) \simeq 1.99 \times 10^{-16}……(13)
f(2;\mu,\sigma^2)=f(5;\mu,\sigma^2) \simeq 2.71 \times 10^{-6}……(14)
f(3;\mu,\sigma^2)=f(4;\mu,\sigma^2) \simeq 0.316……(15)

と見積もられらた。*5従って、(1)(12)式より、nが十分大きいとき、求める確率P(n)は有効数字3桁で、

P(n) \simeq \frac{3.16}{5n+1}……(16)

に漸近する。このままでいくと、この式はnが十分大きいときに良い近似となるが、nが小さいときは値が大きくずれる可能性がある。そこで、この式にある適当な補正の因子をかけてみることにした。*6具体的には、3つのパラメータa,\kappa(\gt 0),\gammaを導入して、

P(n) \simeq \frac{3.16[ 1 + a \exp (-\kappa (n-1)^\gamma) ]}{5n+1}……(17)

と近似してみた。これらのパラメータをnが小さいときの値、つまり、P(1),P(2),P(3)から決定し、最終的な近似式とする。

nが小さいときの確率P(n)

(i)n=1の場合、明らかにどの目が出ても1の倍数なので、

P(1)=1……(18)

(ii)n=2の場合、具体的に全ての場合を書き出すことで、

P(2)=\frac{1}{2}……(19)

とわかる。

(iii)n=3の場合、サイコロを2つふった時点での目の総和が、①3の倍数になるとき、②3で割ったとき1余るとき、③3で割ったとき2余る確率はそれぞれ等しいことが、具体的に全ての場合を書き下すことでわかる。①の場合、3つ目のサイコロの目が3か6ならば、3つのサイコロの目は3の倍数になる。同様にして、②③の場合も、3つ目のサイコロをふったとき、それぞれ2,5と1,4がでれば、目の総和が3の倍数になる。以上をまとめると、

P(3)=\frac{1}{3} \cdot \frac{1}{3}+\frac{1}{3} \cdot \frac{1}{3}+\frac{1}{3} \cdot \frac{1}{3}=\frac{1}{3}……(20)

であることがわかる。

近似パラメータの決定

 (17)式のパラメータをn=1,2,3の時の値を求めて決定する。
 n=1のとき、(17)(18)式より、数値的に見積もって、

a \simeq 0.899……(21)

となる。次に、n=2のとき、(17)(19)式より、数値的に見積もって、

\kappa \simeq 0.194……(22)

となる。最後に、n=3のとき、(17)(20)式より、数値的に見積もって、

\gamma \simeq 0.465……(23)

となる。以上をまとめて、最終的な近似式を書き出すと、

P(n)\simeq \frac{3.16 [ 1+0.899 \exp( - 0.194 (n-1)^{0.465} ) ] }{5n+1}……(24)

となる。

P(6)の見積もりとP(n)の漸近的なふるまい

 得られた近似式を用いて、試しにn=6の場合を見積もってみると、

P(6) \simeq 0.163……(25)

となり、\frac{1}{6}\simeq 0.167に非常に近いことがわかる。
 これより、P(n)は、nが小さいときは、おおよそ\frac{1}{n}となるが、nが大きくなるにつれて\frac{1}{n}からずれていき、(16)式の値に近づいていくことがわかる。n \rightarrow \inftyでの漸近的なふるまいを評価すると、

P(n) \simeq \frac{0.632}{n}……(26)

となる。言い換えると、次のような1より小さい極限値が存在することが予想される。

 \lim_{n \rightarrow \infty } nP(n) \simeq 0.632 (\lt 1) ……(27)

*1:本来の問題が組み合わせ論の問題であったのに対して、今回の近似的な解き方では統計学的な取扱いになっている。

*2:任意性が残るときに、なるべく単純なものを選ぶというのは1つの良い選択肢だと言えるだろう。

*3:この記事では、数値は有効数字3桁まで計算することにする。

*4:全て1が出た時の総和がnで、全て6が出た時の総和が6nとなる。

*5:正規分布で近似しているので、当然、平均値に近い方が値が大きくなる。

*6:この補正の仕方は、単なる思い付きであり、ある種の力技であるが、全くのでたらめというわけではなく、nが大きくなるとき補正でかけた因子がきちんと1に近づくようにしてある。

【問Ⅰ】ぺるせんたげさんの数学コンテスト

追記(2018/10/31):解法は本質的に変わりませんが、(5)式まわりを微妙に修正しました(lの値のとる範囲等を変更しました)。

ぺるせんたげさんの数学コンテストの問Ⅰ

 今回も前回に引き続き、ぺるせんたげさん(@percentage011)の数学コンテストの問Ⅰを解いていきたいと思います。

もし、間違っていたり、解答が不十分だったら優しく教えてくださると助かります。

一回微分からはじまってわかること

 問題は、次の関数f_n(x)k微分せよというもの。

f_n(x)=e^{nx+e^x}……(1)

 いきなり、k微分するのはきびしそうなので、1回微分してみる。合成関数の微分公式より、

\frac{\mathrm{d}}{\mathrm{d}x} f_n(x) = [ \frac{\mathrm{d}}{\mathrm{d}x} (nx +e^x) ] e^{nx+e^x}
=  (n +e^x) e^{nx+e^x}
= n e^{nx+e^x}+ e^{(n+1)x+e^x}
= n f_n (x) + f_{n+1}(x)……(2)

この式を観察すると、n多項式を係数とするf_j(x)の和の形(jは整数)になっていることがわかる。この構造は、繰り返し微分しても明らかに変わらない。また、k微分したときに出てくるf_j(x)j=n,\cdots,n+kであることもただちにわかる。k微分したときのf_{n+k}(x)の係数が1であることもすぐにわかる。
 したがって、k微分した時の関数形は、ある多項式A_{l,k}(n)を用いて、次のように表せることがわかる。*1

\frac{\mathrm{d}^k}{\mathrm{d}x^k} f_n(x) = \sum_{l=0}^{k} A_{l,k}(n) f_{n+l}……(3)

 以下では、このA_{l,k}(n)を求めていく。これが求まれば、f_n(x)k微分が求まり、問題が解けたことになる。

A_{l,k}(n)の漸化式

 まず、A_{l,k}(n)に関する漸化式を求めるために、(3)式をもう1回微分してみる。

 \frac{\mathrm{d}^{k+1}}{\mathrm{d}x^{k+1}} f_n(x) = \frac{\mathrm{d}}{\mathrm{d}x}\sum_{l=0}^{k} A_{l,k}(n) f_{n+l}
= \sum_{l=0}^{k} A_{l,k}(n) \frac{\mathrm{d}}{\mathrm{d}x} f_{n+l}……(4)

 ここで、(2)式より、1回微分の箇所は計算できて、

 \frac{\mathrm{d}^{k+1}}{\mathrm{d}x^{k+1}} f_n(x) = \sum_{l=0}^{k} A_{l,k}(n) [ (n+l)f_{n+l} +f_{n+l+1} ]
= nA_{0,k}(n) f_n + \sum_{l=1}^{k} [ (n+l) A_{l,k}(n) + A_{l-1,k}(n)] f_{n+l}  +f_{n+k+1} ……(5)

と求まった。上式と、(3)式でk+1と置き換えたものと比較することで、

A_{l,k+1}(n)=(n+l) A_{l,k}(n) +A_{l-1,k}(n)……(6)

を得る。ここで、kk-1に置き換えて、l=1,\cdots,k自然数kに対して、

A_{l,k}(n) =(n+l) A_{l,k-1}(n) + A_{l-1,k-1}(n)……(7)

というA_{l,k}に対する漸化式を得た。

数列A_{l,k}(n)から母関数へ

 A_{l,k}(n)の漸化式(7)を直接解くのは難しそうなので、次のように母関数B_{l,n}(x)を定義する。

B_{l,n}(x)=\sum_{k=0}^{\infty} A_{l,k}(n) x^k……(8)

 漸化式(7)を用いて、母関数を求めていきたい。準備として、(n+l)xB_{l,n}(x)xB_{l-1,n}(x)を計算しておく。

(n+l)xB_{l,n}(x)=\sum_{k=0}^{\infty}(n+l) A_{l,k}(n) x^{k+1}
=\sum_{k=1}^{\infty}(n+l) A_{l,k-1}(n) x^{k}……(9)

xB_{l-1,n}(x)=\sum_{k=0}^{\infty} A_{l-1,k}(n) x^{k+1}
=\sum_{k=1}^{\infty} A_{l-1,k-1}(n) x^{k}……(10)

 次は、(8),(9),(10)式を使って母関数を求めていく。漸化式(7)を用いることで、

B_{l,n}(x)-(n+l)xB_{l,n}(x)-xB_{l-1,n}(x)
=A_{l,0}(n) + \sum_{k=1}^{\infty} [ A_{l,k}(n) - (n+l) A_{l,k-1}(n) + A_{l-1,k-1}(n)]x^k
=A_{l,0}(n)……(11)

ここで、(3)式でk=0とおくことで、すなわち0回微分を考えることで、

A_{l,0}(n)=\delta_{l,0}……(12)

とわかる。ここで、\delta_{l,0}クロネッカーのデルタである。(11)式を整理して、

 [1-(n+l)x ]B_{l,n}(x) = \delta_{l,0} + x B_{l-1,n}(x) ……(13)

 ところで、(2)(3)(5)式から、明らかに

A_{0,k}(n)=n^k……(14)

であるので、B_{0,n}(x)は、

B_{0,n}(x)=\sum_{k=0}^{\infty} A_{0,k}(n)x^k
=\sum_{k=0}^{\infty} (nx)^k
=\frac{1}{1-nx}……(15)

と求まる。(13)式でl \ge 1 の場合を考えると、\delta_{l,0}=0となるので、

B_{l,n}(x) = \frac{x}{1-(n+l)x} B_{l-1,n}(x)
=\frac{x}{1-(n+l)x}\frac{x}{1-(n+l-1)x}B_{l-2,n}(x)
=\cdots
=\frac{x}{1-(n+l)x}\frac{x}{1-(n+l-1)x} \cdots \frac{x}{1-(n+1)x} B_{0,n}(x)……(16)

 よって、(15)式より、母関数B_{l,n}(x)が、

B_{l,n}(x) = x^l \frac{1}{1-(n+l)x}\frac{1}{1-(n+l-1)x} \cdots \frac{1}{1-nx} ……(17)

 と求まった。

母関数B_{l,n}(x)から係数A_{l,k}(n)を取り出す

  後は、母関数を再びxの冪で展開してその係数を調べれば、母関数の定義式(8)よりA_{l,k}が求まる。そのために、任意の自然数aに対して成り立つ次の展開式を用いる。

\frac{1}{1-ax}=\sum_{k=0}^{\infty} a^k x^k……(18)

 この展開を(17)式に適用すると、

B_{l,n}(x) = x^l [ \sum_{{k_l}=0}^\infty (n+l)^{k_l} x^{k_l}  ] \cdots [ \sum_{{k_0}=0}^\infty n^{k_0} x^{k_0}  ]……(19)

 ここで、無限和をとる変数はk_0からk_lまでのl+1個ある。(19)式のx^lを除く因子をかけ合わせて、x^kの係数にまとめ直すと、

B_{l,n}(x) = x^l \sum_{k=0}^{\infty} [ \sum_{k_0 + \cdots + k_l = k} n^{k_0}\cdots(n+l)^{k_l}  ] x^k
=  \sum_{k=l}^{\infty} [ \sum_{k_0 + \cdots + k_l = k-l} n^{k_0}\cdots(n+l)^{k_l}  ] x^k ……(20)

となる。ここで、内側の総和は、k_0からk_lまでの和が条件を満たすような0以上の整数の組み合わせ全てについてとる。この式を母関数の定義式(8)と比較して、

A_{l,k}(n) = \{ 0 \qquad ( l \gt k)
                  \{ \sum_{k_0 + \cdots + k_l = k-l} n^{k_0}\cdots(n+l)^{k_l} \qquad (0 \le l \le k)……(21)

が得られた。上式に現れる総和は、l+1個の整数k_0,\cdots,k_l(\ge 0)に対して、条件を満たすもの全てについてとる。
 以上より、(3)式に戻るとその展開係数が全て求まったことになるので、f_n(x)k微分が求まった。

*1:以降、f_n(x)を単にf_nと変数を省略して書く場合があることを注意しておく。

【問Ⅲ】ぺるせんたげさんの数学コンテスト

ぺるせんたげさんの数学コンテスト

 ぺるせんたげさん(@percentage011)がtwitterで数学コンテストを開催していたので、昨日の夜から問題に挑戦していた。問題の詳細は以下のツイートをチェックしてください。

 

 現在、問Ⅰと問Ⅲは解けたが、問Ⅱは解けていない状態です。
今回の記事では問Ⅲの解説をします。

手計算としらみつぶし

 問題Ⅲに取り組んだ最初の頃は、手計算で解がないこと証明しようとしていた。解の候補をどんどん絞り込んでいって、最後には解の候補が無くなるところまで行こうという戦略をとっていた。ところが、解の候補として20パターンの場合分けが必要になり、さらにその1つ1つについて調べなければいけないというところで、心がポキッと折れた。そこで、数値計算を使ってしらみつぶしで解を探すことにした。つまり、1から9の自然数を並び変える順列の全パターンを、逐一a,b,c,d,e,f,g,h,iに代入して、下の式を満たすかどうかをチェックした。

\frac{a}{bc}+\frac{d}{ef}+\frac{g}{hi}=1……(1)

 結論から先に言うと、(1)式を満たす解が見つかってしまっそこで、実際に解があることを、以下で紹介します。

ただ1つの解

 (1)式では、各項の順序や、分母の因子をかけ合わせる順序を適当に変えることで、一般性を失うことなく、

a \lt d \lt g……(2)
b \lt c……(3)
e \lt f……(4)
h \lt i……(5)

という条件を課すことができる。Pythonでの数値計算により、この条件の下での解はただ1つだけ存在することが判明した。その解は、

a=1,b=3,c=6……(6)
d=5,e=8,f=9……(7)
g=7,h=2,i=4……(8)

である。実際、

\frac{1}{3 \cdot 6} +\frac{5}{8 \cdot 9} + \frac{7}{2 \cdot 4}
=\frac{1}{18} +\frac{5}{72} + \frac{7}{8}
=\frac{4+5+7 \cdot 9}{72}
=\frac{4+5+63}{72}
=\frac{72}{72}=1……(9)

という具合に解になっている。正直この解が見つかったときは、結構驚いたし感動もした。というわけで、この問題を作ってくれた人に「ありがとうございました」を言っておきます。

ランダムな係数の多項式の零点分布

多項式の零点~指数・対数関数からランダム多項式へ~

 以前、指数関数と対数関数を多項式で近似したときの零点の分布を調べたことがある。

sky-time-math.hatenablog.jp

sky-time-math.hatenablog.jp

  そこで使ったプログラムを使って何か面白いことができないかと思っていたので、今回はランダムな係数の多項式の零点を調べてみることにした。以下のサイトに、とあるランダム多項式の零点は円形に分布すると書いてある。

www.singularpoint.org

 さて、本当にこんなに綺麗に円形に分布しているのか、非常に気になったのでやってみることにした。ただし、簡単のため、今回使う乱数は一様乱数にしたし、最高次の係数も1に規格化してあるという違いはある。特別な乱数の分布の場合だけ円形になるのかなどの疑問も残る。さて、どうなるだろうか?

ランダムな係数の多項式~パラメータ等の定義~

 今回は、N多項式で、最高次の係数が1のものを考える。

x^N+a_1 x^{N-1} + \cdots +a_N……(1)

 ここで、a_j (j=1,\cdots,N)はそれぞれ独立な乱数で、 (-W/2,W/2)の間に一様に分布しているものとする。今回は、Wの値としてW=0.02,2,20を取った。乱数を変えたときの零点の分布をまとめてみるために、N次方程式を繰り返し生成し、零点を計算した。N毎に、零点の総数が3000個になるように多項式の個数を用意し、計算を行った。

W=2の場合

 まずは、乱数の値がとる幅をW=2にして、2次方程式から5次方程式の零点の分布を見てみよう。画像中の「Seeds=」は方程式の総数を表していて、「N=」は多項式の次数を表している。

f:id:FoxQ:20181022123050p:plain

f:id:FoxQ:20181022123137p:plain

f:id:FoxQ:20181022123156p:plain

f:id:FoxQ:20181022123220p:plain

 実軸上を走る零点はどの次数の場合にもあるが、その他の零点は次数が上がるにつれて、円形に近づいて行っているように見える。そこで、次数をN=10,20,30まで上げてみた。

f:id:FoxQ:20181022123425p:plain

f:id:FoxQ:20181022123440p:plain

f:id:FoxQ:20181022123458p:plain

 以上のように多項式の次数を上げていくと、円形に分布する零点が現れることがわかった。

W=0.02の場合、

 次にやったのは、乱数の幅を小さくとることだった。W=0.2の場合もやったのだが、W=0.02の場合の方が零点の分布がより特徴的で最も面白かったので、そちらを載せることにする。W=2の場合との大きな違いは、多項式の次数が小さい場合に顕著に表れる。

 それでは、次数がN=2,3,4,5の場合を見ていこう。

f:id:FoxQ:20181022123909p:plain

f:id:FoxQ:20181022123922p:plain

f:id:FoxQ:20181022123932p:plain

f:id:FoxQ:20181022123941p:plain

 明らかに、W=2の場合とは異なる零点の分布が見て取れる。N多項式の場合、2N-2個の塊が零点の分布として現れるのだ。そして、中央の「空白」は円形になっているように見える。なんとも不思議である。このような分布が現れる理由を説明できる方がいたらぜひ教えてください。
 N=10,20,30の場合も見ていく。すぐ下にあるように、N=10の場合は、なんとか、上下に9個ずつの塊があるのが見える。

f:id:FoxQ:20181022124137p:plain

f:id:FoxQ:20181022124146p:plain

f:id:FoxQ:20181022124157p:plain

 そして、次数を増やしていくと、零点の分布はやはり円形に近づいていく。

W=20の場合

 最後に、おまけとして、W=20の場合もやってみた。N=2次から5次までをまずはみてみたが、この場合は最初から円形に近い分布が見て取れた。(縦横の比が違うので、楕円っぽくみえますが値はおおよそ円です。)

f:id:FoxQ:20181022124514p:plain

f:id:FoxQ:20181022124525p:plain

f:id:FoxQ:20181022124534p:plain

f:id:FoxQ:20181022124543p:plain

 最後に、次数をさらに大きくした場合を見て終わりにする。

f:id:FoxQ:20181022124738p:plain

f:id:FoxQ:20181022124748p:plain

f:id:FoxQ:20181022124801p:plain

次数を上げると、だんだん零点の分布が原点に絞り込まれていく傾向があることがわかった。また、ここまでを通して、Wが小さいときに見られた円の内部の空白が、W=20では埋まってしまうこともわかった。

【修正版】3が合同数でないことの初等的証明

追記(2018/10/13):(10)式で計算ミスをしており、a^2ではなく、a^4でした。そのため、原始ピタゴラス数の3つ組が(a^2,6b,c)となり、不等式の評価で矛盾を導けなくなってしまいました。いつかうまくいったら、更新します。

前回の失敗から

前回の記事では、3が合同数でないことを初等的に証明しようとしたが、証明に飛躍があり見事に失敗した。(自然数nが合同数であるとは、全ての辺の長さが有理数で、面積がnの直角三角形が存在することと定義されている。)
sky-time-math.hatenablog.jp

 今回の記事では、証明の修正に成功できたと思うので、3が合同数でないことを改めて証明することに挑戦する。

背理法からピタゴラス数へ

 前回と同様に、背理法を用いる。復習のため、あえてもう一度書く。
 3が合同数であると仮定する。すなわち、ある有理数x,y,zが存在して、

x^2+y^2=z^2……(1)
xy=6……(2)

と仮定する。*1上式をx^2で割って、

1+(\frac{y}{x})^2=(\frac{z}{x})^2……(3)
\frac{y}{x}=\frac{6}{x^2}……(4)

(4)式を(3)式に代入して、

1+\frac{6^2}{x^4}=\frac{z^2}{x^2}
\Rightarrow x^2 z^2 = x^4+6^2……(5)

ここで、有理数x,z

x=\frac{a}{b}……(7)
z=\frac{c}{d}……(8)

と表す。ただし、abcdは互いに素であるとする*2
 (5)式に、b^4 d^2をかけて、

(abc)^2 = d^2 (a^4+6^2 b^4)……(9)

を得る。両辺をa^2 b^2で割ると、

c^2 = (\frac{d}{ab})^2 ( a^2+(6b)^2)……(10)

を得る。ここで、cdが互いに素であるためには、abdの倍数でなければいけない。すなわち、ある自然数kが存在して、

dk=ab……(11)

と表せる。(10)式の両辺にk^2をかけて整理すると、

(ck)^2=a^2+6b^2……(12)

を得る。
 ところで、(12)式より、abが互いに素であるためには、kakbがそれぞれ互いに素でなければいけない。このことと、(11)式より、dabの倍数であることがわかるので、結局、

k=1……(14)
d=ab……(15)
c^2=a^2+(6b)^2……(16)

 が言える。まずは、このように合同数の問題が、ピタゴラス数の問題へと書き換わる。

ピタゴラス数から原始ピタゴラス数へ

 次のステップに進むために、(16)式に現れるa,6b,cがそれぞれ互いに素であること、つまり原始ピタゴラス数になっていることを示す。caと互いに素でないと仮定すると、(16)式より明らかに、abが互いに素でなくなり、矛盾する。よって、

caは互いに素である。

同様に、

cbも互いに素である。

 次に、(16)式に現れる数『6』の取扱いがポイントとなる。
 そこで、\gcd (a,6) \neq 1と仮定してみると、aは2,3,6のいずれかの倍数となる。これに対応して、(16)式より、cも2,3,6のいずれかの倍数になる。ところが、(15)式より、daの倍数なので、dも2,3,6のいずれかの対応する倍数になる。このことは、cdが互いに素であったことに矛盾する。従って、

 \gcd (a,6) = 1……(17)

が言えた。さらに、\gcd (c,6) \neq 1と仮定すると、(16)式より、\gcd(a,6) \neq 1となり矛盾するので、

 \gcd (c,6) =1……(18)

も従う。以上より、(16)式に現れる

a,6b,cは原始ピタゴラス数である

ことがわかった。

原子ピタゴラス数の表し方

 任意の原始ピタゴラス数を表す公式があり、(16)式は、互いに素な奇数s,t (s \gt t \ge 1)を用いて、

a=st……(19)
6b=\frac{s^2-t^2}{2}……(20)
c=\frac{s^2+t^2}{2}……(21)

と表せる。この公式の導出は、例えば、『はじめての数論(原著第3版)』の定理2.1等に書いてある。また、英語版のwikipediaPythagorean triple - Wikipediaの2.3にも書いてある。(ちなみに、上記wikiの2の冒頭の公式から(19)~(21)式を導くには、s=m+n,t=m-nとおけばいい。)

はじめての数論 原著第3版 発見と証明の大航海‐ピタゴラスの定理から楕円曲線まで

はじめての数論 原著第3版 発見と証明の大航海‐ピタゴラスの定理から楕円曲線まで

 

  ここで、原子ピタゴラス数の表式から次のステップに進む前に、ある不等式の評価をしておく。

合同数が3の場合に成り立つ不等式

 (10)式に相加・相乗平均の式を用いて、次のように不等式を得る。

c^2=d^2 ( (\frac{a}{b})^2 + (6 \frac{b}{a})^2 )
\ge 2 d^2 \sqrt { (\frac{a}{b})^2  (6 \frac{b}{a})^2}
\ge 12d^2……(22)

 c,dは正の数なので、目的の不等式

c \ge 2 \sqrt {3} d \gt 3 d……(23)

を得る。これが、今回の証明の鍵となる。

証明の完了へ

 (20)式より、

b=\frac{s^2-t^2}{12}……(24)

また、(15)(19)式より、

d=ab=\frac{st(s^2-t^2)}{12}
=\frac{1}{12} st (s+t)(s-t)……(25)

ところで、(17)式に書いたようにaと6は互いに素なので、a=stと表されていることを考慮すると、

s,t \equiv \pm 1 (\mod 6)……(26)

と表されることがわかる。このことから、s-tが最小になるのは、ある整数lによって、s,t

s=6l+1……(27)
t=6l-1……(28)

と表される場合であることが言える。すなわち、

s-t \ge 2……(29)

である。これを(25)式にあてはめると、

d \ge \frac{1}{6} st (s+t)
 \iff 3d \ge \frac{s^2 t+t^2 s}{2}……(30)

を得る。ところが、(26)式を満たす奇数s,tは、

s \ge 5……(31)
t \ge 1……(32)

以上となる。したがって、(30)式は、

3d \ge \frac{s^2 + 5t^2}{2}……(31)

となる。ここで、(21)式より

c=\frac{s^2+t^2}{2}……(21)

 であったので、明らかに、

3d \gt c……(32)

が導かれる。ところが、(23)式より、c \gt 3dであったので、これは矛盾である。
 以上より、一番最初の背理法の仮定が偽であることが言える。すなわち、3は合同数ではないことが証明できた。

*1:zが直角三角形の斜辺の長さになる。

*2:つまり、x,zが既約分数とする。