調べてみると,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)
0 件のコメント:
コメントを投稿