JupyterLab

「JupyterLabのすすめ」という記事があったので，JupyterLabをさっそくインストールしてみた。

退職後の４月から非常勤講師としてこれまでと同じような授業を担当する。とりあえず，現在使用しているコンピュータ（MacProやMacBookProなど）を明け渡さなければならないので，MacBookAirに新しい環境を構築して移動作業をしている。その結果，現在の自分の環境は，次のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
※mbaの設定
(1) homebrewの導入
　　　　xcode-select —install （コマンドライン開発ツールの導入）
　　　　sudo xcodebuild -license
　　　　/usr/bin/ruby -e "＄(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
(2) homebrewの練習と基本環境の設定
　　　　brew update
　　　　brew upgrade
　　　　brew doctor
　　　　brew install wget
　　　　brew install nkf
　　　　brew cask reinstall xquartz
　　　　brew install lynx
　　　　brew install w3m
　　　　brew install xxxx （xxxxをインストール）以下適当に繰り返した。
(3) python3環境（jupyterを含む）の設定
　　　　brew install python3
　　　　python3 -V
　　　　　Python 3.7.2
　　　　pip3 -V
　　　　　pip 19.0.2 from /usr/local/lib/python3.7/site-packages/pip (python 3.7)
　　　　python -V
　　　　　Python 2.7.10
　　　　・これで，python3がインストールできた。
　　　　which python
　　　　　/usr/bin/python
　　　　which python3
　　　　　/usr/local/bin/python3
　　　　・pip3のアップデート方法
　　　　pip3 install --upgrade setuptools
　　　　　Requirement already up-to-date:...
　　　　pip3 install --upgrade pip
　　　　　Collecting pip
　　　　...
　　　　　Successfully installed pip-19.0.3
　　　　・環境設定してその中で作業
　　　　python3 -m venv env
　　　　. env/bin/activate
　　　　(env) … python -V
　　　　　Python 3.7.2
　　　　(env) … pip -V
　　　　　pip 18.1 from /Users/admin/env/lib/python3.7/site-packages/pip (python 3.7)
　　　　(env) ... pip install --upgrade pip
　　　　(env) ... pip install --upgrade setuptools
　　　　(env) ... pip install numpy
　　　　(env) ... pip install scipy
　　　　(env) ... pip install scikit-learn
　　　　(env) ... pip install matplotlib
　　　　(env) ... pip install Pillow
　　　　(env) ... pip install sympy
　　　　(env) ... pip install jupyter
　　　　(env) ... deactivate
(4) julia環境の設定
　　　　https://julialang.org/downloads/
　　　　からダウンロードすればコマンドラインで実行可能に
　　　　echo 'export PATH=$PATH:/Applications/
Julia-1.1.app/Contents/Resources/julia/bin' >> ~/.bashrc
　　　　source ~/.bashrc
　　　　julia
　　　　julia> using Pkg
　　　　julia> Pkg.add(“IJulia”)
　　　　export PATH=＄PATH:/Users/admin/env/bin
　　　　jupyter notebook
　　　　これで，jupyter が起動して，julia 1.1 が使えるようになった。
　　　　versioninfo()
　　　　using Pkg
　　　　Pkg.update()
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

(5) jupyterlabの導入と起動は簡単であった。
　　　　pip3 install jupyterlab
　　　　jupyter lab
（付）JupyterLab環境での設定
Settings -> Advanced Settings Editor -> Notebook
codeCellConfigのlineNumbersをTrueにする。
Settings -> Advanced Settings Editor -> Extension Manager -> "enabled": trueにする。

Mathematicaノートブックを長らく使ってきたので，Jupyter Notebook（この記事もWikipedia中文はあるが日本語はない・・・）はすんなりと馴染むことができた。JupyterLabはJupyter Notebookの後継で，今後はこれに移行するのだそうだ。まだ違いがよくわからないが，Jupyter Notebook のことだってそれほどわかっていたわけではない。ぼちぼちやっていくことにしよう。

金沢泉丘高等学校のプール

大乗寺丘陵公園から円光寺の方に下り，泉丘高校に向って歩いた。前身の一中創立70周年から80周年の間で，理数科ができて２年目の1969年に入学した。1984年には新校舎が落成しており，我々の時代の旧校舎はなくなっている。正門のあたりの古い木々や碑などはそのまま残っているが，校舎も制服も生徒も変わってしまった（昔よりみなさんよくできるのである）。グラウンド横の部室裏のプールは当時のままだと思ったが，沿革によると1992年にプールが完成とあるので同じ場所で新しくしたのだろう。校庭の周辺の桜の木とか，相撲場はどうだっただろうか。

写真　金沢泉丘高等学校のプール（2019.3.22）

小学校５年までは泳げなかったので，水泳の授業（小学校にはまだプールがなかった）では，金沢市営の25mプールに集められた数名の泳げない組で練習していた。というか，先生はあまり教えてくれないのである。同級生の畦地勉君に浮き方をおそわって，ようやく潜った状態でばた足で前進できるようになった。

中学校の体育では25mを泳ぐことが要求された。３年間でようやく25mをなんとか泳げたがほぼ息継ぎのないクロールだったかもしれない。背泳ぎも少しできるようになった。中学の体育の先生は親切だったので，泳ぎの苦手な子に飛び込みを教えてくださった。プールサイドに腰掛けた状態から始まって，徐々に慣らして行き，最後にはスタート台から飛び込むことができるようになった。

高校の体育では200ｍを泳ぐことが要求された。これはちょっと難しかったので，夏休み前に10名程度泳ぎの苦手な連中がプールで特別指導された。というかここでも指導された記憶はあまりないのだが，梅雨も明けきらないやや肌寒いプールでばしゃばしゃしていた。授業本番では背泳ぎで200mをなんとか泳ぐことができた。これ以降は夏休みに市営プールに一人で通って1000mを平泳ぎで泳ぐことができるようになった。

大乗寺丘陵公園

お彼岸に金沢に墓参に帰省したとき，大乗寺から岡部病院を抜けると，その横が大乗寺丘陵公園として整備されていた。泉野小学校や野田中学校のときに，毎年花見遠足やスキーに来たあたりである。丘陵公園の下あたりが，金沢泉丘高等学校時代の体育の時間のマラソンコースになっていた。ここからは，金沢平野が一望できる。まだ肌寒くて花見の季節ではなかった。

写真　大乗寺丘陵公園案内図（2019.3.22）

写真　大乗寺丘陵から金沢市街地の方向（2019.3.22）

MathJaxで数式表示（２）

MathJaxの数式表示で行番号がでないのがどうなのかを調べて試しているうちに，調子が悪くなってしまった。古い設定がどこかに残っているためなのか，Safariで数式が小さく表示される現象が生じている。再起動で回復しても，そのページを閲覧した後にはどこでも発症してしまう。どうしたものか。FirefoxやChromeでは問題ない。

今のヘッダーの設定では，標準で連番あり，"¥begin{equation*}"はだめ，"¥[ ¥]"で行番号なし，"¥tag{n}" で番号指定可になっている（¥はバックスラッシュを表す）。

\begin{equation}
\begin{aligned}
\sqrt{x}+\cos{y}\\
\dfrac{\log(x)}{e^x}
\end{aligned}
\end{equation}
\[
\sqrt{\dfrac{z^2+y^2}{x^2}}
\]
\begin{equation}
y= a x^2 + b x + c \tag{4}
\end{equation}
\begin{equation}
m \ddot{\bm{r}} = \bm{F}
\end{equation}

これで"¥bm{}"も使えるようになったのは次のヘッダーである。
参考にしたのは以下のサイトである。
・Mathjaxを使ってBloggerで数式を書く
・おもちゃ箱：Bloggerで数式を表示する
・MathJaxの使い方（黒木玄）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
<script async='async' src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.2/MathJax.js?config=TeX-MML-AM_CHTML' type='text/javascript'/>
<script type='text/x-mathjax-config'>
    MathJax.Hub.Config({
        HTML: [&quot;input/TeX&quot;,&quot;output/HTML-CSS&quot;],
        TeX: {
               Macros: {
                        bm: [&quot;\\boldsymbol{#1}&quot;, 1],
                        argmax: [&quot;\\mathop{\\rm arg\\,max}\\limits&quot;],
                        argmin: [&quot;\\mathop{\\rm arg\\,min}\\limits&quot;]},
               extensions: [&quot;AMSmath.js&quot;,&quot;AMSsymbols.js&quot;],
               equationNumbers: { autoNumber: &quot;AMS&quot; } },
        extensions: [&quot;tex2jax.js&quot;],
        jax: [&quot;input/TeX&quot;,&quot;output/HTML-CSS&quot;],
        tex2jax: { inlineMath: [ [&#39;$&#39;,&#39;$&#39;], [&quot;\\(&quot;,&quot;\\)&quot;] ],
                   displayMath: [ [&#39;$$&#39;,&#39;$$&#39;], [&quot;\\[&quot;,&quot;\\]&quot;] ],
                   processEscapes: true },
        &quot;HTML-CSS&quot;: { availableFonts: [&quot;TeX&quot;],
                      linebreaks: { automatic: true } }
    });
</script>
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

国立大学法人法と学長の任期

以下は国立大学法人法からの学長の任期に係わる部分の抜粋である。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
国立大学法人法

第十条　各国立大学法人に、役員として、その長である学長及び監事二人を置く。
２　略

第十二条　学長の任命は、国立大学法人の申出に基づいて、文部科学大臣が行う。
２　前項の申出は、第一号に掲げる委員及び第二号に掲げる委員各同数をもって構成する会議（以下「学長選考会議」という。）の選考により行うものとする。
一　第二十条第二項第三号に掲げる者（注１）の中から同条第一項に規定する経営協議会において選出された者
二　第二十一条第二項第三号又は第四号に掲げる者（注２）の中から同条第一項に規定する教育研究評議会において選出された者
３　前項各号に掲げる者のほか、学長選考会議の定めるところにより、学長又は理事を学長選考会議の委員に加えることができる。ただし、その数は、学長選考会議の委員の総数の三分の一を超えてはならない。
４　学長選考会議に議長を置き、委員の互選によってこれを定める。
５　議長は、学長選考会議を主宰する。

第十五条　学長の任期は、二年以上六年を超えない範囲内において、学長選考会議の議を経て、各国立大学法人の規則で定める。
２　略
３　略
４　役員は、再任されることができる。以下略

（注１）経営協議会構成員から，学長と学長が指名する理事及び職員を除いた「三　当該国立大学法人の役員又は職員以外の者で大学に関し広くかつ高い識見を有するもののうちから，次条第一項に規定する教育研究評議会の意見を聴いて学長が任命するもの」

（注２）教育評議会構成員から，学長と学長が指名する理事を除いた「三　学部、研究科、大学附置の研究所その他の教育研究上の重要な組織の長のうち、教育研究評議会が定める者，四　その他教育研究評議会が定めるところにより学長が指名する職員」

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

したがって，国立大学法人法は学長が無期限に再任されることを排除しておらず，それは各国立大学法人の学長選考に関る規程に委ねられているようにみえる。

退職挨拶の練習（２）

懺悔の時間がやってまいりました。理科教育講座の越桐です。

1982年に本学に教務員として着任して37年になります。最初の10年ほどは天王寺分校に研
究室がありました。そのころは附属学校のグラウンドの裏にあった池田宿舎に住んでいました。その当時いっしょにすんでいた方は栗林学長と赤松先生だけになりました。

小専理科の実験の授業を池田分校の物理の伊藤太郎先生とともに担当していましたが，授業は３限目からなので，昼頃に池田分校に行き，夕方前のまだ明るいうちに帰宅するという夢のような日々を過ごしていました。ごめんなさい。

先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん力を合わせてこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。

運動物体からの斜方投射

昨日届いた「大学の物理教育（2019 Vol.25 No.1）」の講義室に，茨城大学工学部の坪井一洋さんが，「質点の投射角を再考する 大学の物理教育25(2019)34-37」という記事を書いていた。とてもおもしろかったので，その計算をフォローしてみることにした。


坪井さんの記事の趣旨は次のとおりである。一様重力場での放物運動では，斜方投射の投射角を45度としたときに質点を最も遠くまで投げることができる。しかし現実の現象（ゴルフ，野球，走り幅跳び，砲丸投げ）ではそうなっていない。空気抵抗で説明できる部分もあるがそうでない部分もある。スポーツ科学では初速度の大きさはが初速角に依存するという議論もある。そこで，これをモデル化して物理教育（初等力学）の題材にしようというものである。

図　運動物体からの斜方投射モデル

原点をOとする慣性系Sを考える。Sのx軸上を一様な速度$u$で正方向に運動している物体がある。この物体がS系の原点Oを通過した瞬間（$t=0$）に，質量$m$の質点を，物体に固定した座標系Mからみて図のように初速度$v$，投射角$\phi$でx-y平面内に投射する。この質点には一様な重力$mg$が y 軸負方向に働いている。放物運動した質点は原点からの距離が$L$であるx軸上の点まで到達したとする。なお，S系からみた質点の速度ベクトル$\boldsymbol{w}$は図のように与えられる。

S系における原点からの距離$L$を最も大きくするにはどうすればよいか。原点における物体の初速度を$(V_x,V_y)$とすれば，高等学校の物理では次の結果を得ている。$L=\dfrac{2V_xV_y}{g}$

ところでこの$V_x$と$V_y$は今のモデルでは次のように与えられる。ただし，$0 < \phi < \pi/2$，$0 < \theta < \pi/2$とする。

\begin{equation}
\begin{aligned}
V_x = w \cos \theta = u + v \cos \phi \\
V_y = w \sin \theta = v \sin \phi
\end{aligned}
\end{equation}

つまり，$u,v$が与えられたときに，$L$を最大にする$\phi$または$\theta$を求めるのがここでの問題となる。なお，$w$は，ベクトルの図に余弦定理をあてはめた次の式に従うため，$\theta$の関数であると考えることにする（$\phi$の関数とすることもできる）。
\begin{equation}
w^2+u^2-2w u \cos\theta = v^2
\end{equation}
$w$についての２次方程式を解いて具体的な$\theta$の関数形を求め，さらに$\dfrac{d w}{d \theta}$を計算しておく。
\begin{equation}
\begin{aligned}
w &= u \cos\theta + \sqrt{v^2-u^2\sin^2\theta}\\
\frac{dw}{d\theta} &= -u\sin\theta + \frac{-u^2\sin\theta \cos\theta}{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -u\sin\theta \ \frac{\sqrt{v^2-u^2\sin^2\theta} + u \cos\theta }{\sqrt{v^2-u^2\sin^2\theta}}\\
&= -\frac{u\sin\theta\ w }{\sqrt{v^2-u^2\sin^2\theta}}
\end{aligned}
\end{equation}
以下では，$\gamma=\dfrac{u}{v}$とおく。

(1) $L(\phi)$を最大にする$\phi$に対する条件を求める。
\begin{equation}
L(\phi)=\frac{2\ (u+v\cos\phi)\ v\sin\phi}{g}\\
\frac{dL}{d\phi}=\frac{2v}{g}\{-v \sin^2\phi +(u + v\cos\phi) \cos\phi \}=0\\
2v \cos^2\phi +u\cos\phi -v=0\\
2 \cos^2\phi + \gamma \cos\phi + 1 = 0\\
\therefore \cos\phi = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
(2) $L(\theta)$を最大にする$\theta$に対する条件を求める。
\begin{equation}
L(\theta)=\frac{2\ w\cos\theta\ w\sin\theta}{g}=\frac{w^2 \sin 2\theta}{g}\\
\frac{dL}{d\theta}=\frac{2w}{g}\,\frac{d w}{d\theta}\,\sin 2\theta + \frac{2w^2}{g}\,\cos 2\theta =0\\
\frac{-u\sin\theta \ \sin 2\theta}{\sqrt{v^2-u^2\sin^2\theta}}+ \cos 2\theta =0\\
4 u^2 \sin^4\theta\,(1-\sin^2\theta) = (v^2-u^2\sin^2\theta)(1-2\sin^2\theta)^2\\
v(1-2\sin^2\theta) = u\sin\theta\\
2 \sin^2\theta + \gamma \sin\theta + 1 = 0\\
\therefore \sin \theta = \frac{\sqrt{\gamma^2+8}-\gamma}{4}
\end{equation}
従って，最大の$L$を与える$\phi$と$\theta$は，$\phi+\theta=\frac{\pi}{2}$を満足している。

バラカ：桐野夏生

まわりからどんどん書店が消えていく。昔は，毎日の仕事の帰りには本屋に立ち寄って，新刊をチェックすることができたが，最近はその習慣が失われてしまった。先日，久しぶりに上本町の近鉄の書店（JUNKUDO）に立ち寄ったときに見つけたのが，桐野夏生の「バラカ（上・下）集英社文庫」だった。

桐野夏生の本は，昔「ＯＵＴ」を図書館で借りて読んだ。その他にいくつか文庫本を買っていたはずだと思って調べると，昨年の部屋の片づけの際に処分していた。確か「柔らかな頬」と「グロテスク」だった。スパーク・ジョイがなかったのである。どんなストーリーだったかを思い出そうとしても印象に残っている部分がない。その点，「ＯＵＴ」にははっきりした読後記憶がある。桐野夏生は1951年に金沢に生まれているけれど，３歳で引っ越しているので，自分の世界線との交わりはほとんどない。

さて，「バラカ」は，福島の原発事故が起こったもう一つの虚構の日本での出来事を描写している。部分的には納得できるのだが，全体として詰めが十分ではないと感じた。自分が期待していたのは，現実の日本を照射するもう一つの世界の精密な描写だった。著者の関心はその物語世界の中の人物造形や関係の方にあって，それはそれでおもしろかったが・・・後半はご都合主義的なまとめ方のような・・・

P. S. 作中でバラカ／おじいちゃんが愛読するのがカズオ・イシグロの「わたしを離さないで」だったので，これは読んでみようと思った。

素数日

ほぼ毎日（多少前後して編集しているけれど）ブログを書き始めて，今日で100記事になった。いつまで続くのだろう？何のために続けているのだろう？惰性なのだろうか？

西暦（グレゴリオ暦）のyyyy年mm月dd日を8桁の整数 yyyymmdd で表したものが素数であるものを素数日と言い習わしているようだ。あちらこちらで目にします。で，Juliaで素数日を求める関数を書いてみた。あらかじめ，Primesパッケージを読み出すので，自分の仕事はほとんどない。

using Pkg
Pkg.add("Primes")
using Primes
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function primeday(y)
  u=28
  if(y%4==0 && y%100!=0 || y%400==0)
    u=29
  end
#  println(u)
  mx=(31,u,31,30,31,30,31,31,30,31,30,31)
  for m in 1:12, d in 1:mx[m]
    if(isprime(y*10000+m*100+d)) 
      println(y*10000+m*100+d) 
    end
  end
end

primeday(2019)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20190221
20190227
20190301
20190319
20190323
20190421
20190523
20190529
20190601
20190613
20190719
20190811
20190823
20190913
20191009
20191027
20191109
20191117
20191231

機械学習と物理

長らく物理学会には足を向けていない。今年の第74回年会は九州大学の伊都キャンパスで開催されている。ちょうど今の時間に，シンポジウム「機械学習と物理」が行われている。こんな感じ。

１（一般シンポジウム講演）はじめに
　阪大理・物理，橋本幸士
２（一般シンポジウム講演）物性物理のグランドチャレンジに対する重回帰分析と機械学習
　東大院工，今田正俊
３（一般シンポジウム講演）強い相互作用の最難問 — 中性子星の状態方程式
　東大理・物理・原子核理論，福嶋健二
４（一般シンポジウム講演）データ駆動手法による相関物質の予測と理解
　産総研 CD-FMat，三宅隆
５（一般シンポジウム講演）機械学習によるマルコフ連鎖モンテカルロ法の高速化へ向けて
　理研(AIP/iTHEMS), 慶應大・数理，田中章詞
６（一般シンポジウム講演）機械学習による特徴抽出と，繰り込み群や熱力学との関係
　OIST，船井正太郎
７（一般シンポジウム講演）広域撮像宇宙サーベイによるビッグデータ宇宙論
　東大理・物理・宇宙理論，吉田直紀
８（一般シンポジウム講演）量子力学と機械学習の数理
　東北大院情報科学，大関真之

ちなみに，昨年の2018年度日本物理学会科学セミナーのテーマも「ＡＩ（人工知能）と物理学（東京大学駒場キャンパス　数理科学研究棟　大講義室）」であった。そのプログラムは以下の通り（ちょっと被っている）。

8月11日（土・祝）10:00-16:30
１　はじめに
　日本物理学会会長　川村　光
２　情報処理技術としてのAI
　中島 秀之（札幌市立大学 学長）
３　思考力を競うゲームの人工知能技術発展の歴史と現状
　保木 邦仁（電気通信大学大学院情報理工学研究科 准教授）
４　広域宇宙撮像データを用いたビッグデータ宇宙論
　吉田 直紀（東京大学大学院理学系研究科 教授）
５　量子コンピュータが人工知能を加速する
　大関 真之（東北大学大学院情報科学研究科 准教授）
６　深層学習と時空
　橋本 幸士 (大阪大学大学院理学研究科　教授)
8月12日（日）10:00-16:40
７　深層学習とは何か、そしてどんなことが出来るようになっているのか
　瀧 雅人（理化学研究所数理創造プログラム(iTHEMS) 上級研究員）
８　多層畳み込みニューラルネットワークで求めた量子相転移の相図
　大槻 東巳（上智大学理工学部機能創造理工学科 教授）
９　人工知能と脳科学
　甘利 俊一（理化学研究所脳神経科学研究センター　特別顧問）
10　量子力学の問題をニューラルネットワークで解く
　斎藤 弘樹（電気通信大学大学院情報理工学研究科 教授）
11　AI は物理において何の役に立つか？
　寺倉 清之（物質・材料研究機構 名誉フェロー・エグゼクティブアドバイザー）
12　おわりに
　日本物理学会科学セミナー担当理事　迫田 和彰

日本物理学会誌の方でも「シリーズ　人工知能と物理学」という特集が開始され，産総研の神嶌敏弘さんの「変わりゆく機械学習と変わらない機械学習」が読める。そのうち日本物理教育学会誌にも，「機械学習と物理教育」とか「人工知能と物理教育」などの特集が組まれる時代がくるのだろうか。５年後？１０年後？

参考：Quantum Machine Learning（Wikipedia）
List of Acronyms
ANN: Artificial neural network
BM: Boltzmann machine
BN: Bayesian network
CDL: Classical deep learning
CML: Classical machine learning
HMM: Hidden Markov model
HQMM: Hidden quantum Markov model
k-NN: k-nearest neighbours
NMR: Nuclear magnetic resonance
PCA: Principal component analysis
QBN: Quantum Bayesian network
QDL: Quantum deep learning
QML: Quantum machine learning
QPCA: Quantum principal component analysis
QRAM: Quantum random access memory 
RAM: Random access memory
SQW: Stochastic quantum walk 
SVM: Support vector machine
WNN: Weightless neural network

円周率の分数近似

円周率の日の余震が続いている。3月14日，Googleは円周率の計算の世界記録を更新したことを発表したとの報道があった。31.4兆桁（31,415,926,535,897桁）までの計算を111.8日（start to end で121.1日）かけて計算したそうだ。このプロジェクトを主導したのは，筑波大学出身の技術者 Emma Haruka Iwao （岩尾・エマ・はるか）さん。今年の1月には終わっていたが，このよき日に発表された。つぎに記録に挑戦するグループは，314兆桁を越えないといけないのか。

さて，だいぶ昔（2011.6.30）のWolfram Blogで，Jon McLooneが，"All Rational Approximations of Pi Are Useless"という記事を書いている。

どういうことかというと・・・円周率のよくある近似値22/7を例にとってみる。22/7 = 3.1426... は小数点以下2桁まで正しいので，小数点も含め，3 .  1 4 の4文字が表現されたことになる。このために必要な数は 2 2 / 7 の4文字である。これじゃ文字数が節約できていない。どちらを覚えてもいっしょじゃないか，ということらしい。

もっとやってみると，355/113 = 3.141592... これは小数点以下6桁すなわち8文字まで正しい円周率に対応している。一方，355/113は7文字，さきほどよりちょっとだけましだ。

McLooneの記事にあるMathematicaの関数を用いて，小数点以下50桁のπとその分数表示を求めると，

\begin{equation}
\frac{26151465932107044561886949}{8324270144388272579650158} = \\
　\\
3.14159265358979323846264338327950288419716939937510...
\end{equation}
ということで共に52文字必要としている。でこぼこはあるが，おおむねこの比が1で推移することから，分数表示意味ないじゃんということになったようだ。

必要なMathematica関数は以下のようになっている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
In[1]:= characterCount[expr_] := StringLength[ 
  StringReplace[ToString[Unevaluated[expr], InputForm], 
   " " -> ""]]; SetAttribute[characterCount, HoldAll];

In[2]:= matchingCharacters[expr_, target_] := 

 Ceiling[-Log10[Abs[expr - target]]] + 1

In[3]:= $MaxExtraPrecision = Infinity;

In[4]:= candidates = Table[Rationalize[Pi, 10^(-n)], {n, 1001}];

In[5]:= a = Table[candidates[[k]], {k, 1, 10}]

Out[5]= {22/7, 22/7, 201/64, 333/106, 355/113, 355/113, 
75948/24175, 100798/32085, 103993/33102, 312689/99532}

In[6]:= b = 
 Table[characterCount[candidates[[k]]], {k, 1, 1001, 50}]

Out[6]= {4, 54, 103, 154, 203, 253, 303, 353, 404, 455, 504,
 553, 603, 652, 704, 753, 803, 853, 904, 953, 1003}

In[7]:= c = 
 Table[matchingCharacters[candidates[[k]], Pi], {k, 1, 1001, 50}]

Out[7]= {4, 53, 103, 153, 204, 253, 303, 353, 404, 453, 503, 
553, 604, 653, 704, 753, 803, 853, 903, 953, 1003}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

ディオファントス方程式 $x^3+y^3+z^3=n$（２）

（ディオファントス方程式 $x^3+y^3+z^3=n$（１）からの続き）

というわけで，今日，３月14日は円周率の日なのかホワイトデイなのかアインシュタインの誕生日なのかホーキングの命日なのか巷では協議が続いている。

x^3+y^3+z^3=n が n=9k+4, 9k+5 では解を持たないことは簡単に示すことができる。xを9k+mとして９で割った余りをmとする。m=0,1,... 8である。このときx^3を９で割った余りは，m^3を９で割った余りになるので，(0,1,8,27,64,125,216,343,512) ≡ (0,1,8,0,1,8,0,1,8)。

そこで，この剰余(0,1,8)の３つの和のすべての組み合わせを９で割った余りの可能性を調べれば良い。(0,0,0)≡0, (0,0,1)≡1, (0,0,8)≡8, (0,1,1)≡2, (0,1,8)≡0, (0,8,8)≡7, (1,1,1)≡3, (1,1,8)≡1, (1,8,8)≡8, (8,8,8)≡6であり，ここから4と5は出てこない。

少しだけ前回のプログラムを高速化してみる。k^3の配列の次元をmとして，前回のものはO(m^3)，今回のものはO(m^2)である。

（注）Juliaは0の配列添字を許さないので，この度もx,y,zから0は除いて考えている。0を含めれば，91=0^3+3^3+4^3 など簡単なものがあると思うけれど。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function pre(m,t)
  for l in 1:m
    t[l]=l^3
  end
end

function cinv(x)
  return Int(ceil(cbrt(abs(x))))
end

function dpe(n,m,t)
  l=0
  if(n%9 ==4 || n%9 ==5)
    return (n, n%9, 0)
  end
  for i in 1:m
    d=cinv(n-t[i])
    if(n>=t[i])
      for j in 1:d
        k=cinv(n-t[i]-t[j])
        l=l+1
        if k !=0 && (t[i]+t[j]+t[k]==n || t[i]+t[j]+t[k]==-n)
          return (n,l)
        end
      end
      for j in d+1:i
        k=cinv(n-t[i]-t[j])
        l=l+1
        if k !=0 && (t[i]+t[j]-t[k]==n || t[i]+t[j]-t[k]==-n)
          return (n,l)
        end
      end
    else
      for j in d+1:i
        k=cinv(t[i]-n-t[j])
        l=l+1
        if k !=0 && (t[i]-t[j]+t[k]==n || t[i]-t[j]+t[k]==-n)
          return (n,l)
        end
      end
      for j in 1:d
        k=cinv(t[i]-n-t[j])
        l=l+1
        if k !=0 && (t[i]-t[j]-t[k]==n || t[i]-t[j]-t[k]==-n)
          return (n,l)
        end
      end
    end
  end
  return (m,l)
end

function pst(l,m,t)
  for n in 1:l
    println(dpe(n,m,t))
  end
end

m=10000
t=ones(Int64,m)
pre(m,t)
@time pst(99,m,t)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1, 1)
(2, 26)
(3, 1)
(4, 4, 0)
(5, 5, 0)
(6, 1)
(7, 5491)
(8, 1)
(9, 23487)
(10, 1)
(11, 4)
(12, 61)
(13, 4, 0)
(14, 5, 0)
(15, 1057)
(16, 100)
(17, 1)
(18, 4)
(19, 185)
(20, 1561)
(21, 131)
(22, 4, 0)
(23, 5, 0)
(24, 2)
(25, 2)
(26, 48674)
(27, 1)
(28, 146)
(29, 1)
(10000, 50004997)
(31, 4, 0)
(32, 5, 0)
(10000, 50004997)
(34, 10)
(35, 96)
(36, 1)
(37, 1574)
(38, 364)
(10000, 50004997)
(40, 4, 0)
(41, 5, 0)
(10000, 50004997)
(43, 2)
(44, 30)
(45, 152329)
(46, 422)
(47, 31)
(48, 6)
(49, 4, 0)
(50, 5, 0)
(51, 317009)
(10000, 50004997)
(53, 9)
(54, 70)
(55, 1)
(56, 240)
(57, 726)
(58, 4, 0)
(59, 5, 0)
(60, 9)
(61, 466761)
(62, 2)
(63, 23)
(64, 5)
(65, 6188)
(66, 5)
(67, 4, 0)
(68, 5, 0)
(69, 342)
(70, 219)
(71, 12)
(72, 50)
(73, 5)
(10000, 50004998)
(10000, 50004998)
(76, 4, 0)
(77, 5, 0)
(78, 1509)
(79, 612)
(80, 6)
(81, 3)
(82, 96)
(83, 11)
(10000, 50004994)
(85, 4, 0)
(86, 5, 0)
(87, 9120551)
(88, 13)
(89, 884141)
(90, 7)
(91, 72576)
(92, 1)
(93, 20)
(94, 4, 0)
(95, 5, 0)
(96, 239)
(97, 6)
(98, 108)
(99, 2)
 14.962059 seconds (103.76 k allocations: 4.814 MiB)

ディオファントス方程式 $x^3+y^3+z^3=n$（１）

昨日，ディオファントス方程式，$x^3+y^3+z^3=n$ の $n=33$ の解が見つかったというニュースが飛び込んできた。ここで，$n$は自然数，$x, y, z$ は整数である。例えば，$(x,y,z)=(1,1,-1)$ ならば　$n=1$であり，$(x,y,z)=(7,-6,-5)$ならば　$n=2$であり，$(x,y,z)=(1, 1, 1)$ならば　$n=3$である。
\begin{equation}
33= (8,866,128,975,287,528)^3 \\
+(-8,778,405,442,862,239)^3 \\
+(-2,736,111,468,807,040)^3
\end{equation}
20年ほどまえの，KOYAMA, TSURUOKA & SEKIGAWA の論文では，nが1000以下では，

n=  30, 33, 42, 52, 74, 75, 110, 114, 156, 165, 195, 290, 318, 366, 390, 420, 435, 444, 452, 501, 530, 534, 564, 579, 588, 600, 606, 609, 618, 627, 633, 732, 735, 758, 767, 786, 789, 795, 830, 834, 861, 894, 903, 906, 912, 921, 933, 948, 964, 969, 975.

が未発見となっている。

簡単なJuliaのプログラムをつくってみた。$1^3$から$1000^3$までの組み合わせでつくれない２桁以下の$n$は，30,33,39,42,52,74,75,84,87の９つである。『数学者の密室（三島久典さん）』のページにはさらに詳しい情報があり，$0 \le n \le 99$, not equal 4 or 5 (mod 9) , $0 \le |x| \le |y| \le |z| \le 10^{10}$の範囲で，未解決なものはその記事の執筆時点で 33, 42, 74となっていた。

つまり，
30=(-283059965)^3+(-2218888517)^3+(2220422932)^3
33= ? → 上述のとおり解決
39=(117367)^3+(134476)^3+(-159380)^3
42= ?
52=(23961292454)^3+(60702901317)^3+(-61922712865)^3
74= ?
75=(4381159)^3+(435203083)^3+(-435203231)^3
84=(-8241191)^3+(-41531726)^3+(41639611)^3
87=(-1972)^3+(-4126)^3+(4271)^3
である。

（注）オンライン整数列大事典では，このあたり（A173515）。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c=zeros(Int64,100)

function dpt(n,c)
  a=ones(Int128,n)
  for m in 1:n
   a[m]=m^3
  end
  for i in 1:n, j in -n:n, k in j:n
    if(i*j*k != 0 && j+k !=0 && i+j !=0 && i+k !=0)
      b=a[i]+sign(j)*a[abs(j)]+sign(k)*a[abs(k)]
      if(b>0 && b<100)
        c[b]=b
      end
    end
  end
end

for m in 10:10
  println(m)
  @time dft(c,100*m)
  println(c,count(x->x==0,c))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
10

  8.256924 seconds (38.49 k allocations: 1.883 MiB)
[1, 2, 3, 0, 0, 6, 7, 8, 9, 10, 11, 12, 0, 0, 15, 16, 17, 18, 19, 20, 21, 0, 0, 24, 25, 26, 27, 28, 29, 0, 0, 0, 0, 34, 35, 36, 37, 38, 0, 0, 0, 0, 43, 44, 45, 46, 47, 48, 0, 0, 51, 0, 53, 54, 55, 56, 57, 0, 0, 60, 61, 62, 63, 64, 65, 66, 0, 0, 69, 70, 71, 72, 73, 0, 0, 0, 0, 78, 79, 80, 81, 82, 83, 0, 0, 0, 0, 88, 89, 90, 91, 92, 93, 0, 0, 96, 97, 98, 99, 0]32

（ディオファントス方程式 $x^3+y^3+z^3=n$（２）に続く）

World Wide Web

本日2019年3月12日のグーグルは，ワールド・ワイド・ウェブ(WWW)誕生30周年のデザインになっていた。ワールド・ワイド・ウェブの歴史が平成の歴史とほぼ重なっているのか。1989年の3月12日に欧州原子核研究機構（CERN=Conseil Europeen pour la Recherche Nucleaire）のティム・バーナーズ＝リーが，情報管理システムの提案をしたのが出発点となっている。1990年の12月には，NeXTコンピュータ上にシステムを構築し，1991年には，世界最初のウェブサイト http://info.cern.ch/  が公開されている。

その2年後の1993年に，画像とテキストが統合されたブラウザ，NCSAモザイクが登場する。そのシステムを開発したマーク・アンドリーセンはNetscape社を立ち上げ，しばらくは，ネットスケープ・ナビゲータが標準ブラウザの位置を占めていた。まもなく，それはMicrosoftのインターネット・エクスプローラに取って代わられてしまった。

ウェブサーバのソフトとしては，ティム・バーナーズ＝リーによって最初に作られたCERN-httpd，NCSAモザイクに対応した NCSA HTTPdなどがあった。20年ほど前に，最初にサーバにインストールしたのはNCSA HTTPdだったが，これが，Apache HTTP Serverに引き継がれたので切り替えた。10年前ごろには忙しくなったのか歳を取ったせいか，サーバーのお守りもできなくなってしまった。

このあたりの歴史は，TheWebKANZAKIのウェブブラウザ小史2001に詳しい。

写真：NCSA MOSAIC beta（ウェブアーカイブから）

第二種ソフトウェア危機とマクロユーザサービス

標題の文書が発掘された。たぶん，1993年前後のものではないかと思われる。大阪大学大型計算機センターの何かによせて寄稿することを想定した文書だが（利用相談員の自己紹介かな），どこかの段階でボツになったのかもしれない。とりあえず，参考のためにテキストを起こして再現してみる。第二種というのは，研究テーマだった原子核の弱い相互作用に存在するがどうかが問題となっていたＧパリティ非保存の「第二種カレント」にかけているのだが，そんなもの誰にも分からない。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
第二種ソフトウェア危機とマクロユーザサービス
−ユーザのオタク化による情報環境破壊の進行と対策−
大阪教育大学　物理学教室　越桐國雄　(メールアドレス略)

　自己紹介：　私の専門分野は原子核理論（電弱相互作用と核構造）です。センターのACOSをTSSで使い始めたのは1980年ごろで，それまでは阪大核物理研究センターのTOSBACを使っていました。最近は主にワークステーションを使っています。スーパーコンピュータでUNIXが使えれば，また変わるかもしれません。相談分野は次のとおりです。
機種：ACOS，FACOM，NeXTstation，SPARCstation，Macintosh，PC-9801，FM-TOWNS．言語：FORTRAN，Mathematica．
（#ただ，どれもほとんどあてになりませんのでよろしく。）

　大阪教育大の紹介：　大教大では柏原市への移転統合を機会に，情報処理センターが新たに設置されました。ホストはFACOM-M770で，ワークステーション30台，パーソナルコンピュータ50台等の構成です。各研究室からはTCP/IPで学内LANにアクセスします。N1ネットワークへは64kbpsの専用回線で接続し，SINETによる学外へのIP接続も予定しています。しかしあまりの急激な変化（これまではホスト無し，TSS端末数台とMacII40台等があっただけ）に戸惑っているような次第です。（#ほんとにもうたいへんなのです。）

　センターへの希望：　ダウンサイジングとネットワークング及びフリーソフトウェアの浸透により，スタッフのそろわない中小規模大学や情報系以外の研究室では，ユーザのシステム管理やソフトウェア管理の負担が年々重くなっています。また，ネットワークの発展に伴うコンピュータのメディア化が我々の情報処理能力の限界を越えて進もうとしています（ニュースとメールを読み，ファイル転送してmakeすると１日が終わってしまう!?）。こうしてあふれた情報ゴミがユーザのオタク化を促進しています。いわゆる「ソフトウェア危機」＝ソフトウェア生産場面での需給ギャップに対して，これらはソフトウェア消費場面での需給ギャップといえるかもしれません。これを，「第二種ソフトウェア危機」と呼ぶことにします。第二種ソフトウェア危機は人間の基本的な欲求に密接に関係しているため，より本質的で深刻です。これを回避するための一助として，研究室，教室あるいは小規模大学のユーザグループ（＝マクロユーザ）に対する適切な情報提供，事例紹介や講習等のサポートをセンターにお願いしたいと考えています。（#でも難しいですよね。中途半端な結論で終わってしまった・・・）
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

(注) FM-TOWNSは，日本で最初に発売されたCD-ROMドライブ付きパソコンだった。情報処理センターの副センター長だった定金晃三先生らと，これいいよね（デジタル教材メディア活用のはしりだぁ），とかいいながらネットワーク整備予算を使って一教室分導入したものである。そのうちビジネス用のPCにもCD-ROMがつくのが当たり前になってしまった。そして，これからはCD-ROMドライブが付かないのものが普通になるようだ。

JuliaのDoubleFloats

Juliaの任意精度実数としてはBigFloatがあるが，拡張精度として，DoubleFloatsのモジュールを見つけた。Double64はBigFloatより高速であるという説明があったので，ためしてみると，sqrtにバグがあって動かないような気がする。どうしたものか。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
using DoubleFloats
z=Double64(0.5)
#println((sqrt(z),cbrt(z)))
println((sin(z),cos(z),tan(z)))
#println((asin(z),acos(z),atan(z)))
println((csc(z),sec(z),cot(z)))
#println((acsc(1/z),asec(1/z),acot(1/z)))
println((sinh(z),cosh(z),tanh(z)))
#println((asinh(z),acosh(z),atanh(z)))
println((log(z),exp(z),log2(z),log10(z)))
sqrt(z)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(0.479425538604203, 0.8775825618903728, 0.5463024898437905)
(2.085829642933488, 1.139493927324549, 1.830487721712452)
(0.5210953054937474, 1.1276259652063807, 0.46211715726000974)
(-0.6931471805599453, 1.6487212707001282, -1.0, -0.3010299956639812)

UndefVarError: v not defined

Stacktrace:
 [1] mul_ddfp_dd at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_ddfp_dd.jl:16 [inlined]
 [2] sqrt_dd_dd(::Tuple{Float64,Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_dd_dd.jl:74
 [3] sqrt_db_db at /Users/koshi/.julia/packages/DoubleFloats/obU7N/src/math/ops/op_db_db.jl:26 [inlined]
 [4] sqrt(::DoubleFloat{Float64}) at /Users/xxxxx/.julia/packages/DoubleFloats/obU7N/src/math/ops/arith.jl:9
 [5] top-level scope at In[1]:11

Juliaの配列添字範囲

Juliaでは配列の添字が１から始まっている。Cなどでは配列の添字が０から始まっている。Fortranでは，配列の添字範囲を指定できた。そこで，Juliaでも配列の添字範囲を指定したい。

調べてみると，OffsetArrays.jl というモジュールがあった。
using Pkg
Pkg.add("OffsetArrays")
using OffsetArrays
とすればよい。

使用例は以下の通り。昔作ったFortranプログラムをJuliaにしたもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
setprecision(BigFloat,2^6)
o=BigFloat(0)
e=BigFloat(1)
d=BigFloat(2)
s=e/d

N = 300
fc = OffsetArray{BigFloat}(undef,-N:N);
fd = OffsetArray{BigFloat}(undef,-N:N);
fj = OffsetArray{BigFloat}(undef,-N:N);
bn = OffsetArray{BigFloat}(undef,-N:N,-N:N);

function factorials(fc,fd,fj,bn,n)
# k! = k*(k-1)*...*2*1
# fc[k] = sqrt(k!); fc[-k] = k!
# k!! = k*(k-2)*...*3*1 or k*(k-2)*...*4*2
# fd[k] = sqrt((2*k-1)!!/2^k); fd[-k] = fd[k]^2
# fj[k] = sqrt(k+1); fj[-k] = 1/sqrt(k+1)
#
  fc[0]=e
  fd[0]=e
  fj[0]=e
  for k in 1:n
    kk=BigFloat(k)
    fc[k]=sqrt(kk)*fc[k-1]
    fc[-k]=kk*fc[-1+k]
    fd[k]=sqrt(kk-s)*fd[k-1]
    fd[-k]=(kk-s)*fd[-1+k]
    fj[k]=sqrt(kk+e)
    fj[-k]=e/sqrt(kk+e)
  end
  bn[0,0]=e
  for k in 1:n
    bn[k,0]=e
    bn[0,k]=e
    bn[k,k]=e 
    pd=e
      for l in 1:div(k,2)
       j=k-l
       pd=pd*BigFloat(j+1)/BigFloat(l) 
       sd=sqrt(pd)
       bn[k,l]=sd 
       bn[k,j]=sd 
       bn[l,k]=pd 
       bn[j,k]=pd
     end
   end
  println(fc[n])
  println(fc[-10*div(n,10)])
  println((bn[45,90]))
end

@time factorials(fc,fd,fj,bn,300)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1.74944994845934542383e+307

3.03013619603034000008e+308

1.0382742128755341137e+26

  0.732669 seconds (1.29 M allocations: 58.421 MiB, 10.08% gc time)

退職挨拶の練習（１）

理科教育講座の越桐です。1982年に本学に着任してこの３月末で37年になります。

本学に在席していた時間を２つに分けると，20世紀が19年，21世紀が18年ということで，20世紀のほうが長い旧世紀型人間のひとりなので，最近の社会や大学の変化の速さにはとてもついていけない今日この頃です。

さて，退職の挨拶とは，これまで言えなかったことを告白する懺悔の時間のようなものなので，１つだけ紹介させてください。昔の天王寺分校にはデータステーションという，情報処理センターの前身がありました。雨の日には雨漏りするのでコンピュータにシートをかけていました。２階の私の研究室から徒歩30秒ところです。

インターネットの普及の直前の1991年に，本学ははじめてインターネットにつながりました。当時，キャンパスネットワークはないので，利用するにはデータステーションまで行って，コンピュータを使う必要がありました。

そこで，データステーションの垣本先生と結託して，キャンパスネットワークの実験をしようということで，日本橋（にっぽんばし）までいって，私費で10BASE2のケーブルを30メートル買い，研究室のNeXTステーションにヤミで接続したのは私です。ごめんなさい。これが，本学におけるキャンパスＬＡＮのはしりでした。

今では無線LANがどこでも使える夢の世界が実現したのですが，いいのか悪いのか・・・

そんなわけで，先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん仲たがいせずにこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。

Juliaの数学的定数と精度

（ネイピア数の計算からの続き）

円周率やネイピア数の計算で，JuliaのBigFloatの精度が$10^{-77}$であったが，これは，
$2^{-255}$に対応していた。すなわち，BigFloatの精度は次のようになる。
\begin{equation}
{\rm eps(BigFloat)}=2^{-255}=1.7272337110188889250772703725600\\
79914223200072887256277004740694033718360632485\ {\rm e}^{-77}
\end{equation}
円周率πやネイピア数eの有効桁数がこれによって決まるが，BigFloatの精度は制御することができる。このためには，setprecision(BigFloat,n) とすればよい。デフォルト値はn=256に対応しており，m = log10(BigInt(2)^(-n+1))  から，10進法における有効桁数mが求まることになる。

その様子は次のプログラムでわかる。Juliaの数学的定数として，円周率πやネイピア数ℯに加えて，カタラン数，オイラーのガンマγ，黄金比φなどが Base.MathConstants に準備されている。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cnst(n)
  setprecision(BigFloat,2^n)
  println(" ")
  println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
  println(2^n,"　eps=",eps(BigFloat))
  println("- - - - - - - - - - - - - - - - - - - - - - - - - - - -")
  println("　π=",BigFloat(Base.MathConstants.π))
  println("　 =",BigFloat(Base.MathConstants.ℯ))
  println("　 c=",BigFloat(Base.MathConstants.catalan))
  println("　γ=",BigFloat(Base.MathConstants.γ))
  println("　φ=",BigFloat(Base.MathConstants.φ))
end

for k in 4:2:10
 @time cnst(k)
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

- - - - - - - - - -
16　eps=3.05176e-05
- - - - - - - - - -
　π=3.1416
　 ℯ=2.71826
　 c=0.91597
　γ=0.577209
　φ=1.61804
  0.000738 seconds (475 allocations: 36.148 KiB)

- - - - - - - - - -
64　eps=1.08420217248550443401e-19
- - - - - - - - - -
　π=3.14159265358979323851
　 ℯ=2.71828182845904523543
　 c=0.915965594177219015048
　γ=0.577215664901532860616
　φ=1.61803398874989484821
  0.000528 seconds (471 allocations: 37.184 KiB)

- - - - - - - - - -
256　eps=1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77
- - - - - - - - - -

π=3.141592653589793238462643383279502884197169399375105820974944592307816406286198

ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353555

c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762584

γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776685

φ=1.61803398874989484820458683436563811772030917980576286213544862270526046281891
  0.000546 seconds (471 allocations: 40.613 KiB)

- - - - - - - - - -
1024　eps=1.112536929253600691545116358666202032109607990231165915276663708443602217406959097927141579506255510282033669865517905502576217080776730054428006192688859410565388996766001165239805073721291818035960782523471251867104187625403325308329079474360245589984295819824250317954385059152437399890443876874974725790226e-308
- - - - - - - - - -

π=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724586997

ℯ=2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082264800168477411853742345442437107539077744992069551702766

c=0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062485744226191961995790358988033258590594315947374811584069953320287733194605190387274781640878659090247064841521630002287276409423882599577415088163974702524820115607076448838078733704899008647751132259971343386

γ=0.5772156649015328606065120900824024310421593359399235988057672348848677267776646709369470632917467495146314472498070824809605040144865428362241739976449235362535003337429373377376739427925952582470949160087352039481656708532331517766115286211995015079847937450857057400299213547861466940296043254215190587755369

φ=1.618033988749894848204586834365638117720309179805762862135448622705260462818902449707207204189391137484754088075386891752126633862223536931793180060766726354433389086595939582905638322661319928290267880675208766892501711696207032221043216269548626296313614438149758701220340805887954454749246185695364864449242
  0.000710 seconds (481 allocations: 56.924 KiB)

ネイピア数の計算

（円周率の計算からの続き）

ネイピア数は自然対数の底 $e$ のこと。欧米では Euler's Number とよばれることが多いとあるが，Eulerが関わる数はたくさんあってややこしいので，ここではネイピア数とする。

円周率の計算方法というか表現方法はたくさんあって，その中には与えられた級数の少数の項で良い精度の近似を与えるものがあった。ところで，ネイピア数の表式をざっと調べても，なんだかぱっとしないのである。
\begin{equation}
\begin{aligned}
e^x &= 1 + x + \frac{x^2}{2!}+\frac{x^3}{3!}+\frac{x^4}{4!}+\cdots\\
e &= 1 + 1 +  \frac{1}{2!}+\frac{1}{3!}+\frac{1}{4!}+\cdots
\end{aligned}
\end{equation}
を多少変形したものか，連分数で表現したものはあるけれど，ラマヌジャンの公式のようにはっとするものが見当たらない。単に調査不足だけかもしれないし，そもそも指数関数を微分や積分しても元に戻るだけなのが原因かもしれないが，よくわからないのだった。

５項くらいで精度が70桁程度になる公式はないものなのだろうか？Improving the Convergence of Newton's Series Approximation for $e$ とかあったけど，目的のものではない。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)

function fact(f,n)
# calculation of factorial 1!..n!
  for i in 2:n
    f[i]=BigFloat(i)*f[i-1]
  end
# println(f[n])
end

function cf(x,a,b,n)
# continued fraction
  for k = n:-1:1
    x = a[k] + b[k]/x
  end
  return x
end

function napier1(n)
# by series expansion
  o=BigFloat(1)
  sum=o
  for k in 1:n
    sum=sum+o/f[k]
  end
  println(sum-exp(o))
end

function napier2(n)
# by series expansion
  o=BigFloat(1)
  sum=o
  for k in 1:n
    sum=sum+o/((BigFloat(4)^k)*f[k])
  end
  println(sum*sum*sum*sum-exp(o))
end

function napier3(n)
# by continued fraction
  o=BigFloat(1)
  z=BigFloat(2)
  a=ones(BigFloat,1000)
  for k in 1:1000
    a[k]=z*(z*k-z-o)
  end
  a[1]=o
  a[2]=o
  b=ones(BigFloat,1000)
  b[1]=z
  println(cf(z,a,b,n)-exp(o))
end

function napier4(n)
# sqrt(e) by series expansion
  o=BigFloat(1)
  z=BigFloat(2)
  sum=BigFloat(0)
  for k in 0:n
    sum=sum+BigFloat(4*k+3)/(z^(2*k+1)*f[2*k+1])
  end
  println(sum*sum-exp(o))
end

fact(f,160)
@time napier1(57)
@time napier2(41)
@time napier3(26)
@time napier4(24)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
  0.018460 seconds (30.73 k allocations: 1.488 MiB)
-2.763573937630222280123632596096127862757120116619610043207585110453949377011977e-76
  0.029598 seconds (49.70 k allocations: 2.358 MiB)
-1.381786968815111140061816298048063931378560058309805021603792555226974688505988e-76
  0.054272 seconds (100.99 k allocations: 5.090 MiB)
1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-76
  0.032685 seconds (72.55 k allocations: 3.431 MiB)

（Juliaの数学的定数と精度へ続く）

円周率の計算

（電子計算機と人間からの続き）

現代科学シリーズの「電子計算機と人間」では，円周率の計算のためのFORTRANプログラムが書かれていたが，高校１年生の自分にはなぜそれで円周率が求まるのかがさっぱり分からなかった。50年ぶりに，図書館でその本を手にとると，たしかに分かりやすくはなかったが，アークタンジェントのマクローリン展開に１を代入した単純な公式なのであった。
\begin{equation}
\begin{aligned}
\arctan(x) & = x - \frac{x^3}{3!} +  \frac{x^5}{5!} - \frac{x^7}{7!} + \cdots \\
\frac{\pi}{4} & =  1 - \frac{1}{3!} +  \frac{1}{5!} - \frac{1}{7!} + \cdots
\end{aligned}
\end{equation}
この公式は収束がとても悪い見本のような式なので，1000項でようやく3.14が0.1%精度で求まる始末である。

そこで，いくつかの代表的な公式をJuliaプログラムに落としてみた。BigFloatの精度がeps(BigFloat) = $10^{-77}$程度なので，有効数字は77桁弱である。そのBigFloatの使い方が，いまひとつわかっておらず適当に入れているが，たぶん冗長だと思われる。

円周率の計算で，有効数字77桁程度の精度を出すには，昔よくみかけたマチンの公式で 54項，なんだかすごいのだけれど今回初めて使った ラマヌジャンの公式で ９項，スマートなガウス＝ルジャンドルの方法で５項が必要であった。もっと楽しそうなJuliaでのπの話はこちら。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
f=ones(BigFloat,1000)

function fact(f,n)
# calculation of factorial 1!..n!
  for i in 2:n
    f[i]=BigFloat(i)*f[i-1]
  end
# println(f[n])
end

function pie1(n)
# Gregory-Leibniz Formula
# pi=4*arctan(1)
  sum=BigFloat(0)
  for i in 1:n
    sum=sum+(-1)^(i-1)/BigFloat(2*i-1)
  end
  println(4*sum-pi)
end

function pie2(n)
# Machin's Formula
# pi = 16*arctan(1/5)-4*arctan(1/239)
  sum=BigFloat(0)
  sun=BigFloat(0)
  for i in 1:n
    sum=sum+(-1)^(i-1)/((2*i-1)*BigFloat(5)^(2*i-1))
    sun=sun+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
  end
  println(16*sum-4*sun-pi)
end

function pie3(n)
# Takano-Kikuo's Formula
# pi = 48*arctan(1/49)+128*arctan(1/57)-30*arctan(1/239)+48*arctan(1/110443)
  sum1=BigFloat(0)
  sum2=BigFloat(0)
  sum3=BigFloat(0)
  sum4=BigFloat(0)
  for i in 1:n
    sum1=sum1+(-1)^(i-1)/((2*i-1)*BigFloat(49)^(2*i-1))
    sum2=sum2+(-1)^(i-1)/((2*i-1)*BigFloat(57)^(2*i-1))
    sum3=sum3+(-1)^(i-1)/((2*i-1)*BigFloat(239)^(2*i-1))
    sum4=sum4+(-1)^(i-1)/((2*i-1)*BigFloat(110443)^(2*i-1))        
  end
  println(48*sum1+128*sum2-20*sum3+48*sum4-pi)
end

function pie4(n)
# Ramanujan Formula
  sum=BigFloat(1103)
  for k in 1:n
    sum=sum+f[4*k]/(f[k]^4*BigFloat(396)^BigFloat(4*k))*BigFloat(1103+26390*k)
  end
  println(BigFloat(9801)/(sqrt(BigFloat(8))*sum)-pi)
end

function pie5(n)
# Chudnovsky formula 
  sum=BigFloat(13591409)
  for k in 1:n
    sum=sum+f[6*k]/(f[3*k]*f[k]^3*BigFloat(-640320)^BigFloat(3*k))*BigFloat(13591409+545140134*k) 
  end
  println(BigFloat(640320)^1.5/(BigFloat(12)*sum)-pi)
end

function pie6(n)
# Gauss-Legendre algorithm
  z=BigFloat(2)
  a=BigFloat(1)
  b=1/sqrt(z)
  t=1/(z*z)
  p=BigFloat(1)
  for k in 1:n
    o=a
    a=(o+b)/z
    b=sqrt(o*b)
    t=t-p*(a-o)*(a-o)
    p=z*p
  end
  println((a+b)*(a+b)/(4*t)-pi)
end

fact(f,160)
@time pie1(10000000)
@time pie2(54)
@time pie3(22)
@time pie4(9)
@time pie5(5)
@time pie6(5)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

-9.999999999999975000000000000312499999999990468750000000541015625007904590936847e-08
  5.580711 seconds (60.04 M allocations: 3.131 GiB, 4.69% gc time)
3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
  0.055656 seconds (76.18 k allocations: 3.601 MiB, 7.32% gc time)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
  0.071662 seconds (123.15 k allocations: 5.740 MiB)
-6.908934844075555700309081490240319656892800291549025108018962776134873442529942e-77
  0.057192 seconds (93.11 k allocations: 4.357 MiB)
-3.454467422037777850154540745120159828446400145774512554009481388067436721264971e-77
  0.068310 seconds (106.84 k allocations: 5.084 MiB)
0.0
  0.036871 seconds (53.77 k allocations: 2.544 MiB)

（ネイピア数の計算に続く）

若冲の樹花鳥獣図屏風

（旧三井家下鴨別邸からの続き）

旧三井家下鴨別邸の望楼特別公開にあわせて，２階では伊藤若冲の樹花鳥獣図屏風の高精細複写品（キャノン）が展示されていた。本物は静岡県立美術館にあって，これは実際にみたことがある。たぶん平面にして展示してあったと思うが，今回は屏風立てで，床に置いたロウソク風のLED照明をあてている。以前の印象より小振りに見え，表現もすっきりしているように感じた。

写真：若冲の樹花鳥獣図屏風の高精細複写（2019.3.4）

それは，プライスコレクション「鳥獣花木図屏風」のイメージと比較したからだろうか。2006年に「若冲と江戸絵画」という展覧会があって，愛知県美術館で「鳥獣花木図屏風」をはじめてみた。白象の背中に敷物があるほうだ。これがこの屏風を見た初めての機会であったがとても印象深かった。最終日に近く，ちょうどジョー・プライス夫妻も来ていてので図録にサインをもらうことができた。静岡県立美術館の方はその数年後に訪れたが，設置の位置のせいか，プライスコレクションのそれより小さく感じたのだった。

「樹花鳥獣図屏風（静岡）」と「鳥獣花木図屏風（プライス）」については，佐藤康宏と辻惟雄の間に真贋論争がある。

旧三井家下鴨別邸

京都の下鴨神社の南にある旧三井家下鴨別邸で，平成31年２月７日から３月19日まで，主屋の三階望楼の特別公開をしていたので夫婦で行ってみた。毎年この時期に一ヶ月程度公開しているらしい。明治の末に三井家の祖霊社がこのあたりに移されたので，大正の末に木屋町にあった別邸を休憩所として移築したそうだ。戦後は国に譲渡され，家庭裁判所の所長宿舎として使われた後，近代和風建築として重要文化財に指定された。一帯の敷地は八千坪にもなるが，一部は京都家庭裁判所となっている。

あいにくの雨模様だったが，到着すると順番に十名程度ずつ登ってもらいますとのことだった。二階の待合室で説明を聞いてから，それほど待たずに登ることができた。昔は五山の送り火もすべて一望できたようだが，いまは，左大文字と妙法の一部がみえるにとどまる。南側の木立を整理すれば，鴨川と高野川の合流地点がきれいに見えるはずなのに惜しい。

三井家から広岡家に嫁入りした広岡あさをモデルにした「あさが来た」が2015年度下半期のNHKの朝ドラとして放映され，その縁の品なども展示されていた。この建物が当時の撮影にも登場したそうで，放送当時は表まで行列ができるほどのにぎわいだったそうだ。

二階から急な階段の脇に張った太い組み紐につかまりながら，恐る恐る登っていくと三畳ほどの四方が開けた望楼になる。雨戸が下から引き出す方式が珍しいという説明を受けた。プライバシーの問題があるので，三階の望楼からの撮影は禁止である。

十分ほど説明があって，降りた二階には，伊藤若冲の「樹花鳥獣図屏風」の高精彩複製品（キャノン）が展示されていた。一階にはしゃれたデザインの洗面所や風呂場などがあり，瓢箪池のあり庭園にでて散策したが，そこからは建物の全体が一望できた。

 

写真：三井下鴨別邸（2019.3.3）

（若冲の樹花鳥獣図屏風に続く）