三次方程式の解（２）

 三次方程式の解（１）からの続き

大学入試問題で扱われるような三次方程式$\ x^3 + p x + q =0\ $の解についての別のアプローチがあった。まず，$\ x=y+\frac{a}{y}\ $とおいて， $ y $の6次方程式に変換すると，$\ y^3+3a (y+\frac{a}{y}) +\frac{a^3}{y^3} + p (y+\frac{a}{y}) + q = 0\ $が得られる。これが$\ t=y^3\ $の2次方程式になるように，$\ 3a+p=0\ $という条件をつければ，$ a=-\frac{0}{3} $となる。このとき，$\ t^2 + q t -\frac{p^3}{27}=0\ $が成り立つ。

この2次方程式の解は，$\ t_1=-\frac{q}{2} + \sqrt{q^2/4 + p^3/27}, \ t_2=-\frac{q}{2} - \sqrt{q^2/4 + p^3/27} \ $となる。また，6次方程式の解は，$y_1^3=t_1, \ y_2^3=t_2$を満足している。ここで，$x^3=1$の解を，$ \{ 1, \ \omega=\frac{-1+\sqrt{3} i}{2}, \ \omega^2=\frac{-1-\sqrt{3} i}{2} \} $とする。

(1) $\ t_1,\ t_2$が実数の場合：

$\alpha=t_1^{1/3},\ \beta=t_2^{1/3}$を実数とすると，$y_1 = \{ \alpha,\ \alpha \omega,\ \alpha \omega^2 \},\quad  y_2 = \{ \beta,\ \beta \omega,\ \beta \omega^2 \}\ $となる。この6つの解からそれぞれ，$\ x=y -\frac{p}{3 y}\ $によって，もとの三次方程式の解がえられるので，解の組が過剰に存在するようにみえるが，$y_2$のセットは，$y_1$のセットと同じになることが確かめられる。

$ \alpha^3 \beta^3 = t_1 t_2 =  -\frac{p^3}{27}$で，$\alpha, \beta$がともに実数であることから，$\alpha \beta = -\frac{p}{3}$である。

これから，$x_2 = y_2 - \frac{p}{3 y_2} = \{ \beta-\frac{p}{3\beta},\ \beta \omega - \frac{p}{3 \beta \omega},\ \beta \omega^2 - \frac{p}{3 \beta \omega^2} \} $

$ = \{ -\frac{p}{3\alpha} + \alpha,\ -\frac{p}{3 \alpha \omega}+ \alpha \omega ,\ - \frac{p}{3 \alpha  \omega^2} +\alpha \omega^2  \} =  y_1 - \frac{p}{3 y_1} = x_1$

つまり，$\beta $からくる解はすべて$\alpha$からくる解の三つ組に帰着する。

三次方程式の解（３）に続く。

平方完成

 平方完成は，入試問題を解くときなど，条件設定の場面でたいへん重宝する技法だ。ちょっと手計算が面倒な式がでてきたので，Mathematicaに任せようと思った。

ところが，探してみてもMathematicaで平方完成する関数が組み込まれていないようなのだ。もしかしたら調べ方が足りないのかもしれないが，普通に考えるとイの一番に出てきても良さそうな機能なのだが。

それらしいユーザ定義関数がいくつか見つかったけれど，2変数の整式を代入しても思ったような変形ができず，望みのものではなかった。しかたがないので，自分で関数パーツを考えることにした。これを一般化するには，Mathematicaプログラミングにおける文法の知識が足りなさすぎる。

ここで考えたのは，ある変数の二次式を与えたときに，平方完成された部分と残余部分のリストを返すユーザ定義関数 sq[式, 変数]だ。変数がn個ある場合は，n回繰り返して使う必要があるという残念なコード素片だ。

sq[f_, v_] :=
　Module[{a, b},
　　a = Coefficient[f, v^2];
　　b = Coefficient[f, v];
　　{a (v + b /(2 a))^2,
　　c = f - a (v + b /(2 a))^2}] // Simplify

これを使って次のような計算ができる。

In[1]:= sq[2 x^2 - 4 x y + 2 y^2 + 24 x - 24 y + 288 + 3 x y, x] 

Out[1]= {1/8 (-24 - 4 x + y)^2, 216 - 18 y + (15 y^2)/8}

In[2]:= sq[%[[2]], y] 

Out[2]= {3/40 (24 - 5 y)^2, 864/5}

三次方程式の解

 二次方程式の解は，与えられた２次式を平方完成すればよいので，公式を忘れても導ける。まあ，平方完成の手順を理解して導出できるくらいなら，かつて曾野綾子に「二次方程式の解の公式を学んだことは，人生において何の役にもたたなかった」とボロクソに腐された解の公式もすぐに出てくるだろうから心配する必要はない。

準備として，$x^3=1$の解を，$\{ 1,\ \omega=\frac{-1+\sqrt{3} i}{2}, \ \omega^2=\frac{-1-\sqrt{3} i}{2} \}$としておく。

三次方程式$a x^3 + b x^2+ c x + d = 0$は，$x^3 +p x + q =0$ の形にすることができる。次に，因数分解の公式，$x^3+y^3+z^3-3xyz = (x+y+z)(x^2+y^2+z^2-x y -y z -z x)$を用いる。つまり，$p = - 3 y z$，$q = y^3 + z^3$とすれば，もとの三次方程式は因数分解できることになり，すなわち，解が求まることになる。

ここで，$y^3$と$z^3$の対称式を考えるのがポイントである。$p^3=-27 y^3 z^3$から，$y^3$と$z^3$は，$t^2-q t -(p/3)^3=0$の解である。$t = (q/2) \pm \sqrt{(q/2)^2+(p/3)^3}$

因数分解された右辺の第2項を$x$の2次式と考えてさらに因数分解するため，$x^2-(y+z)x +y^2 -yz + z^2 = 0$とおいて，2次方程式の解の公式を使うと，

$x=\frac{1}{2} \bigl( y + z \pm \sqrt{(y+z)^2-4(y^2+z^2-yz)} \bigr) = \frac{1}{2} \bigl( y + z \pm \sqrt{-3y^2+6yz-3z^2} \bigr)$

$\quad= \frac{1}{2} \bigl( y + z \pm (y - z ) \sqrt{3}i  \bigr) = -y \frac{-1 \mp \sqrt{3}i}{2} -z \frac{-1 \pm \sqrt{3}i}{2}$

したがって，$x^3 +p x + q =0$の解は，$\{ -y -z, \ -\omega^2 y -\omega z, \ -\omega y - \omega^2 z \}$，ただし，$\{ y , z \}= \{ \bigl( q/2 + \sqrt{(q/2)^2+(p/3)^3} \bigr)^{1/3}, \ \bigl ( q/2 - \sqrt{(q/2)^2+(p/3)^3} \bigr)^{1/3} \} $である。

二項分布と正規分布

統計物理学のための準備シリーズが続く。ここでは，$N$が大きいときの二項分布を正規分布で近似する方法を確かめる。

アボガドロ数$N$個の粒子を，左右2つの箱に確率$p$と$q$（$p+q=1$）で入れるとき，分配される粒子の個数の確率分布は二項分布に従う。すなわち，左の箱に入る粒子の数を$n$，その場合の確率を$r(n)$とすると，$r(n)={}_N C_{N-n} p^n q^{N-n}=\frac{N!}{n!(N-n)!} p^n q^{N-n}, \quad \sum_{n=0}^N {}_N C_{N-n} p^n q^{N-n} =(p+q)^N = 1$

ここに，スターリングの公式，$n! \simeq \sqrt{2\pi n} (\frac{n}{e})^n $ 等を当てはめると，

$r(n) \simeq \sqrt{\frac{N}{2\pi n(N-n)}} \frac{N^N p^n q^{N-n}}{n^n (N-n)^{N-n}} =  \sqrt{\frac{N}{2\pi n(N-n)}} \bigl( \frac{Np}{n}\bigr)^n \bigl(\frac{Nq}{N-n} \bigr)^{N-n}$

$\therefore \log r(n) \simeq n \log \frac{Np}{n} + (N-n) \log \frac{Nq}{N-n} $

 ただし，$O(\{n,N\}^{-1/2})$である初項はおとす。極値を求めるため，$\log r(n)$を$n$で微分して，

$\log Np -1 -\log n -\log Nq +\log(N-n) +1 =0, \quad \log \frac{Np}{n} = \log \frac{Nq}{N-n}$

極値を与えるのは$n=Np$であり，このとき$r(n)=\sqrt{\frac{1}{2\pi N p q}}$ となる。

次に，$n=Np+x$とおき，$r(n)$を$n=Np$のまわりに展開して$x$の2次近似式を求める。ただし，$x \ll Np$であり，$\log (1\pm x) \simeq \pm x + \frac{x^2}{2}$を用いる。

$r(n )= \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ -n \log \frac{n}{Np} - (N-n) \log \frac{N-n}{Nq} \}$

$\quad\quad = \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ -(Np+x) \log (1+ \frac{x}{Np}) - (Nq-x) \log (1-\frac{x}{Nq}) \}$

$\quad\quad = \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ -(Np+x)  ( \frac{x}{Np}+ \frac{x^2}{2 (Np)^2} ) - (Nq-x)  (-\frac{x}{Nq} + \frac{x^2}{2 (Nq)^2} ) \}$

$\quad\quad = \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ -(x + \frac{x^2}{2 Np}) - (-x + \frac{x^2}{2 Nq}) \}$

$\quad\quad = \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ - \frac{x^2}{2 p q N} \} = \sqrt{\frac{1}{2\pi p q N}} \ \exp \{ - \bigl(\frac{n-Np}{\sqrt{2 p q N }}\bigr)^2 \} $

このとき，次の規格化条件が満たされる。$\sigma = p q N$とおいて，$\int_{-\infty}^{\infty} \sqrt{\frac{1}{2 \pi \sigma}} \ \exp \{ - \bigl(\frac{n-Np}{\sqrt{2 \sigma }}\bigr)^2 \} dn = 1$

［１］De Moivre - Laplace Theorem

スターリングの公式

 d次元球の体積からの続き

統計力学をはじめるには，スターリングの公式が必須なので菊池誠のテキストで復習する。

ガンマ関数$\Gamma(x)$は，階乗の実数への拡張と考えることができる関数である。$x$が自然数$n$のときに$\Gamma(n+1)=n!$を満足し，次の積分で定義されている。

$\Gamma(x)=\int_0^\infty t^{x-1} e^{-t} dt ,\quad \Gamma(1)=1,\quad  \Gamma(1/2)=\int_{-\infty}^\infty e^{-x^2} dx = \sqrt{\pi}$

ここで，階乗の雰囲気をただよわせる漸化式が成り立つ。$\Gamma(x+1) = \int_0^\infty t^x e^{-x} dt = \bigl[ (-1) t^x e^{-x}\bigr]_0^\infty + \int_0^\infty x t^{x-1} e^{-x} dt = x \Gamma(x)$

引数$x$に，後にアボガドロ数オーダーの非常に大きな数になる予定の$N$を入れてみる。

$N!= \Gamma(N+1) = \int_0^\infty t^{N} e^{-t} dt = \int_0^\infty t^{N} e^{N\log t -t} dt $

指数関数の肩の項$N\log t -t$を微分して極値を探すと，$t=N$で極大値が$N\log N -N$となる。そこで，この$t=N$の周りで$N\log t-t$をテイラー展開して2次の項まで残す。

$N\log t-t \simeq N\log N - N -\frac{1}{2N}(t-N)^2$

指数関数の肩にある$t$の2次の項をほとんど寄与がない積分領域である$\int_{-\infty}^0$まで拡張した上でガウス積分すると，$N!= \int_0^\infty t^{N} e^{N\log N -N -\frac{1}{2N}(t-N)^2} dt = N^N e^{-N} \int_{-\infty}^\infty e^{-\frac{1}{2N}(t-N)^2} = \sqrt{2\pi N} \bigl( \dfrac{N}{e}\bigr)^N$

これから，スターリングの公式，$\log N! = N \log N -N + \log  \sqrt{2\pi N}$ が得られる。

d次元球の体積

 来年度の統計物理学の準備をしている。

先日購入したばかりだが，阪大の湯川諭さんの統計力学の教科書はなかなか読みやすいと思ったが，その師匠の菊池誠「統計力学のはじめの一歩」も結構よかった。これまでテキストに指定されてきた，講談社基礎物理学シリーズの北原・杉山の統計力学は取っつきにくかったが改めて見直すと丁寧でやさしそうだ。高橋康さんの統計力学入門は好感が持てるけれど・・・，あれこれ目移りするのでなかなか講義ノートがまとまらない。

とりあえず，ミクロカノニカル集団で出会う最初の関門が，半径$R$の$d$次元球の体積${\cal V}_d(R)$の導出だ。これは次のように$d-1$次元球面を動径方向に積分して求めることができる。ただし，半径$R$の$d-1$次元球面積は，その次元を考慮して，${\cal S}_d(R) = c_d \ R^{d-1}$で与えられるとする。$c_d$は幾何学的な係数である。

${\cal V}_d( R) = \int_0^R {\cal S}_d(r) dr =\int_0^R c_d\ r^{d-1} dr = c_d \dfrac{R^d}{d}$

したがって，幾何学的な係数である$c_d$が求まれば良いわけだ。これを計算するために，$r^2= x_1^2 + \cdots + x_d^2$として，次のガウス積分を利用する。

$\int_{-\infty}^{\infty} dx_1 \cdots \int_{-\infty}^{\infty}dx_d \ e^{- (x_1^2 + \cdots + x_d^2)} =    \int_0^\infty c_d\ r^{d-1}\ e^{- r^2} dr$

$l. h. s. = \Bigl( \int_{-\infty}^{\infty} e^{- x^2} dx \Bigr)^d = \sqrt{\pi}^{\ d}= \pi ^{d/2} $

$r. h. s. = c_d \int_0^\infty r^{d-1} e^{- r^2} dr = c_d \int_0^\infty \frac{1}{2} t^{d/2 -1} e^{-t} dt = \dfrac{c_d}{2} \Gamma(\frac{d}{2})$

$\therefore c_d = 2 \pi ^{d/2} / \Gamma(\frac{d}{2}),\quad  \therefore \ {\cal V}_d(R) = \dfrac{\pi^{d/2} R^d}{\Gamma(\frac{d}{2}+1)}$

［１］超球の体積

乗法的積分

新型コロナウィルス感染症のオミクロン株の感染者数を見ている。かつては，実行再生産数が話題になることが多かった。いまは，前週との感染者数の比較によって増加や減少の傾向が議論されている。その場合でも，指数関数的（等比級数的）な振る舞いかどうかをみるには，曜日による変動を除くために，前週の感染者数との比率が重要になる。

例えば，東京と大阪の各曜日の前週に対する感染者数比率の相乗平均を週ごとに計算してプロットすると次のようになる。増加率のピークで6倍になったのは1月の初旬である。この比率は1.0に近づいているので，感染者数のピークアウトがそろそろ視野に入ってきた。

図：東京・大阪の感染者数前週比の推移

この比率関数のモデル化が容易ならば，逆算して一日当たり感染者数関数を求めることができる。もっとも，大阪維新によって破壊されている行政統計は出鱈目なので，そもそも意味があるかという問題は残るのだが。

さて，増加比率$r(t)$からもとの関数$f(t)$を求める問題を$f(t+h)/f(t) =r(t)^h$と定式化した。$\lim h \rightarrow 0$で両辺は1になる。この両辺の対数をとると，$\log f(t+h) - \log f(t) = h \log r(t)$ であるから，$\lim_{h \to 0} \dfrac{\log f(t+h) - \log f(t) }{h} = \dfrac{d}{dt} \log f(t) = \log r(t) $ が得られる。

これから，$\log f(t) = \int_c^t \log r(t') dt'$ となり，$f(t) = e^{\int_c^t \log r(t') dt'}$ が得られる。適当な関数を使って試してみたが，比率のピーク関数を時間方向に対称にとれば，感染者数の関数は対称ではなく，後方に広がるといった程度の定性的なイメージしか得られなかった。

ところが，こうした関数の比率についての解析的な議論は，乗法的積分のうちの幾何積分としてすでに知られていることがわかった。数列の和から普通の積分が導かれるように，数列の積から乗法的積分が得られるのだった。なるほど，数学の奥は深い。

大学入学共通テスト（２）

 大学入学共通テスト（１）からの続き

数学ＩＡの問題について批判ばかりしていてもしょうがないので，時間を決めて解こうとした。問題と解答は公表されており，次のようなものだ。

70分で4問（必答2問，選択2問）
第１問　必答：［１］対称式，［２］三角関数，［３］円と三角形
第２問　必答：［１］2次方程式の解と集合・論理，［２］統計処理
第３問　選択（2/3）［１］場合の数と確率
第４問　選択（2/3）［１］不定方程式と剰余
第５問　選択（2/3）［１］三角形の重心

スタート早々につまづいた。第１問の［１］の問題の条件が(1)だけでなく，(2)に及ぶことが読み取れずに全く進まない。なんとかクリアしたら［３］の問題が解けない。うーん，円周角の定理と三角関数はどういう関係になっていたっけ・・・，とりあえずおいて第２問に進むと，簡単な2次方程式の話なのだが頭に入ってこない。

というわけで30分経っても1/3もできずに脱落してしまった。チーン。不定方程式はなおさらだし，統計や幾何の問題は頭に入ってこないのであった。時間制限にあわてると平均点もとれないかもしれない。むしろこれで40点も取れるほうがすごい。

共通一次試験やセンター試験のころは，監督の時間に問題をパラパラと見て，これなら解けそうだと思ったものだ。いまでは新聞に掲載されている問題の字が小さすぎて読むことすらできない。ネットに掲載されている問題ならば字が大きいので大丈夫かと思ったら，共通テストになって読解力重視になったせいで，文章が多すぎてこれまた最後まで読み通せないのだった。

大学入学共通テスト（１）

 今年の大学入学共通テストはたいへんだった。オミクロン株の感染者数が急上昇するなか，1月15日（土）には東大の試験会場前で傷害事件があり，1月16日（日）にはトンガの海底火山の爆発にともなう津波警報で試験ができなかった会場がでた。

大学入試センター試験が，いろいろあって，大学入学共通テストに模様替えしたときに「思考力，判断力，表現力等を発揮して解くことが求められる問題を重視した」結果，今年の問題は読解力を必要とする部分が増えて，難易度があがったらしい。

そんなわけで，数学の試験を眺めてみたら，あら，防衛省の不祥事をヒントとしたようにみえる問題があるではないか。2019年にイージス・アショア配備の調査にかかわって露呈した事件だ。配備地の設定根拠を説明する資料において，Google Earthを使って水平方向と鉛直方向の縮尺比を考慮せずに山の俯角を導いて，判断のための重要なデータにしたというお粗末な（悪質な）ものだった。

今回の問題では，水平10万分の1と鉛直2万5千分の1の縮尺比を考慮せずに図を読み取った結果，キャンプ地から山頂を見込む俯角が16度になったのを検証するというものだ。三角関数表からこの16度に対応する正接を読み取り，これを1/4した値を解答し，さらに，この正接値から逆引して真の俯角（4度〜5度）を導くというものだ。

はい，いい問題ですね。しかし，平均点は数学I・数学Aで20点，数学II・数学Bで18点も下がった。

図：朝日新聞2019年6月記事と大学入学共通テスト数学I・A問題の比較

P. S. Twitterの感想によれば，第2問のデータ処理，第4問の整数問題はかなり評判悪い。やはり，その出発点が利権がらみ（後に消失か）だった共通テストは失敗で，センター試験のほうが良質だったようだ。

四項間漸化式

 三項間漸化式からの続き

前回の問題を次のように修正する。2次方程式$\ x^2-p x + q = 0\ $の解を$\alpha, \beta$として，$a_n = \alpha^n +\beta^n$と定義する。ただし，$a_0=2,\  a_1 =\alpha + \beta = p,\  a_{-1}=\frac{1}{\alpha}+\frac{1}{\beta} = \frac{\alpha + \beta}{\alpha \beta} = \frac{p}{q}$である。

このとき，次の漸化式，$a_{n+1}=p\ a_{n}-q\ a_{n-1}, \ (n=0,1,2 \cdots)$が成り立つ。

また，数列の一般項は $a_n = \alpha^n + \beta^n = \Bigl( \dfrac{p + \sqrt{p^2-4q}}{2} \Bigr)^n +  \Bigl(\dfrac{p - \sqrt{p^2-4q}}{2} \Bigr)^n$ で与えられる。

これを3次方程式に拡張すると次のようになる。3次方程式$\ x^3-p x^2 + q x -r = 0\ $の解を$\alpha, \beta, \gamma$として，$b_n = \alpha^n + \beta^n + \gamma^n$と定義する。ただし，$b_0=3,\  b_1 =\alpha + \beta + \gamma = p,\  b_{-1}=\frac{1}{\alpha}+\frac{1}{\beta} +\frac{1}{\gamma} = \frac{\alpha \beta + \beta \gamma + \gamma \alpha}{\alpha \beta \gamma} = \frac{q}{r}$である。

このとき，漸化式は，$b_{n+2}=p\ b_{n+1}-q\ b_{n}+r\ b_{n-1}, \ (n=0,1,2 \cdots)$となる。

さらに，数列の一般項は $b_n = \alpha^n + \beta^n + \gamma^n$なので，基本対称式の値，$p=\alpha+\beta+\gamma,\ q=\alpha \beta + \beta \gamma + \gamma \alpha,\ r=\alpha \beta \gamma$の多項式になる。2次方程式の場合のように，3次方程式の一般解を代入してもよいが，非常に煩雑な表現になってしまい，Mathematicaを持ってしてもExpandとSimplifyを組み合わせて n=8で行き詰まってしまった。もちろん，漸化式を使えば大丈夫なのだが。

a = x /. Solve[x^3 - p x^2 + q x - r == 0, x];

b[n_] := a[[1]]^n + a[[2]]^n + a[[3]]^n

Table[Expand[b[i]] // Simplify, {i, 1, 8}]

{p, p^2 - 2 q, p^3 - 3 p q + 3 r, p^4 - 4 p^2 q + 2 q^2 + 4 p r, p^5 - 5 p^3 q + 5 p q^2 + 5 p^2 r - 5 q r, p^6 - 6 p^4 q + 9 p^2 q^2 - 2 q^3 + 6 p^3 r - 12 p q r + 3 r^2, p^7 - 7 p^5 q + 14 p^3 q^2 - 7 p q^3 + 7 p^4 r - 21 p^2 q r + 7 q^2 r + 7 p r^2, p^8 - 8 p^6 q + 20 p^4 q^2 + 8 p^5 r - 32 p^3 q r + 24 p q^2 r + 2 q (q^3 - 4 r^2) - 4 p^2 (4 q^3 - 3 r^2)}

三項間漸化式

鈴木貫太郎のチャンネルで， $p = a+\dfrac{1}{a}\ $が与えられたときに，$a^6+\dfrac{1}{a^6}\ $を求める問題があった。これを一般化して，$b_n \equiv a^n+\dfrac{1}{a^n}\ $を求める問題を考えた。

ここで，$b_{n+1}=p\ b_n- b_{n-1}$ が成り立つ。これは三項間漸化式なので，このあたりに一般解法がある。漸化式の特性方程式 $x^2 - p\ x + 1 = 0$の解を$\alpha=\frac{p + \sqrt{p^2-4}}{2},  \beta=\frac{p - \sqrt{p^2-4}}{2}$とする。このとき，$b_n = A \alpha^n + B \beta^n \ $で一般解が与えられることにより，$b_0=2, b_1=p\ $の初期条件から$A, B$を決定すると $A=B=1$となる。

つまり，$b_n = \Bigl( \frac{p + \sqrt{p^2-4}}{2} \Bigr)^n +  \Bigl(  \frac{p + \sqrt{p^2-4}}{2} \Bigr)^n$ となる。例えば，$b_{10} = -2 + 25 p^2 - 50 p^4 + 35 p^6 - 10 p^8 + p^{10}$ だ。

これをMathematicaでプロットしたら，こんな感じのグラフが出てきた。ただし，$a$が実数ならば，$| p | \ge 2$ なのだ。

図：Plot[b[10], {p, -2.4, 2.4}]のグラフ

三角関数の積

その鈴木貫太郎の問題で次の式の値を求めよというのがあった。

$\cos{\dfrac{\pi}{33}}\cos{\dfrac{2\pi}{33}}\cos{\dfrac{4\pi}{33}}\cos{\dfrac{8\pi}{33}}\cos{\dfrac{16\pi}{33}}$

三角関数の積和の式を繰り返し使うのか，それにしても面倒だろう，どうするのかなあと解答を見ると，$\sin{\dfrac{\pi}{33}}$をかけて倍角の公式からドミノ倒しのようにしてあっという間に解けてしまった。なるほどね。

そこでMathematicaで一般化してみた。

f1[n_, m_] := Product[Cos[2^k Pi/(2^n + 1)], {k, 0, m}] // Simplify
g1[j_] := Table[f1[j, i], {i, j - 1, 20, j}]
g1[5]
{1/32, -(1/1024), 1/32768, -(1/1048576)}

f2[n_, m_] := Product[Cos[2^k Pi/(2^n - 1)], {k, 0, m}] // Simplify
g2[j_] := Table[f2[j, i], {i, j - 1, 20, j}]
g2[5]
{-(1/32), -(1/1024), -(1/32768), -(1/1048576)}

フェルマーの小定理

フェルマーの小定理は，素数を$p$とし，$p$とは互いに素な整数を$a$として，$a^{p-1} \equiv 1 \quad (\mod p \ ) $というものである。

頭の体操のためによく見ている鈴木貫太郎の YouTubeチャンネルでは，ちょっと目には解きにくそうな整数問題に，縦横無尽にmodを使って条件を絞り込んでいくというパターンがよく見られる。自分が，高校生のときには，こんなふうにしてmodを活用するということがなかったので，新鮮な感じがしている。ここまで自由に使えれば便利には違いない。

そんなわけで，フェルマーの小定理あるいは，その他のmod問題の勘を養成するための１行コードをMathematicaで書いてみた。

f[n_] := Table[ Table[Mod[i^k, Prime[n]], {i, 1, Prime[n]}], {k, 1, Prime[n] - 1}]
g[n_] := Table[Table[Mod[i^k, n], {i, 1, n}], {k, 1, n - 1}]

Do[{Print[g[i], " "]}, {i, 2, 7}]
{{1,0}}
{{1,2,0},{1,1,0}}
{{1,2,3,0},{1,0,1,0},{1,0,3,0}}
{{1,2,3,4,0},{1,4,4,1,0},{1,3,2,4,0},{1,1,1,1,0}}
{{1,2,3,4,5,0},{1,4,3,4,1,0},{1,2,3,4,5,0},{1,4,3,4,1,0},{1,2,3,4,5,0}}
{{1,2,3,4,5,6,0},{1,4,2,2,4,1,0},{1,1,6,1,6,6,0},{1,2,4,4,2,1,0},{1,4,5,2,3,6,0},{1,1,1,1,1,1,0}} 

お化け煙突

 NHKの映像の世紀プレミアムの再放送，第15回の「東京 夢と幻想の1964年」では，オリンピック前後の東京の様子が映し出されていた。もちろんオリンピックの話題が中心なのが，その前後には意外な都市の表情があった。

１つは，当時の東京が非常に汚かったということだ。家庭のゴミはその辺の道路や川に無造作に捨てられて，街中に腐臭が漂っているようだ。また，血液銀行という名の売血制度があった時代でもあり，日銭を稼ぐための建設労働者が行列していた。その中で，首都高速道路が建設され，東海道新幹線が開業している。

エピソードの一つにお化け煙突の解体というのがあった。隅田川沿いの足立区にあった東京電力の千住火力発電所の４本の煙突のことだ。もちろん現物を見たことはないが，子供のころから知っていた。

それは，偕成社の少年少女ものがたり百科（10）「数のふしぎ・形のなぞ」で取り上げられていたからだ。お化け煙突は，4本の煙突が2本づつ平行に配置されていて，見る方向によって，1本，2本，3本，4本に見えるというものだ。

この他にも，曽呂利新左衛門の褒美の話（等比数列の和），地球の鉢巻の話（円周と半径の関係），一億まで数えるのにかかる時間の話，ガウスのレンガ積みの話（等差数列の和），アルキメデスの原理の話，パスカルの三角形の内角の和の話などが印象深かった。

1962年の出版で，小学校低学年で買ってもらったお気に入りの本の一つで何度も何度も繰り返して読んでいた。

写真：「数のふしぎ・形のなぞ」の書影（日本の古本屋から引用）

立体角

 平均自由行程を考えるために色々試行錯誤していたら，立体角の計算が必要になった。球の外側にある点から球を見込む立体角である。図のように，半径$a$の球の中心${\rm Q}$から，$\overline{\rm PQ}=\ell (>a)$の距離に点${\rm P}$をとる。

${\rm P}$から球を見込んだ時の球との接円が$x-y$平面にできるとする。接円と$y$軸の交点を${\rm A, B}$とすると，$\overline{\rm PA} = \overline{\rm PB} = \sqrt{\ell^2-a^2}$となる。

立体角$\Omega$は次の積分で与えられる。$\Omega = \int_0^{2\pi} \int_{\pi/2}^{\pi/2+\alpha} \sin\theta d\theta d\phi = 2\pi [-\cos \theta]_{\pi/2}^{\pi/2 + \alpha} = 2\pi \sin \alpha$ 。ただし，$ \alpha = \angle {\rm QPA}$であり，$\sin \alpha = \frac{a}{\ell}$。

図：球を見込む立体角

ルジャンドル変換（４）

 ルジャンドル変換（３）からの続き

ルジャンドル変換の例題を考えていて，$f(x=0)=f^*(p=0)=0$ならば，積分の図とうまく整合するのだが，そうでない場合にどうなるのかちょっと困った。そこで，簡単な２次関数の例で試してみた。

図：定数分の補正で説明できるルジャンドル変換の例

図左のように元の関数を，$f(x)=(x-a)^2+b$とする。このとき$x$が増加して$f(x)$が減少する部分があるので，正の面積だけでは表現できない。実際，図中で$f'(x)$は$0<x<a$で負になる。そこで，$f(x) = \int_0^x f(y) dy + C$として，初期値を$f(0)=(a^2+b)=C$，$x$軸以下の部分の面積を負の量として考える。

このとき，$f^{*}(p)= \int_{-2a}^p {f^{*}}'(q) dq - C$として，$f(x) + f^{*}(p) = x\,p$が成り立つことになる。図右には，結果としての$f^{*}(p)$のグラフを示している。このように積分領域や定数を適当に調整することで，面積によるルジャンドル変換の解釈はそのまま使えると思われる。

この例では，$x$の範囲は，$0<x$としているが，$x_\min < x < x_\max$の範囲で$f(x)$を定義し，これに対応する $p_\min < p < p_\max$ の範囲の$f*(p)$を考えて，定数分の補正をすれば良い。

ルジャンドル変換（３）

 ルジャンドル変換（２）からの続き

田崎さんの熱力学の付録HのLegendre変換によれば（変数$\alpha \rightarrow p$として），$f(x)+f^*(p) \ge x p$ というYoungの不等式が成り立つとある。ルジャンドル変換を一般化した凸共役性のところでも，フェンシェル=ヤングの不等式として示され，これらはルジャンドル変換の定義（min/maxを用いるもの）から直ちに導かれるとしている。

ところが，積分で表示した場合はどうなのかちょっと困った。まあ，不等号でmin/maxの条件に当てはまらないところを考えるということならば，単調増加する連続関数の積分における不等式で説明できたということにすれば良いのかな。

（例1）$f(x)=a x^2 + b x$のルジャンドル変換。$f'(x)=2 a x+ b =p$から$x=\frac{p-b}{2a}$が最小値を与える$x$である。これを用いて，$f^*(p) = x p - f(x) = x p - (a x^2 + b x) |_{x=\frac{p-b}{2a}} = \frac{(p-b)^2}{4a}$

（例2）$\alpha, \beta > 1, \frac{1}{\alpha}+\frac{1}{\beta}=1$として，$f(x)=\frac{1}{\alpha}x^\alpha$のルジャンドル変換。$f'(x)=x^{\alpha-1}=p$から$x=p^{1/(\alpha-1)}$が最小値を与える$x$である。これを用いて，$f^*(p) = x p - f(x) = x p -\frac{1}{\alpha}x^\alpha |_{x=p^{1/(\alpha-1)}} = p^{\alpha/(\alpha-1)} - \frac{1}{\alpha}p^{\alpha/(\alpha-1)} = \frac{\alpha - 1 }{\alpha} p^{\alpha/(\alpha-1)} = \frac{1}{\beta} p^\beta$となる。すなわち，$\dfrac{x^\alpha}{\alpha} + \dfrac{p^\beta}{\beta} \ge x p$ が成り立つ。

（例3）対応する変数を$(\dot{x}, p)$として，$f(\dot{x})=L(\dot{x},x)=\frac{1}{2}\dot{x}^2-\frac{1}{2}x^2$のルジャンドル変換。$\frac{\partial L(\dot{x},x)}{\partial \dot{x}}=\dot{x} = p$となる。これを用いて，$f^*(p)=H(p,x)=\dot{x} p - L(\dot{x},x) = p^2 - L(p,x) = \frac{1}{2}p^2 +\frac{1}{2} x^2$となる。

ルジャンドル変換（２）

 ルジャンドル変換（１）からの続き

ルジャンドル変換にかかわる関数として，２階微分が正であるC1級関数ではなくて，凸関数という表現を使っているのは，相転移点で微分可能でなくなるような熱力学的関数に対してもルジャンドル変換をしたいがためだとあった。

ということで，ある関数$f(x)$の１階微分$f'(x)$が不連続になるような場合を図示してみると次のようになる。この場合も$f(x)$は連続になっている。

通常，このような場合の関数$f(x)$のルジャンドル変換$f^*(p)$は，次の式で表現されている。

$f^*(p) = - \underset{x}{\min} \{ f(x) - p\, x \} =  \underset{x}{\max} \{ p\, x - f(x) \}$

一階微分できる点の場合はカッコ内を$x$で微分すれば，$p$と$x$が一意的に対応する。そうでない点の場合は，その点の左微分から右微分の値の範囲の$p$に対して，上記の条件から$f^*(p)$を定めることになる。

図：一階微分が不連続な場合のルジャンドル変換のイメージ

ルジャンドル変換（１）

 ルジャンドル変換は，積分を使って表現するとわかりやすいという説と通常の説明の対応について考えてみたい。この過程でtikzにおける塗りつぶし方法を練習した。

図：ルジャンドル変換のための説明図

原点を通る単調増加関数上の点C $(x,p)$があって，図のように矩形の領域をとる。矩形領域内で関数の下の部分の面積を$f(x)$，上の部分の面積を$f^*(p)$とすると，$f(x)+f^*(p)=xp$が成り立つ。$f^*(p)$が$f(x)$のルジャンドル変換になる。同様に，$f(x)$は$f^*(p)$のルジャンドル変換である。$f^*(p)$のルジャンドル変換は$f^{**}(x)$とも書けるから，$f^{**}(x)=f(x)$となって，ルジャンドル変換を2回繰り返すと元の関数に戻る。

この単調増加関数は，$x$の関数とすると$f'(x)$であり，$p$の関数と見れば${f^{*}}'(p)$となる。そこで，$f^*(p)$を求めるには，$p=f'(x)$を$x$について解いて，$x=\varphi(p)$を求めてから，$f^*(p) = x p - f(x) = \varphi(p) \cdot p - f(\varphi(p))$として求めることができる。

いいかえれば，$f(x) = x p - f^*(p) |_{p=f'(x)}$ として，元の関数$f(x)$が，傾き$p$と切片$ - f^*(p)$でも表現されるということになる。

／begin{tikzpicture}
\tikzstyle{every node}=[font = \Large];
\filldraw (0,0) circle(1pt) node[below left]{O};
\draw[->] (-2,0) -- (8,0) node[right]{$x$};
\draw[->] (0,-2) -- (0,8) node[above]{$p$};
\draw[step=1.0, dotted] (-2,-2) grid (8,8);
\draw [dotted, pattern=north west lines, pattern color=blue] (0,0) -- (1,1.3) -- (2,2) -- (3,2.4) -- (4,3)-- (5,4) -- (6,5.4) -- (6,0);
\draw [dotted, pattern=north west lines, pattern color=red] (0,0) -- (1,1.3) -- (2,2) -- (3,2.4) -- (4,3) -- (5,4) -- (6,5.4) -- (0,5.4);
\draw[domain=0:3, thick] plot(\x,{-0.2*(\x-3.5)^2+2.45});
\draw[domain=3:7, thick] plot(\x,{0.2*(\x-2)^2+2.2});
\filldraw (6,0) circle(1pt) node[below]{$x$};
\filldraw (3,2.4) circle(1pt);
\filldraw (6,5.4) circle(1pt) node[right]{C};
\node[below right] at (6.5,5.4) {$f^{'}(x)$};
\node[above left] at (6,5.4) {$f^{*'}(p)$};
\filldraw (0,5.4) circle(1pt) node[left] {$p$};
\draw[blue] (4,1) node{$f(x)$};
\draw[red] (2,4) node{$f^*(p)$};
\end{tikzpicture}

凸関数

 ルジャンドル変換について調べようと，田崎さんの熱力学や谷村さんの資料を見ていたら，事前知識として凸関数（Convex Function）が要求された。それは次のようなものだった。

(1) 定義：ある区間で定義された実数値関数 $f(x)$ において，区間内の任意の点，$x_1 < x_2$に対して，$0 < \lambda < 1$ として，$f((1 - \lambda) x_1 + \lambda x_2) \le (1 - \lambda) f(x_1) + \lambda f(x_2)$ を満足する$f(x)$は凸関数であるという。

(2) 凸関数は連続である。

(3) 凸関数が２回微分可能な点$x$では，$f''(x) >0$である（２回微分できない点もある）。

(4) 任意の点$x_0$ で右微分$f'_{-}(x_0)$と左微分$f'_{+}(x_0)$が存在し，次の条件，$f'_{-}(x_0) < \alpha < f'_{+}(x_0)$を満足する定数$\alpha$に対して，$f(x) \ge f(x_0) + \alpha (x - x_0)$となる。

(5) 一階微分可能な点$x_0$では，$f(x) \ge f(x_0) + f'(x_0) (x - x_0)$ が成り立つ。

図：凸関数の定義のための補助図