円周率の日の余震が続いている。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}
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2019年3月14日木曜日
ディオファントス方程式 $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 など簡単なものがあると思うけれど。
（注）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)
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)
2019年3月13日水曜日
ディオファントス方程式 $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
（ディオファントス方程式 $x^3+y^3+z^3=n$（２）に続く）
\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
（ディオファントス方程式 $x^3+y^3+z^3=n$（２）に続く）
2019年3月12日火曜日
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に詳しい。
その2年後の1993年に，画像とテキストが統合されたブラウザ，NCSAモザイクが登場する。そのシステムを開発したマーク・アンドリーセンはNetscape社を立ち上げ，しばらくは，ネットスケープ・ナビゲータが標準ブラウザの位置を占めていた。まもなく，それはMicrosoftのインターネット・エクスプローラに取って代わられてしまった。
ウェブサーバのソフトとしては，ティム・バーナーズ＝リーによって最初に作られたCERN-httpd，NCSAモザイクに対応した NCSA HTTPdなどがあった。20年ほど前に，最初にサーバにインストールしたのはNCSA HTTPdだったが，これが，Apache HTTP Serverに引き継がれたので切り替えた。10年前ごろには忙しくなったのか歳を取ったせいか，サーバーのお守りもできなくなってしまった。
このあたりの歴史は，TheWebKANZAKIのウェブブラウザ小史2001に詳しい。
写真：NCSA MOSAIC beta（ウェブアーカイブから）
2019年3月11日月曜日
第二種ソフトウェア危機とマクロユーザサービス
標題の文書が発掘された。たぶん，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ドライブが付かないのものが普通になるようだ。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
第二種ソフトウェア危機とマクロユーザサービス
−ユーザのオタク化による情報環境破壊の進行と対策−
大阪教育大学 物理学教室 越桐國雄 (メールアドレス略)
自己紹介： 私の専門分野は原子核理論（電弱相互作用と核構造）です。センターの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ドライブが付かないのものが普通になるようだ。
2019年3月10日日曜日
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
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
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
2019年3月9日土曜日
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)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
調べてみると，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)
2019年3月8日金曜日
退職挨拶の練習（１）
理科教育講座の越桐です。1982年に本学に着任してこの３月末で37年になります。
本学に在席していた時間を２つに分けると，20世紀が19年，21世紀が18年ということで，20世紀のほうが長い旧世紀型人間のひとりなので，最近の社会や大学の変化の速さにはとてもついていけない今日この頃です。
さて，退職の挨拶とは，これまで言えなかったことを告白する懺悔の時間のようなものなので，１つだけ紹介させてください。昔の天王寺分校にはデータステーションという，情報処理センターの前身がありました。雨の日には雨漏りするのでコンピュータにシートをかけていました。２階の私の研究室から徒歩30秒ところです。
インターネットの普及の直前の1991年に，本学ははじめてインターネットにつながりました。当時，キャンパスネットワークはないので，利用するにはデータステーションまで行って，コンピュータを使う必要がありました。
そこで，データステーションの垣本先生と結託して，キャンパスネットワークの実験をしようということで，日本橋（にっぽんばし）までいって，私費で10BASE2のケーブルを30メートル買い，研究室のNeXTステーションにヤミで接続したのは私です。ごめんなさい。これが，本学におけるキャンパスＬＡＮのはしりでした。
今では無線LANがどこでも使える夢の世界が実現したのですが，いいのか悪いのか・・・
そんなわけで，先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん仲たがいせずにこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
本学に在席していた時間を２つに分けると，20世紀が19年，21世紀が18年ということで，20世紀のほうが長い旧世紀型人間のひとりなので，最近の社会や大学の変化の速さにはとてもついていけない今日この頃です。
さて，退職の挨拶とは，これまで言えなかったことを告白する懺悔の時間のようなものなので，１つだけ紹介させてください。昔の天王寺分校にはデータステーションという，情報処理センターの前身がありました。雨の日には雨漏りするのでコンピュータにシートをかけていました。２階の私の研究室から徒歩30秒ところです。
インターネットの普及の直前の1991年に，本学ははじめてインターネットにつながりました。当時，キャンパスネットワークはないので，利用するにはデータステーションまで行って，コンピュータを使う必要がありました。
そこで，データステーションの垣本先生と結託して，キャンパスネットワークの実験をしようということで，日本橋（にっぽんばし）までいって，私費で10BASE2のケーブルを30メートル買い，研究室のNeXTステーションにヤミで接続したのは私です。ごめんなさい。これが，本学におけるキャンパスＬＡＮのはしりでした。
今では無線LANがどこでも使える夢の世界が実現したのですが，いいのか悪いのか・・・
そんなわけで，先生方，職員の皆様方には長らく大変お世話になりました。大学はこれから大変な時代を迎えますが，どうかお身体にお気をつけ，皆さん仲たがいせずにこの難関に立ち向かっていただければと思います。本日はどうもありがとうございました。
2019年3月7日木曜日
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)
円周率やネイピア数の計算で，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)
2019年3月6日水曜日
ネイピア数の計算
（円周率の計算からの続き）
ネイピア数は自然対数の底 $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の数学的定数と精度へ続く）
ネイピア数は自然対数の底 $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の数学的定数と精度へ続く）
2019年3月5日火曜日
円周率の計算
（電子計算機と人間からの続き）
現代科学シリーズの「電子計算機と人間」では，円周率の計算のための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)
（ネイピア数の計算に続く）
現代科学シリーズの「電子計算機と人間」では，円周率の計算のための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)
（ネイピア数の計算に続く）
2019年3月4日月曜日
若冲の樹花鳥獣図屏風
（旧三井家下鴨別邸からの続き）
旧三井家下鴨別邸の望楼特別公開にあわせて，２階では伊藤若冲の樹花鳥獣図屏風の高精細複写品（キャノン）が展示されていた。本物は静岡県立美術館にあって，これは実際にみたことがある。たぶん平面にして展示してあったと思うが，今回は屏風立てで，床に置いたロウソク風のLED照明をあてている。以前の印象より小振りに見え，表現もすっきりしているように感じた。
それは，プライスコレクション「鳥獣花木図屏風」のイメージと比較したからだろうか。2006年に「若冲と江戸絵画」という展覧会があって，愛知県美術館で「鳥獣花木図屏風」をはじめてみた。白象の背中に敷物があるほうだ。これがこの屏風を見た初めての機会であったがとても印象深かった。最終日に近く，ちょうどジョー・プライス夫妻も来ていてので図録にサインをもらうことができた。静岡県立美術館の方はその数年後に訪れたが，設置の位置のせいか，プライスコレクションのそれより小さく感じたのだった。
「樹花鳥獣図屏風（静岡）」と「鳥獣花木図屏風（プライス）」については，佐藤康宏と辻惟雄の間に真贋論争がある。
旧三井家下鴨別邸の望楼特別公開にあわせて，２階では伊藤若冲の樹花鳥獣図屏風の高精細複写品（キャノン）が展示されていた。本物は静岡県立美術館にあって，これは実際にみたことがある。たぶん平面にして展示してあったと思うが，今回は屏風立てで，床に置いたロウソク風のLED照明をあてている。以前の印象より小振りに見え，表現もすっきりしているように感じた。
写真：若冲の樹花鳥獣図屏風の高精細複写（2019.3.4）
それは，プライスコレクション「鳥獣花木図屏風」のイメージと比較したからだろうか。2006年に「若冲と江戸絵画」という展覧会があって，愛知県美術館で「鳥獣花木図屏風」をはじめてみた。白象の背中に敷物があるほうだ。これがこの屏風を見た初めての機会であったがとても印象深かった。最終日に近く，ちょうどジョー・プライス夫妻も来ていてので図録にサインをもらうことができた。静岡県立美術館の方はその数年後に訪れたが，設置の位置のせいか，プライスコレクションのそれより小さく感じたのだった。
「樹花鳥獣図屏風（静岡）」と「鳥獣花木図屏風（プライス）」については，佐藤康宏と辻惟雄の間に真贋論争がある。
2019年3月3日日曜日
旧三井家下鴨別邸
京都の下鴨神社の南にある旧三井家下鴨別邸で，平成31年２月７日から３月19日まで，主屋の三階望楼の特別公開をしていたので夫婦で行ってみた。毎年この時期に一ヶ月程度公開しているらしい。明治の末に三井家の祖霊社がこのあたりに移されたので，大正の末に木屋町にあった別邸を休憩所として移築したそうだ。戦後は国に譲渡され，家庭裁判所の所長宿舎として使われた後，近代和風建築として重要文化財に指定された。一帯の敷地は八千坪にもなるが，一部は京都家庭裁判所となっている。
あいにくの雨模様だったが，到着すると順番に十名程度ずつ登ってもらいますとのことだった。二階の待合室で説明を聞いてから，それほど待たずに登ることができた。昔は五山の送り火もすべて一望できたようだが，いまは，左大文字と妙法の一部がみえるにとどまる。南側の木立を整理すれば，鴨川と高野川の合流地点がきれいに見えるはずなのに惜しい。
三井家から広岡家に嫁入りした広岡あさをモデルにした「あさが来た」が2015年度下半期のNHKの朝ドラとして放映され，その縁の品なども展示されていた。この建物が当時の撮影にも登場したそうで，放送当時は表まで行列ができるほどのにぎわいだったそうだ。
二階から急な階段の脇に張った太い組み紐につかまりながら，恐る恐る登っていくと三畳ほどの四方が開けた望楼になる。雨戸が下から引き出す方式が珍しいという説明を受けた。プライバシーの問題があるので，三階の望楼からの撮影は禁止である。
十分ほど説明があって，降りた二階には，伊藤若冲の「樹花鳥獣図屏風」の高精彩複製品（キャノン）が展示されていた。一階にはしゃれたデザインの洗面所や風呂場などがあり，瓢箪池のあり庭園にでて散策したが，そこからは建物の全体が一望できた。
（若冲の樹花鳥獣図屏風に続く）
あいにくの雨模様だったが，到着すると順番に十名程度ずつ登ってもらいますとのことだった。二階の待合室で説明を聞いてから，それほど待たずに登ることができた。昔は五山の送り火もすべて一望できたようだが，いまは，左大文字と妙法の一部がみえるにとどまる。南側の木立を整理すれば，鴨川と高野川の合流地点がきれいに見えるはずなのに惜しい。
三井家から広岡家に嫁入りした広岡あさをモデルにした「あさが来た」が2015年度下半期のNHKの朝ドラとして放映され，その縁の品なども展示されていた。この建物が当時の撮影にも登場したそうで，放送当時は表まで行列ができるほどのにぎわいだったそうだ。
二階から急な階段の脇に張った太い組み紐につかまりながら，恐る恐る登っていくと三畳ほどの四方が開けた望楼になる。雨戸が下から引き出す方式が珍しいという説明を受けた。プライバシーの問題があるので，三階の望楼からの撮影は禁止である。
十分ほど説明があって，降りた二階には，伊藤若冲の「樹花鳥獣図屏風」の高精彩複製品（キャノン）が展示されていた。一階にはしゃれたデザインの洗面所や風呂場などがあり，瓢箪池のあり庭園にでて散策したが，そこからは建物の全体が一望できた。
写真：三井下鴨別邸（2019.3.3）
（若冲の樹花鳥獣図屏風に続く）
2019年3月2日土曜日
電子計算機と人間
SSS 現代の科学シリーズというのが，1967年から1980年にかけて，河出書房新社から出版されていた。ちょうど中学の終わりから高校にかけて何冊かそろえていた。講談社のブルーバックスよりも少し上といった感じのシリーズである。最初に買ったのが，アイザック・アシモフの「化学の歴史」であり，中学の時に理科クラブで化学に取り組んでいたことに関係があるかもしれない。この本は今では，ちくま学芸文庫に収録されていているので入手しやすい。
もう一冊よく憶えているのがドナルド・フィンクの「電子計算機と人間」である。あと何冊かあったような気がするが，タイトルをみてもどれもピンとこないのである。もしかすると２冊だけしか持っていなかったのかもしれない。
「電子計算機と人間」は，石田晴久先生が最初に手がけられた本のようである。情報処理 Vol.50 No.7（2009年）に「あの時代」に想いをはせて−証言者たちからのメッセージ−として，2009年の３月に亡くなられた石田晴久の追悼特集が載っている。その中に，元bit誌編集長の小山透さんが石田先生の著述物の一覧を整理されている。その第１号が，1969年に出版された，高橋秀俊先生との共訳によるこの本なのだった。
写真：電子計算機と人間（アマゾンより引用）
金沢泉丘高等学校の理数科の１年でプログラミングのさわりを学ぼうとしていた自分にとっては，まさに時宜を得た読書体験であった。ところで，この本には，FORTRANのプログラミングで円周率を求めるという話が延々と書いてあったのだが，どうして円周率が求まるのかがさっぱりわからず，狐につままれたような思いが残っているのだった。いったい，どんなアルゴリズムを用いようとしていたのだろうか。
（円周率の計算に続く）
（円周率の計算に続く）
2019年3月1日金曜日
NHK-FM放送開始50周年
NHKが1969年の３月１日にFM放送を開始してから50周年を迎え，３月１日から３日間にわたって，特別番組を放送している。1969年といえば，高校１年生だった。その頃，親にねだって買ってもらったSONYのテープレコーダサーボマチックＦ１（TC-222）と家にあった古いトランジスタラジオでAM放送のエアチェック（洋楽ポップス＋日本のフォーク少々）をするのが趣味となっていた。その後，半裸ケーブルでイヤホンジャックから結線していたラジオから，SONYのトランジスタラジオ（TFM-110F）に更新してもらったように思う。
ノートに楽曲名や演奏者を記録しながら，SONYの直径５インチ（テープ幅1/4インチ）の家庭用オープンリールテープ10数枚に楽曲がたまっていった。サーボマチックＦ１は，テープスピードが 4.8 cm/sと 9.6cm/s に切り替えられるようになっていたが，テープ代を節約したい高校生は，音質を無視して 低速で録音していた。この場合の１本のテープの録音時間は１時間くらいだったと思う。したがって１本に15-20曲程度入ることになる。10本で200曲か。いまだと，iPhoneに楽々4000曲保存できるのであった。
1972年に大学に入ってしばらくすると，もうオープンリールの時代ではないということで，SONYのラジカセ（CF-1600）に乗り換えた。なんと景気の良い時代であったこと。NHK-FMでは大学の同級生の佐藤秀明君の影響でクラシックを聴くようになり，やはり同級生の佐藤正泰君にならって，FM-fanを何度か購入したこともあった。
写真：SONY TC-222 ヤフオクより引用
写真：SONY TFM-110F ゴールデン横丁より引用
写真：SONY-TAPE 100 オープンリール復刻図鑑より引用
写真：SONY CF-1600 ヤフオクより引用
1972年に大学に入ってしばらくすると，もうオープンリールの時代ではないということで，SONYのラジカセ（CF-1600）に乗り換えた。なんと景気の良い時代であったこと。NHK-FMでは大学の同級生の佐藤秀明君の影響でクラシックを聴くようになり，やはり同級生の佐藤正泰君にならって，FM-fanを何度か購入したこともあった。
2019年2月28日木曜日
小山内宏
ハノイでの第二回米朝首脳会談であるが，ベトナム戦争で米国に勝利して，その後の経済発展が目覚ましいベトナム社会主義共和国が朝鮮民主主義人民共和国のモデルとなるかどうか，ということで話題になっているのか，どうだろう。
1972年にジェーン・フォンダが北ベトナムのハノイを訪れて物議をかもしていたころ，大学に入学したばかりだったが，ベトナム反戦はキャンパスの学生運動における主要なテーマではなかったような気がする。それでも，外部から識者をよんでベトナム戦争の現状について学習するための講演会があった。ベトナム戦争が実質的に終結するのは，1975年4月のサイゴン陥落と南ベトナム政府崩壊の時点である。その２年ほど前だろうか，すでに戦況が北ベトナム有利となりつつあるころに，ベトナム戦争の状況を説明してくれた講師が，軍事評論家の小山内宏先生（小山内薫の次男）だった。
阪大豊中キャンパスのロ号館の中教室に集まった学生はそれほど多くはなく，当時阪大キャンパスの学生運動の実権を握っていた民学同系の団体の主催だった。小山内さんは，黒板にチョークで絵を描きながら，ベトコンの移動・輸送路の網目が如何に北爆の影響を回避しながら活動しているかをわかりやすく説明してくれた。
1972年にジェーン・フォンダが北ベトナムのハノイを訪れて物議をかもしていたころ，大学に入学したばかりだったが，ベトナム反戦はキャンパスの学生運動における主要なテーマではなかったような気がする。それでも，外部から識者をよんでベトナム戦争の現状について学習するための講演会があった。ベトナム戦争が実質的に終結するのは，1975年4月のサイゴン陥落と南ベトナム政府崩壊の時点である。その２年ほど前だろうか，すでに戦況が北ベトナム有利となりつつあるころに，ベトナム戦争の状況を説明してくれた講師が，軍事評論家の小山内宏先生（小山内薫の次男）だった。
阪大豊中キャンパスのロ号館の中教室に集まった学生はそれほど多くはなく，当時阪大キャンパスの学生運動の実権を握っていた民学同系の団体の主催だった。小山内さんは，黒板にチョークで絵を描きながら，ベトコンの移動・輸送路の網目が如何に北爆の影響を回避しながら活動しているかをわかりやすく説明してくれた。
写真：秋田書店「拳銃百科」奥付けより引用
P. S. ここには1924年生まれとあるが，経歴からして，Wikipediaの1916年が正しいのではないか。
2019年2月27日水曜日
Sofitel Legend Metropole Hanoi
今日，明日とベトナムのハノイにある，ソフィテル・レジェンド・メトロポール・ハノイで，第二回の米朝首脳会談が行われている。あいかわらず，Wikipedia日本語版には記事がなく，英語版のSofitel Legend Metropole Hanoiはこちら。
2016年の5月に家族旅行でベトナムのハノイを訪問した。ハロン湾とハノイ市内の観光である。そのとき，後半の２泊ほどをこのホテルに宿泊した。そのときの写真。
いろいろな著名人も宿泊している歴史のある重厚なホテルである。ベトナム戦争が終結する少し前に，ジェーン・フォンダが泊まったときの資料が展示されていた。ベトナム戦争のキーワードによって，脳内でジョーン・バエズに変換して記憶されていたのを今回訂正した。ジェーン・フォンダはＳＦマガジンの解説でドキドキしたあのバーバレラ（ロジェ・ヴァディム監督作品）であった。
2016年の5月に家族旅行でベトナムのハノイを訪問した。ハロン湾とハノイ市内の観光である。そのとき，後半の２泊ほどをこのホテルに宿泊した。そのときの写真。
写真：ソフィテル・レジェンド・メトロポール・ハノイの中庭（2016.5）
2019年2月26日火曜日
折田先生像
昨日今日は，国立大学の個別学力検査，前期試験の日。
ネットでは例年のように京大の折田彦市先生の像が話題になっていた。
折田彦市先生は，旧制第三高等学校の初代校長で，1950年に製作された銅像が今に伝わっている。1991年ごろから落書きの連鎖が始まり，これを避けるために1997年には像が総合人間学部図書館の地下に収納された。しかし，このムーブメントは治まることはなく，台座と像のパロディが受験シーズンなどに毎年繰り返して作られることになった。
ちょうどその1991年から1997年にかけて，総合人間学部の情報科学実習を担当していたので，折田先生の像の変化はしばしば目にしていたのだった。例えば，こんなページや，そんなページや，あんなページ（からのリンク）が見つかる。そして，そのリンクの先には1996年の日本のWWWサーバ一覧（なつかしい）などもある。
P. S. 折田先生の胸像は，現在は，京都大学の百周年時計台記念館１Ｆの歴史展示室に安置され，やすらかな余生を過ごしていらっしゃるようだ。
ネットでは例年のように京大の折田彦市先生の像が話題になっていた。
写真：折田彦市（Wikipediaより）
折田彦市先生は，旧制第三高等学校の初代校長で，1950年に製作された銅像が今に伝わっている。1991年ごろから落書きの連鎖が始まり，これを避けるために1997年には像が総合人間学部図書館の地下に収納された。しかし，このムーブメントは治まることはなく，台座と像のパロディが受験シーズンなどに毎年繰り返して作られることになった。
ちょうどその1991年から1997年にかけて，総合人間学部の情報科学実習を担当していたので，折田先生の像の変化はしばしば目にしていたのだった。例えば，こんなページや，そんなページや，あんなページ（からのリンク）が見つかる。そして，そのリンクの先には1996年の日本のWWWサーバ一覧（なつかしい）などもある。
P. S. 折田先生の胸像は，現在は，京都大学の百周年時計台記念館１Ｆの歴史展示室に安置され，やすらかな余生を過ごしていらっしゃるようだ。
2019年2月25日月曜日
辺野古米軍基地建設のための埋め立ての賛否（２）
辺野古米軍基地建設のための埋め立ての賛否（１）からの続き
昨日の沖縄県民投票の結果が出た。ざっくりまとめると，
有権者数[A]： 115万人
投票数[B]（投票率=B/A）： 60.5万票（52.8%）
有効投票数[C]（有効投票率=C/A）： 60.2万票（52.2%）
賛成票[D]（賛成比率=D/C｜賛成割合=D/A）： 11.5万票（19.1%｜10.0%）
反対票[E]（反対比率=E/C｜反対割合=E/A）： 43.4万票（72.2%｜37.6%）
どちらでもない票[F]（中間比率=F/C ｜中間割合=F/A ）：5.3万票（8.8%｜4.6%）
・有権者数に占める割合が最も多いものは反対の37.6%で，判定基準の25%を越えた。
・沖縄県民の７割以上が反対している（と普通は表現するであろう）。
・玉城デニー知事の県知事選の得票数 39.7万票を43.4万票で上回っている。
・NHKや産経新聞や読売新聞などは，結果を極力過小に示すような報道を続けた。
・政府はこの結果を無視して引き続き埋め立て工事を進めるとした。
参考： 元山仁士郎さんのTwitter
昨日の沖縄県民投票の結果が出た。ざっくりまとめると，
有権者数[A]： 115万人
投票数[B]（投票率=B/A）： 60.5万票（52.8%）
有効投票数[C]（有効投票率=C/A）： 60.2万票（52.2%）
賛成票[D]（賛成比率=D/C｜賛成割合=D/A）： 11.5万票（19.1%｜10.0%）
反対票[E]（反対比率=E/C｜反対割合=E/A）： 43.4万票（72.2%｜37.6%）
どちらでもない票[F]（中間比率=F/C ｜中間割合=F/A ）：5.3万票（8.8%｜4.6%）
・有権者数に占める割合が最も多いものは反対の37.6%で，判定基準の25%を越えた。
・沖縄県民の７割以上が反対している（と普通は表現するであろう）。
・玉城デニー知事の県知事選の得票数 39.7万票を43.4万票で上回っている。
・NHKや産経新聞や読売新聞などは，結果を極力過小に示すような報道を続けた。
・政府はこの結果を無視して引き続き埋め立て工事を進めるとした。
参考： 元山仁士郎さんのTwitter
写真：沖縄県名護市辺野古（CC-BY-SA Kugel~commonswiki）
2019年2月24日日曜日
辺野古米軍基地建設のための埋め立ての賛否（１）
沖縄県のウェブサイト(2/24/2019)ではこうなっている。
沖縄県民投票、午後11時ごろ大勢判明 期日前に23万人、辺野古に直接審判
沖縄タイムズ(2/24 5:00)ではこうなっている。
きょう沖縄県民投票 辺野古埋め立て賛否を判断 得票数・率に注目
日本経済新聞朝刊(2/24)は５面の５段の小さな記事でこれ。
沖縄きょう県民投票
米軍普天間基地（沖縄県宜野湾市）の名護市辺野古移設を巡る県民投票が...
（いきなりテーマが誘導されている）
TBSの報道特集(2/23)では，金平さんの現地取材を含めて長時間取り上げていた。
NHKのNW9(2/23)では，棄権派とどちらでもない派だけを取り上げていた。2/21の特集では，22年前の移設賛否投票に関する防衛施設庁・防衛庁関係者の意見とりあげて，県民分断や沖縄の苦悩という感情的フレーズにまぶして投票忌避を誘導していた。
（辺野古米軍基地建設のための埋め立ての賛否（２）に続く）
琉球新報(2/24 5:00)ではこうなっている。平成31年２月24日（日）は県民投票です。辺野古米軍基地建設のための埋立ての賛否を問う県民投票について県民投票とは、通常の選挙とは異なり、特定の候補者に投票するものではありません。『普天間飛行場の代替施設として国が名護市辺野古に計画している米軍基地建設のための埋立て』について、県民の意思を示すための投票です。
沖縄県民投票、午後11時ごろ大勢判明 期日前に23万人、辺野古に直接審判
沖縄タイムズ(2/24 5:00)ではこうなっている。
きょう沖縄県民投票 辺野古埋め立て賛否を判断 得票数・率に注目
日本経済新聞朝刊(2/24)は５面の５段の小さな記事でこれ。
沖縄きょう県民投票
米軍普天間基地（沖縄県宜野湾市）の名護市辺野古移設を巡る県民投票が...
（いきなりテーマが誘導されている）
TBSの報道特集(2/23)では，金平さんの現地取材を含めて長時間取り上げていた。
NHKのNW9(2/23)では，棄権派とどちらでもない派だけを取り上げていた。2/21の特集では，22年前の移設賛否投票に関する防衛施設庁・防衛庁関係者の意見とりあげて，県民分断や沖縄の苦悩という感情的フレーズにまぶして投票忌避を誘導していた。
（辺野古米軍基地建設のための埋め立ての賛否（２）に続く）
2019年2月23日土曜日
はやぶさ２
2019年2月22日の午前7時半ごろ，小惑星探査機はやぶさ２が小惑星リュウグウへの最初の着陸（タッチダウン）に成功した。2010年の６月，はやぶさがイトカワからサンプルを持ち帰った。これに関して，当時大阪大学理学研究科宇宙地球科学専攻惑星物質学グループ教授だった土山明先生（2012年から京都大学大学院理学研究科地球惑星科学専攻鉱物学講座教授）に，日本物理教育学会近畿支部の講演でお話をうかがったことを思い出した。
イトカワは直径390 m（体積0.018 km^3），質量3.5 × 10^10 kgであり，その密度は1.9 g/cm^3 である。地球表面に持ってくればグズグズにくずれそうな物体である，と聞いたような気がする。
リュウグウはどうかと思い Wikipediaをみると，直径 700 m（この辺が日本のWikipediaの限界か，英語版では 865±15 mとして参考文献も示されている），質量が 4.5 × 10^11 kgであり，イトカワより一回り大きいことがわかる。体積や密度は記載されていないので，体積を直径の３乗の半分だと仮定すると（立方体の周囲を半分だけ削ったような形だと想像）0.7 g/cm^3となる。おかしいな，水より軽い。やはり体積の式が間違っているので，その結果が報告されるのを待つほうがよろしい。
リュウグウは炭素質コンドライトからなるC型小惑星であり，それらの平均の密度は1.38 g/cm^3 程度のようだ。M型小惑星はニッケルや鉄などの金属からできていて，その平均の密度は 5.32 g/cm^3 になる（Wikipedia のStandard asteroid physical characteristics より，これも中文版はあるが日本語版はまだない）。
イトカワは直径390 m（体積0.018 km^3），質量3.5 × 10^10 kgであり，その密度は1.9 g/cm^3 である。地球表面に持ってくればグズグズにくずれそうな物体である，と聞いたような気がする。
リュウグウはどうかと思い Wikipediaをみると，直径 700 m（この辺が日本のWikipediaの限界か，英語版では 865±15 mとして参考文献も示されている），質量が 4.5 × 10^11 kgであり，イトカワより一回り大きいことがわかる。体積や密度は記載されていないので，体積を直径の３乗の半分だと仮定すると（立方体の周囲を半分だけ削ったような形だと想像）0.7 g/cm^3となる。おかしいな，水より軽い。やはり体積の式が間違っているので，その結果が報告されるのを待つほうがよろしい。
リュウグウは炭素質コンドライトからなるC型小惑星であり，それらの平均の密度は1.38 g/cm^3 程度のようだ。M型小惑星はニッケルや鉄などの金属からできていて，その平均の密度は 5.32 g/cm^3 になる（Wikipedia のStandard asteroid physical characteristics より，これも中文版はあるが日本語版はまだない）。
写真：小惑星リュウグウ（Wikipedia 英語版より）
2019年2月22日金曜日
NeXTの話
1985年にアップルから追いやられたスティーブ・ジョブズが立ち上げたのが，NeXTである。1990年の6月にMITで，素粒子と原子核（の中間領域）の国際会議PANIC XIIが開催された。森田正人先生に連れられて，一年先輩の佐藤透さんといっしょに，ワシントン経由でボストンに着いた。
会場のMITでは，コンピュータルームが参加者のオープン利用に供されていた。そこにあったのが，初代のNeXT Cubeであった。1990年の３月に天王寺分校のデータステーションにS4/330が導入されて，ワークステーションにはなじんでおり，NeXTの感触を確かめることができた。Mathematicaも使えてこれはすごいということで，翌1991年に，大阪教育大学物理１研の予算を工面してもらって，NeXT Station（CPU 68040 25Hz, 主記憶 8MB, ディスク120MB?）を１台購入することができた。天王寺分校は移転を控えてバタバタしていたので，そのどさくさである。
天王寺分校のデータステーションでは，1990年にサンのワークステーションS/4-330を導入すると同時に室内にはイーサネット（10BASE5 イエローケーブル）を張って室内の端末等と接続していた。垣本さんと相談して，学内ネットワークの練習のためにその一部を延長して，同じ２階で2-30m離れた私の研究室までシン・イーサネット（10BASE2）の同軸ケーブルを引き回した（ケーブルは日本橋で自費で調達し，コネクタは垣本さんにもらったのだった）。彼は技術教育専攻の出身で，ハンダ付けのプロだったので，コネクタ部分はお願いして作ってもらった。
これにより，晴れて（ヤミか？），研究室のNeXT Stationはインターネットにつながることになった。次の年の1992年にはデータステーションは情報処理センターとして教養学科とともに柏原キャンパスに移転し，引き続いて，我々教員養成課程も1993年から新キャンパスでの生活が始まった（柏原キャンパスでは，各研究室への情報コンセントを整備していた）。
若くして亡くなられた山口英（阪大基礎工情報科学）さんが，1990年から1991年にかけて情報処理教育センター（豊中）の助手をされていた。そのとき，阪大の情報教育用のシステムとしてNeXT stationを400台導入するというので，豊中地区に行ったついでに旧豊中データステーションに立ち寄り，資料をもらうことができた。自分の記憶の中では，気軽に資料をくれた当時の山口さんはヒゲを生やしていたが，すずきひろのぶ氏と混線しているのかもしれない。
会場のMITでは，コンピュータルームが参加者のオープン利用に供されていた。そこにあったのが，初代のNeXT Cubeであった。1990年の３月に天王寺分校のデータステーションにS4/330が導入されて，ワークステーションにはなじんでおり，NeXTの感触を確かめることができた。Mathematicaも使えてこれはすごいということで，翌1991年に，大阪教育大学物理１研の予算を工面してもらって，NeXT Station（CPU 68040 25Hz, 主記憶 8MB, ディスク120MB?）を１台購入することができた。天王寺分校は移転を控えてバタバタしていたので，そのどさくさである。
写真：NeXTstation（Wikipediaより）
天王寺分校のデータステーションでは，1990年にサンのワークステーションS/4-330を導入すると同時に室内にはイーサネット（10BASE5 イエローケーブル）を張って室内の端末等と接続していた。垣本さんと相談して，学内ネットワークの練習のためにその一部を延長して，同じ２階で2-30m離れた私の研究室までシン・イーサネット（10BASE2）の同軸ケーブルを引き回した（ケーブルは日本橋で自費で調達し，コネクタは垣本さんにもらったのだった）。彼は技術教育専攻の出身で，ハンダ付けのプロだったので，コネクタ部分はお願いして作ってもらった。
これにより，晴れて（ヤミか？），研究室のNeXT Stationはインターネットにつながることになった。次の年の1992年にはデータステーションは情報処理センターとして教養学科とともに柏原キャンパスに移転し，引き続いて，我々教員養成課程も1993年から新キャンパスでの生活が始まった（柏原キャンパスでは，各研究室への情報コンセントを整備していた）。
若くして亡くなられた山口英（阪大基礎工情報科学）さんが，1990年から1991年にかけて情報処理教育センター（豊中）の助手をされていた。そのとき，阪大の情報教育用のシステムとしてNeXT stationを400台導入するというので，豊中地区に行ったついでに旧豊中データステーションに立ち寄り，資料をもらうことができた。自分の記憶の中では，気軽に資料をくれた当時の山口さんはヒゲを生やしていたが，すずきひろのぶ氏と混線しているのかもしれない。
2019年2月21日木曜日
スーパーコンピュータ登場前夜
大阪大学大型計算機センターニュース（センター25周年のあゆみ：出口弘）によると，センターシステムは1980年代の前後には次のように変遷している。我々がRCNPから大型計算機センターに乗換え始めたのは，1980年に入出力棟が新設されたころからである。
1969.05 大型計算機センター開設
1976.09 ACOS700導入
1977.12 ACOS800導入
1978.10 ACOS900導入
1979.05 センター設置10周年
1980.03 入出力棟新設
1981.12 ACOS1000導入
1982.05 ACOS1000増設・IAPサービス開始・FORTRAN77
1985.01 HFPサービス開始
1986.06 SX-1サービス開始
1987.04 学術情報ネットワーク・大学間ネットワーク
1987.11 ACOS2020サービス開始
1988.01 SX-2Nサービス開始
1993.02 SX-3Rサービス開始
1994,01 ACOS3900サービス開始
1994.05 センター設置25周年
2019.01 大型計算機センター法制化50周年（記念シンポジウムのページが保護中...orz）
1970年代の半ばに，CRAY-1が開発されたのに続いて，国産の大型汎用機メーカーがスーパーコンピュータの開発を進めていた。日本電気の開発したベクトル型スーパーコンピュータのSXシリーズが登場する前，汎用機であるACOS1000に統合アレイプロセッサ（IAP）とよばれるベクトル計算機構が導入されて，1982年からサービスが始まった。ちょうど，大阪教育大学に就職したときで，天王寺のデータステーションから阪大大型計算機センターのTSS利用が可能になったころである。
1983年ごろの研究ノートをみると，FORTRAN66からFORTRAN77へのプログラム書き換えや，ACOS1000のTSSのジョブ処理手順への変更などで四苦八苦している。FORTANのコンパイル時の最適化オプションで実行時間が1/3になり，IAPオプションでさらにその1/2になるというのが，最も効果のある場合の例だった。
1985年には，ACOS1000のバックエンドシステムとして高速FORTRANプロセッサ（HFP）が導入され，ここでもIAPが使えた。これが，SX-1の導入イメージのための練習用システムとなっていたようだ。翌1986年にはSX-1が使えるようになった。
1986年に核物理センターの計算費をもらって，大型計算機センターでSX-1を利用した話（スーパーコンピュータと原子核殻模型計算 1987）は，Juliaでパズル（５）で紹介したとおりである。
1969.05 大型計算機センター開設
1976.09 ACOS700導入
1977.12 ACOS800導入
1978.10 ACOS900導入
1979.05 センター設置10周年
1980.03 入出力棟新設
1981.12 ACOS1000導入
1982.05 ACOS1000増設・IAPサービス開始・FORTRAN77
1985.01 HFPサービス開始
1986.06 SX-1サービス開始
1987.04 学術情報ネットワーク・大学間ネットワーク
1987.11 ACOS2020サービス開始
1988.01 SX-2Nサービス開始
1993.02 SX-3Rサービス開始
1994,01 ACOS3900サービス開始
1994.05 センター設置25周年
2019.01 大型計算機センター法制化50周年（記念シンポジウムのページが保護中...orz）
1970年代の半ばに，CRAY-1が開発されたのに続いて，国産の大型汎用機メーカーがスーパーコンピュータの開発を進めていた。日本電気の開発したベクトル型スーパーコンピュータのSXシリーズが登場する前，汎用機であるACOS1000に統合アレイプロセッサ（IAP）とよばれるベクトル計算機構が導入されて，1982年からサービスが始まった。ちょうど，大阪教育大学に就職したときで，天王寺のデータステーションから阪大大型計算機センターのTSS利用が可能になったころである。
1983年ごろの研究ノートをみると，FORTRAN66からFORTRAN77へのプログラム書き換えや，ACOS1000のTSSのジョブ処理手順への変更などで四苦八苦している。FORTANのコンパイル時の最適化オプションで実行時間が1/3になり，IAPオプションでさらにその1/2になるというのが，最も効果のある場合の例だった。
1985年には，ACOS1000のバックエンドシステムとして高速FORTRANプロセッサ（HFP）が導入され，ここでもIAPが使えた。これが，SX-1の導入イメージのための練習用システムとなっていたようだ。翌1986年にはSX-1が使えるようになった。
1986年に核物理センターの計算費をもらって，大型計算機センターでSX-1を利用した話（スーパーコンピュータと原子核殻模型計算 1987）は，Juliaでパズル（５）で紹介したとおりである。
登録:
投稿 (Atom)