　　　　On a Thread of the Web

2019年1月29日火曜日

Juliaでパズル（５）

（Juliaでパズル（４）から続く）

結論：失敗でした。７時間以上かかるプログラムを高速化するために，重複計算部分を配列に落とし込もうとしたが，プログラムの見通しが非常に悪くなって，しかも時間も短縮できないという残念な結果に終わってしまった。

30年くらい前，スーパーコンピュータを使って原子核の殻模型計算を高速化するための試みについて，大阪大学の大型計算機センターニュースに投稿したことがある（スーパーコンピュータと原子核殻模型計算 1987）。導入されて間もなかったNECのベクトル型計算機の最初のモデルのSX-1を使い，$sd$殻の11粒子系である${}^{29}{\rm Si}$の基底状態（$J=\frac{1}{2}, T=\frac{1}{2}$，2360次元）を求める計算を高速化しようというものである（西村道明さんとの共同研究のための予備計算，大坪久夫先生にRCNPの計算費を確保していただいた）。

結果は，ベクトル化よりも重複計算を間接アドレッシングやメモリに蓄えて呼び出すことで避ける効果が高いというなさけないものであったが，その時から重複計算を避けようというのがスローガンとしてたたき込まれていたのであった（しかし，ベクトルマシンそのものが陳腐化して消えてしまうわけで，世の中の栄枯盛衰，諸行無常なり）。

今回の結果は，そもそも重複計算がそんなにあったのかというところからして問題であり，コードが複雑になって保守性にも欠けるという結末で終わってしまった。もともと重複だと思っている関数呼び出しも大したコストがかかっておらず，多次元の配列呼び出しとあまりかわらなかったのかもしれない。とりあえず，悪戦苦闘した（残念な）プログラムは以下の通り

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
return sqrt(abs(x))
end

function cbr(x)
return cbrt(abs(x))
end

function rog(x)
return log(abs(x))
end

function cc(a,b)
if (a==1 && b==1) || (a==11 && b==1)
return 10*a+b
else
return 0
end
end

function ne(op, op1)
expr = Expr(:call, op, op1)
return expr
end

function me(op, op1, op2)
expr = Expr(:call, op, op1, op2)
return expr
end

op1=[:+,:sqr,:cbr,:exp,:rog]
op2=[:+,:-,:*,:/,:cc]
ne1=Array{Any}(undef, 5)
ne2=Array{Any}(undef, 5, 5, 5)
ne11=Array{Any}(undef, 5, 5, 5, 5)
nex1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
ne1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5)
mexy11=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey1x1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
mey11x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1yx1=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)
me1y1x=Array{Any}(undef, 5, 5, 5, 5, 5, 5, 5, 5, 5)

function stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)
for i in 1:5
ne1[i]=Expr(:call,op1[i],1)
end
for i in 1:5, j in 1:5, k in 1:5
ne2[i,j,k]=Expr(:call,op2[i],ne1[j],ne1[k])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5
ne11[i,j,k,l]=Expr(:call,op1[i],ne2[j,k,l])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
nex1[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne11[k,l,m,n],ne1[o]))
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5
ne1x[i,j,k,l,m,n,o]=Expr(:call,op1[i],Expr(:call,op2[j],ne1[k],ne11[l,m,n,o]))
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mexy11[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne11[j,k,l,m],ne11[n,o,p,q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mey1x1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],nex1[j,k,l,m,n,o,p],ne1[q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
mey11x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1x[j,k,l,m,n,o,p],ne1[q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
me1yx1[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],nex1[k,l,m,n,o,p,q])
end
for i in 1:5, j in 1:5, k in 1:5, l in 1:5, m in 1:5, n in 1:5, o in 1:5, p in 1:5, q in 1:5
me1y1x[i,j,k,l,m,n,o,p,q]=Expr(:call,op2[i],ne1[j],ne1x[k,l,m,n,o,p,q])
end
end

@time stor(ne1,ne11,ne1x,nex1,mexy11,mey1x1,mey11x,me1yx1,me1y1x,op1,op2)

function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
# x=ne11[d[5],a,d[1],d[2]]
# y=ne11[d[6],c,d[3],d[4]]
z=ne(op1[d[7]],mexy11[b,d[5],a,d[1],d[2],d[6],c,d[3],d[4]])
if b==5
z=:0
end
return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
# x=ne11[d[3],a,d[1],d[2]]
# y=nex1[d[5],b,d[3],a,d[1],d[2],d[4]]
z=ne(op1[d[7]],mey1x1[c,d[5],b,d[3],a,d[1],d[2],d[4],d[6]])
if (a!=5 && b==5) || ((a!=5 || b!=5) && c==5)
z=:0
end
return (z,eval(z))
end

function doe3(a,b,c,d)
# ((a,(b,c)),d)
# x=me(b,ne(d[1],1),ne(d[2],1))
# y=ne1x[d[5],a,d[3],d[4],b,d[1],d[2]]
z=ne(op1[d[7]],mey11x[c,d[5],a,d[3],d[4],b,d[1],d[2],d[6]])
if a==5 || c==5
z=:0
end
return (z,eval(z))
end

function doe4(a,b,c,d)
# (a,((b,c),d))
# x=me(b,ne(d[1],1),ne(d[2],1))
# y=nex1[d[5],c,d[3],b,d[1],d[2],d[4]]
z=ne(op1[d[7]],me1yx1[a,d[6],d[5],c,d[3],b,d[1],d[2],d[4]])
if a==5 || (b!=5 && c==5)
z=:0
end
return (z,eval(z))
end

function doe5(a,b,c,d)
# (a,(b,(c,d)))
# x=me(c,ne(d[1],1),ne(d[2],1))
# y=ne1x[d[5],b,d[3],d[4],c,d[1],d[2]]
z=ne(op1[d[7]],me1y1x[a,d[6],d[5],b,d[3],d[4],c,d[1],d[2]])
if a==5 || b==5
z=:0
end
return (z,eval(z))
end

function main(pz)
d=[1,1,1,1,1,1,1]

for a in 1:5, b in 1:5, c in 1:5
for d[1] in 1:5, d[2] in 1:5, d[3] in 1:5, d[4] in 1:5, d[5] in 1:1, d[6] in 1:1, d[7] in 1:5
e,f = doe1(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe2(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe3(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe4(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe5(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
end
end
end

function fine(pz)
println(length(pz))
qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(+1, +1)"=>"11")
rz=replace(rz,"cc(+(11), +1)"=>"111")
rz=replace(rz,"+1"=>"1")
println(i," -> ",rz," = ",pz[j][2])
end
end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18.116161 seconds (20.53 M allocations: 1.368 GiB, 85.01% gc time)
420.117500 seconds (115.60 M allocations: 6.676 GiB, 4.81% gc time)
740539
1 -> (exp(1 + +((rog(1) - rog(exp(1) + exp(1)))))) = 0.5000000000000001
2 -> (cbr(1 + +(1 / exp(11)))) = 1.000005567202603
3 -> (cbr(1 - +(exp(1) * sqr(11)))) = 2.0012925727345556
4 -> (+((1 + +((exp(1) - cbr(1 / exp(1))))))) = 3.001750517885256
5 -> (exp(1 + +((exp(1) - sqr(exp(1) + exp(1)))))) = 4.0013741789600505
6 -> (+(1 * +(exp(1) / rog(1 - exp(1))))) = 5.021535230269061
7 -> (exp(1 + +(1 / cbr(1 + 1)))) = 6.011657650957131
8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402
9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874
10 -> (exp(1 + +((rog(1) + cbr(1 - exp(1)))))) = 9.004695764350682
11 -> (rog(1 + +(exp(11) / exp(1)))) = 10.000045398899218
12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002
13 -> (rog(1 + +(exp(1) * exp(11)))) = 12.000006144193478
14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708
15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938
16 -> (sqr(1 + +((exp(1) - exp(exp(1) + exp(1)))))) = 15.031080541832813
17 -> (+((1 + +(exp(1) / exp(1 - exp(1)))))) = 16.15426224147926
18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262
19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822
20 -> (exp(1 + +((rog(1) + cbr(exp(1) * exp(1)))))) = 19.062709434872296
21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664
22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717
23 -> (exp(1 + +((exp(1) - sqr(1 / exp(1)))))) = 22.46034183113636
24 -> (exp(1 + +(exp(1) / cbr(1 + 1)))) = 23.511783406204675
25 -> (exp(1 / +((rog(1 + exp(1)) - 1)))) = 24.34239023161335
26 -> (sqr(1 + +(exp(exp(1) + exp(1)) * exp(1)))) = 25.005158374895853
27 -> (exp(1 / +((1 - rog(1 + 1))))) = 26.020673411168517
28 -> (exp(1 + +((1 + sqr(1 - exp(1)))))) = 27.40793292860369
29 -> (cbr(1 - +(exp(11) / exp(1)))) = 28.031200676839138
30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598
31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514
32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058
33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581
35 -> (exp(1 / +((1 - cbr(1 / exp(1)))))) = 34.046473983702036
36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533
37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296
38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464
39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416
40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508
41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533
42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161
43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116
44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116
45 -> (exp(1 + +((exp(1) + exp(rog(1) - exp(1)))))) = 44.00353027860415
47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535
48 -> (exp(1 + +((exp(exp(1) - 1) - exp(1))))) = 47.30706718593521
50 -> (exp(1 + +((exp(1) + exp(1 - exp(1)))))) = 49.286780801520734
51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701
52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723
54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407
55 -> (cbr(1 - +(exp(1) * exp(11)))) = 54.598038212039256
56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455
57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197
58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067
60 -> (exp(1 + +((exp(1) + exp(rog(1) - 1))))) = 59.511005963978846
62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258
63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041
65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827
67 -> (sqr(1 - +(exp(exp(1) * exp(1)) * exp(1)))) = 66.31488392984271
68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707
69 -> (exp(1 + +((1 + cbr(11))))) = 68.30480329758417
70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886
71 -> (exp(1 + +(cbr(1 - exp(1)) * exp(1)))) = 70.51403099154743
72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819
74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825
75 -> (exp(1 + +((rog(1) + sqr(11))))) = 74.93527870314952
76 -> (exp(1 + +((exp(1) + sqr(1 / exp(1)))))) = 75.55134476063684
77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655
79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831
82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715
83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322
84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863
85 -> (exp(1 + +((exp(1) + cbr(1 / exp(1)))))) = 84.33636424122227
86 -> (+((1 + +(exp(exp(1) + exp(1)) / exp(1))))) = 85.48412584713863
88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563
89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601
96 -> (exp(1 + +(sqr(1 - exp(1)) * exp(1)))) = 95.89100192175613
97 -> (exp(1 + +(rog(1 + exp(1)) * exp(1)))) = 96.52628755961122
98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673
100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271
105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627
107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063
110 -> (+((+(111) - 1))) = 110
111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968
112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675
113 -> (+((1 + +(exp(1 + exp(1)) * exp(1))))) = 112.9756938401968
120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112
121 -> (exp(1 + +(cbr(exp(1) + rog(1)) * exp(1)))) = 120.74343163832631
10.825410 seconds (63.96 M allocations: 993.083 MiB, 36.17% gc time)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

2019年1月28日月曜日

地獄からの復活劇

新聞によりますと，岡山の小学校５年生が，地獄−極楽の移動時間を求めたとのこと。移動時間は12.7年とのことで，移動に用いられる蜘蛛の糸の強度も割り出したそうだ。

これは，是非とも調べてみなければいけません。一般財団法人理数教育研究所（代表理事岡本和夫）が主催している「算数・数学の自由研究」作品コンクールの審査員特別賞の受賞作の１つであった。「地獄からの復活劇〜御釈迦様からの試練〜岡山県岡山大学教育学部附属小学校５年末光由季」

地獄と極楽の距離は，何万里という芥川の本文の表現の下限から１万里≒４万ｋｍとしている。静止軌道半径に近いので，極楽は静止軌道上にあるのかもしれない。ということは蜘蛛の糸は宇宙エレベータの暗示だったのか・・・

すばらしいのは，蜘蛛の糸を登る時間を自分で実験しているところだ。近所の公園の登り棒を３ｍ昇る時間の平均から，５秒で１ｍ（0.2ｍ／秒）を得ている。最近は，危険だということで公園から次々と遊具が撤去されている。うちの近くの公園でも子どもたちが遊んでいたブランコや回転遊具は何年か前に消えてしまい，もうすぐ鉄棒までなくなるそうだ。組み体操をやめるのは賛成だが，公園の遊具は残してほしかった。

話がそれた。結局，４万ｋｍ÷ 0.2 ｍ／秒＝２億秒＝ 6.4 年・・・あれ？ファクター２違うな。２億秒まではあっている。で，１年≒ π × 10^７秒ですよね。２億秒から時間に換算する計算が２倍になっているようだ。１日の半分は休憩とすればつじつまがあうということにしておこう。

さらに蜘蛛の糸の強度まで計算してピアノ線と比較しているところがすごい。最後のコメントも含蓄がある。宇宙エレベータ協会の皆様は是非，彼女を顧問に迎えられたい。

軌道エレベータの基部の想像図（Wikipediaより）

2019年1月27日日曜日

Juliaでパズル（４）

（Juliaでパズル（３）からの続き）

目的はJuliaプログラミングの習得ということなのでこの話題がだらだたと続いている。高速が話題のJuliaで初めて実行時間が１時間以上のプログラムを作ったが，アルゴリズムがポンコツだといくらでも遅いプログラムができるのであった・・・orz

昨日の２項演算だけの組み合わせで整数を作るものから，単項演算を加えた組み合わせに拡張する。$5 \times 5^3 = 625$回から$5 \times 5^3 \times 5^7 ≒ 5000$万回に増えて７時間で終わらなかったので，次は時間短縮の方法を考えなければならない。

単項演算は１つのエントリーポイントで１回だけ作用するものとしている。
以下は，７重ループを５重ループに落として７分程度で実行したもの。
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function sqr(x)
return sqrt(abs(x))
end

function cbr(x)
return cbrt(abs(x))
end

function rog(x)
return log(abs(x))
end

function cc(a,b)
if (a==1 && b==1) || (a==11 && b==1)
return 10*a+b
else
return 0
end
end

function ne(op, op1)
expr = Expr(:call,Symbol(op), op1)
return expr
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function doe1(a,b,c,d)
# ((1,1)_a,(1,1)_c)_b
x=me(a,ne(d[1],1),ne(d[2],1))
y=me(c,ne(d[3],1),ne(d[4],1))
z=ne(d[7],me(b,ne(d[5],x),ne(d[6],y)))
if b=="cc"
z=:0
end
return (z,eval(z))
end

function doe2(a,b,c,d)
# (((1,1)_a,1)_b,1)_c
x=me(a,ne(d[1],1),ne(d[2],1))
y=me(b,ne(d[3],x),ne(d[4],1))
z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
if (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
z=:0
end
return (z,eval(z))
end

function doe3(a,b,c,d)
# ((1,(1,1)_b)_a,1)_c
x=me(b,ne(d[1],1),ne(d[2],1))
y=me(a,ne(d[3],1),ne(d[4],x))
z=ne(d[7],me(c,ne(d[5],y),ne(d[6],1)))
if a=="cc" || c=="cc"
z=:0
end
return (z,eval(z))
end

function doe4(a,b,c,d)
# (1,((1,1)_b,1)_c)_a
x=me(b,ne(d[1],1),ne(d[2],1))
y=me(c,ne(d[3],x),ne(d[4],1))
z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
if a=="cc" || (b!="cc" && c=="cc")
z=:0
end
return (z,eval(z))
end

function doe5(a,b,c,d)
# (1,(1,(1,1)_c)_b)_a
x=me(c,ne(d[1],1),ne(d[2],1))
y=me(b,ne(d[3],1),ne(d[4],x))
z=ne(d[7],me(a,ne(d[6],1),ne(d[5],y)))
if a=="cc" || b=="cc"
z=:0
end
return (z,eval(z))
end

function main(pz)
uno=["+","-","*","/","cc"]
dos=["+","sqr","cbr","exp","rog"]
d=["+","+","+","+","+","+","+"]
for a in uno, b in uno, c in uno
#a="*"; b="-"; c="/";
for d[1] in dos, d[2] in dos, d[3] in dos, d[4] in dos, d[7] in dos
#, d[5] in dos, d[6] in dos, d[7] in dos
e,f = doe1(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe2(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe3(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe4(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
e,f = doe5(a,b,c,d)
if f>0.5 && f<200
push!(pz,(e,f))
end
end
end
end

function fine(pz)
println(length(pz))
qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(+1, +1)"=>"11")
rz=replace(rz,"cc(+(11), +1)"=>"111")
rz=replace(rz,"+1"=>"1")
println(i," -> ",rz," = ",pz[j][2])
end
end
end

pz=[]
@time main(pz)
@time fine(pz)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

408.936226 seconds (144.22 M allocations: 8.591 GiB, 5.96% gc time)

751191

1 -> (exp(+((rog(1) - rog(exp(1) + exp(1)))) + 1)) = 0.5000000000000001

2 -> (cbr(+(1 / exp(11)) + 1)) = 1.000005567202603

3 -> (cbr(+(exp(1) * sqr(11)) - 1)) = 2.0012925727345556

4 -> (+((+((exp(1) - cbr(1 / exp(1)))) + 1))) = 3.001750517885256

5 -> (exp(+((exp(1) - sqr(exp(1) + exp(1)))) + 1)) = 4.0013741789600505

6 -> (+(+(exp(1) / rog(1 - exp(1))) * 1)) = 5.021535230269061

7 -> (exp(+(1 / cbr(1 + 1)) + 1)) = 6.011657650957131

8 -> (exp(+((1 + cbr(exp(1) * exp(1)))) - 1)) = 7.012778894114402

9 -> (+((+(exp(1) * sqr(11)) - 1))) = 8.015520899439874

10 -> (exp(+((rog(1) + cbr(1 - exp(1)))) + 1)) = 9.004695764350682

11 -> (rog(+(exp(11) / exp(1)) + 1)) = 10.000045398899218

12 -> (exp(+((1 + rog(11))) - 1)) = 11.000000000000002

13 -> (rog(+(exp(1) * exp(11)) + 1)) = 12.000006144193478

14 -> (exp(+(rog(1 + exp(1)) * exp(1)) - 1)) = 13.063412466658708

15 -> (exp(+((sqr(exp(1) + rog(1)) + rog(1))) + 1)) = 14.135951028495938

16 -> (sqr(+((exp(1) - exp(exp(1) + exp(1)))) + 1)) = 15.031080541832813

17 -> (+((+(exp(1) / exp(1 - exp(1))) + 1))) = 16.15426224147926

18 -> (+((+((1 + exp(exp(1) + rog(1)))) + 1))) = 17.154262241479262

19 -> (exp(+((exp(1) + exp(1 - exp(1)))) * 1)) = 18.131593378402822

20 -> (exp(+((rog(1) + cbr(exp(1) * exp(1)))) + 1)) = 19.062709434872296

21 -> (+(+((exp(1) + rog(1))) * +(exp(1) * exp(1)))) = 20.085536923187664

22 -> (exp(+(+(11) / exp(1)) - 1)) = 21.045228353448717

23 -> (exp(+((exp(1) - sqr(1 / exp(1)))) + 1)) = 22.46034183113636

24 -> (exp(+(exp(1) / cbr(1 + 1)) + 1)) = 23.511783406204675

25 -> (sqr(+(exp(exp(1) * exp(1)) / exp(1)) - 1)) = 24.37815447036041

26 -> (sqr(+(exp(exp(1) + exp(1)) * exp(1)) + 1)) = 25.005158374895853

27 -> (exp(+((exp(1) + rog(1 - exp(1)))) * 1)) = 26.039293433236857

28 -> (exp(+((1 + sqr(1 - exp(1)))) + 1)) = 27.40793292860369

29 -> (cbr(+(exp(11) / exp(1)) - 1)) = 28.031200676839138

30 -> (+(+((exp(1) + exp(1))) * +((exp(1) + exp(1))))) = 29.556224395722598

31 -> (exp(+((exp(1) + rog(exp(1) + exp(1)))) - 1)) = 30.308524482958514

32 -> (exp(+((exp(1) + cbr(1 / exp(1)))) * 1)) = 31.025614547492058

33 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) - 1)) = 32.350947201581

35 -> (exp(+((exp(1) - exp(1 - exp(1)))) + 1)) = 34.42929324111139

36 -> (exp(+(sqr(1 - exp(1)) * exp(1)) * 1)) = 35.27632820034533

37 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) - 1)) = 36.68805458104296

38 -> (exp(+(11) - +(exp(1) * exp(1)))) = 37.00096158422464

39 -> (exp(+((sqr(1 + exp(1)) + exp(1))) - 1)) = 38.34279034771416

40 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) - 1)) = 39.09583226076508

41 -> (+(+((exp(1) + exp(1))) * +(exp(1) * exp(1)))) = 40.17107384637533

42 -> (exp(+((+((1 + 1)) + exp(1))) - 1)) = 41.1935556747161

43 -> (+((+((exp(1 + exp(1)) + rog(1))) + 1))) = 42.193555674716116

44 -> (+((+((1 + exp(1 + exp(1)))) + 1))) = 43.193555674716116

45 -> (exp(+((exp(1) + exp(rog(1) - exp(1)))) + 1)) = 44.00353027860415

47 -> (exp(+(sqr(1 + 1) * exp(1)) * 1)) = 46.72274206040535

48 -> (exp(+((exp(exp(1) - 1) - exp(1))) + 1)) = 47.30706718593521

50 -> (exp(+((exp(1) + exp(1 - exp(1)))) + 1)) = 49.286780801520734

51 -> (exp(+((exp(1) + cbr(1 - exp(1)))) * 1)) = 50.20065233451701

52 -> (exp(+((exp(1) + cbr(11))) - 1)) = 51.5350376483723

54 -> (exp(+((cbr(1 + 1) + exp(1))) * 1)) = 53.42094397564407

55 -> (cbr(+(exp(1) * exp(11)) - 1)) = 54.598038212039256

56 -> (exp(+(exp(1) / rog(1 - exp(1))) - 1)) = 55.78668552594455

57 -> (exp(+((exp(1) + sqr(1 - exp(1)))) * 1)) = 56.21110430560197

58 -> (exp(+(1 / exp(1)) * +(11))) = 57.20686180895067

60 -> (exp(+((exp(1) + exp(rog(1) - 1))) + 1)) = 59.511005963978846

62 -> (exp(+((cbr(exp(1) + rog(1)) + exp(1))) * 1)) = 61.18452227596258

63 -> (exp(+((sqr(1 + 1) + exp(1))) * 1)) = 62.33327490494041

65 -> (exp(+((exp(1) + exp(1 / exp(1)))) * 1)) = 64.2607927021827

67 -> (sqr(+(exp(exp(1) * exp(1)) * exp(1)) - 1)) = 66.31488392984271

68 -> (exp(+(cbr(1 + exp(1)) * exp(1)) * 1)) = 67.43919204479707

69 -> (exp(+((1 + cbr(11))) + 1)) = 68.30480329758417

70 -> (exp(+(sqr(1 + exp(1)) * exp(1)) - 1)) = 69.52046855397886

71 -> (exp(+(cbr(1 - exp(1)) * exp(1)) + 1)) = 70.51403099154743

72 -> (exp(+((cbr(1 + exp(1)) + exp(1))) * 1)) = 71.34344142876819

74 -> (exp(+(cbr(exp(1) * exp(1)) * exp(1)) - 1)) = 73.2948282500825

75 -> (exp(+((rog(1) + sqr(11))) + 1)) = 74.93527870314952

76 -> (exp(+((exp(1) + sqr(1 / exp(1)))) + 1)) = 75.55134476063684

77 -> (exp(+((1 + sqr(exp(1) + exp(1)))) + 1)) = 76.06924026135655

79 -> (exp(+((sqr(exp(1) + rog(1)) + exp(1))) * 1)) = 78.80710038074831

82 -> (+(+(exp(1) * exp(1)) * +(11))) = 81.27961708823715

83 -> (exp(+((exp(1) + rog(1 + 1))) + 1)) = 82.3871113494322

84 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) - 1))) = 83.48412584713863

85 -> (exp(+((exp(1) + cbr(1 / exp(1)))) + 1)) = 84.33636424122227

86 -> (+((+(exp(exp(1) + exp(1)) / exp(1)) + 1))) = 85.48412584713863

88 -> (exp(+((cbr(exp(1) + exp(1)) + exp(1))) * 1)) = 87.93899191149563

89 -> (exp(+(sqr(exp(1) + rog(1)) * exp(1)) * 1)) = 88.38383317988601

96 -> (exp(+(sqr(1 - exp(1)) * exp(1)) + 1)) = 95.89100192175613

97 -> (exp(+(rog(1 + exp(1)) * exp(1)) + 1)) = 96.52628755961122

98 -> (exp(+(1 / exp(1 - exp(1))) - 1)) = 97.02236556502673

100 -> (exp(+(rog(exp(1) + exp(1)) * exp(1)) * 1)) = 99.72847208916271

105 -> (exp(+((sqr(1 + exp(1)) + exp(1))) * 1)) = 104.22651025460627

107 -> (exp(+((exp(1) + cbr(exp(1) * exp(1)))) * 1)) = 106.27349040292063

110 -> (+((+(cc(+(11), 1)) - 1))) = 110

111 -> (+((+(exp(1 + exp(1)) * exp(1)) - 1))) = 110.9756938401968

112 -> (exp(+((1 + 1)) + +((exp(1) + rog(1))))) = 111.97569384019675

113 -> (+((+(exp(1 + exp(1)) * exp(1)) + 1))) = 112.9756938401968

120 -> (exp(+(cbr(exp(1) + exp(1)) * exp(1)) * 1)) = 119.07124802191112

121 -> (exp(+(cbr(exp(1) + rog(1)) * exp(1)) + 1)) = 120.74343163832631

6.644247 seconds (64.72 M allocations: 1004.683 MiB, 1.93% gc time)

2019年1月26日土曜日

Juliaでパズル（３）

（Juliaでパズル（２）からの続き）

世の中で流布している数学パズルだと，４つの数字と加減乗除を組み合わせて10をつくる Make 10 が有名らしい。もちろん，この手の問題では計算結果はすべて整数となるのが普通であるが，今回は，実数の範囲で計算して最後に整数化できるものを探している。

とりあえず，２項演算だけを組み合わせた結果をすべて探すプログラムを作ってみた。

(1) findfirst()のように，ある目的の数に一致する数をコレクション（配列またはタプル）の中から探すとともに，その位置の添字を返すことができれば有り難い。いまのところ「タプルの配列に対して１コマンドである数字を見つける方法」がわからない。

(2) ２次元配列ならば融通がききそうなのであるが，「タプルの配列を１コマンドで多次元配列に展開する方法」がわからない。もちろん地道にやればできることなのだが。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cc(a,b)
return 10*a+b
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
x=me(a,1,1)
y=me(c,1,1)
z=me(b,x,y)
if b=="cc"
z=:0
end
return (z,eval(z))
end

function dop2(a,b,c)
# (((a,b),c),d)
x=me(a,1,1)
y=me(b,x,1)
z=me(c,y,1)
if (a!="cc" && b=="cc") || ((a!="cc" || b!="cc") && c=="cc")
z=:0
end
return (z,eval(z))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
x=me(b,1,1)
y=me(a,1,x)
z=me(c,y,1)
if a=="cc" || c=="cc"
z=:0
end
return (z,eval(z))
end

function dop4(a,b,c)
# (a,((b,c),d))
x=me(b,1,1)
y=me(c,x,1)
z=me(a,1,y)
if a=="cc" || (b!="cc" && c=="cc")
z=:0
end
return (z,eval(z))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
x=me(c,1,1)
y=me(b,1,x)
z=me(a,1,y)
if a=="cc" || b=="cc"
z=:0
end
return (z,eval(z))
end

pz=[]
uno=["+","-","*","/","cc"]
for a in uno, b in uno, c in uno
d,e = dop1(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop2(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop3(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop4(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
d,e = dop5(a,b,c)
if e>0.5 && e<Inf
push!(pz,(d,e))
end
end

qz=[]
sort!(pz, by=x->x[2])
pl=length(pz)
for i in 1:pl
push!(qz,Int(ceil(pz[i][2])))
end

for i in 1:121
if in(i,qz)==true
j=findfirst(isequal(i),qz)
rz=repr(pz[j][1])
rz=replace(rz,":"=>"")
rz=replace(rz,"cc(1, 1)"=>"11")
rz=replace(rz,"cc(11, 1)"=>"111")
println(i," -> ",rz," = ",pz[j][2])
end
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

1 -> (1 - 1 / 11) = 0.9090909090909091

2 -> (1 + 1 / 11) = 1.0909090909090908

3 -> ((1 + 1) + 1 * 1) = 3

4 -> ((1 + 1) + (1 + 1)) = 4

6 -> (11 / (1 + 1)) = 5.5

9 -> (11 - (1 + 1)) = 9

10 -> (1 * 11 - 1) = 10

11 -> ((1 + 11) - 1) = 11

12 -> (1 + 1 * 11) = 12

13 -> ((1 + 1) + 11) = 13

22 -> ((1 + 1) * 11) = 22

110 -> (111 - 1) = 110

111 -> (1 * 111) = 111

112 -> (1 + 111) = 112

121 -> (11 * 11) = 121

「冬麗の不思議をにぎる赤ン坊」（野澤節子　1920-1995）

（Juliaでパズル（４）に続く）

2019年1月25日金曜日

Juliaでパズル（２）

（Juliaでパズル（１）からの続き）

４つの数字から整数をつくるパズルをコンピュータで処理させるためには，式に対応する文字列を，プログラム内で式として評価する必要がある。これは，JuliaドキュメントのMetaprogrammingのところに詳しいが，ここでのポイントは，Juliaでは，Juliaプログラムの中で，Juliaコードを生成したり操作したりすることができるということである。具体例としては次のような例文が示されている。

julia> function math_expr(op, op1, op2)
           expr = Expr(:call, op, op1, op2)
           return expr
       end
math_expr (generic function with 1 method)

julia>  ex = math_expr(:+, 1, Expr(:call, :*, 4, 5))
:(1 + 4 * 5)

julia> eval(ex)
21

そこで（1, 1, 1, 1）パズルで上記を用い，２項演算子を３つ与えた場合の５つのパターンの構文とその値を出力する関数を書いてみる。cc(a, b) = 10*a + bは２つの数字を並べて２桁の数（例えば a=1, b=1 なら cc(1,1) = 11である）を作る関数である。funcion me() においてSymbol(op)は文字列として与えた２項演算子を Expr型に変換する関数である（注：Ver 0.6以前の資料では symbol と小文字で始まっているので，Ver. 1.0 以降では大文字にすること）。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
function cc(a,b)
return 10*a+b
end

function me(op, op1, op2)
expr = Expr(:call,Symbol(op), op1, op2)
return expr
end

function dop1(a,b,c)
# ((a,b),(cd))
x=me(a,1,1)
y=me(c,1,1)
z=me(b,x,y)
println((z,eval(z)))
end

function dop2(a,b,c)
# (((a,b),c),d)
x=me(a,1,1)
y=me(b,x,1)
z=me(c,y,1)
println((z,eval(z)))
end

function dop3(a,b,c)
# ((a,(b,c)),d)
x=me(b,1,1)
y=me(a,1,x)
z=me(c,y,1)
println((z,eval(z)))
end

function dop4(a,b,c)
# (a,((b,c),d))
x=me(b,1,1)
y=me(c,x,1)
z=me(a,1,y)
println((z,eval(z)))
end

function dop5(a,b,c)
# (a,(b,(c,d)))
x=me(c,1,1)
y=me(b,1,x)
z=me(a,1,y)
println((z,eval(z)))
end
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

dop1("cc","/","+")
dop2("cc","+","*")
dop3("-","*","/")
dop4("+","\\","^")
dop5("+","-","*")

(:(cc(1, 1) / (1 + 1)), 5.5)
(:((cc(1, 1) + 1) * 1), 12)
(:((1 - 1 * 1) / 1), 0.0)
(:(1 + (1 \ 1) ^ 1), 2.0)
(:(1 + (1 - 1 * 1)), 1)

（Juliaでパズル（３）に続く）

2019年1月24日木曜日

Juliaでパズル（１）

４つの同じ数字に演算をほどこして，１から順に整数を作るというパズルがある。
例えば，6÷6+6-6 = 1,　 6÷6+6÷6=2,　 (6+6+6)/6 = 3,　などなど。
これを，Juliaの基本演算子や基本関数だけを使って実行してみると，

1　1+1-1*1
2　1+1-1+1
3　1+1*1+1
4　1+1+1+1
5　Int(floor(11/(1+1)))
6　Int(ceil(11/(1+1)))
7　Int(floor(exp(1*1+1*1)))
8　Int(ceil(exp(1*1+1*1)))
9　11-1-1
10　11-1*1
11　11+1-1
12　11+1*1
13　11+1+1
14　Int(ceil(sqrt(sqrt(exp(11)))-1-1))
15　Int(floor(exp(1+sqrt(1+1+1))))
16　Int(ceil(exp(1+sqrt(1+1+1))))
17　Int(floor(exp((1+1)*sqrt(1+1))))
18　Int(ceil(exp((1+1)*sqrt(1+1))))
19　Int(floor(exp(1+1+1)-1))
20　Int(floor(exp(1+1*1+1)))
21　Int(floor(exp(1+1+1)+1))
22　Int(ceil(exp(1+1+1)+1))
23　Int(floor(exp(sqrt(1+1)*cbrt(11))))
24　Int(ceil(exp(sqrt(1+1)*cbrt(11))))
25　Int(floor(exp(sqrt(11))-1-1))
26　Int(floor(exp(sqrt(11))-1*1))
27　Int(floor(exp(sqrt(11))+1-1))
28　Int(ceil(exp(sqrt(11))+1-1))
29　Int(ceil(exp(sqrt(11))+1*1))
30　Int(ceil(exp(sqrt(11))+1+1))
31　Int(floor(sqrt(sqrt(exp(11)))*(1+1)))
32　Int(floor(cbrt(exp(sqrt(111)))-1))
33　Int(floor(cbrt(exp(sqrt(111)))*1))
34　Int(floor(cbrt(exp(sqrt(111)))+1))
35　Int(ceil(cbrt(exp(sqrt(111)))+1))
36　Int(floor(exp(sqrt(11+1+1))))
37　Int(floor(cbrt(exp(11))-1-1))
38　Int(ceil(cbrt(exp(11))-1-1))
39　Int(ceil(cbrt(exp(11))-1*1))
40　Int(ceil(cbrt(exp(11))+1-1))
41　Int(ceil(cbrt(exp(11))+1*1))
42　Int(ceil(cbrt(exp(11))+1+1))
43　Int(floor(cbrt(sqrt(1+1)*exp(11))))
44　Int(ceil(cbrt(sqrt(1+1)*exp(11))))
45　Int(ceil(exp(cbrt(111)-1)))
46　Int(ceil(exp(sqrt(sqrt(11))+1+1)))
47　Int(ceil(sqrt(exp(11-sqrt(11)))))
などなど

（１）もし，ユーザ定義関数が使えるならば，任意の n に対して，すべての整数を作り出すことが簡単にできるのであった・・・

f(n1,n2,n3,n4)=Int(n1/n2+n3-n4)
f(n)=Int(ceil(n+1/n))

f(1,1,1,1)
1
f(f(1,1,1,1))
2
f(f(f(1,1,1,1)))
3

（２）juliaドキュメントの算術的演算と基本関数から主な演算と関数を選んで使うことにする。例えば，２項演算ならば，+, -, *, / の４種とする（ここでは，整数ゼロ除算をさけるために÷と%は含めない）。数字を２〜４桁並べる操作が必要なのでこれに対応する２項演算を１つ付け加えると５種になる。この２項演算は cc(a,b)=10*a+b で実現できる（例えば３桁並べるには，cc(cc(a, b), c)とすればよい），ただし適応できる位置に制限があることに注意する。次に単項演算を選ぶ。ここでは，+, sqrt(abs()), qbrt(abs()), exp(), log(abs()) の５種とする（+は恒等演算子）。

（３）４つの数字（a, b, c, d, a = b = c = d）に対応する２項演算のパターンは５つある。
２項演算を（a, b）などで表し，ふつうのカッコの優先順位で計算するものとして，
（（a, b）,（c, d）），（（（a, b）, c）, d），（a,（b,（c, d））），
（（a,（b, c））, d），（a,（（b, c）, d））の５パターンである。１パターンあたりの２項演算の数は３であるから，パターンあたりの総数は$5^3$通りである。

（４）上記の各数字や２項演算の結果に対して５種の単項演算を任意の階数だけ作用することができる。単項演算のエントリーポイントは１パターンにつき$n$個あるとする。単項演算子の種類を$m$個とし，エントリーポイントに$k$個まで入るならば，パターンあたりの配位の数は，$\sum_{j=0}^k {}_{n+j-1}C_j m^k$となる。

したがって，この他に条件をつけなければ，総数で $5 \times 5^3 \times (5^k \dfrac{n^k}{k!}+\cdots) $個の式が得られる。これをしらみつぶしに調べるのだろうか・・・

（Juliaでパズル（２）に続く）

2019年1月23日水曜日

円周率が22/7より小さいことの証明

というタイトルのWikipedia記事がある（円周率が22/7より小さいことの証明）。

証明は次の通りである。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^4(1-x)^4}{1+x^2} dx \\
& = \int_0^1\dfrac{x^4-4x^5+6x^6-4x^7+x^8}{1+x^2} dx\\
&= \int_0^1x^6 -4x^5 + 5x^4 -4x^2 +4 -\dfrac{4}{1+x^2} dx\\
&=\frac{x^7}{7}-\frac{2 x^6}{3}+x^5-\frac{4 x^3}{3}+4 x-4 \tan^{-1} x\ \Big|_0^1\\
&= \frac{1}{7}-\frac{2}{3}+1-\frac{4}{3}+4-\pi\\
&=\frac{22}{7}-\pi
\end{aligned}
\end{equation}

同様にして円周率が355/113より小さいことも証明される。
\begin{equation}
\begin{aligned}
0 &< \int_0^1 \dfrac{x^8(1-x)^8(25+816x^2)}{3164(1+x^2)} dx \\
&= \cdots\\
&= \frac{12 x^{17}}{791}
-\frac{102 x^{16}}{791}
+\frac{3151 x^{15}}{6780}
-\frac{703 x^{14}}{791}
+\frac{393 x^{13}}{452}
-\frac{23 x^{12}}{113}\\
&-\frac{145 x^{11}}{452}
-\frac{5 x^{10}}{791}
+\frac{1409 x^9}{3164}
-\frac{4 x^7}{7}
+\frac{4 x^5}{5}
-\frac{4 x^3}{3}
+4 x−4 \tan^{-1}x \Bigr|_0^1 \\
&= \cdots\\
&=\frac{355}{113}-\pi
\end{aligned}
\end{equation}

えーっ，3164と816と25はどこから連れてきたのだ・・・という方は上記記事から参考文献をたどってみよう。この世界にはまだまだ自分が知らないことが沢山埋もれている。

2019年1月22日火曜日

角運動量の合成への道（４）

（角運動量の合成への道（３）からの続き）

次に，スピン１の粒子の軌道角運動量$\boldsymbol{L}$とスピン角運動量$\boldsymbol{T}$の合成を考える。$\boldsymbol{J}=\boldsymbol{L}+\boldsymbol{T}$であり，前回までの手順と同様に，$(\boldsymbol{J}^2, J_z) | J M \rangle = (J(J+1)\hbar^2, M\hbar ) | J M \rangle$ となる状態は，$| \ell m \rangle | t \rangle$ の線形結合で表されるとする。

（１）$J_z$を$| \ell m \rangle | t \rangle$に作用させる。
$J_z | \ell m \rangle | t \rangle = (m + t) \hbar | \ell m \rangle | t \rangle = M | \ell m \rangle | t \rangle$。したがって，$m=M-t$となる。

（２）$\boldsymbol{J}^2$を$| \ell M-t \rangle | t \rangle$に作用させる。
\begin{equation}
\begin{aligned}
(\boldsymbol{L}^2+\boldsymbol{T}^2+2L_zT_z &+L_{+}T_{-}+L_{-}T_{+})| \ell M-t \rangle | t \rangle \\
= & \sqrt{\ell(\ell+1)-m(m-1)} \sqrt{2-t(t+1)}\ \hbar^2 | \ell m-1 \rangle | t+1 \rangle\\
+ & (\ell(\ell+1)+2+2 m t )\ \hbar^2 | \ell m \rangle | t \rangle \\
+ & \sqrt{\ell(\ell+1)-m(m+1)}\sqrt{2-t(t-1)}\ \hbar^2 | \ell m+1 \rangle | t-1 \rangle
\end{aligned}
\end{equation}
ここで，$\boldsymbol{J}^2$の固有値を決める固有値方程式は次のようになる。
\begin{equation}
\boldsymbol{J}^2 \sum_{t=-1}^1 \alpha_t | \ell m-t \rangle | t \rangle
= \lambda \sum_{t=-1}^1 \alpha_t | \ell m-t \rangle | t \rangle
\end{equation}
ここで$\ell(\ell+1)=j$, $m(m\pm 1)=\mu_{\pm}$と略記する。
\begin{equation}
\begin{pmatrix}j+2m & \sqrt{2(j-\mu_{-})} & 0 \\
\sqrt{2(2(j-\mu_{-})} & j+2 & \sqrt{2(j-\mu_{+})}\\
0& \sqrt{2(j-\mu_{+})}&j-2m \end{pmatrix}
\begin{pmatrix} \alpha_{+1} \\ \alpha_0 \\ \alpha_{-1}\end{pmatrix}
= \lambda \begin{pmatrix} \alpha_{+1} \\ \alpha_0 \\ \alpha_{-1} \end{pmatrix}
\end{equation}
この固有値方程式を解くと，$\lambda=(\ell+1)(\ell+2), \ell(\ell+1), \ell(\ell-1)$となる。
また，固有ベクトルの係数は，
\begin{array}{c|ccc}
\lambda & \alpha_{+1} & \alpha_{0} & \alpha_{-1} \\
\hline
\ell(\ell-1)
& \sqrt{\frac{(\ell-m)(\ell-m+1)}{2\ell(2\ell+1)}}
& -\sqrt{\frac{(\ell+m)(\ell-m)}{\ell(2\ell+1)}}
& \sqrt{\frac{(\ell+m)(\ell+m+1)}{2\ell(2\ell+1)}} \\
\hline
\ell(\ell+1)
& -\sqrt{\frac{(\ell+m)(\ell-m+1)}{2\ell(\ell+1)}}
& \sqrt{\frac{m}{\ell(\ell+1)}}
&\sqrt{\frac{(\ell-m)(\ell+m+1)}{2\ell(\ell+1)}} \\
\hline
(\ell+1)(\ell+2)
& \sqrt{\frac{(\ell+m)(\ell+m+1)}{2(\ell+1)(2\ell+1)}}
& \sqrt{\frac{(\ell+m+1)(\ell-m+1)}{(\ell+1)(2\ell+1)}}
&\sqrt{\frac{(\ell-m)(\ell-m+1)}{2(\ell+1)(2\ell+1)}} \\
\hline
\end{array}

2019年1月21日月曜日

角運動量の合成への道（３）

（角運動量の合成への道（２）からの続き）

２つの独立な自由度からなる系の別の例として，１つの粒子が軌道角運動量$\boldsymbol{L}$とスピン角運動量$\boldsymbol{S}$を持つ場合を考える。これらを合成した全角運動量を $\boldsymbol{J}=\boldsymbol{L}+\boldsymbol{S}$とする（２粒子系の場合と同様，より正確には通常の３次元空間とスピン空間に作用する演算子のテンソル積 $\boldsymbol{J}=\boldsymbol{L}\otimes \boldsymbol{1}_S+\boldsymbol{1}_L\otimes\boldsymbol{S}$）。

軌道角運動量$\boldsymbol{L}$の固有状態を$|\ell m \rangle$，スピン角運動量$\boldsymbol{S}$の固有状態を$| s \rangle$とし，その積$|\ell m \rangle | s \rangle$を考える。このとき，次の関係が成り立っている。
\begin{equation}
\begin{aligned}
\boldsymbol{L}^2 |\ell m\rangle = \ell(\ell+1)\hbar^2 |\ell m\rangle, \quad &
\boldsymbol{S}^2 | s\rangle = \frac{3}{4} \hbar^2 |s \rangle, \\
L_z |\ell m\rangle = m \hbar |\ell m\rangle, \quad &
S_z | s\rangle = s \hbar |s \rangle, \\
L_\pm |\ell m\rangle = \sqrt{\ell(\ell+1)-m(m\pm 1)} \ \hbar |\ell m \pm 1 \rangle, \quad &
S_\pm | s\rangle = \sqrt{\frac{3}{4}-s(s\pm 1)} \ \hbar |s \pm 1 \rangle
\end{aligned}
\end{equation}
また，合成角運動量$\boldsymbol{J}$も一般の角運動量の交換関係，$[J_i, J_j]=i\hbar \epsilon_{ijk} J_k$を満足するので，$\boldsymbol{J}^2$と$J_z$の同時固有状態を$|j\mu \rangle$とすると，次の関係が成り立つ。
\begin{equation}
\boldsymbol{J}^2 |j \mu \rangle = j(j+1)\hbar^2 |j \mu \rangle, \quad
J_z |j \mu\rangle = \mu \hbar |j \mu\rangle
\end{equation}
$[\boldsymbol{J}^2,\boldsymbol{L}^2]=0, [\boldsymbol{J}^2,\boldsymbol{S}^2]=0, [J_z,L_z]=0,[J_z,S_z]=0$などが成立するので，$|j\mu\rangle$を$|\ell m \rangle | s \rangle$から構成することができる（$[\boldsymbol{J}^2,L_z]\neq 0$や$[\boldsymbol{J}^2,S_z]\neq 0$なので，一般には単独の$|\ell m \rangle | s \rangle$は$\boldsymbol{J}^2$の固有状態ではない）。

（１）$J_z$を$|\ell m \rangle | s \rangle$に作用させる。
\begin{equation}
J_z |\ell m \rangle | s \rangle = (L_z+S_z) |\ell m \rangle | s \rangle = (m+s)\hbar |\ell m \rangle | s \rangle
\end{equation}
したがって，$|\ell m \rangle | s \rangle$は，$J_z$の固有値$\mu\hbar = (m+s)\hbar$の固有状態である。以下，$m=\mu-s$と表記することもある。

（２）$\boldsymbol{J}^2$を$|\ell \mu-s \rangle | s \rangle$に作用させる。
\begin{equation}
\begin{aligned}
\boldsymbol{J}^2 &= \boldsymbol{L}^2+ \boldsymbol{S}^2 + 2 \boldsymbol{L}\cdot\boldsymbol{S} = \boldsymbol{L}^2+ \boldsymbol{S}^2 + 2 L_z S_z + L_{+} S_{-} + L_{-} S_{+} \\
\boldsymbol{J}^2 |\ell \mu-s \rangle |s \rangle &= \{ \ell(\ell+1)+\frac{3}{4} + 2 (\mu-s) s \} \hbar^2 |\ell \mu-s \rangle |s \rangle \\
&+ \sqrt{\ell(\ell+1)-(\mu-s)(\mu+s)} \ \hbar^2 |\ell \mu+s \rangle | -s \rangle\\
&= \{ j^2+2s\mu \}\ \hbar^2 |\ell \mu-s \rangle |s \rangle
+ \sqrt{j^2-\mu^2}\ \hbar^2 |\ell \mu+s \rangle | -s \rangle\\
\end{aligned}
\end{equation}
ここで，$j=\ell+\frac{1}{2}$とし，$s^2=\frac{1}{4}$を用いた。これから，$\boldsymbol{J}^2$の固有状態を構成するために，$|\mathscr{j} \mu\rangle = \alpha |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle + \beta |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle$ とする。ただし規格化条件より$|\alpha|^2+|\beta|^2 = 1$である。この２つの独立な状態に対する固有値方程式は $\boldsymbol{J}^2 | \mathscr{j} \mu\rangle = \lambda \hbar^2 | \mathscr{j} \mu\rangle$ であり，$\alpha, \beta$を用いると次のように行列形式で表される。
\begin{equation}
\begin{pmatrix} j^2+\mu & \sqrt{j^2-\mu^2} \\ \sqrt{j^2-\mu^2} & j^2-\mu \end{pmatrix}
\begin{pmatrix} \alpha \\ \beta \end{pmatrix}
= \lambda \begin{pmatrix} \alpha \\ \beta \end{pmatrix}
\end{equation}
この$\alpha, \beta$に対する連立方程式が自明でない解を持つための条件は，右辺を移行して0にしたときの左辺の行列式が0になることであり，$(j^2+\mu -\lambda)(j^2-\mu -\lambda)-(j^2-\mu^2)=0$という$\lambda$の２次方程式になる。$\therefore (\lambda -j(j+1))(\lambda -j(j-1)) = 0$から$\boldsymbol{J}^2$の固有値$\lambda \hbar^2$は，$j(j+1)\hbar^2$と$j(j-1)\hbar^2$である（したがって，$\mathscr{j}=j$ と $j-1$）。

（３）固有状態の構成
上記の$\alpha$と$\beta$の連立方程式に，固有値$\lambda$を代入して，規格化条件とともに$\alpha$と$\beta$を求めると次のようになる。
\begin{array}{c|cc||cc}
\lambda & \alpha & \beta & \alpha & \beta\\
\hline
j(j+1) & \sqrt{\dfrac{j+\mu}{2j}} & \sqrt{\dfrac{j-\mu}{2j}}
&\sqrt{\dfrac{\ell+m}{2\ell+1}} & \sqrt{\dfrac{\ell-m+1}{2\ell+1}} \\
\hline
j(j-1) & \sqrt{\dfrac{j-\mu}{2j}} & -\sqrt{\dfrac{j+\mu}{2j}}
&\sqrt{\dfrac{\ell-m+1}{2\ell+1}} & -\sqrt{\dfrac{\ell+m}{2\ell+1}}
\end{array}
まとめると，
\begin{equation}
\begin{array}{l}
\begin{array}{|l|}
\hline

\quad |\ j\quad\ \mu\rangle = \sqrt{\dfrac{j+\mu}{2j}} |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle+\sqrt{\dfrac{j-\mu}{2j}} |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle \quad (\lambda = j(j+1)) \quad \\
\quad |j-1\mu\rangle = \sqrt{\dfrac{j-\mu}{2j}} |\ell \mu-\frac{1}{2} \rangle |\frac{1}{2} \rangle-\sqrt{\dfrac{j+\mu}{2j}} |\ell \mu+\frac{1}{2} \rangle | -\frac{1}{2} \rangle \quad (\lambda = j(j-1) ) \quad \\

\hline
\end{array} \\
\end{array}
\end{equation}

（角運動量の合成へ道（４）に続く）

2019年1月20日日曜日

角運動量の合成への道（２）

（角運動量の合成への道（１）の続き）

２つの独立な自由度がある系の例として，スピン$\frac{1}{2}$を持った２粒子の系を考え，この系の合成スピン$\boldsymbol{S}=\boldsymbol{S}_1+\boldsymbol{S}_2$を考える（より正確にいうと，２つの状態ベクトルのテンソル積の状態に作用する演算子のテンソル積を考えるので，$\boldsymbol{S}=\boldsymbol{S}_1 \otimes \boldsymbol{1}_2+\boldsymbol{1}_1 \otimes \boldsymbol{S}_2$が $|s_1\rangle_1 \otimes |s_2\rangle_2$ に作用する。ここでは，演算子は $\boldsymbol{S}=\boldsymbol{S}_1+\boldsymbol{S}_2$とし，状態の方も，$|s_1\rangle |s_2\rangle$と表記する）。

スピン$\frac{1}{2}$演算子の固有状態を$|\frac{1}{2}s\rangle$ として$\frac{1}{2}$を省略すると，２粒子状態は，$|s_1\rangle |s_2\rangle$と表わされ，次の関係が成立する（$ |\pm \frac{3}{2} \rangle = 0$ とする）。
\begin{equation}
\begin{aligned}
\boldsymbol{S}_1^2 |s_1\rangle |s_2 \rangle = \frac{3}{4}\hbar^2 |s_1\rangle |s_2 \rangle,
\quad & \boldsymbol{S}_2^2 |s_1\rangle |s_2 \rangle = \frac{3}{4}\hbar^2 |s_1\rangle |s_2 \rangle \\
S_{1z} |s_1\rangle |s_2 \rangle = s_1\hbar |s_1\rangle |s_2 \rangle,\
\quad & S_{2z} |s_1\rangle |s_2 \rangle = s_2 \hbar |s_1\rangle |s_2 \rangle \\
S_{1\pm} |s_1 \rangle |s_2 \rangle = \hbar |s_1 \pm 1 \rangle |s_2 \rangle,\quad & S_{2\pm} |s_1\rangle |s_2 \rangle = \hbar |s_1\rangle |s_2 \pm 1 \rangle \\
\end{aligned}
\end{equation}

$\boldsymbol{S}=\boldsymbol{S}_1+\boldsymbol{S}_2$については，一般の角運動量の交換関係$[S_i, S_j]=i\hbar\epsilon_{ijk}S_k$が成立するので，合成スピン角運動量の大きさとz成分の同時固有状態$|S M\rangle$とすると，次の関係が成り立つ。
\begin{equation}
\begin{aligned}
\boldsymbol{S}^2 |S M\rangle &= S(S+1) \hbar^2 |S M\rangle \\
S_z |S M\rangle &= M \hbar |S M\rangle
\end{aligned}
\end{equation}
一方，$[\boldsymbol{S}^2, \boldsymbol{S}_i^2]=0$, $[S_z, S_{iz}]=0 $などが成り立つので，$|S M\rangle$を２粒子状態から構成することができる（ただし，$[\boldsymbol{S}^2, S_i]\neq0$ なので，一般の $|s_1\rangle |s_2\rangle$ は $\boldsymbol{S}^2$ の固有状態とは限らない）。

（１）$S_z$ を $|s_1\rangle |s_2\rangle$ に作用させる。
\begin{equation}
S_z |s_1\rangle |s_2\rangle = (S_{1z}+S{2z}) |s_1\rangle |s_2\rangle = (s_1+s_2)\hbar |s_1\rangle |s_2\rangle
\end{equation}
したがって， $|s_1\rangle |s_2\rangle$ は，$S_z$ の固有値 $M \hbar =(s_1+s_2)\hbar$ の固有状態である。

（２）$\boldsymbol{S}^2$ を $|s_1\rangle |s_2\rangle$ に作用させる。
　　（注: $s_i \neq \pm\frac{1}{2}$の場合は状態ベクトル$ |s_i\rangle$ は0とする）
\begin{equation}
\begin{aligned}
\boldsymbol{S}^2 &= \boldsymbol{S}_1^2+ \boldsymbol{S}_2^2 + 2 \boldsymbol{S}_1\cdot\boldsymbol{S}_2 = \boldsymbol{S}_1^2+ \boldsymbol{S}_2^2 + 2 S_{1z}S_{2z} + S_{1+} S_{2-} + S_{1-} S_{2+} \\
\boldsymbol{S}^2 |s_1\rangle |s_2\rangle &= (\frac{3}{2}\hbar^2 + 2 s_1 s_2 \hbar^2) |s_1\rangle |s_2\rangle + \hbar^2 |s_1+1\rangle |s_2-1\rangle + \hbar^2 |s_1-1\rangle |s_2+1\rangle
\end{aligned}
\end{equation}
簡単のために，$|s \rangle \rightarrow |+\rangle$ ($s=\frac{1}{2}$の場合), $|-\rangle$ ($s=-\frac{1}{2}$の場合) と表記すると，
\begin{equation}
\begin{aligned}
\boldsymbol{S}^2 | + \rangle | + \rangle &= 2\hbar^2 | + \rangle | + \rangle\\
\boldsymbol{S}^2 | + \rangle | - \rangle &= \hbar^2 | + \rangle | - \rangle + \hbar^2 | - \rangle | + \rangle \\
\boldsymbol{S}^2 | - \rangle | + \rangle &= \hbar^2 | - \rangle | + \rangle + \hbar^2 | + \rangle | - \rangle \\
\boldsymbol{S}^2 | - \rangle | - \rangle &= 2\hbar^2 | - \rangle | - \rangle
\end{aligned}
\end{equation}
したがって，次のようにして$\boldsymbol{S}^2$の固有状態を構成できる。
\begin{equation}
\begin{aligned}
\boldsymbol{S}^2 ( | + \rangle | - \rangle + | - \rangle | + \rangle) &= 2\hbar^2 ( | + \rangle | - \rangle + | - \rangle | + \rangle)\\
\boldsymbol{S}^2 ( | + \rangle | - \rangle - | - \rangle | + \rangle) &= 0 ( | + \rangle | - \rangle - | - \rangle | + \rangle)\\
\end{aligned}
\end{equation}
まとめると，$\boldsymbol{S}^2, S_z$の固有値 $S(S+1)\hbar^2, M\hbar$の同時固有状態 $ | S M \rangle $に対して，
\begin{equation}
\begin{array}{l}
\begin{array}{|l|}
\hline

\quad | 1 +1 \rangle = |+ \rangle |+ \rangle \quad \\
\quad | 1 \quad 0 \rangle = \frac{1}{\sqrt{2}} ( |+ \rangle |- \rangle + |- \rangle |+ \rangle ) \quad \\
\quad | 1 -1 \rangle = |- \rangle |- \rangle \quad \\
\quad | 0 \quad 0 \rangle = \frac{1}{\sqrt{2}} ( |+ \rangle |- \rangle - |- \rangle |+ \rangle ) \quad\\

\hline
\end{array} \\
\end{array}
\end{equation}
$ | 1 M \rangle $がスピン３重項（トリプレット）状態，$| 0 0 \rangle $がスピン１重項（シングレット）状態である。

「人の世の窓打ちにけり冬の雨」（西嶋あさ子　1938-）

（角運動量の合成への道（３）に続く）

2019年1月19日土曜日

角運動量の合成への道（１）

（MathJaxの記号一覧は，このEasy Copy MathJax が便利である。）

　角運動量演算子については，次の交換関係から出発して，固有値と固有状態の議論ができる。一般化された角運動量演算子$\boldsymbol{J}$はエルミート演算子であり，その成分は交換関係$[J_i, J_j] = i \hbar \epsilon_{ijk}J_k$を満足すると仮定する。このとき，$\boldsymbol{J}^2=\sum_{i} J_i^2=\sum_i J_i^\dagger J_i$も固有値が正のエルミート演算子となり，$[\boldsymbol{J}^2, J_i]=0$であることから，$\boldsymbol{J}^2$と$J_z$の同時固有状態$|\lambda \mu \rangle$が存在し，その固有値$\lambda \hbar^2, \mu \hbar $は実数となる（$\lambda \ge 0$）。
\begin{equation}
\begin{aligned}
\boldsymbol{J}^2 |\lambda \mu \rangle &= \lambda \hbar^2 |\lambda \mu \rangle\\
J_z |\lambda \mu \rangle &= \mu \hbar |\lambda \mu \rangle
\end{aligned}
\end{equation}
$J_x^\dagger J_x + J_y^\dagger J_y = \boldsymbol{J}^2- J_z^2$の左辺の $|\lambda \mu \rangle$による期待値が正であることから，固有値 $\lambda - \mu^2 \ge 0$が成り立つ。すなわち，$-\sqrt{\lambda} \le \mu \le \sqrt{\lambda}$である。

また，昇降演算子，$J_{\pm}\equiv J_x \pm i J_y$を定義すると，$J_{\pm}\dagger=J_{\mp}$, $[\boldsymbol{J}^2, J_{\pm}]=0$,$ [J_z, J_{\pm}]=\pm \hbar J_{\pm}$, $ [J_+, J_-]=2\hbar J_z$などが成り立つ。

ここで，$J_{\pm} |\lambda \mu \rangle$ がどんな状態かを調べると次のことがわかる。
\begin{equation}
\begin{aligned}
\boldsymbol{J}^2 J_{\pm} |\lambda \mu \rangle &= \lambda \hbar^2 J_{\pm} |\lambda \mu \rangle\\
J_z J_{\pm}|\lambda \mu \rangle &= (\mu\pm1) \hbar J_{\pm} |\lambda \mu \rangle
\end{aligned}
\end{equation}
すなわち，$J_{\pm} |\lambda \mu \rangle = C_{\lambda\mu}^{\pm}\hbar |\lambda \mu\pm 1 \rangle$。ここで，$C_{\lambda\mu}^{\pm}$は比例定数であり，規格化条件を用いて，$\langle \lambda \mu| J_{\mp} J_{\pm} |\lambda \mu \rangle = |C_{\lambda\mu}^{\pm}|^2 \langle \lambda \mu\pm 1 |\lambda \mu\pm 1 \rangle$ より，$C_{\lambda\mu}^{\pm} = \sqrt{\lambda -\mu(\mu\pm 1)}$
ただし，$J_\mp J_\pm = \boldsymbol{J}^2-J_z^2 \mp \hbar J_z$であることに注意する。

ある固有状態から出発して，昇降演算子$J_{\pm}$を繰り返して作用すると，$\boldsymbol{J}^2$の固有値を共有し，$J_z$の固有値が離散的に変化する一連の固有状態のシリーズが得られるが，これは，固有値$\mu$に対する条件と矛盾することから，$\mu$の上限$\mu_max$と下限$\mu_min$においては比例定数$C_{\lambda\mu}^{\pm}$が0となって，シリーズが中断される必要がある。すなわち，
\begin{equation}
\begin{aligned}
J_{+} |\lambda \mu_{max} \rangle = 0, \quad &C_{\lambda\mu_{max}}^+ = \sqrt{\lambda-\mu_{max}(\mu_{max}+1)} = 0\\
J_{-} |\lambda \mu_{min} \rangle = 0, \quad &C_{\lambda\mu_{min}}^- = \sqrt{\lambda-\mu_{min}(\mu_{min}-1)} = 0
\end{aligned}
\end{equation}
これから，$\mu_{max}(\mu_{max}+1)= \mu_{min}(\mu_{min}-1)$が成り立ち，因数分解すると，$(\mu_{max}-\mu_{min}+1)(\mu_{max}+\mu_{min})=0$となる。$\mu_{max}$と$\mu_{min}$の差は整数$n=0,1,2,\cdots$となることから，$\frac{n}{2}=j$とおいて，この半整数 $j=0, \frac{1}{2}, 1, \frac{3}{2}, 2, \cdots$に対して，$\mu_{max}=j, \mu_{min}=-j, \lambda=j(j+1)$となる。そこでこの場合の$\mu$を$m=-j, -j+1, \cdots j-1, j$と書くことにする。

そこで，固有状態と固有値を表すシンボルを$\lambda, \mu$から$j, m$に変えてまとめると，
\begin{equation}
\begin{array}{l}
\begin{array}{|l|}
\hline

\quad \boldsymbol{J}^2 | j m \rangle &= j(j+1) \hbar^2 | j m \rangle \quad \\
\quad J_z | j m \rangle &= m \hbar | j m \rangle \quad \\
\quad J_{\pm} | j m \rangle &= \sqrt{j(j+1)-m(m \pm 1)}\ \hbar | j m \pm1 \rangle \quad \\

\hline
\end{array} \\
\end{array}
\end{equation}

（角運動量の合成への道（２）に続く）

2019年1月18日金曜日

日本古典籍データセット

（Yahoo!知恵袋データ(第3版) からの続き）

国立情報学研究所（NII）の情報研究データリポジトリの末尾には，提供を完了したデータセットとして，国文研データセットがあげられており，情報・システム研究機構人文学オープンデータ共同利用センターに移管されたとされている。

それがこちら，人文学オープンデータ共同利用センターの日本古典籍データセットである。このセンターにはこれ以外にも下記のような７種のデータセット等が提供されている。日本古典籍データセットの方は，2017年12月現在で（最近は滞っているのだろうか・・・）1767点（329,702コマ）のデータが集積されている。すべてのzipデータを合計すると665GB以上とのことで取り扱い注意だ。その内容は以下のとおり。

古典籍画像データ　
　各作品の画像データをJPEG形式で保存したものです。国文学分野のほか、国文学研究資料館が所蔵する医学や理学、産業など多分野の古典籍、さらに味の素食の文化センターが所蔵する料理本等で、国文学研究資料館が撮影した古典籍を含みます。

書誌データ　
　各作品の書誌データをCSV形式でまとめたものです。国文学研究資料館で公開している「新日本古典籍総合データベース」より、書誌ID/書名/著者名/巻数/刊写の別/出版事項/形態/注記などを抽出したものとなっています。なお一部の作品には国文学研究資料館にて付与した略解題も含まれています。

本文テキストデータ　
　翻刻した本文テキストデータをプレーンテキストまたはDOCX形式で保存したものです。一部の作品に限ります。

タグデータ　
　国文学研究資料館で付与作業を行っている、1枚1枚の画像に対する文中の固有名詞のタグ情報をCSV形式でまとめたものです。一部の作品に限ります。

こちらはライセンスが，

なので，ちょっと安心である。
本文テキストデータを充実させてほしいかな。浄瑠璃の床本データとか。

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1) 江戸料理レシピデータセット
(2) 日本古典籍くずし字データセット
(3) KMNISTデータセット
(4) 武鑑全集
(5) 顔貌コレクション
(6) 近代雑誌データセット
(7) Geoshapeリポジトリ

P. S. 「国文研の古典籍オープンデータセットであそぶ」

2019年1月17日木曜日

Yahoo!知恵袋データ（第3版）

国立情報学研究所（NII）が，情報学研究データリポジトリとして，様々なデータセットを公開している（下記リスト参照）。最近公開されたのが，Yahoo!知恵袋データ（第3版）だ。Yahoo!知恵袋は，Yahoo株式会社が2004年から公開している利用者参加によるＱ＆Ａ型のコミュニケーションサイトで，現時点の質問総数が 2.0 億件，回答総数が 4.9 億件，利用登録数が 0.47 億件である。ここには，研究機関への研究データ提供に関する注意も掲載されている。

「Yahoo! JAPANでは投稿者のYahoo! JAPAN IDを暗号化するなど、個人を特定することができない情報に処理したうえで投稿内容、投稿日時などの投稿に関する情報を大学、独立行政法人などの研究機関に提供します。Yahoo! JAPANが提供する情報によって、当該大学、独立行政法人などが投稿者が誰であるかを知ることはありません。」

一般に，Yahoo!知恵袋の回答の信頼度はあまり高くないといわれているような気がするが（実証データを見たことはない），ある分野でどんな質問が出されており，何が問題とされているのかについては，一定の知見を得ることができる。これをある種のｅラーニングシステムや，自学自習システムのベースとして利用することは可能かもしれない。

１０年ほど前に教育学部理科の卒業研究で「小学校教員のための理科質問検索システムの構築について」というテーマを取り上げたときに，質問のデータベース構築のために，インターネット上のＱ＆Ａシステムでよく取り上げられる質問を参考にしたり，附属学校で子どもたちにアンケートをとったりしたことがあったのを思い出した。

ところで，Yahoo!知恵袋データ（第3版）の利用条件はかなりきつい。基本的に大学等の研究機関の研究者に限定されているため，企業との共同研究なども難しそうである。

「申請の単位は大学の研究室等とし，申請者は研究室等を代表する常勤の職員（大学の場合は教員等，研究機関の場合は研究員等）の方として下さい。また，契約締結者は法人を代表してデータ利用契約書を締結する権限のある方で，データ利用に関する同意書に押印いただく公印をお持ちの方（通常は学部長相当以上の方）としてください。
　本データを利用できるのは，申請者自身および申請者と同一機関に所属（職員，在学生など正規の在籍に限る）し，かつ同一研究室等において申請者の管理の下で研究活動を行う方のみです。たとえ共同研究で使用する場合あっても，他機関や他の独立した研究室の方が使用する場合は，別途に申請して下さい。」となっている。

Yahoo!知恵袋の場合は年齢情報はないが，学習者の抱く疑問（Ｑ＆ＡのＱのみで良いかもしれない）のデータベースが，完全なオープン・リソースとして，学校教育の各分野や領域で整備されるとおもしろいことができそうな気がするのだが・・・

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
(1) Yahoo!データセット
　・Yahoo!知恵袋データ（第3版）
(2) 楽天データセット
　・楽天市場の全商品データ，レビューデータ
　・楽天トラベルの施設データ，レビューデータ
　・楽天GORAのゴルフ場データ，レビューデータ
　・楽天レシピのレシピ情報，レシピ画像
　・PriceMinisterのユーザレビュー，レビュー有効性情報アノテーション付きデータ
(3) ニコニコデータセット
　・ニコニコ動画コメント等データ
　・ニコニコ大百科データ
(4) リクルートデータセット
　・ホットペッパービューティーデータ
(5) クックパッドデータセット
　・レシピデータ，献立データ
(6) LIFULL HOME'Sデータセット
　・賃貸物件スナップショットデータ（賃貸物件データ＋画像データ）
　・高精細度間取り図画像データ
　・賃貸・売買物件月次データ
(7) 不満調査データセット
　・投稿された不満データ，ユーザ情報
　・カテゴリ別不満特徴語辞書
(8) Sansanデータセット
　・サンプル名刺データ
(9) インテージデータセット
　・i-SSPデータ
　・みんレポデータ

（日本古典籍データセットに続く）

2019年1月16日水曜日

Juliaの四捨五入

Juliaの実数を整数に丸めたいので，a=1.5 b=Int(a) とするとダメなのであった，ここでつまづくのは何回目（カプレカ数の項をみる）。documentをちゃんと読むと，実数に対して丸め関数を作用させた後に型変換せよとある。丸め関数に型指定を含むものもある。V1.03ドキュメントの数学的操作と基本関数を参照。

round(x) round x to the nearest integer typeof(x)
round(T, x) round x to the nearest integer T
floor(x) round x towards -Inf typeof(x)
floor(T, x) round x towards -Inf T
ceil(x) round x towards +Inf typeof(x)
ceil(T, x) round x towards +Inf T
trunc(x) round x towards zero typeof(x)
trunc(T, x) round x towards zero T

こんな感じ。

2019年1月15日火曜日

Juliaで積分

（１）数式処理の場合
Juliaのバージョンが古いが，たけしの忘備録「Julia+SymPyの勝手なまとめ」に設定やその他の数式処理の例についても書いてある。最低限の部分は以下のとおりである。

import PyCall
import SymPy

@syms x
(x,)
f(x)=1/x
f (generic function with 1 method)
integrate(f)
log(x)

（２）数値計算の場合
Juliaの具体的な活用例については，黒木さんのページが最も充実していると思う。

import QuadGK

f(x)=1/x
f (generic function with 1 method)
quadgk(f, 1.0, 2.0)
(0.6931471805599453, 1.9931833961095435e-11)

2019年1月14日月曜日

Juliaにおける変数のクリア

Jupyter 上の Juliaでいろいろな作業をしていると，使用している変数や関数の名前が重なってしまうことがある。Ver. 0.6には，workspace() というコマンドがあって，変数のバインディングをクリアできたようだが，Ver. 1.0ではこの機能がなくなっている。対処法を調べてみたのだが・・・

(1) Revise.jl を用いる。ちょっと違うような感じ。

(2) a=nothing とする。ちょっと違うような感じ。

(3) 以下のコードを実行する（Startup.jl に導入），ちょっと違うような感じ。

function workspace()
　atexit() do
　　run(`$(Base.julia_cmd())`)
　end
　exit()
end

これならば，Kernelの再起動でよいのか。

2019年1月13日日曜日

献灯使：多和田葉子

講談社文庫の献灯使（多和田葉子）を読了。中編小説の「献灯使」と３つの短編小説１つの戯曲からなる短編集。「韋駄天どこまでも」に稲垣足穂のフレーバーを感じる。斎藤美奈子の岩波新書の「日本の同時代小説」の2010年代の項では，ディストピア小説に向う純文学として，佐藤友哉の「ベッドサイド・マーダーケース（未読）」，吉村萬壱の「ボラード病（未読）」と並んで多和田葉子の上記短編集の中から「不死の島」が取り上げられている。最近の純文学は半分くらいＳＦフレーバーな作品によって占められているのか・・・。

多和田葉子の献灯使は書評でかなり好意的に取り上げられていたので，期待が大きかった分だけ少し残念かもしれない。もう少し他の本も読んでみよう。

http://bkmkn.s3-website-ap-northeast-1.amazonaws.com/9784062937283/9784062937283_w.jpg

P. S. 「日本の同時代小説」に登場しなかった作家で気になったのは，三枝和子。たぶん，一定の広がりを持った作品を中心に構成しているためだろう。

P. P. S. 「日本の同時代小説」では，主に1960年代から2010年代に渡る350名程度の作家が取り上げられているが，自分が読んだことがある作家は，そのうち約1/4の90名程度か。知らない世界が多すぎる。

2019年1月12日土曜日

Juliaにおける変数のスコープ

２つのソートされた配列を比較する練習問題プログラムを実行しようとしてひっかかった。

#3 function main()
a=sort(rand(1:100,30))
b=sort(rand(1:100,30))
j=1
k=1
while j<=30 && k<=30
#2 global j
if a[j] > b[k]
#1 j=j+0
k=k+1
elseif a[j] < b[k]
j=j+1
else
println((j,a[j],k,b[k]))
j=j+1
k=k+1
end
end
#3 end
#3 main()

UndefVarError: j not defined

Stacktrace:
[1] top-level scope at ./In[25]:8

Whileの中の変数 j が未定義だというエラーである。調べてみると，これが Julia 1.0の仕様であるらしい。そこで，上記のプログラムにおいて，コメント#1，#2，#3をのいずれかをはずすとエラーは回避された。

このあたりの Scope of Variables をちゃんと読んでおけということですね。

2019年1月11日金曜日

元号の予想（４）

（元号の予想（３）の続き）

こんなことをしていていいのだろうか？新元号の予想を，過去の日本の元号に用いられた既出漢字の組み合わせで考えるという課題。可能な組み合わせから，既存の各国の元号を消去するコードをJuliaで書いてみた（既存の元号と重なる組み合わせは出現率を０にリセット）。

日本の元号との文字列比較をすることから，中国の元号は簡体字を含んだ中国語版ではなく，日本語版Wikipediaの元号一覧（中国）に変更した（満州国の２つの元号を含む）。元号の数が約1200から約840に縮小した。両言語版でこれだけ違うため，採用されている範囲や基準が統一されていないと思われる。また，信頼性も不十分であることに留意しよう。これに，元号一覧（朝鮮），元号一覧（ベトナム），元号一覧（台湾）を加えたセットを考える。必要に応じて，漢字２文字セットで表現される他の制約条件を簡単に追加することができる。

P. S. 昭和の「昭」や平成の「成」は日本の元号使用漢字としては初出で採用されているので，頻度リストにこだわることにはほとんど意味がないことを改めて注意する。

function gfrq(c,n,r,m)
#
# c input array of kanji char
# n size of c
# r output array of sorted (kanji, count)
# m size of r
#
sort!(c)
k=0
for j=2:n
k=k+1
if c[j]!=c[j-1]
push!(r,(c[j-1],k))
k=0
end
end
sort!(r, by=x->x[2], rev=true)
m=length(r)
# println((m,r))
end

#japan
g1="大化,白雉,朱鳥,大宝,慶雲,和銅,霊亀,養老,神亀,天平,感宝,勝宝,宝字,神護,景雲,宝亀,天応,延暦,大同,弘仁,天長,承和,嘉祥,仁寿,斉衡,天安,貞観,元慶,仁和,寛平,昌泰,延喜,延長,承平,天慶,天暦,天徳,応和,康保,安和,天禄,天延,貞元,天元,永観,寛和,永延,永祚,正暦,長徳,長保,寛弘,長和,寛仁,治安,万寿,長元,長暦,長久,寛徳,永承,天喜,康平,治暦,延久,承保,承暦,永保,応徳,寛治,嘉保,永長,承徳,康和,長治,嘉承,天仁,天永,永久,元永,保安,天治,大治,天承,長承,保延,永治,康治,天養,久安,仁平,久寿,保元,平治,永暦,応保,長寛,永万,仁安,嘉応,承安,安元,治承,養和,寿永,元暦,文治,建久,正治,建仁,元久,建永,承元,建暦,建保,承久,承応,元仁,嘉禄,安貞,寛喜,貞永,天福,文暦,嘉禎,暦仁,延応,仁治,寛元,宝治,建長,康元,正嘉,正元,文応,弘長,文永,建治,弘安,文永,建治,弘安,正応,永仁,正安,乾元,嘉元,徳治,延慶,応長,正和,文保,元応,元亨,正中,嘉暦,元徳,元弘,正慶,建武,延元,興国,正平,建徳,文中,天授,弘和,元中,暦応,康永,貞和,観応,文和,延文,康安,貞治,応安,永和,康暦,永徳,至徳,嘉慶,康応,明徳,応永,正長,永享,嘉吉,文安,宝徳,享徳,康正,長禄,寛正,文正,応仁,文明,長享,延徳,明応,文亀,永正,大永,享禄,天文,弘治,永禄,元亀,天正,文禄,慶長,元和,寛永,正保,慶安,承応,明暦,万治,寛文,延宝,天和,貞享,元禄,宝永,正徳,享保,元文,寛保,延享,寛延,宝暦,明和,安永,天明,寛政,享和,文化,文政,天保,弘化,嘉永,安政,万延,文久,元治,慶応,明治,大正,昭和,平成,"
#china
g2="建元,元光,元朔,元狩,元鼎,元封,太初,天漢,太始,征和,後元,始元,元鳳,元平,本始,地節,元康,神爵,五鳳,甘露,黄龍,初元,永光,建昭,竟寧,建始,河平,陽朔,鴻嘉,永始,元延,綏和,建平,太初,元将,元寿,元始,居摂,初始,始建,天鳳,地皇,更始,建武,建武,中元,永平,建初,元和,章和,永元,元興,延平,永初,元初,永寧,建光,延光,永建,陽嘉,永和,漢安,建康,永憙,本初,建和,和平,元嘉,永興,永寿,延熹,永康,建寧,熹平,光和,中平,光熹,昭寧,永漢,初平,興平,建安,延康,黄初,太和,青龍,景初,正始,嘉平,正元,甘露,景元,咸熙,章武,建興,延熙,景耀,炎興,黄武,黄龍,嘉禾,赤烏,太元,神鳳,建興,五鳳,太平,永安,元興,甘露,宝鼎,建衡,鳳凰,天冊,天璽,天紀,泰始,咸寧,太康,太熙,永熙,永平,元康,永康,永寧,太安,永安,建武,永興,光熙,永嘉,建興,建武,大興,永昌,太寧,咸和,咸康,建元,永和,升平,隆和,興寧,太和,咸安,寧康,太元,隆安,大亨,元興,義熙,元熙,桓玄,天康,永鳳,河瑞,光興,嘉平,建元,麟嘉,漢昌,光初,太和,建平,延熙,建武,太寧,青龍,永寧,永興,元璽,光寿,建熙,建初,建興,晏平,玉衡,玉恒,漢興,太和,嘉寧,建興,和平,升平,皇始,寿光,永興,甘露,建元,太安,太初,延初,元光,燕元,建興,永康,建始,延平,青龍,建平,長楽,光始,建始,燕興,更始,昌平,建明,建平,建武,中興,建平,太上,正始,太平,太興,白雀,建初,皇初,弘始,永和,建義,太初,更始,永康,建弘,永弘,龍昇,鳳翔,昌武,真興,承光,勝光,太安,麟嘉,龍飛,承康,咸寧,神鼎,太初,建和,弘昌,嘉平,庚子,建初,嘉興,永建,神璽,天璽,永安,玄始,承玄,義和,承和,承平,真興,承陽,縁禾,太縁,永初,景平,元嘉,太初,孝建,大明,永光,景和,泰始,泰豫,元徽,昇明,建元,永明,隆昌,延興,建武,永泰,永元,中興,天監,普通,大通,大通,大同,大同,太清,大宝,天正,太始,承聖,天成,紹泰,太平,天啓,大定,天保,広運,永定,天嘉,天康,光大,太建,至徳,禎明,登国,皇始,天興,天賜,永興,神瑞,泰常,始光,神䴥,延和,太延,太平,真君,正平,承平,興安,興光,太安,和平,天安,皇興,延興,承明,太和,景明,正始,永平,延昌,熙平,神亀,正光,孝昌,武泰,建義,永安,建明,普泰,中興,太昌,永興,永熙,天平,元象,興和,武定,大統,天保,乾明,皇建,太寧,河清,天統,武平,隆化,徳昌,承光,武成,保定,天和,建徳,宣政,大成,大象,大定,開皇,仁寿,大業,義寧,皇泰,武徳,貞観,永徽,顕慶,龍朔,麟徳,乾封,総章,咸亨,上元,儀鳳,調露,永隆,開耀,永淳,弘道,嗣聖,文明,光宅,垂拱,永昌,載初,神龍,景龍,唐隆,景雲,太極,延和,先天,開元,天宝,至徳,乾元,上元,宝応,広徳,永泰,大暦,建中,興元,貞元,永貞,元和,長慶,宝暦,大和,開成,会昌,大中,咸通,乾符,広明,中和,光啓,文徳,龍紀,大順,景福,乾寧,光化,天復,天祐,天授,如意,長寿,延載,証聖,天冊,万歳,万歳,登封,万歳,通天,神功,聖暦,久視,大足,長安,仁安,大興,宝暦,中興,正暦,永徳,朱雀,太始,建興,咸和,建初,承平,義熙,甘露,章和,永平,和平,建昌,延昌,延和,義和,重光,延寿,賛普,長寿,見龍,上元,元封,応道,龍興,全義,大豊,保和,天啓,建極,法尭,貞明,大同,嵯耶,中興,安国,始元,天瑞,景星,安和,貞祐,初歴,孝治,天応,尊聖,興聖,大明,鼎新,光聖,文徳,神武,文経,至治,明徳,広徳,順徳,明政,広明,明応,明統,明聖,明治,明啓,乾興,明通,正治,聖明,天明,保安,正安,正徳,保徳,明侯,上徳,広安,上明,保定,建安,天祐,上治,天授,開明,天政,文安,日新,文治,永嘉,保天,広運,永貞,大宝,龍興,盛明,建徳,利貞,盛徳,嘉会,元亨,安定,鳳歴,元寿,天開,天輔,仁寿,道隆,利正,興正,天定,彝泰,開平,乾化,鳳歴,貞明,龍徳,同光,天成,長興,応順,清泰,天福,開運,天福,乾祐,広順,顕徳,天復,天祐,武義,順義,乾貞,大和,天祚,昇元,保大,中興,交泰,天宝,宝大,宝正,龍啓,永和,通文,永隆,天徳,乾亨,白龍,大有,光天,応乾,乾和,大宝,天復,武成,永平,通正,天漢,光天,乾徳,咸康,明徳,広政,乾祐,天会,広運,応天,同慶,天尊,中興,天興,天寿,建隆,乾徳,開宝,太平,興国,雍熙,端拱,淳化,至道,咸平,景徳,大中,祥符,天禧,乾興,天聖,明道,景祐,宝元,康定,慶暦,皇祐,至和,嘉祐,治平,熙寧,元豊,元祐,紹聖,元符,建中,靖国,崇寧,大観,政和,重和,宣和,靖康,建炎,紹興,隆興,乾道,淳熙,紹熙,慶元,嘉泰,開禧,嘉定,宝慶,紹定,端平,嘉熙,淳祐,宝祐,開慶,景定,咸淳,徳祐,景炎,祥興,神冊,天賛,天顕,会同,大同,天禄,応暦,保寧,乾亨,統和,開泰,太平,景福,重熙,清寧,咸雍,太康,太安,寿昌,乾統,天慶,保大,甘露,建福,徳興,神暦,延慶,康国,咸清,紹興,崇福,天禧,収国,天輔,天会,天眷,皇統,天徳,貞元,正隆,大定,明昌,承安,泰和,大安,崇慶,至寧,貞祐,興定,元光,正大,開興,天興,顕道,開運,広運,大慶,天授,礼法,延祚,延嗣,寧国,天祐,垂聖,福聖,承道,奲都,拱化,乾道,天賜,礼盛,国慶,大安,天安,礼定,天儀,治平,天祐,民安,永安,貞観,雍寧,元徳,正徳,大徳,大慶,人慶,天盛,乾祐,天慶,応天,皇建,光定,乾定,宝義,中統,至元,元貞,大徳,至大,皇慶,延祐,至治,泰定,致和,天順,天暦,至順,元統,至元,至正,至正,宣光,天元,洪武,建文,永楽,洪熙,宣徳,正統,景泰,天順,成化,弘治,正徳,嘉靖,隆慶,万暦,泰昌,天啓,崇禎,弘光,隆武,紹武,永暦,天命,天聡,崇徳,順治,康熙,雍正,乾隆,嘉慶,道光,咸豊,祺祥,同治,光緒,保慶,宣統,民国,洪憲,大同,康徳,"
#korea
g3="永楽,延寿,建元,開国,大昌,鴻済,建福,仁平,太和,慶雲,正開,武泰,聖冊,水徳,万歳,政開,天授,光徳,天開,開国,建陽,光武,光武,隆熙,大韓,民国,"
#vietnam
g4="天徳,太平,天福,興統,応天,景瑞,順天,天成,通瑞,乾符,有道,明道,天感,聖武,崇興,大宝,龍瑞,太平,彰聖,嘉慶,龍彰,天嗣,天貺,宝象,神武,太寧,英武,昭勝,広祐,会豊,龍符,会祥,大慶,天符,睿武,天符,慶寿,天順,天彰,宝嗣,紹明,大定,政隆,宝応,天感,至宝,貞符,天資,嘉瑞,天嘉,宝祐,治平,龍応,建嘉,天彰,有道,建中,天応,政平,元豊,紹隆,宝符,紹宝,重興,興隆,大慶,開泰,開祐,紹豊,大治,大定,紹慶,隆慶,昌符,光泰,建新,聖元,紹成,開大,興慶,重光,天慶,順天,紹平,大宝,大和,延寧,天興,光順,洪徳,景統,泰貞,端慶,洪順,光紹,統元,明徳,大正,広和,永定,景暦,光宝,淳福,崇康,延成,端泰,興治,洪寧,武安,宝定,康佑,乾統,隆泰,順徳,元和,順平,天祐,正治,洪福,嘉泰,光興,慎徳,弘定,永祚,徳隆,陽和,福泰,慶徳,盛徳,永寿,万慶,景治,陽徳,徳元,永治,正和,永盛,保泰,永慶,龍徳,永佑,景興,昭統,泰徳,光中,景盛,宝興,嘉隆,明命,紹治,嗣徳,協和,建福,咸宜,同慶,成泰,維新,啓定,保大,"
#taiwan
g5="永暦,康熙,永和,雍正,乾隆,順天,天運,嘉慶,光明,道光,咸豊,祺祥,同治,光緒,永清,大靖,民国,"
#east-asia
g=g1*g2*g3*g4*g5
println((div(length(g1),3),div(length(g2),3),div(length(g3),3),div(length(g4),3),div(length(g5),3)))
m=div(length(g),3)
c=fill((" ",0),m)
for i=1:7:m*7
c[div(i,7)+1]=(g[i]*g[i+3],div(i,7)+1)
end
#
n=div(length(g1),3)
a=fill(' ',n)
b=fill(' ',n)
for i=1:7:n*7
a[div(i,7)+1]=g1[i]
b[div(i,7)+1]=g1[i+3]
end
p=[ ]
q=[ ]
mp=0
mq=0
gfrq(a,n,p,mp)
gfrq(b,n,q,mq)
mp=length(p)
mq=length(q)
#
r=[ ]
for j in 1:mp
for k in 1:mq
if p[j][1]!=q[k][1]
push!(r,[string(p[j][1],q[k][1]),p[j][2]*q[k][2]/(mp*mq)])
end
end
end
mr=length(r)
#
for i in 1:mr, j in 1:m
if r[i][1] == c[j][1]
r[i][2]=0
end
end
sort!(r, by=x->x[2],rev=true)
for k in 1:10
　println((r[k][1],r[k][2]))
end

(250, 839, 26, 159, 17)
("嘉治", 0.0860082304526749)
("延治", 0.0860082304526749)
("永応", 0.08559670781893004)
("嘉和", 0.08148148148148149)
("寛暦", 0.08065843621399177)
("元安", 0.08024691358024691)
("文元", 0.07901234567901234)
("承治", 0.07818930041152264)
("寛安", 0.07489711934156379)
("寛応", 0.07489711934156379)
("正永", 0.07489711934156379)
("元保", 0.07407407407407407)
("天久", 0.0662551440329218)
("嘉徳", 0.06337448559670782)
("延永", 0.06337448559670782)
("長永", 0.06337448559670782)
("嘉安", 0.0588477366255144)
("延安", 0.0588477366255144)
("長応", 0.0588477366255144)
("承永", 0.05761316872427984)
("応治", 0.05473251028806585)
("延保", 0.05432098765432099)
("建応", 0.053497942386831275)
("文仁", 0.05267489711934156)
("文長", 0.05267489711934156)
("元長", 0.04938271604938271)
("天亀", 0.047325102880658436)
("寛長", 0.04609053497942387)
("宝和", 0.044444444444444446)
("正仁", 0.04279835390946502)

登録: 投稿 (Atom)