ヒルベルト行列

Julia で書くと以下のようになる。

一行で書けるが決してこの書き方をすすめるわけではない。普通に2重ループで書くほうがよい。

演算子に // を使ったのは，分数型で答えを出力するためである。

julia> h(n) = [1 // (i + j - 1) for i in 1:n, j in 1:n]

h (generic function with 1 method)

julia> h(5)

5×5 Matrix{Rational{Int64}}:

 1//1  1//2  1//3  1//4  1//5

 1//2  1//3  1//4  1//5  1//6

 1//3  1//4  1//5  1//6  1//7

 1//4  1//5  1//6  1//7  1//8

 1//5  1//6  1//7  1//8  1//9

実数型に変換すると。

julia> float.(h(5))

5×5 Matrix{Float64}:

 1.0       0.5       0.333333  0.25      0.2

 0.5       0.333333  0.25      0.2       0.166667

 0.333333  0.25      0.2       0.166667  0.142857

 0.25      0.2       0.166667  0.142857  0.125

 0.2       0.166667  0.142857  0.125     0.111111

Python だと

>>> a = [[1/(i+j-1) for i in range(1, n+1)] for j in range(1, n+1)]

のようになるが，range(1, n+1) が鬱陶しい。

>>> a = [[1/(i+j+1) for i in range(n)] for j in range(n)]

とすればよい。

しかし，n = 5000 のとき，Julia で // を / にして直接実数型の結果を返す場合に比べて，Python は，ほぼ 24 倍遅い。

追記：

Nemo パッケージに hilbert 関数があった。

 

julia> using Nemo

julia> M = MatrixSpace(QQ, 5, 5)

Matrix Space of 5 rows and 5 columns over Rational Field

 

julia> h = hilbert(M)

[   1   1//2   1//3   1//4   1//5]

[1//2   1//3   1//4   1//5   1//6]

[1//3   1//4   1//5   1//6   1//7]

[1//4   1//5   1//6   1//7   1//8]

[1//5   1//6   1//7   1//8   1//9]