Julia の SymPy は，Python の sympy を利用している。Python の整数型は必要に応じて長精度整数になる。
求める式には実数が含まれるが，この式を 4 倍することで整式になる。
Rump の例題
https://www.tuhh.de/ti3/paper/rump/Ru88a.pdf
a = 77617, b = 33096 として，以下の計算を行う
R
f = function(a, b) 333.75*b^6 + a^2 *(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b)
library(gmp)
a = as.bigq(77617)
b = as.bigq(33096)
ans = f(a, b)
as.numeric(ans) # -0.8273960599468214
Julia
f(a, b) = 333.75b^6 + a^2*(11a^2*b^2 - b^6 - 121b^4 - 2) + 5.5b^8 + a/2b
a = big(77617)
b = big(33096)
f(a, b) # -0.8273960599468213681411650954798162919990331157843848199178148416727096930142628
Julia 1.8.1 が公開されました
Version 1.8.1 (2022-09-06)
Official https://julialang.org/ release
一般解は
https://stackoverflow.com/questions/69385888/parse-multi-line-result-text-file-with-awk-bash-or-r
にあった。AWK スクリプトを平易に書き変えた。
★ 以下のような定式化された入力ファイル melting.txt
The MELTING results are :
Enthalpy : -181,100 cal/mol ( -756,998 J /mol)
Entropy : -467.3 cal/mol-K ( -1,953.31 J /mol-K)
Melting temperature : 75.13 degrees C.
The MELTING results are :
Enthalpy : -170,800 cal/mol ( -713,944 J /mol)
Entropy : -444 cal/mol-K ( -1,855.92 J /mol-K)
Melting temperature : 70.6 degrees C.
からデータを取り出し，R で
Enthalpy Entropy MeltingTemperature
1 -181100 -467.3 75.13
2 -170800 -444.0 70.60
のようなデータフレームを作る。
★ 用意する AWK スクリプト melting.awk
/Enthalpy/ || /Entropy/ { # "Enthalpy" または "Entropy" を含む行について
n = split($0, a, " ") # 空白でフィールドに区切る
gsub(",", "", a[3]) # a[3] が求めるデータである。数値区切りの「,」を消去する
printf "%g, ", a[3] # 書き出す
}
/Melting temperature/ { # "Melting temperature" を含む行について
n = split($0, a, " ") # 空白でフィールドに区切る
print a[4] # 求めるデータ a[4] を書き出す
}
★ R コンソールで以下を実行する。
CSV ファイルに結果を書き出す。
以下で |> を使っているが，tidyverse ではなく Base R のパイプである。
最初の文字列がコンソールで入力されるコマンドである。（脚注参照）
"awk -f melting.awk melting.txt" |>
pipe() |>
readLines() |>
write("results.csv")
後は，ありふれた処理。CSV ファイルをデータフレームに読み込むだけである。
df = read.csv("results.csv", header=FALSE,
col.names=c("Enthalpy", "Entropy", "MeltingTemperature"))
出来上がり。
Enthalpy Entropy MeltingTemperature
1 -181100 -467.3 75.13
2 -170800 -444.0 70.60
なお，入力ファイルがちゃんと同じフォーマットであることを仮定している。
その仮定ができない場合の処理を AWK スクリプトに含めるかどうかは，程度問題。
＝＝＝＝＝＝＝＝＝＝
脚注
以下の一行をコンソールで入力すれば
$ > awk -f melting.awk melting.txt
以下が表示される
-181100, -467.3, 75.13
-170800, -444, 70.6
結果をファイルに入れるには出力を results.csv へリダイレクトすればよい。
$ > awk -f melting.awk melting.txt > results.csv
=================
なお，R だけでやろうとすると以下のようにもなるだろうか。
txt = readLines("melting.txt")
df = NULL
x = NULL
for (i in 1:length(txt)) {
if (grepl("Enthalpy", txt[i]) || grepl("Entropy", txt[i])) {
target = sub(",", "", unlist(strsplit(txt[i], " "))[4])
x = append(x, target)
} else if (grepl("Melting temperature", txt[i])) {
x = append(x, unlist(strsplit(txt[i], " "))[5])
df = rbind(df, x)
x = NULL
}
}
df = as.data.frame(df, col.names=c("Enthalpy", "Entropy", "MeltingTemperature"))
print(df)
Julia でやろうとすると以下のようになる（書式に関して別の仮定を置く）
using DataFrames
txt = readlines("melting.txt")
x = Float64[]
for s in txt
pos = findfirst(r"[-+]?[0-9,]+\.?[0-9]*", s)
isnothing(pos) == false && append!(x, parse(Float64, replace(s[pos], "," => "")))
end
df = DataFrame(Matrix(reshape(x, 3, :)'), ["Enthalpy", "Entropy", "MeltingTemperature"])
以下のような入力ファイル CodeFormat
Code
Format
Example
0
dd-mmm-yyyy HH:MM:SS
07-Sep-2000 15:38:09
1
dd-mmm-yyyy
07-Sep-2000
2
mm/dd/yy
09/07/00
3
mmm
Sep
4
m
S
からデータを読み込み
Code | Format | Example
:--|:--|:--
0 | dd-mmm-yyyy HH:MM:SS | 07-Sep-2000 15:38:09
1 | dd-mmm-yyyy | 07-Sep-2000
2 | mm/dd/yy | 09/07/00
3 | mmm | Sep
4 | m | S
のように整形して出力ファイル table に書き込む。
つまり，最終形の表は3列からなり，各列は | で区切られている。元のデータは3行で表の1行になる。
awk で書くとこんなスクリプト
{
getline a
getline b
printf "%s | %s | %s\n", $0, a, b
}
FNR == 3 {
print ":--|:--|:--"
}
こんなふうに実行する。
awk -f temp.awk < CodeFormat >table
前の記事で ```{awk} ... ``` が動かないと書いた。
https://gedevan-aleksizde.github.io/rmarkdown-cookbook/custom-engine.html
の「15.1 カスタム言語エンジンを登録する」を読んで，以下のように解決した（のかな?）。
最初に一回だけ，以下の R スクリプトを実行する。
knitr::knit_engines$set(awk = function(options) {
code <- paste(options$code, collapse = '\n')
out <- system2('awk', shQuote(code), stdout = TRUE)
knitr::engine_output(options, code, out)
})
この後，awk によるチャンクを書く。
```{awk}
BEGIN {
for(i = 1; i <= 10; i++) {
print i, sqrt(i)
}
}
```
実行すると，以下の結果を得る。
## 1 1 ## 2 1.41421 ## 3 1.73205 ## 4 2 ## 5 2.23607 ## 6 2.44949 ## 7 2.64575 ## 8 2.82843 ## 9 3 ## 10 3.16228 ```{awk} BEGIN {print "Hello, world!"} ```
## Hello, world!
【数学のひっかけ問題】三角形の面積を求めよ
https://nazesuugaku.com/fake_triangle/
で，「底辺が10cm、高さ6cmの直角三角形の面積を求めなさい。」
に対して，「そのような直角三角形は存在しない」とある。
確かに記事中のように底辺をとると「そのような直角三角形は存在しない」であるが，明らかに別の直角三角形もあるじゃないか？
短い問題文の中に，この三角形を排除する記述はないと思うのだが。
元々は「マイクロソフトの入社試験」だそうなので，原文を見ないとどうとも言えないが。
どうも，ご丁寧に図が描いてあったらしい。これじゃ「そんな三角形は存在しない」わなぁ。
真面目な解答もありました。（すごいなあ）
これ解けますか？（後編：マイクロソフトの三角形ありました）
https://blog.tech-monex.com/entry/2021/08/06/130000
sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。
"id","label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"
これを入力とし，以下のように変更を加え出力する。
"id",label"
1,"aaa"
1,"あああ"
2,"aaa"
2,"いいい"
3,"bbb"
3,"ううう"
4,"bbb"
4,"えええ"
5,"ccc"
5,"おおお"
6,"ccc"
6,"かかか"
# tidyr
あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 125.581, 122.979, 121.734 であった。
平均値は 123.431 秒。
しかし，システムが使用する時間が同程度あり，結局経過時間（計算開始から結果が出るまでの時間）は実に388.996, 283.905, 297.495 で，平均値は 323.465 秒，実に 5 分超である。
library(dplyr)
library(tidyr)
system.time({
sample = read.csv("sample.csv")
sample_row <- sample %>% separate_rows(label,sep="_")
write.csv(sample_row, "tidyr.csv", row.names=FALSE)
})
# Base R
最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 4 行のプログラムでの実行時間は 25.035, 25.480, 25.446 であった。
平均値は 25.320 秒。
経過時間は 24.190 秒である。tidyr によるものより 30 倍速い。
system.time({
sample = read.csv("sample.csv")
a = unlist(strsplit(sample$label, "_"))
sample = data.frame(id=rep(sample$id, each=2), label=a)
write.csv(sample, "BaseR.csv", row.names=FALSE)
})
# AWK
AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 5 行のプログラムを書く。実行時間は 21.380, 22.024, 21.848 であった。
平均値は 21.751 秒。
BaseR によるよりも 16% 速い。tidyr に比べると 5 倍速い。
BEGIN {
FS = ","
getline
print $0
}
{
split($2, t, "_")
printf "%s,%s\"\n%s,\"%s\n", $1, t[1], $1, t[2]
}
結果
またまた，AWK が最速という結果になった。
R はデータをすべてメモリ中に保存している。しかも，場合によってはそのコピーを作ったりすることもある。このプログラムが動いているとき R は 2GB ものメモリを消費している。
AWK はデータを「一行読んで，処理して，書き出して」を繰り返すのでメモリは殆ど使わない。
そのあたりがこういう結果になっているのであろう。
特に tidyr は，BaseR に比べてメモリ操作が下手くそ（余計なことをしている）なのであろう。
sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。
"id","label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"
これを入力とし，以下のように変更を加え出力する。
"id","new_label1","new_label2"
1,"aaa","あああ"
2,"aaa","いいい"
3,"bbb","ううう"
4,"bbb","えええ"
5,"ccc","おおお"
6,"ccc","かかか"
# tidyr
あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 111.195, 101.132, 100.276 であった。
平均値は 104.201 秒。
しかし，システムが使用する時間が同程度あり，結局経過時間（計算開始から結果が出るまでの時間）は実に751.059, 686.955, 827.347 で，平均値は 755.120 秒，実に 12 分超である。
library(dplyr)
library(tidyr)
system.time({
sample = read.csv("sample.csv")
sample_col <- sample %>% separate(label,into=c("new_label1","new_label2"),sep="_")
write.csv(sample_col, "tidyr.csv", row.names=FALSE)
})
# Base R
最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 6 行のプログラムでの実行時間は 23.509, 23.268, 22.892 であった。平均値は 23.223 秒。
経過時間は 24.190 秒である。tidyr によるものより 30 倍速い。
system.time({
sample = read.csv("sample.csv")
a = unlist(strsplit(sample$label, "_"))
sample$new_label1 = a[seq(1, 12, by=2)]
sample$new_label2 = a[seq(2, 12, by=2)]
sample = sample[c("id","new_label1","new_label2")]
write.csv(sample, "BaseR.csv", row.names=FALSE)
})
# AWK
AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 5 行のプログラムを書く。実行時間は 20.328, 20.503, 20.126 であった。
平均値は 20.319 秒。
BaseR によるよりも二割方速い。tidyr に比べると 5 倍速い。
BEGIN {
FS = ","
getline
print "\"id\",\"new_label1\",\"new_label2\""
}
{
split($2, t, "_")
printf "%s,%s\",\"%s\n", $1, t[1], t[2]
}
結果
AWK が最速という結果になった。
sample.csv は 17035201 行の CSV ファイルである。最初の 7 行は以下のようになっている。
"id","label","label2"
1,"aaa","あああ"
2,"aaa","いいい"
3,"bbb","ううう"
4,"bbb","えええ"
5,"ccc","おおお"
6,"ccc","かかか"
これを入力とし，以下のように変更を加え出力する。
"id","new_label"
1,"aaa_あああ"
2,"aaa_いいい"
3,"bbb_ううう"
4,"bbb_えええ"
5,"ccc_おおお"
6,"ccc_かかか"
# tidyr
あまり 'tidy' なプログラムでもないが，この 3 行のプログラムでの実行時間は 11.633, 11.909, 12.028 であった。
平均値は 11.857 秒。
library(dplyr)
library(tidyr)
system.time({
sample = read.csv("sample.csv")
sample_col <- sample %>% unite(new_label,label,label2,remove=TRUE,sep="_")
write.csv(sample_col, "tidyr.csv", row.names=FALSE)
})
# Base R
最近では tidyr は知っているが，Base R は知らない（書けない）という人も多いようであるが，この 4 行のプログラムでの実行時間は 12.017, 11.987, 12.020 であった。
平均値は 12.008 秒。
system.time({
sample = read.csv("sample.csv")
sample$new_label = paste(sample$label, sample$label2, sep="_")
sample = sample[c("id", "new_label")]
write.csv(sample, "BaseR.csv", row.names=FALSE)
})
# AWK
AWK を知らない人も多いと思うが，この例のような簡単なデータ前処理には十分使える。
この CSV ファイルの構造をフルに利用して，実質 4 行のプログラムを書く。
これくらい短ければ，使い捨てのプログラムをすぐ書ける。
実行時間は 11.026, 11.087, 10.952 であった。
平均値は 11.0217 秒。
BEGIN {
getline
print "\"id\",\"new_label\""
}
{
sub("\",\"", "_", $0)
print $0
}
結果
AWK が最速という結果になった。
階乗の限界と計算速度を Python, Julia, R, Octave で比較してみた
追記: 2022/08/18 R での結果を更新。x86-64 版でやったら早々にクラッシュしていた。M1 チップ対応版で Julia と同等以上の速度である。
# Python 3.10.6
from math import factorial
from time import time
s = time()
a = factorial(100000)
print(time() - s)
0.1885082721710205 sec.
# Julia 1.7.3
@time a = factorial(big(100000));
0.011688 seconds (302 allocations: 3.076 MiB)
# R 4.2.1 (M1)
library(gmp)
> system.time({a = factorialZ(100000)})
ユーザ システム 経過
0.009999999999999787 0.000000000000000000 0.009999999999990905
# Octave 6.2.0
tic
a = factorial(18); # 整数値解を表示するにはこれが限界
toc
Elapsed time is 0.00777102 seconds.
ということで，Julia と R の圧勝...かな?
Julia も R も，gmp を使っているので，Python も gmpy を使えば同じ速度になるのかも知れない。が，現在の環境では pip install gmpy はできるのだが，いざ import gmpy すると，エラーになるので，確認できない。
ヒルベルト行列
ヒルベルト行列とその逆行列を求め，両者の行列乗算で単位行列になることを，4 種の言語でやってみた。
Octave と Julia はさすが。
R はちょっと誤差が大きい?
Python は R より誤差が大きい?
Octave
hilb
Return the Hilbert matrix of order N.
invhilb
Return the inverse of the Hilbert matrix of order N.
hilb(5)
ans =
1.0000 0.5000 0.3333 0.2500 0.2000
0.5000 0.3333 0.2500 0.2000 0.1667
0.3333 0.2500 0.2000 0.1667 0.1429
0.2500 0.2000 0.1667 0.1429 0.1250
0.2000 0.1667 0.1429 0.1250 0.1111
invhilb(5)
ans =
25 -300 1050 -1400 630
-300 4800 -18900 26880 -12600
1050 -18900 79380 -117600 56700
-1400 26880 -117600 179200 -88200
630 -12600 56700 -88200 44100
hilb(5) * invhilb(5)
ans =
1.0000 0 0 0 0
0 1.0000 0 0 0
0 0 1.0000 -0.0000 0
0 0 0 1.0000 0
0 0 0 0 1.0000
Julia
# import Pkg; Pkg.add("Nemo")
using Nemo
M = MatrixSpace(QQ, 5, 5)
h = hilbert(M)
Welcome to Nemo version 0.32.2
Nemo comes with absolutely no warranty whatsoever
[ 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]
inv(h)
[ 25 -300 1050 -1400 630]
[ -300 4800 -18900 26880 -12600]
[ 1050 -18900 79380 -117600 56700]
[-1400 26880 -117600 179200 -88200]
[ 630 -12600 56700 -88200 44100]
h * inv(h)
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
[0 0 0 1 0]
[0 0 0 0 1]
R
# install.packages("matrixcalc")
library(matrixcalc)
h = hilbert.matrix(5)
h
A matrix: 5 × 5 of type dbl
1.0000000 0.5000000 0.3333333 0.2500000 0.2000000
0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
0.3333333 0.2500000 0.2000000 0.1666667 0.1428571
0.2500000 0.2000000 0.1666667 0.1428571 0.1250000
0.2000000 0.1666667 0.1428571 0.1250000 0.1111111
solve(h)
A matrix: 5 × 5 of type dbl
25 -300 1050 -1400 630
-300 4800 -18900 26880 -12600
1050 -18900 79380 -117600 56700
-1400 26880 -117600 179200 -88200
630 -12600 56700 -88200 44100
h %*% solve(h)
A matrix: 5 × 5 of type dbl
1.000000e+00 -4.547474e-13 -1.818989e-12 0.000000e+00 0.000000e+00
-1.421085e-14 1.000000e+00 0.000000e+00 -5.456968e-12 -2.728484e-12
0.000000e+00 -4.547474e-13 1.000000e+00 0.000000e+00 -9.094947e-13
0.000000e+00 2.273737e-13 0.000000e+00 1.000000e+00 -9.094947e-13
1.421085e-14 -2.273737e-13 0.000000e+00 -1.818989e-12 1.000000e+00
round(h %*% solve(h), 10)
A matrix: 5 × 5 of type dbl
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1
Python
from scipy.linalg import hilbert, invhilbert
hilbert(5)
array([[1. , 0.5 , 0.33333333, 0.25 , 0.2 ],
[0.5 , 0.33333333, 0.25 , 0.2 , 0.16666667],
[0.33333333, 0.25 , 0.2 , 0.16666667, 0.14285714],
[0.25 , 0.2 , 0.16666667, 0.14285714, 0.125 ],
[0.2 , 0.16666667, 0.14285714, 0.125 , 0.11111111]])
invhilbert(5)
array([[ 2.500e+01, -3.000e+02, 1.050e+03, -1.400e+03, 6.300e+02],
[-3.000e+02, 4.800e+03, -1.890e+04, 2.688e+04, -1.260e+04],
[ 1.050e+03, -1.890e+04, 7.938e+04, -1.176e+05, 5.670e+04],
[-1.400e+03, 2.688e+04, -1.176e+05, 1.792e+05, -8.820e+04],
[ 6.300e+02, -1.260e+04, 5.670e+04, -8.820e+04, 4.410e+04]])
hilbert(5) @ invhilbert(5)
array([[ 1.00000000e+00, -1.39888101e-13, 6.29496455e-13, 2.65876210e-12, -1.32938105e-12],
[-2.00395256e-14, 1.00000000e+00, -2.34356978e-12, 2.63500333e-12, -1.31750166e-12],
[ 2.34257058e-14, -1.27453603e-13, 1.00000000e+00, -2.93853830e-12, 5.59774449e-13],
[ 0.00000000e+00, -2.27373675e-13, 9.09494702e-13, 1.00000000e+00, 0.00000000e+00],
[-1.80966353e-14, 3.05089287e-13, -3.49720253e-13, 2.36299869e-12, 1.00000000e+00]])
import numpy as np
np.round(hilbert(5) @ invhilbert(5), 10)
array([[ 1., -0., 0., 0., -0.],
[-0., 1., -0., 0., -0.],
[ 0., -0., 1., -0., 0.],
[ 0., -0., 0., 1., 0.],
[-0., 0., -0., 0., 1.]])
以下の例を上げてきたが，これはもはや「早期リターンではない」
class Student:
def __init__(self, student_id, name):
self.student_id = student_id
self.name = name
def student_check(self):
if self.student_id != "" or self.name != "":
if len(self.student_id) > 12:
flg = False
else:
flg = True
else:
flg = False
return flg
student_exe = Student("101122222212","山田花子")
student_exe.student_check()
「if のネストを避ける」というなら，以下のようにすればよい。これで「早期リターン」だ！！
def student_check(self):
if self.student_id != "" or self.name != "":
return len(self.student_id) <= 12
else:
return False
蛇足だが，flg なんて変数名はダメダメ。flag にしよう。
早期リターンとガード節
ガード節を使わないプログラムとして挙げられているのが prog1
if value === "hoge"
result = "hoge"
else
if value === "fuga"
result = "fuga"
else
if value === "hogefuga"
result = "hogefuga"
end
return result
ガード節を使うプログラムとして挙げられているのが prog2
if value === "hoge"
return "hoge"
end
if value === "huga"
return "fuga"
end
return "hogefuga"
この言語は何かわからないが，やっていることはわかる。
疑問は，「この言語は if-else しかないのか?if-elseif-else はないのか?」ということ
prog2 は，以下のように書くとわかりやすいのではないか? prog3
if value === "hoge"
return "hoge"
else if value === "huga"
return "fuga"
else
return "hogefuga"
end
prog1 では，value が "hoge", "huga", "hogefuga" 以外の場合，未定義の result を返そうとする。言語によっては単に何も返さないと解釈されることもあるが，いずれにせよ，論理エラーである。
prog2, prog3 では，value が "hoge", "huga" 以外の場合，無条件に "hogefuga" を返す。これは論理エラーである。
prog3 はもう一段 else if があるべきだ。 prog4
これは prog1 を else if を使って，更に条件のチェック漏れを補完したプログラムだ。
if value === "hoge"
return "hoge"
else if value === "huga"
return "fuga"
else if value === "hogefuga"
return "hogefuga"
else
return "ERROR!!"
end
最後に，このバカバカしいプログラムに限ってだが，以下のように書くのがベストだ。prog5
if value === "hoge" || value === "huga" || value === "hogefuga"
return value
else
return "ERROR!!"
end
あるプログラミングスタイルを推奨する記事を書くとき，悪い例と良い例を提示するが，提示されたプログラムがクソだと，プログラミングスタイル自体がクソに見える。
早期リターンは必要なら使えばよい。
いずれにしても，プログラムをよく解析して，最良なアルゴリズムを適用することだ。
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]
PVアクセスランキング にほんブログ村
