軸対称電荷分布（５）

軸対称電荷分布（４）からの続き　　

クーロンポテンシャルを場合分けによって書くことができたが，実際に計算する場合には，いちいち分けるのは面倒だ。といっても高々六類型なのでどうってことはないといえばない。

軸対称電荷分布（２）でやったように，ヘヴィサイドの階段関数を使えば，$r_<,\ r_> \ $が表現できるので，クーロンポテンシャルのルジャンドル成分$\ \widetilde{V}_n(r)\ $を一つの関数を使って表わせる。そこで，$R(t)=\dfrac{1}{\sqrt{1+t^2}}, \ a=1.0,\ c= 0.707107$の場合について具体的に計算してみた。

rs[s_, r_] := s* HeavisideTheta[r - s] + r*HeavisideTheta[s - r];

rl[s_, r_] := r* HeavisideTheta[r - s] + s*HeavisideTheta[s - r];

R[t_] := 1/Sqrt[1 + t^2];

v[n_, r_] := 3*NIntegrate[ 1 + 0 * LegendreP[n, t] 

    *rs[s, r]^n * s^2/rl[s, r]^(n + 1), {t, 0, 1}, {s, 0, R[t]}]

v0 = v[0, 0]

2.64412

Plot[{v0, v[0, r]}, {r, 0, 3}]

Plot[{v0, v[2, r]}, {r, 0, 3}, PlotRange -> {2.55, 2.65}]

Plot[{v0, v[4, r]}, {r, 0, 3}, PlotRange -> {2.64, 2.665}]

Plot[{v0, v[6, r]}, {r, 0, 3}, PlotRange -> {2.638, 2.648}]

ここで，被積分関数の中に 1 を加えたのは，積分値が 0 になって計算が不安定になるのを防ぐためである。その部分の定数$\ v_0\ $を取り除くために，ルジャンドル関数の前に 0 をかけて積分することで$\ v_0\ $を求め，実際の計算ではこれを取り除く。

その結果は次のようになった。$\widetilde{V}_4(r),\ \widetilde{V}_6(r)$が，$r<c\ $でゼロになるのは，$\int ds\ $からでてくる$\ R(t)^{n-2}\ $がルジャンドル関数の$n$次未満の線形結合になるためである。

図：クーロンポテンシャルのルジャンドル成分 $\widetilde{V}_n(r)\  (n=0,2,4,6)$

世界はＳＦになる？

「10年後，世界はSFになると思う」というAIのレビューをみつけた。著者は@bioshok3さん。

その結論は，ほぼ確実に訪れるであろう未来は一般化知能レベルのAIが2020年代後半から今後社会のあらゆる領域を変えていくであろうこと，である。丁度 1990年後半からのインターネット革命がイノベーションをもたらしたように，30年経って新しい波がきたということらしい。

最近の風潮は，何でもかんでもAIというキャッチコピーでまぶして，利権に結びつけるというものであり，話題も玉石混交でみんな辟易している。知に足のついた[1]丁寧な分析や考察が必要とされている。この論説では，予測集計サイト metaculcus や，合理的推論と意志決定のコミュニティ LessWrong なども参照しながら，未来予測をしている。

少し詳しく見ると，つぎのようにまとめている。

① ANI（狭いAI）～2025年まで

② GI   （一般化知能）　指示待ち人間または猫レベルの知能　2026-2030年頃実現

③ AGI（汎用人工知能）　2035-2040年頃実現

④ ASI （超知能）　2038~2045年頃実現

今，GPT-3などの大規模言語モデルから派生している，テキスト指示からの画像生成AI（text2img），img2img（画像の編集），text2music（音楽生成），text2motion（人間の動き），text2 3D（3次元物体生成），text 2video（映像生成），text2code，text2text（自然言語処理AI）などは，①の範疇のできごとだ。

著者は，そこそこな汎用性と実用性は持つが能動性や好奇心は本質的に持たないようなAIを一般化知能と呼んでいる。問題は，能動性や好奇心を獲得した汎用人工知能の段階がいつくるかということなのだが，そこに至れば超知能まではあっという間だという予想だ。

［１］誤変換だったのが，これはこれでありかと思ったところ，既出のフレーズだった。

［２］Preventing an AI-related catastrophe

軸対称電荷分布（４）

軸対称電荷分布（３）からの続き　　

軸対称電荷分布を持つ変形核の作るクーロン場のルジャンドル成分は次式で与えられる。

$\displaystyle V_n (r, s) = -\dfrac{Ze^2}{4\pi\varepsilon_0}\dfrac{4\pi}{v}P_n(s) \int_0^1 P_n(t)  \int_0^{R(t) } \dfrac{1}{r_>} \Bigl(\dfrac{r_<}{r_>} \Bigr)^n \ r'^2 dr' \ dt$

ただし，$Ze = \rho_0 v$，$\displaystyle v=4\pi \int_0^1 \int_0^{R(t)} r'^2 dr' dt$ である。

そこで，図のように観測点の座標$r$と電荷分布の関係で3つの場合に分けることができる。

図：軸対称電荷分布の断面と観測点の関係

観測点の半径$r$が$a>r>c$を満足する場合の交点（$z$軸に垂直な円）は，図O$_2$における交点の座標を$(x,z)$，$z$軸からの方位角のコサインを$t=\tau$とすると，$(x=r \sqrt{1-\tau^2},\  z=r \tau)$である。そこで，電荷分布断面の曲線と観測点の半径$r$の円の交点$(x,z)$を求めれば，それらから$\tau$が定まる。例えば，楕円の場合は$\tau=\dfrac{a}{r}\sqrt{\dfrac{r^2-c^2}{a^2-c^2}}$である。

$V_n(r, s)= -\dfrac{Ze^2}{4\pi\varepsilon_0}\dfrac{4\pi}{v}P_n(s) \widetilde{V}_n(r)$とおくと

(1) ${\rm O_1}\ (r>a>c, r'<r) $の場合

$\displaystyle \widetilde{V}_n(r)= \dfrac{1}{r^{n+1}}\int_0^1  P_n(t) \int_0^{R(t)}r' ^{\ n+2}dr'\ dt = \dfrac{1}{r^{n+1}}\int_0^1  P_n(t) \dfrac{R(t)^{n+3}}{n+3} dt$

(2) ${\rm O_2}\ (a>r>c) $の場合

　(2-1) $(1>t>\tau,\ r'<r)$の場合

$\displaystyle \widetilde{V}_n(r)= \dfrac{1}{r^{n+1}}\int_\tau^1  P_n(t) \int_0^{R(t)}r' ^{\ n+2}dr'\ dt = \dfrac{1}{r^{n+1}}\int_\tau^1  P_n(t) \dfrac{R(t)^{n+3}}{n+3} dt$

　(2-2) $(\tau>t>0, \ r'<r)$の場合

$\displaystyle \widetilde{V}_n(r)= \dfrac{1}{r^{n+1}}\int_0^\tau  P_n(t) \int_0^r r' ^{\ n+2}dr'\ dt = \dfrac{1}{r^{n+1}}\int_0^\tau  P_n(t) \dfrac{r^{n+3}}{n+3} dt = \frac{r^2}{n+3}\int_0^\tau  P_n(t) dt $

　(2-3) $(\tau>t>0,\  r<r')$の場合

$\displaystyle \widetilde{V}_n(r)= r^n \int_0^\tau  P_n(t) \int_r^{R(t)} r' ^{\ 1-n}dr'\ dt = r^n \int_0^\tau  P_n(t) \Bigl[ \dfrac{r' ^{2-n}}{2-n} \Bigr]_r^{R(t)} dt \cdot (1-\delta_{n2}) $

$\displaystyle + r^n \int_0^\tau  P_n(t) \Bigl[ \log r' \Bigr]_r^{R(t)} dt \cdot \delta_{n2}$

(3) ${\rm O_3}\ (a>c>r) $の場合

　(3-1) $(r'<r)$の場合

$\displaystyle \widetilde{V}_n(r)= \dfrac{1}{r^{n+1}}\int_0^1  P_n(t) \int_0^r r' ^{\ n+2}dr'\ dt = \dfrac{1}{r^{n+1}}\int_0^1 P_n(t) \dfrac{r^{n+3}}{n+3} dt = \frac{r^2}{n+3}\delta_{n0}$

　(3-2) $(r<r')$の場合

$\displaystyle \widetilde{V}_n(r)= r^n \int_0^1  P_n(t) \int_r^{R(t)} r' ^{\ 1-n}dr'\ dt = r^n \int_0^1 P_n(t) \Bigl[ \dfrac{r' ^{2-n}}{2-n} \Bigr]_r^{R(t)} dt \cdot (1-\delta_{n2}) $

$\displaystyle + r^n \int_0^1  P_n(t)  \log R(t) dt \cdot \delta_{n2}\quad$（注）上式の後半は$\Bigl( -\dfrac{r^2}{2-n}\delta_{n0}\Bigr)$とまとまる。

軸対称電荷分布（３）

軸対称電荷分布（２）からの続き　

今から25年前，1997年度の出世・松井・渡会の卒業論文では，変形核の軸対称電荷分布の形として原点からの距離を$R(\cos\theta)$として，2種類の関数形を採用した。なお，$\theta$は対称軸である$z$軸正方向からの方位角を表わし，$x-y$平面に対して対称な分布を考えるので，$R(\cos\theta)$は$\cos\theta$の偶関数に限定される。$t=\cos\theta$とおいて，電荷密度$\rho_0$は，$Ze = 4 \pi \rho_0 \int_0^1 \int_0^{R(t)} s^2 ds\ dt $で与えられる。

(1) 回転楕円体：$R(\cos\theta) = \dfrac{ac}{\sqrt{c^2+(a^2-c^2)\cos^2\theta}}$である。

この場合の電荷密度は，$\rho_0=\dfrac{3Ze}{4\pi a c^2}$となる。

(2) 2次ルジャンドル関数：$R(\cos\theta) = a+(c-a) \cos^2\theta = \alpha + \beta P_2(\cos\theta)$である。

ただし，$\alpha = \dfrac{2a+c}{3c},\quad \beta = \dfrac{2(c-a)}{3}$。

この場合の電荷密度は，$\rho_0=\dfrac{105 Ze}{4 \pi (16 a^3+8a^2c + 6 a c^2 + 5 c^3)}$となる。

ところで，ルジャンドル関数は直交多項式であり，回転楕円体の断面の楕円はルジャンドル関数による展開で近似できる。この近似がどのくらい妥当するかを，長軸と短軸の比が1:2や2:1の場合について確かめてみる。楕円の離心率は，$e=\sqrt{\dfrac{a^2-c^2}{a^2}}$で与えられる，$\varepsilon=\dfrac{e^2 a^2}{c^2}$とおくと，$R(t) = \dfrac{a}{\sqrt{1+ \varepsilon t^2}}$となる。まず，$n=4$までの展開係数を求めておく。

$\dfrac{1}{\sqrt{1+ \varepsilon t^2}} \approx 1 - \dfrac{\varepsilon}{2}t^2 + \dfrac{3}{8}\varepsilon^2 t^4 = \alpha P_0 (t)+ \beta P_2(t) + \gamma P_4(t)$

$= \alpha + \beta \dfrac{3 t^2-1}{2} + \gamma \dfrac{35t^4-30t^2+3}{8}=1-\dfrac{\beta}{2}+\dfrac{3\gamma}{8} +\Bigl( \dfrac{3\beta}{2}-\dfrac{15\gamma}{4}\Bigr) t^2+\dfrac{35\gamma}{8}t^4$

これを解いて，$\alpha = 1 -\dfrac{\varepsilon}{6} + \dfrac{3\varepsilon^2}{40}$，$\beta = -\dfrac{\varepsilon}{3} + \dfrac{3 \varepsilon^2}{14}$，$\gamma = \dfrac{3\varepsilon^2}{35}$

a = 0.5; c = 1.0;

(* a = 1.0; c = 0.5; *)

d[a_,c_,t_] := a/Sqrt[1 + (a^2-c^2)/c^2 * t^2]

4 Pi * Integrate[Integrate[r^2, {r,0,d[t,a,c]}], {t,0,1}]

g0 = Plot[d[t,a,c], {t,0,1}, PlotStyle->Hue[0]]

f[n_,a_,c_] := Integrate[LegendreP[n,t]*d[a,c,t], {t,0,1}]

b = Table[f[n,a,c], {n,0,16,2}]

g[m_] := Plot[Sum[b[[n/2+1]]*LegendreP[n,t]*(2*n+1), {n,0,m,2}],

         {t,0,1}, PlotStyle-> Hue[(m+1)/10]]

Show[Table[{g0, g[k]}, {k, 4, 6, 2}], PlotRange -> {0.4,1}]

図：軸対称電荷分布（左 oblate，右 prolate，横軸 $t=\cos\theta$，縦軸 $R(t)$）

扁平率 $f=\dfrac{a-c}{a}\ $が-1〜1/2の範囲では，4次のルジャンドル展開で回転楕円体をほぼ表現できる。6次ならば十分だ。

タイパ

「神への長い道」とツイッターからの続き

アインシュタインの特殊相対性理論では，異なった座標系（慣性系）の観測者は，それぞれの時間軸を持っていて，時間の進み方が必ずしも共通とはならない。しかし，我々の日常生活の世界では，相対速度が光の速度に近い観測者はいないから，時間の進み方は共通のはずだった。

NHKの朝のニュースでタイパ（タイム・パフォーマンス）を取り上げていた。具体的には，動画を2倍速でみるような視聴形態が若者には普通に見られるが，年寄りではそうでもないことを結論づける文脈だった。

コスパ（コスト・パフォーマンス）という言葉は，費用対効果という意味で完全に社会に定着し，企業だけでなく行政や公共分野でもあたりまえに使われている。タイパもコスパと結びつけばそうなるだろう。なんといっても時は金なりという，E = m c^2に匹敵するゴールデンルールで結ばれるのだから。

教育系YouTubeでいえば，ヨビノリ（1993-）では，板書を時間短縮する動画編集があたりまえのように導入されていたが，鈴木貫太郎（1966-）の場合は，リアルタイムだけれど，できるだけ高速な論理展開がなされている。コロナ禍で全国的に導入されたオンライン授業の非同期型コンテンツは，学習者が視聴速度を自由にコントロールできる。若者にとっては速度可変視聴があたりまえの光景なのだろう。

ビデオテープの時代だと，早送りはできるとはいってもそれほど便利ではなかった。スマートフォンの中でデジタル動画が自由に操作編集できる時代には，時間の制御はとても簡単になった。その結果，意識や思考における時間の進み方は個人ごとに相対的であるという相対論の世界が実現する。

いや，もともとそうだったのかもしれない。子供の時間と若者の時間と大人の時間と老人の時間は，そもそもズレていたのだ。あるいは，個々人によって理解の速さも表現の速さも異なっている。そのズレをコミュニケーションで同期しながら社会は回ってきた。

新しいテクノロジーは，そのギャップ（＝時間の進み方のズレ）を埋める方に作用するのか，あるいは拡大する方向を推し進めるのか。

一例をあげてみる。大学在職時代に，聴覚障害を持った学生のためのノートテイク支援というボランティアによる活動があった。板書の間に教員が話す説明を聞き取ることができないため，補助者がそれをパソコンで文字起しして，当該学生のノートパソコンに伝えるというものだった。今では，補助者なしのAI音声認識がその水準に達している。これによって，時間のズレは解消されそうだ。

もう一例。放送大学の松井哲男さんの授業で，「初歩からの物理学」や「物理の世界」という講義をときどき聞くことがある。彼は昔からそうだったけれど，頭の回転がはやくて早口で適確に圧縮された説明をする。確かに，教科書の内容はカバーできているしエピソードもふんだんに加えられているのだけれど，分かりやすいかというと，なぜか微妙に未消化感が残ってしまう。講師と受講生に流れる時間のズレのせいではないかと想像している。これはテクノロジーとは関係なかった。

物理的な時間と意識の時間を混同しながら議論してきたのであたまはぐるぐる回って話はまとまらない。いずれにせよ気になっている問題なので，忘れないようにしよう。

図：stopwatch time performance fine style of dali , sketch（DiffusionBeeによる）

軸対称電荷分布（２）

軸対称電荷分布（１）からの続き

前回の$z$軸対称で，$x-y$平面にも対称な電荷分布のつくるクーロンポテンシャルは，

$\displaystyle V(r, s) = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  P_n(s)  \int_{0}^{1}  \int_0^{R(t)} \rho(r', t) P_n (t) \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  4\pi r'^2 dr' d t\ $であった（ただし，$s=\cos\theta$，$t=\cos\theta'$）。

最初に球対称の場合で，上記の式を確認する。このとき$R(t) =R$， $\rho(r', t) = \rho(r') = H(R-r') \rho_0$である。$H(x)\ $はヘヴィサイドの階段関数，$\dfrac{4\pi}{3} R^3 \rho_0 = Ze $である。

このとき$\ n=0\ $の項だけが残り，$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r_>}  4\pi r'^2 dr' $

(1) $r>R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  \int_0^{R} \rho_0  \dfrac{1}{r}  4\pi r'^2 dr'  =  \dfrac{-e}{4 \pi \varepsilon_0}  4\pi  \rho_0  \dfrac{1}{r} \int_0^{R}  r'^2 dr'  =   \dfrac{-Ze^2}{4 \pi \varepsilon_0} \dfrac{1}{r} $

(2) $r<R\ $の場合：

$\displaystyle V(r) =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \int_0^r   \dfrac{1}{r}  r'^2 dr'  + \int_r^R  \dfrac{1}{r'} r'^2 dr'  \Bigr\} $

$\displaystyle  =  \dfrac{-e}{4 \pi \varepsilon_0}  4 \pi \rho_0 \Bigl\{ \dfrac{1}{r} \dfrac{r^3}{3}  +   \left [ \dfrac{r'^2}{2} \right ]_r^R  \Bigr\} =  \dfrac{-Ze^2}{4 \pi \varepsilon_0} \Bigl\{\dfrac{r^2}{R^3} + \dfrac{3(R^2-r^2)}{2R^3}\Bigl\}$

$\displaystyle  =  \dfrac{-Ze^2}{4 \pi \varepsilon_0 }\dfrac{1}{R} \Bigl\{\dfrac{3}{2} - \dfrac{r^2}{2R^2}\Bigl\}$

これで，球対称一様電荷分布によるクーロンポテンシャルが再現できた。

Mathematicaでは次のように計算する。

v[r_] := -3＊NIntegrate[s^2/

       (r*HeavisideTheta[r-s] + s*HeavisideTheta[s-r]), {s, 0, 1}]

Plot[v[r], {r, 0, 5}, PlotRange -> {0, -1.6}]

図：球対称一様電荷分布のクーロンポテンシャル(R=1.0, V(R)=-1.0)

軸対称電荷分布（１）

ドーナツ地球（３）からの続き

ドーナツ地球の重力場のことを考えていたら，昔，卒業論文のテーマに変形核のクーロン場中のミューオン原子波動関数を求めるという課題を出していたことを思い出した。ちょっと再現してみようと試みたけれど，まったく頭が働かない。

過去の卒論は，退職前にあらかた pdf にして保存したはずなのに見つからない。最近ありがちなパターンである。朝早くうとうとしていたら，ふと探し方に気がついた。一度探していた外付けハードディスクを再度検索して，無事に目的のフォルダを発見することができた。

1997年度に研究室配属された出世・松井・渡会の3人チームにテーマを与えて，手間のかかる変形核のクーロンポテンシャルの積分はMathematicaで計算するように指示したものだった。自分ではすぐ導出できたものが，4回生のチームが追いかけてくるのに，かなり時間がかかったのを憶えている。ところがいまでは，その導出方法を自分では簡単に再現できないような始末だ。

さいわい，彼らの卒業論文をみると手順を理解することができそうな・・・。とりあえず，砂川さんの理論電磁気学で学んだように，任意の電荷分布がつくるクーロン場$V(\bm{r})$をルジャンドル関数で多重極展開すればよい。変形核の軸対称電荷分布が$\ell=2$の球面調和関数$Y_{20}(\theta, \phi)$の形の一様分布とすると，回転楕円体の場合よりも計算は簡単になる。

$V(\bm{r}) = \dfrac{-e}{4 \pi \varepsilon_0} \int_V \dfrac{\rho(\bm{r'})}{|\bm{r}-\bm{r'}|}  d \bm{r'} = \dfrac{-e}{4 \pi \varepsilon_0} \int_V \dfrac{\rho(\bm{r'})}{\sqrt{r^2+r'^2-2r r' \cos \omega}}  d \bm{r'} $

$ \displaystyle = \dfrac{-e}{4 \pi \varepsilon_0} \int_V\rho(\bm{r'}) \Bigl\{ \sum_{n=0}^\infty  P_n(\cos\omega) \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  d \bm{r'} $

ここで，$(r_<, r_>)$は $\ r<r' \rightarrow (r, r')\ $または$\ r'<r \rightarrow (r',r) \ $である。さらに，$\cos\omega = \hat{\bm{r}}\cdot \hat{\bm{r'}} = \sin\theta \sin\theta' \cos(\phi-\phi') + \cos \theta \cos \theta'$であり，これはルジャンドル加法定理でさらに次のように展開できる。

$\displaystyle P_n(\cos\omega) = \dfrac{4 \pi}{2 n + 1}\sum_{m=-n}^n Y_{nm}(\theta, \phi) Y_{nm}^*(\theta', \phi') $

したがって， $\displaystyle V(\bm{r}) = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  \sum_{m=-n}^{n} \dfrac{4 \pi}{2n+1} Y_{nm}(\theta, \phi)  \int_V \rho(\bm{r'}) Y_{nm}^*(\theta', \phi') \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  d \bm{r'} $

軸対称かつ$x-y$平面に対称な電荷分布に限定すれば，$m=0 \ $の項だけが残り，$Y_{n0}(\theta, \phi) = \sqrt{\dfrac{2n+1}{4\pi}} P_n (\cos\theta)$ を使うと，クーロンポテンシャルは，$\displaystyle V(\bm{r}) = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  P_n(\cos \theta)  \int_V \rho(r', \cos \theta' ) P_n (\cos \theta') \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  2\pi r'^2 dr' \sin \theta ' d \theta' $

$\displaystyle = \dfrac{-e}{4 \pi \varepsilon_0} \sum_{n=0}^\infty  P_n(s)  \int_{0}^{1}  \int_0^{R(t)} \rho(r', t ) P_n (t) \Bigl\{   \dfrac{1}{r_>} \bigl (\dfrac{r_<}{r_>}\bigr)^n \Bigr\}  4\pi r'^2 dr' d t$

ただし，$s=\cos\theta, t=\cos\theta'$と置いた。軸対称な一様電荷分布は，$\int_0^{R(t)} \rho(r', t) $ で表現され，$R(t)\ $が方向に依存する一様電荷分布の境界までの距離を表わしている。

図：ルジャンドル加法定理による展開

Mathics（２）

Mathics（１）からの続き

ちょうど2年前の今ごろ，Mathicsについての記事をかいていたが，完全に記憶の彼方に跳んでいた。別の記事を復習していてたまたま目に付いたのである。

Mathicsは，Pythonベースのフリーでオープンソースの汎用計算機代数システムであり，Mathematica®と互換性のある構文と関数を特徴としている。早速，MacBook Air M1で動かしてみると動かない。pip3 install Mathics で簡単に復活した。

Mathics 5.0.2

on CPython 3.10.8 (main, Oct 13 2022, 09:48:40) [Clang 14.0.0 (clang-1400.0.29.102)]

using SymPy 1.10.1, mpmath 1.2.1, numpy 1.23.4

Copyright (C) 2011-2022 The Mathics Team.

This program comes with ABSOLUTELY NO WARRANTY.

This is free software, and you are welcome to redistribute it

under certain conditions.

See the documentation for the full license.

前回より，かなりバージョンがあがっている。 GUIの画面も実装しているようなので試してみたが，これがなかなか難渋する。どうやら，Three.jsというJavascriptによる三次元グラフックシステムを利用しているらしい。これも7-8年前に流行ったようだが全く知らなかった。MathicsのGUI版の方は，Mathics-Djangoという名前で，PythonのウェブフレームワークのDjangoの上で動いていた。

pip3 install Mathics-Django 

cd mathics-django-master/mathics_django

*You have 21 unapplied migration(s). 

python3 manage.py migrate

python3 manage.py runserver 

*error AttributeError: module 'scipy.linalg' has no attribute 'inv'

pip3 install scipy-stack

python3 manage.py runserver

*error sqlite3.OperationalError: unable to open database file

*error django.db.utils.OperationalError: unable to open database file

mkdir ~/.local/var/mathics 

djangoは，作業ディレクトリ（~/mathics-django-master/mathics_django）で，manage.py runserver を実行すると，ローカルのウェブサーバー（ http://127.0.0.1:8000/）が立ち上がる仕組みになっている（エラーが出たら指示に従ってディレクトリに必要なファイルを持ってくる）。最初のエラーは，scipy-stackのインストールで，次のエラーは，データベースの場所の問題だったが，~/.local/var/mathicsというデフォルトのディレクトリを作成することでようやく解決することができた。

早速実行してみたが，相変わらず，計算量の制限は大きく，1500!（10進 4114桁）くらいまでしか計算できない。本物のMathematicaならば，10000000!（10進 65 657 059桁）でも平気で計算できるのだけれど。

なお，List of computer algebra systems には無事に掲載されていた。

図：MacBook Air M1にインストールしたMathics GUI版

マストドン（４）

マストドン（３）からの続き

マストドンについては，5年前にすこし齧っただけだった。そこで，初心者として一から勉強してみることにした。

(1) 入門テキストはいろいろあったけれど，マストドンWikiがよかった。

マストドンはメールとよく似ている。複数のメールサーバは，共通のプロトコル（SMTP）に基づいて，異なったサーバの登録ユーザ同士が互いにメールをやりとりできる仕組みを提供している。一方，マストドンとは個人や組織が運営する複数のサーバの連合（Federation）間の共通のプロトコル（ActivityPub）で記事（ポスト）をやりとりする仕組みのことだ。SNSの特徴であるフォロー関係は自サーバだけでなく，連合に属する他サーバのユーザに及んでいる。

マストドンの記事における文字数上限は，ツイッター（140字）の3.5倍の500字である。ツイート（トゥート）やリツイート（ブースト）はできるが，引用リツイートはできないし，記事内容の全文検索もない。タイムラインには三種ある。(1) ホーム：自分のフォロー関係（連合を含む），(2) ローカル：所属しているマストドンサーバの利用者全体，(3) 連合：所属サーバが連合しているすべてのサーバの利用者の全体。

ツイッターは（他の主なSNSもそうだけれど），一企業（組織）が運営する一つのサーバシステムの元にすべてのユーザが登録されている。このサーバシステムの運営が終了すればそれでおしまい。そんなわけで，イーロン・マスクの動きで不安になったユーザがエクソダスの準備を始めた。マストドンは，一つのサーバの運営が終了しても，そのアカウントはフォロー関係を維持しながらサーバ間で引っ越すことができる（できない部分もある？）。フェデレーション環境は全体として存続しているので，相対的には安心な仕組みだ。

(2) マストドンの公式統計資料は探してもなかなか出てこなかったが，Mastodon instances and network APIを使った，What goes on Mastodon というサイトが見やすかった。そこにある図から読み取ったマストドンサーバの登録者数と投稿数は次のとおり。

サーバ           ユーザ数 ポスト数

mastodon.social 87.7万 42.3M

pawoo.net         77.9万 67.3M

mstdn.jp            28.8万 62.0M

mastodon.cloud   22.8万 4.5M

mastodon.online 16.5万 2.5M

mas.to               12.7万 1.7M

mastodon.world    7.4万 0.2M

図：1700個のマストドンサーバーの分布（横軸: ユーザ数，縦軸: 投稿数）

イーロンマスクは，6兆円を使ってツイッターを買収してCEOになった。その目的の一つが右寄り（≒差別的）言論の開放であり，そのシンボルとしてトランプ元大統領のアカウントを復活させた [1]。しかしトランプ本人はツイッターに戻るつもりはないらしい。

［１］イーロン・マスクがツイッターを買収したわけ（drskyskull@mastodon.social）

［２］The Federation（Fediverseについて）

［３］マストドンサーバー情報 Ver 2

マストドン（３）

マストドン（２）からの続き

そもそも現在のSNS事情はどうなっているかを復習してみた。国内で月当りアクティブユーザ数（MAU）の多い7つのSNSを順に並べると以下のようになる（Social Media Lab 2022.11）。

                国内(MAU) 世界(MAU)

LINE       9,200万 1億9,300万

YouTube     7,000万 20億

note       6,300万  −  

Twitter   4,500万 3億3,300万

Instagram 3,300万 10億

Facebook   2,600万 29億3,000万

TikTok              1,700万 10億

なお，世界的にはこれらの他に，WhatsApp（20億，メタ=米国） ，WeChat（13億…中国国内分，テンセント=中国）などもある。

7つのSNSすべてについて自分のアカウントを作っている。LINEは家族連絡用のみ，YouTubeは情報検索・情報入力用，noteは無料コンテンツに誘導された場合だけ，Twitterは世の中の動きをみるため（エコーチェンバーになっているのであまり意味がないかも），Instagramは散歩スナップの発信用，高齢者ユーザが多いFacebookでは友人・知人の様子がわかる，若者向けTikTokはほとんど使っていない。

マストドンには，分散サーバ（インスタンス）の数が1650ほどある。その合計ユーザ数が700万人だ。各サーバ平均4000人だが，フォローしていればどのインスタンスのユーザでも自分のタイムラインに見える。ただし，Twitterに比べて1桁小さい規模である。Twitterでは900人ほどをフォローしフォローされているが，5年前のマストドンで自分がリンクしていたのは200人ほどであったが多くは休眠状態だ。そんなわけで，タイムラインを眺めていても世界の様子を知ることはまだできない。

図：SNSのイメージ（by DiffusionBee）

弾道ミサイルの軌道（５）

弾道ミサイルの軌道（４）からの続き

NHKの報道によれば，11月18日（金）午前10時14分，北朝鮮のICBMが発射され，午前11時23分に北海道渡島大島の西200kmの海上に落下した。飛行時間は69分，飛行距離は1,000km，最高高度は6,100km，速度はマッハ22とされている。

また，弾頭質量によっては15,000kmを越える射程となり，米国本土を狙えるとあった。しかし，今回の実験もロフテッド軌道だったので射出角は90度にかなり近いはずだ。これをより低い角度にすれば質量を変化させなくても十分到達できると思うけれどどうなのか。とりあえず，自分のモデルで計算して確かめてみよう。

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  45 Degree; T = 5100; 

g = 0.0098; R = 6350; τ = 87; p = 0.75; a = 0.0446; s =  86.2 Degree; T = 4140; 

fr[t_, τ_] := a*Sin[s]*HeavisideTheta[τ - t]

ft[t_, τ_] := a*Cos[s]*r[t]*HeavisideTheta[τ - t]

fm[t_, τ_] := -p/(τ - p*t)*HeavisideTheta[τ - t]

sol = NDSolve[{r''[t] == -fm[t, τ]*r'[t] + h[t]^2/r[t]^3 - g R^2/r[t]^2 + fr[t, τ], 

r[0] == R, r'[0] == 0, h'[t] == -fm[t, τ]*h[t] + ft[t, τ], h[0] == 0}, {r, h}, {t, 0, T}]

f[t_] := r[t] /. sol[[1, 1]]

d[t_] := h[t] /. sol[[1, 2]]

Plot[{6350, f[t]}, {t, 0, T}]

Plot[{f[t + 1] - f[t], d[t]*R/f[t]^2, d[t]/f[t]}, {t, 0, T},  PlotRange -> {-4, 8}]

tyx[T_] := {T, f[T] - R, NIntegrate[R d[t]/f[t]^2 , {t, 0, T}]}

v[T_] := Sqrt[(f[T] - f[T - 1])^2 + (R d[T]/f[T]^2)^2]

{tyx[T], v[T]/0.340}

{{4140, 41.9866, 1036.42}, 22.9628}

g0 = ParametricPlot[{NIntegrate[R d[t]/f[t]^2 , {t, 0, tt}],  f[tt] - R}, {tt, 0, T}]

このモデルの仮定は次のようなものである。(1) 距離の二乗に反比例する万有引力が働く極座標での運動方程式を解く。(2) コリオリの力や空気抵抗は考えない。(3) 推進剤の燃焼噴射によって，投射体の速度ベクトルの方向に等加速度（a=0.0446 km/s^2）を発射時から一定時間加える。これにより投射体の質量は燃料が0になるまで一定の割合で減少する。(4) 多段式ロケットは考えず，一段式として全質量に対する燃料比（p=0.75）を与える。このとき，運動方程式は全質量に依存しない。(5) 飛行時間の観測値を与え，投射角度と推進加速時間を主パラメータとして，飛行時間と最高高度の観測値が再現されるように調整する。

結論としては，弾頭質量（というか全質量に対する燃料比）を変化させなくても，射出角を86.2度から45度に変更すれば，飛行時間5100秒（85分）で到達距離が14,600kmになる。なお，ロフテッド軌道の場合の弾道ミサイルの最高速度はマッハ23となり，報告されたマッハ22をほぼ再現した。

図：同一パラメタで射出角と飛行時間だけを変えた場合の弾道軌道

P. S. 11月19日の北朝鮮の朝鮮中央通信によれば，このロフテッド軌道の弾道ミサイルは，新型大陸間弾道弾「火星17」の試験発射であり，最高高度 6,040.9 km ，飛行距離 999.2 km ，飛行時間 4135 秒とのことだ。なお，同じ火星17は，今年の3月24日には，最高高度 6,248.5 km ，飛行距離 1090 km ，飛行時間 4052 秒で飛行している。

［１］【詳報】北朝鮮ICBM級ミサイル発射 北海道の西 EEZ内落下か（NHK）

［２］北朝鮮のミサイル等関連情報（続報第２報，防衛省・自衛隊）

　　（防衛省からは，わかった部分から順に発表されている。発表時刻を入れるべきでは）

マストドン（２）

マストドン（１）からの続き

イーロン・マスク（1971-）が，2022年10月31日にツイッター社を買収完了してCEOに着任してから混乱が始まっている。なぜかWikipediaの日本語版記事には，ツイッター買収後の話題はほとんどない。

社員の半数の3500人に解雇通告を出したため，モデレーションが機能しなくなり，差別的な投稿内容がほぼ野放し状態になってしまったようだ。また，本人確認アカウントが有料化されて，公式マークと分離されたことで，ニセ公式アカウントが横行しているらしく，この方針はペンディングになっているらしい。さらに，11月18日には，本社オフィスが従業員に対して11月21日まで閉鎖されるという告知がきた。

リバタリアンであるイーロン・マスクが，言論の自由というときには，大きな副作用がともなう。公共セクターが私的所有されることは，メディアだけでなく，エネルギーや交通その他の社会インフラでもこれまでに長い歴史があるのが，マスクは7兆円で買ったツイッターを玩具のようにいじくり回しているようだ。

これに危機感をいだいたユーザや広告主のツイッター離れが進んでいる。その行き先の一つが，マストドンである。2016年にはじまり，2017年の前半にブレイクしたもののあっというまにしぼんでしまった分散型 Twitter 類似SNSだ。数学記号を自由に使えるmathtod.online の方はときどきアクセスしていたが，5年ぶりに日本で最大規模のユーザをかかえるサーバのひとつであるmstdn.jpに帰ってみた。もっかのはやりは，フェディバードらしい。

結城浩や，るまたんなどごく少数のアクティブユーザが帰還しているだけなので，まだどうなるのかわからない。

図：.mastodonの公式イメージかな？

無敵

「無敵」というのは，スタニスラフ・レム（1921-2006，今ではスタニスワム）のSF小説，砂漠の惑星（1964）の原題である。エデン（1959），ソラリスの陽のもとに（1961），砂漠の惑星が，レムのファーストコンタクト三部作とよばれている。フランク・ハーバート（1920-1986）のデューンシリーズの第一作が砂の惑星（1965）と訳されたので，なんだか紛らわしくて困る。早川書房がわるいのか。

そのレムは，サイバネティックスに大きな影響を受けていて，砂漠の惑星では，砂つぶのようなサイズで全体として知能を持ち自己増殖する自動機械の集合体が登場する。それがいよいよ現実の存在になろうとしている。

コーネル大学のグループが，アリの頭よりも小さな電子脳（1000トランジスタのCMOSクロック回路）をマイクロロボット（40 μm × 70 μm）に搭載し，外部からの制御なしに自律的に動くことを可能にした。ロボットの脚は白金ベースのアクチュエーターで，回路と脚の両方が太陽光発電で駆動している。

このマイクロロボットは半導体製造プロセスのような手法を用い，シリコンウェーハで大量生産できる。光センサーや化学物質センサーを搭載して体内に注入すれば，ロボットは折り紙の原理で変形して，その集合体が標的となるガン細胞を取り囲むような将来イメージも描かれている。ほとんど，アイザック・アシモフ（1920-1992）がノベライズした映画のミクロの決死圏（1966）の世界でもある。

写真：コーネル大学のマイクロロボットの顕微鏡写真（2020）

［１］「4つ足で歩く」毛髪よりも小さいマイクロボットが誕生、足の作りは”日本の折り紙”がヒントに

［２］Robots Walk Autonomously Via Microscopic Brains

［３］Brains on board: Smart microrobots walk autonomously（Cornell University）

三宅雪嶺記念館

好きな作家を三人あげなさいといわれたら，最近ならば町田康がその一人に入るような気がする。飯田橋文学会では，日本文学の作家のインタビューシリーズの動画を制作しているが，その最近号に町田康がでてきて，聞き手はちょっといやだったけれど，面白くまじめな話をしていた。

それで，気がつくと三宅雪嶺（1860-1945）を検索していた。別に町田康がその名前を出していたわけではないので，因果関係は特にない。なんだかよくわからない。三宅雪嶺の名前は金沢出身の文筆家ということで意識の片隅にあった。町田康のインパクトで意識空間にさざ波が立ったのかもしれない。

三宅雪嶺記念資料館が，流通経済大学に設けられていて，10年ほど前に10周年を迎えたという記念誌が発刊されていた。雪嶺の孫の三宅立雄名誉教授が雪嶺の資料を寄贈したことで開設に至った。金沢市の新堅町に生まれた三宅雪嶺の資料は，金沢ふるさと偉人館にもある。雪嶺は15歳で金沢を離れ，名古屋の愛知英学校を経て東京開成学校（東京大学）に進学する。

金沢泉丘高等学校の前身の石川県尋常中学校（後の金沢第一中学校）が設立されるのが1893年なので，雪嶺の学生時代とは重ならない。泉丘高校の前庭にあった厳霜碑は明治41年に日露戦争で亡くなった一中卒業生の慰霊碑だ。現在は全物故同窓会員の慰霊の碑となっており，その碑文を三宅雪嶺が書いている。3年間この慰霊碑を目にしており，前でクラス写真を撮ったりしたけれど，それが何かまったく理解していなかった。

三宅雪嶺は，哲学者，国粋主義者，評論家である。その妻の三宅花圃（1869-1943）は，明治時代の初めての女性小説家であり，同門の樋口一葉の誕生に影響を及ぼしている。右翼，ジャーナリスト，政治家の中野正剛は，雪嶺の娘婿である。ナチスに近い右翼政党の東方会を設立した中野正剛は東条英機と激しく対立し，中野正剛事件で自殺に至っている。

写真：金沢泉丘高校の厳霜碑（三宅雪嶺の碑文）

80億人の日

世界人口デーからの続き

2022年11月15日に世界の推計人口が80億人を超えたというニュースが流れている。その内容は，7月11日の世界人口デーに公開された国連人口基金の世界人口推計で示されていた。

11月15日は，「80億人の日」と命名されたらしい。そもそも世界人口をちゃんと数えられているかという問題もある。日本では国勢調査も行っているので，それなりの精度があるのではないかと昔は思っていた。が，昨今の政府の統計偽装か改竄からすれば，日本の公的なデータの信頼性も眉唾でかからないといけないような時代になった。

世界人口時計によれば，1日に34万人生まれ，17万人亡くなっているので，人口の有効桁数だけでなくゆらぎも考慮する必要があるかもしれない。11月15日とピンポイントで指定しているのは，動物園の1000万人目入場者おめでとうみたいな，象徴的な意味だけなのだろう。

国連は，世界人口の場合の○○億人目の節目の赤ちゃんを選んでいるらしい。まあ，エイヤっと決めるのだろうけれど。仮想世界で完全に世界の情報が把握できている場合は，すべての誕生アバターには，その誕生時刻における世界の総アバター数をリンクすることができる。そこで，総アバター数から逆リンクをたどれば，これに対応する誕生アバターの集合が得られる。つまり，○○億人目のアバターはゆらぎのせいで（あるいは単調増加関数でないため）一意的に定まらないということになる。

図：世界人口80億人のスクリーンショット

Intel Fortran

量子物理学の授業は，前期の復習モードが終わりつつある。1次元の有限井戸型ポテンシャル問題の解となるエネルギー固有値と固有関数を求めるというものだった。とりあえず，普通の教科書にあるような，$y=-\frac{x}{\tan x}$と，$x^2+y^2=r^2$のグラフを描いて交点からエネルギー固有値を求めようというお話で終るわけだ。

ところで，現代的な量子力学の教科書ではどうなっているのだろうか。サクライの教科書では，井戸型ポテンシャルは，付録のシュレーディンガー方程式の解の例のところで，結果だけがほんの1pにまとめられていた。堀田さんの教科書では，第11章の粒子の量子的挙動の演習問題にチョロッと顔を出しているだけだった。まあそういうことだ。

気分を変えて1次元ポテンシャルの一般的な数値解法がないかと探してみると，インテルのフォートランを使った解説ページが見つかった。インテルのフォートランは有償だと思っていたけれど，どうやら無償で使えるし，macOSにもかろうじて対応しているようだ。そこで参考サイトにしたがって，早速インストールしてみた。

(1) Intel OneAPI Toolkitsで，Intel oneAPI Base Toolkit とIntel oneAPI HPC Toolkit の２つのパッケージをダウンロードする。その際にユーザ登録が必要となるが無償である。

(2) それぞれのToolkitを展開する（ネットワークインストール版を使った）。

(3) 下記の場所にある環境設定スクリプトに誘導されるのでこれを実行する。

https://www.intel.com/content/www/us/en/develop/documentation/get-started-with-intel-oneapi-base-hpc-macos/top/before-you-begin.html?cid=oth&campid=iags_install&source=installer?cid=oth&campid=iags_install&source=installer

. /opt/intel/oneapi/setvars.sh

 (4) 簡単なプログラム a.f を作って実行してみた。

      implicit complex*16(a-h,o-z)

      write(*,*) "input a,b,c"

      read(*,*) a,b,c

      x1=(-b+sqrt(b**2-4*a*c))/2

      x2=(-b-sqrt(b**2-4*a*c))/2

      write(*,*) " x1= ",x1," x2= ",x2

      end

$ ifort a.f

$ ./a.out

$  input a,b,c

$  1,1,1

$   x1=  (-0.500000000000000,0.866025403784439)  x2= 

 (-0.500000000000000,-0.866025403784439)

できた。

SIN基底とQUADPACKによる1次元時間依存しないシュレーディンガー方程式のページの例題も試してみる。調和振動子のエネルギー準位が再現できた。なお，行列の対角化ライブラリがはいっているインテル・マス・カーネルライブラリのオプション -mkl を -qmkl にせよとの警告がでていた。

$ ifort -qmkl quadpack.f90 main.f90

$ ls

a.out main.f90 quadpack.f90

$ ./a.out

           1  0.499999999999997     

           2   1.50000000000004     

           3   2.50000000000008     

           4   3.50000000000073     

           5   4.50000000003176     

           6   5.50000000018886     

           7   6.50000000502383     

           8   7.50000002123218     

           9   8.50000037227348     

          10   9.50000121353051   

ホール・アース・カタログ

NHKの映像の世紀で先日放映された「世界を変えた“愚か者”フラーとジョブズ」がすごく面白かった。

フラーはアメリカの思想家，建築家であるバックミンスター・フラー（1895-1983）。宇宙船地球号という概念の提唱者であり，ジオデシック・ドームを考案した。1985年に見つかった炭素の同素体のフラーレン（buckminsterfullerene，buckyballs）はフラーにちなんで名付けられている。一方，ジョブズは（Apple）のスティーブ・ジョブズ（1955-2011）である。

この二人の愚か者が，ホール・アース・カタログを経由して強く結ばれていたというのが放送のテーマだった。2005年6月12日にスタンフォード大学の卒業式に招かれたジョブズのスピーチは，簡単な挨拶の後にMy story is about death から始まるが，とても印象深いものだった。その結びの言葉である，Stay Hungry. Stay Foolish. は彼が大きな影響を受けたホール・アース・カタログの最終号からきている。

ホール・アース・カタログは，スチュアート・ブラント（1938-）が創刊したヒッピー文化やハッカー文化などのカウンターカルチャーの雑誌だった。そのブラントが，バックミンスター・フラーの多大な影響を受けていた。カタログの最初には，フラーの4冊の本が取り上げられ，The insights of Buckminster Fuller are what initiated this catalog.とあった。（スチュアート・ブラントが後に原発推進の立場になっていたのはまた別の話，松岡正剛の千夜千冊を参照）

ホール・アース・カタログが1974年に終刊を迎えたとき，余剰金の2万ドルをどうするかが議論となり，それは結局，平和活動家のフレッド・ムーア（1941-1997）に託された。そのフレッド・ムーアは，1975年3月にメンロパークで立ち上げられたホームブリュー・コンピュータ・クラブ（1975-1986）の創設に加わる。ジョブズはアップルの共同創業者となるスティーブ・ウォズニアック（1950-）とともにこのクラブに参加していたのだった。

2011年にジョブズは，MacやiPhoneやiPadを残し膵臓ガンで亡くなった。彼が生前に計画していた，ノーマン・フォスター（1935-）のデザインによるアップル・パーク本社の特徴的なドーナツ型の建物は，2017年にできあがった。ノーマン・フォスターはバックミンスター・フラーの弟子でもあった。

写真：アップルパークのアップル本社（Wikipediaから引用）

一谷嫩軍記

久々の国立文楽劇場で，第二部の一谷嫩軍記を観劇する。

弥陀六内の段（睦太夫・團吾）から，脇ヶ浜宝引の段（織太夫・燕三），熊谷桜の段（希太夫・清丈）まで気持ちよく寝てしまった。織太夫は病休の咲太夫の代演だったようで，ここはちゃり場だとはいうものの，声が大きいだけで粗っぽくあまり聞く気にはならなかった。家人によれば希太夫が良かったということだが，これは残念ながら聞き逃してしまった。

熊谷陣屋の段の前半が竹本錣太夫と竹澤宗助。モチモチした語り癖はあるとはいうものの，うまく語り分けていて，熊谷次郎直実と妻相模と敦盛の母の藤の局のややこしい話が進んでゆく。後半の切は豊竹呂太夫と鶴澤清介。語り始めの三味線の前奏部分は，似たようなメロディだけれど，宗助に比べると清介の音の方が一段と澄んでいる。バチ捌きはともにするどい。呂太夫は相変わらず声が出ていない。瞬間的なバーストがあったものの，なんだか残念だ。

熊谷陣屋の段は何度か見ているはずだけれど，いまだに物語が飲み込めていなかった。義経が出した桜の前の制札にある「一枝を伐らば一指を剪るべし」の後半の「一指」が「一子」にかかっていて，敦盛を討取ったと見せかけて自分（熊谷直実）の子どもの小次郎の首を差し出すことにせよと解釈するところまではよい。そこに絡む藤の局と相模の関係とか，弥陀六の位置づけがわからなかった。今回，肝腎のその説明部分で寝ていたのだけれど，最終盤で義経を挟んで弥陀六，藤の局，相模，熊谷次郎直実の5人が対称的に並んで見得を切っていたのでなんとなく雰囲気はわかった。

今日は珍しく花道が設置されていたが，第三部に弁慶の勧進帳があるからだった。

図：一谷嫩軍記の熊谷次郎直実（Wikipediaから引用）

最小交換硬貨枚数（２）

最小交換硬貨枚数（１）からの続き

前回の結果からいくつかの疑問がでてきた。最小交換硬貨枚数の平均値が最小となる中間硬貨は5円でよいのだろうか。5の倍数以外だと，ドルでは25セント，ユーロでは，0.02ユーロ，0.2ユーロ，2ユーロなどの硬貨がある。中国（元，角）や韓国（ウォン）は日本と同じで5の倍数だけだ。

そこで，2円〜9円までの中間硬貨を設定したときの最小交換硬貨枚数の平均値を求めてみることにする。ついでに，硬貨種類別の交換枚数の合計や情報エントロピーも合わせて計算する。なお，1円から1000円までの一円刻みの商品の支払いを硬貨のみで行い，手持ち硬貨数やおつりの硬貨数は冗長性をはぶいた最小の値の範囲で考える。

function foop(m,y)

    nmin = 100

    mm = div(9,m)

    c=[1,m,10,10*m,100,100*m]

    k=[0,0,0,0,0,0]

    t=[0,0,0,0,0,0]

    for k[1] in -(m-1):(m-1)

    for k[2] in -mm:mm

        for k[3] in -(m-1):(m-1)

        for k[4] in -mm:mm

            for k[5] in -(m-1):(m-1)

            for k[6] in 0:mm+1

                z = 0

                n = 0

                for j in 1:6

                    n = n + abs(k[j])

                    z = z + c[j]*k[j]

                end

                if z==y

#                   println(z," : ",k," : ",n)

                    if n < nmin

                        nmin = n

                        for l in 1:6

                            t[l] = abs(k[l])

                        end

                    end

                end

            end

            end

        end

        end

    end

    end

    return nmin,t

end

function main(m)

  p=[0.,0.,0.,0.,0.,0.]

  s=[0,0,0,0,0,0]

  t=[0,0,0,0,0,0]

  a=zeros(Int,1000)

  sum=0

  for i in 1:1000

    (n,t) = foop(m,i)

    a[i]=n

    sum = sum + n

    for j in 1:6

        s[j] = s[j] + t[j]

    end

#    println(i," ",n)

  end

  inf = 0.0

  for j in 1:6

     p[j] = s[j]/sum

     inf = inf -p[j]*log2(p[j])

  end

  @printf("%.3f : %04d %04d %04d %04d %04d %04d -> %.3f\n",

          inf,s[1],s[2],s[3],s[4],s[5],s[6],sum/1000)

#    println(a)

  plot(a)

end

for i in 2:9

    print(i," ")

    main(i)

end

2 2.301 : 0500 1025 0445 1000 0455 2275 -> 5.700

3 2.473 : 0753 0747 0601 0758 0547 1611 -> 5.017

4 2.552 : 0700 0738 0564 0728 0826 1118 -> 4.674

5 2.487 : 1254 0446 1020 0450 1000 0956 -> 5.126

6 2.466 : 0984 0500 0814 0500 1530 0690 -> 5.018

7 2.491 : 0900 0600 0793 0600 1518 0636 -> 5.047

8 2.458 : 1358 0614 1065 0640 1640 0545 -> 5.862

9 2.314 : 2054 0446 1644 0442 1650 0510 -> 6.746

等比数列の中間値である√10や根拠はないけれどネピア数e 近辺あたりの3円がよいのかと思いきや，4円にすると平均交換硬貨枚数が 4.67枚となって，5円の5.13枚よりも小さくなるのだった。ただし，情報エントロピーの最小値は2円の 2.30なのである。どういう意味かまではわかっていない。

図：1円から1000円までの交換硬貨枚数（中間硬貨は5円）

最小交換硬貨枚数（１）

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

2025年度の大学入学共通テストに新規参入する「情報Ｉ」の試作問題が，日経の朝刊に掲載されていた。

退職前には，監督に当たるのはこれ以上体力的精神的な限界をはるかにこえてもう絶対無理と思っていた大学入学センター試験だ。これが，大学入学共通テストと看板を替えて，教科「情報」を新しい試験科目として追加採用してしまった。情報教育関係の大学や高校の教員の皆さんはお喜びのようであるが，なんだかなぁの案件である。

60分で全問必答の大問が4問出題されている。掲載されていた第3問がプログラミングの問題だったので，さっそくチャレンジしてみた。その問題のテーマは最小交換硬貨枚数だ。最近は，PayPayで支払うことが増えてきたが，財布の硬貨数の制御は老人の認知症防止のために最適である。おつりが増えすぎないように，少し硬貨を加えて切りのよいおつりにするあれね。

例えば，46円の支払いだと，10円×４ + 5円×１ + 1円×１ でもよいし，50円×１ + 1円×１ - 5円×１ もある。ここで，マイナスはお店から戻ってくる硬貨を表わしている。前者の交換硬貨枚数は6枚であり，後者では3枚となり，後者の方が交換硬貨枚数は少ない。

問題の誘導部分を読む前に，さっそくプログラムを作ってみたが，肝腎の多めに払っておつりが返ってくるところに不備がありまくりだった。

function pay(m,y)

    n = 0

    c=[500,100,50,10,5,1]

    d=[y]

    for i in 1:6

        while y-c[i] >= 0

            y = y - c[i]

            n = n + 1

            push!(d,c[i])

        end

    end

    push!(d,-(m+n))

    println(d)

    return m+n

end

function change(y)

    m = 0

    n = 0

    c=[1,5,10,50,100,500]

    d=[y]

    pay(0,y)

    for i in 1:5

        while mod(y,c[i+1])!=0

            y = y + c[i]

            n = n + 1

            push!(d,c[i])

        end

        m1 = m + n

        if m1!=m

            push!(d,-m1)

            print(d)

            pay(m1,y)

            n=0

            deleteat!(d,length(d))

        end

        m = m1

    end

end

change(46)

[46, 10, 10, 10, 10, 5, 1, -6]
[46, 1, 1, 1, 1, -4][50, 50, -5]
[46, 1, 1, 1, 1, 50, -5][100, 100, -6]
[46, 1, 1, 1, 1, 50, 100, 100, 100, 100, -9][500, 500, -10]

そこで，1000円以下なら数もしれているということで，方針を変えて総当たりで確認するアルゴリズムに変更した。全くスマートではないのである。この中から最小値を探すのは目の子でできる。出力される配列 [ ] の中身は，交換される硬貨枚数が並んでいて，マイナスがついたものは店から客へのバック分を表わしている。最後の出力値が交換硬貨枚数だ。

function foop(y)
    c=[1,5,10,50,100,500]
    k=[0,0,0,0,0,0]
    for k[1] in -4:4
    for k[2] in -1:1
        for k[3] in -4:4
        for k[4] in -1:1
            for k[5] in -4:4
            for k[6] in 0:2
                z = 0
                n = 0
                for j in 1:6
                    n = n + abs(k[j])
                    z = z + c[j]*k[j]
                end
                if z==y
                    println(z," : ",k," : ",n)
                end
            end
            end
        end
        end
    end
    end
end

foop(46)

46 : [-4, 0, 0, -1, -4, 1] : 10
46 : [-4, 0, 0, -1, 1, 0] : 6
46 : [-4, 0, 0, 1, 0, 0] : 5
46 : [1, -1, 0, -1, -4, 1] : 8
46 : [1, -1, 0, -1, 1, 0] : 4
46 : [1, -1, 0, 1, 0, 0] : 3
46 : [1, 1, -1, -1, -4, 1] : 9
46 : [1, 1, -1, -1, 1, 0] : 5
46 : [1, 1, -1, 1, 0, 0] : 4
46 : [1, 1, 4, 0, 0, 0] : 6

ここまでくるのに2時間くらいかかったので，とてもじゃないけれど自分の場合60分で大問4問も解けそうにはない。ループ変数を配列にするなんていうのは初めての経験だった。


写真：日本の通常硬貨（造幣局から引用）