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)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -