🍝

単精度・倍精度・四倍精度・八倍精度で固有値問題を解く

2022/12/19に公開

数値計算 Advent Calendar 2022 19日目です🎄

このノートでは単精度・倍精度・(ついでに)疑似四倍精度・四倍精度・八倍精度でフランク行列ヒルベルト行列の固有値を求めます. フランク行列の場合, 1000×1000程のサイズで単精度では最小固有値が1桁も正しく計算できなくなることを確かめました. より条件数の大きいヒルベルト行列の場合, さらに小さい5×5~35×35程で単精度・倍精度・疑似四倍精度・四倍精度・八倍精度では最小固有値が1桁も正しく計算できなくなることを確かめました.

パッケージ・環境

Quadmath.jlの四倍精度・Float128はIEEE 754 binary128準拠です. ついでにDoubleFloats.jlの疑似四倍精度・Double64も追加しました. QD.jlの八倍精度・Float256はインストール, ビルド時にエラーが出たので, 今回は代わりにBigFloatを使いました. BigFloatの仮数は自由に変えられるので237ビットにしてIEEE 754 binary256に合わせましたが, 指数部分はデフォルトから変えられないらしいのでIEEE 754 binary256より大きくなっています. なお, BigFloatのバックエンドはMPFRで, Juliaに最初から入っているので特別なパッケージ等は不要です. 後ほど, setprecision(237)のようにして仮数の桁を変えます. 固有値については, Float32とFloat64までは標準ライブラリのLinearAlgebra.jl(固有値部分はLAPACK)で計算できます. Double64, Float128, BigFloatはGenericLinearAlgebra.jlで計算します.

入力
# using Pkg
# Pkg.add("Quadmath")
# Pkg.add("DoubleFloats")
# Pkg.add("GenericLinearAlgebra")

using Printf               # 表示
using DoubleFloats         # Double64型のサポート
using Quadmath             # Float128型のサポート
using LinearAlgebra        # Float32, Float64までの固有値計算
using GenericLinearAlgebra # Double64, Float128, BigFloatの固有値計算

versioninfo()
出力
Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
  Threads: 1 on 8 virtual cores

BigFloatを八倍精度の代わりに使う

まず, BigFloatが単精度・倍精度・疑似四倍精度・四倍精度の代わりに使えるのか試します. それぞれの型の仮数部分のビット数はそれぞれ下記のようになっています.

入力
for t in [Float32, Float64, Double64, Float128, BigFloat]
    println(t, " \t", precision(t))
end
出力
Float32 	24
Float64 	53
Double64 	106
Float128 	113
BigFloat 	256

例えば, setprecision(53)とするとBigFloatとFloat64の仮数部分のビット数が一致しますので, 近い振る舞いをすると考えられます. 有理数・Rational型の1/3からそれぞれの型に変換すると下記のようになります. ただし, 仮数のビット数が53のBigFloatをBigFloat(53)のように表しました.

入力
for t in [Float32, Float64, Double64, Float128]
    p = precision(t)
    setprecision(p) do
        println("$t \t", big(convert(t,1//3)))
        println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
        println()
    end
end

setprecision(237) do
    println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
    println()
end

println("BigFloat($(precision(BigFloat)))\t", BigFloat(1//3))
出力
Float32 	0.333333343
BigFloat(24)	0.333333343

Float64 	0.33333333333333331
BigFloat(53)	0.33333333333333331

Double64 	0.333333333333333333333333333333335
BigFloat(106)	0.333333333333333333333333333333335

Float128 	0.333333333333333333333333333333333317
BigFloat(113)	0.333333333333333333333333333333333317

BigFloat(237)	0.3333333333333333333333333333333333333333333333333333333333333333333333326

BigFloat(256)	0.3333333333333333333333333333333333333333333333333333333333333333333333333333348

BigFloatの仮数部分を他の型に合わせると, 少なくとも有理数・Rationalから変換した際の桁数は一致することが確かめられました. 完全には一致しないかもしれませんが, BigFloat(237)で概ねFloat256の目安にはなると思われます. 実際, 後ほどの固有値計算でもFloat128とBigFloat(113)は同じ行列サイズまで計算できるという結果を与えています.

フランク行列

フランク行列の実装例を示します. まず整数型で要素を作ってから, datatypeキーワードで指定した型へ変換するようにしてあります. なお, 速度よりもFortran等への移植性を優先しました.

入力
# 速い例:
# Frank(n; datatype=Int64) = convert.(datatype, [n - max(i,j) + 1 for i=1:n, j=1:n])

# 遅い例:
function Frank(n; datatype=Int64)
    A = zeros(datatype, n, n)
    for j in 1:n
        for i in 1:n
            A[i,j] = convert(datatype, n - max(i,j) + 1)
        end
    end
    return A
end

A = Frank(5)
出力
5×5 Matrix{Int64}:
 5  4  3  2  1
 4  4  3  2  1
 3  3  3  2  1
 2  2  2  2  1
 1  1  1  1  1

幸谷智紀『LAPACK/BLAS入門』(森北出版, 2016) p.76 表4.7を参考にn=10,100,1000における最小固有値を計算します. 表4.7と同様にn=1000では単精度で計算した最小固有値が1桁も合わないことが確かめられます.

入力
# using LinearAlgebra
# using GenericLinearAlgebra

for n in [10,100,1000]
    println("\nn = $n  \tλ[1]")
    for datatype in [Float32,Float64]#,Double64,Float128,BigFloat]
        A = Frank(n, datatype=datatype)
        λ = eigvals(A)
        println(datatype, " \t", λ[1])
    end
end
出力
n = 10  	λ[1]
Float32 	0.25567994
Float64 	0.2556795627964355

n = 100  	λ[1]
Float32 	0.2500508
Float64 	0.2500610827206968

n = 1000  	λ[1]
Float32 	0.12658979
Float64 	0.2500006162349238

Float32では1000×1000のフランク行列の最小固有値が1桁も正しく計算できないことが確かめられました. 実際には900×900程度から振動し始め, おかしな挙動を見せます. フランク行列だとメモリや計算時間の限界が先に来ますので, より悪条件なヒルベルト行列で単精度・倍精度・疑似四倍精度・四倍精度・八倍精度の限界を見てみます.

ヒルベルト行列

ヒルベルト行列の実装例を示します. まず有理数型で要素を作り, それをdatatypeキーワードで指定した型に変換します. 1 / (i + j - 1)だと勝手にFloat64と解釈されて丸め誤差が発生してしまいますので, 1 // (i + j - 1)とすることでRationalから他の型に直接変換し, 丸め誤差を小さくできます. なお, 速度よりもFortran等への移植性を優先しました.

入力
# 速い例:
# Hilbert(n; datatype=Rational) = convert.(datatype, [1 // (i + j - 1) for i=1:n, j=1:n])

# 遅い例:
function Hilbert(n; datatype=Rational)
    A = zeros(datatype, n, n)
    for j in 1:n
        for i in 1:n
            A[i,j] = convert(datatype, 1 // (i + j - 1))
        end
    end
    return A
end

H = Hilbert(5)
出力
5×5 Matrix{Rational}:
 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

ヒルベルト行列のサイズを変えてFloat32, Float64, Double64, Float128, BigFloat(113), BigFloat(237), BigFloat(256)で固有値を計算します.

入力
# using Printf
# using DoubleFloats
# using Quadmath
# using LinearAlgebra
# using GenericLinearAlgebra

for n in 2:40
    println("\nn = $n   \tλ[1]")
    for datatype in [Float32, Float64, Double64, Float128]
        A = Hilbert(n, datatype=datatype)
        λ = eigvals(A)
        @printf("%s \t%.9e\n", datatype, λ[1])
    end

    for p in [113,237,256,1024]
        setprecision(p) do
            A = Hilbert(n, datatype=BigFloat)
            λ = eigvals(A)
            @printf("%s \t%.9e\n", "BigFloat($(precision(BigFloat)))", λ[1])
        end
    end

end

# setprecision(BigFloat, 256)
出力
n = 2   	λ[1]
Float32 	6.574144959e-02
Float64 	6.574145409e-02
Double64 	6.574145409e-02
Float128 	6.574145409e-02
BigFloat(113) 	6.574145409e-02
BigFloat(237) 	6.574145409e-02
BigFloat(256) 	6.574145409e-02
BigFloat(1024) 	6.574145409e-02

n = 3   	λ[1]
Float32 	2.687271917e-03
Float64 	2.687340356e-03
Double64 	2.687340356e-03
Float128 	2.687340356e-03
BigFloat(113) 	2.687340356e-03
BigFloat(237) 	2.687340356e-03
BigFloat(256) 	2.687340356e-03
BigFloat(1024) 	2.687340356e-03

n = 4   	λ[1]
Float32 	9.672460874e-05
Float64 	9.670230402e-05
Double64 	9.670230402e-05
Float128 	9.670230402e-05
BigFloat(113) 	9.670230402e-05
BigFloat(237) 	9.670230402e-05
BigFloat(256) 	9.670230402e-05
BigFloat(1024) 	9.670230402e-05

n = 5   	λ[1]
Float32 	3.361120889e-06
Float64 	3.287928772e-06
Double64 	3.287928772e-06
Float128 	3.287928772e-06
BigFloat(113) 	3.287928772e-06
BigFloat(237) 	3.287928772e-06
BigFloat(256) 	3.287928772e-06
BigFloat(1024) 	3.287928772e-06

n = 6   	λ[1]
Float32 	1.306005259e-07
Float64 	1.082799486e-07
Double64 	1.082799485e-07
Float128 	1.082799485e-07
BigFloat(113) 	1.082799485e-07
BigFloat(237) 	1.082799485e-07
BigFloat(256) 	1.082799485e-07
BigFloat(1024) 	1.082799485e-07

n = 7   	λ[1]
Float32 	7.071950847e-08
Float64 	3.493898716e-09
Double64 	3.493898606e-09
Float128 	3.493898606e-09
BigFloat(113) 	3.493898606e-09
BigFloat(237) 	3.493898606e-09
BigFloat(256) 	3.493898606e-09
BigFloat(1024) 	3.493898606e-09

n = 8   	λ[1]
Float32 	-4.826334887e-08
Float64 	1.111539635e-10
Double64 	1.111538966e-10
Float128 	1.111538966e-10
BigFloat(113) 	1.111538966e-10
BigFloat(237) 	1.111538966e-10
BigFloat(256) 	1.111538966e-10
BigFloat(1024) 	1.111538966e-10

n = 9   	λ[1]
Float32 	2.470161853e-09
Float64 	3.499736372e-12
Double64 	3.499676403e-12
Float128 	3.499676403e-12
BigFloat(113) 	3.499676403e-12
BigFloat(237) 	3.499676403e-12
BigFloat(256) 	3.499676403e-12
BigFloat(1024) 	3.499676403e-12

n = 10   	λ[1]
Float32 	-6.898610394e-08
Float64 	1.093880753e-13
Double64 	1.093153819e-13
Float128 	1.093153819e-13
BigFloat(113) 	1.093153819e-13
BigFloat(237) 	1.093153819e-13
BigFloat(256) 	1.093153819e-13
BigFloat(1024) 	1.093153819e-13

n = 11   	λ[1]
Float32 	-8.925326256e-08
Float64 	3.375236254e-15
Double64 	3.393218595e-15
Float128 	3.393218595e-15
BigFloat(113) 	3.393218595e-15
BigFloat(237) 	3.393218595e-15
BigFloat(256) 	3.393218595e-15
BigFloat(1024) 	3.393218595e-15

n = 12   	λ[1]
Float32 	-3.159367168e-08
Float64 	1.630874793e-16
Double64 	1.047946398e-16
Float128 	1.047946398e-16
BigFloat(113) 	1.047946398e-16
BigFloat(237) 	1.047946398e-16
BigFloat(256) 	1.047946398e-16
BigFloat(1024) 	1.047946398e-16

n = 13   	λ[1]
Float32 	-6.428467003e-08
Float64 	-3.723992712e-17
Double64 	3.222901015e-18
Float128 	3.222901015e-18
BigFloat(113) 	3.222901015e-18
BigFloat(237) 	3.222901015e-18
BigFloat(256) 	3.222901015e-18
BigFloat(1024) 	3.222901015e-18

n = 14   	λ[1]
Float32 	-2.169667468e-08
Float64 	6.961897267e-18
Double64 	9.877051735e-20
Float128 	9.877051735e-20
BigFloat(113) 	9.877051735e-20
BigFloat(237) 	9.877051735e-20
BigFloat(256) 	9.877051735e-20
BigFloat(1024) 	9.877051735e-20

n = 15   	λ[1]
Float32 	-5.260721281e-08
Float64 	-1.970813429e-17
Double64 	3.017915296e-21
Float128 	3.017915296e-21
BigFloat(113) 	3.017915296e-21
BigFloat(237) 	3.017915296e-21
BigFloat(256) 	3.017915296e-21
BigFloat(1024) 	3.017915296e-21

n = 16   	λ[1]
Float32 	-2.375249863e-08
Float64 	-7.991212304e-17
Double64 	1.971645939e-22
Float128 	9.197419821e-23
BigFloat(113) 	9.197419821e-23
BigFloat(237) 	9.197419821e-23
BigFloat(256) 	9.197419821e-23
BigFloat(1024) 	9.197419821e-23

n = 17   	λ[1]
Float32 	-7.643087940e-08
Float64 	-3.084917681e-17
Double64 	2.796721645e-24
Float128 	2.796721632e-24
BigFloat(113) 	2.796721632e-24
BigFloat(237) 	2.796721632e-24
BigFloat(256) 	2.796721632e-24
BigFloat(1024) 	2.796721632e-24

n = 18   	λ[1]
Float32 	-3.352907640e-08
Float64 	-2.922712756e-17
Double64 	1.804347127e-25
Float128 	8.487413205e-26
BigFloat(113) 	8.487413251e-26
BigFloat(237) 	8.487413230e-26
BigFloat(256) 	8.487413230e-26
BigFloat(1024) 	8.487413230e-26

n = 19   	λ[1]
Float32 	-2.656284792e-08
Float64 	-1.790409527e-17
Double64 	3.747605947e-27
Float128 	3.747605889e-27
BigFloat(113) 	3.747605889e-27
BigFloat(237) 	2.571239705e-27
BigFloat(256) 	2.571239705e-27
BigFloat(1024) 	2.571239705e-27

n = 20   	λ[1]
Float32 	-2.125183407e-08
Float64 	-2.797264572e-16
Double64 	-1.905490628e-25
Float128 	1.233279673e-28
BigFloat(113) 	1.233279672e-28
BigFloat(237) 	7.777377397e-29
BigFloat(256) 	7.777377397e-29
BigFloat(1024) 	7.777377397e-29

n = 21   	λ[1]
Float32 	-5.624558241e-08
Float64 	-1.348650289e-17
Double64 	-2.714301231e-23
Float128 	3.678871784e-30
BigFloat(113) 	3.678871738e-30
BigFloat(237) 	2.349182538e-30
BigFloat(256) 	2.349182538e-30
BigFloat(1024) 	2.349182538e-30

n = 22   	λ[1]
Float32 	-9.640318410e-08
Float64 	-1.951507152e-16
Double64 	-2.703080487e-30
Float128 	-2.703634429e-30
BigFloat(113) 	-2.703493685e-30
BigFloat(237) 	7.086841582e-32
BigFloat(256) 	7.086841582e-32
BigFloat(1024) 	7.086841582e-32

n = 23   	λ[1]
Float32 	-1.269348502e-07
Float64 	-1.840525074e-16
Double64 	-1.395136546e-26
Float128 	-1.395136676e-26
BigFloat(113) 	-1.395136676e-26
BigFloat(237) 	2.135463438e-33
BigFloat(256) 	2.135463438e-33
BigFloat(1024) 	2.135463438e-33

n = 24   	λ[1]
Float32 	-5.662219849e-08
Float64 	-1.536692857e-16
Double64 	-6.294692986e-25
Float128 	-6.294693199e-25
BigFloat(113) 	-6.294693199e-25
BigFloat(237) 	6.428050325e-35
BigFloat(256) 	6.428050325e-35
BigFloat(1024) 	6.428050325e-35

n = 25   	λ[1]
Float32 	-9.834495529e-08
Float64 	-1.158400702e-16
Double64 	-9.786027902e-24
Float128 	3.475441588e-36
BigFloat(113) 	7.644790285e-36
BigFloat(237) 	1.933092811e-36
BigFloat(256) 	1.933092811e-36
BigFloat(1024) 	1.933092811e-36

n = 26   	λ[1]
Float32 	-8.256382245e-08
Float64 	-1.745151220e-16
Double64 	-9.530727276e-23
Float128 	-3.164126445e-29
BigFloat(113) 	-3.164126445e-29
BigFloat(237) 	5.808263531e-38
BigFloat(256) 	5.808263531e-38
BigFloat(1024) 	5.808263531e-38

n = 27   	λ[1]
Float32 	-1.064803996e-07
Float64 	-3.982034472e-17
Double64 	-2.741452206e-34
Float128 	-1.690538599e-27
BigFloat(113) 	-1.690538566e-27
BigFloat(237) 	1.743773036e-39
BigFloat(256) 	1.743773036e-39
BigFloat(1024) 	1.743773036e-39

n = 28   	λ[1]
Float32 	-8.385376304e-08
Float64 	-3.548348719e-16
Double64 	-1.942851263e-34
Float128 	-2.793018499e-26
BigFloat(113) 	-2.793018543e-26
BigFloat(237) 	5.231305960e-41
BigFloat(256) 	5.231305960e-41
BigFloat(1024) 	5.231305960e-41

n = 29   	λ[1]
Float32 	-1.160353804e-07
Float64 	-6.199549616e-17
Double64 	2.031105710e-35
Float128 	-2.822973889e-25
BigFloat(113) 	-2.822973889e-25
BigFloat(237) 	1.568304495e-42
BigFloat(256) 	1.568304495e-42
BigFloat(1024) 	1.568304495e-42

n = 30   	λ[1]
Float32 	-5.065797382e-08
Float64 	-5.289389957e-17
Double64 	-1.102984167e-34
Float128 	-2.146399554e-24
BigFloat(113) 	-2.146399554e-24
BigFloat(237) 	4.698636565e-44
BigFloat(256) 	4.698636565e-44
BigFloat(1024) 	4.698636565e-44

n = 31   	λ[1]
Float32 	-7.601806118e-08
Float64 	-7.329944390e-17
Double64 	-2.916407324e-34
Float128 	-1.255399562e-26
BigFloat(113) 	-1.255399570e-26
BigFloat(237) 	1.406868380e-45
BigFloat(256) 	1.406868380e-45
BigFloat(1024) 	1.406868380e-45

n = 32   	λ[1]
Float32 	-6.675949749e-08
Float64 	-2.374236517e-17
Double64 	-1.025210430e-34
Float128 	-3.257260908e-36
BigFloat(113) 	-3.257273966e-36
BigFloat(237) 	4.210099032e-47
BigFloat(256) 	4.210099032e-47
BigFloat(1024) 	4.210099032e-47

n = 33   	λ[1]
Float32 	-3.973698597e-07
Float64 	-6.226694532e-16
Double64 	7.189674368e-35
Float128 	-2.193596836e-36
BigFloat(113) 	-2.196233777e-36
BigFloat(237) 	1.259226436e-48
BigFloat(256) 	1.259226436e-48
BigFloat(1024) 	1.259226436e-48

n = 34   	λ[1]
Float32 	-2.554655794e-07
Float64 	-3.335481520e-16
Double64 	-3.129737789e-34
Float128 	-1.976669371e-25
BigFloat(113) 	-1.976656135e-25
BigFloat(237) 	3.764454505e-50
BigFloat(256) 	3.764454505e-50
BigFloat(1024) 	3.764454505e-50

n = 35   	λ[1]
Float32 	-2.991216945e-07
Float64 	-5.137067885e-16
Double64 	-2.859278165e-34
Float128 	-2.330807894e-36
BigFloat(113) 	-4.558199889e-36
BigFloat(237) 	2.318729279e-51
BigFloat(256) 	1.124863268e-51
BigFloat(1024) 	1.124863268e-51

n = 36   	λ[1]
Float32 	-1.958887168e-07
Float64 	-4.311007920e-16
Double64 	-2.255124275e-34
Float128 	-3.238481422e-36
BigFloat(113) 	-3.284703834e-36
BigFloat(237) 	6.919587088e-53
BigFloat(256) 	3.359762135e-53
BigFloat(1024) 	3.359762135e-53

n = 37   	λ[1]
Float32 	-1.831468239e-07
Float64 	-6.535979384e-16
Double64 	-2.116415176e-34
Float128 	-4.635128397e-36
BigFloat(113) 	-4.608945882e-36
BigFloat(237) 	2.064207412e-54
BigFloat(256) 	1.003088162e-54
BigFloat(1024) 	1.003088162e-54

n = 38   	λ[1]
Float32 	-3.336613190e-07
Float64 	-4.512357590e-16
Double64 	-4.432253018e-34
Float128 	-7.250642462e-36
BigFloat(113) 	-4.329796608e-36
BigFloat(237) 	4.617382896e-56
BigFloat(256) 	6.155701539e-56
BigFloat(1024) 	2.993651511e-56

n = 39   	λ[1]
Float32 	-2.744432948e-07
Float64 	-5.247024262e-16
Double64 	-1.290923053e-34
Float128 	-3.133177866e-36
BigFloat(113) 	-3.137829276e-36
BigFloat(237) 	1.376487504e-57
BigFloat(256) 	1.835102354e-57
BigFloat(1024) 	8.931079340e-58

n = 40   	λ[1]
Float32 	-2.591677912e-07
Float64 	-4.441535386e-16
Double64 	-4.783821138e-34
Float128 	-1.174660257e-36
BigFloat(113) 	-1.174234562e-36
BigFloat(237) 	3.420000952e-59
BigFloat(256) 	4.102179409e-59
BigFloat(1024) 	2.663517273e-59

Float32はn=6\rightarrow7で1桁も合わなくなり, Float64はn=12\rightarrow13で1桁も合わなくなりました. Float128とBigFloat(113)は同じサイズで固有値が計算できなくなっているので, BigFloat(237)はFloat256の目安になると思われます.

まとめ

単精度・倍精度・四倍精度(ついでに疑似四倍精度)・八倍精度でフランク行列ヒルベルト行列の固有値を求めました. フランク行列の場合, 1000×1000程のサイズで単精度では固有値が1桁も正しく計算できなくなることを確かめました. より条件数の大きいヒルベルト行列の場合, 下表のような結果になりました.

固有値が1桁も正しく計算できなくなるサイズ
Float32
n=6\rightarrow7
Float64
n=12\rightarrow13
Double64
n=17\rightarrow18
Float128
n=18\rightarrow19
BigFloat(113) ≒ Float128
n=18\rightarrow19
BigFloat(237) ≒ Float256
n=34\rightarrow35
BigFloat(256)
n=37\rightarrow38

数値計算では速度至上主義のような価値観に陥りがちですが, 正しい結果が得られないことには速度を追求しても仕方がありません. 四倍精度, 八倍精度は倍精度よりも格段に遅いとされますが, この記事の執筆を通して, 高精度演算のライブラリ開発のモチベーションが理解できました. Juliaの高精度計算のポテンシャルも感じられたと思います.

付録:メモリ消費量

各要素のメモリ消費量をそのまま要素数にかけているだけですので, 実際にはもう少し大きい値になります. 正確な値はvarinfo()で確認できます.

入力
# using Printf

function KMG(byte)
    if byte > 10^9
        return @sprintf("%5.1f GB", byte/10^9)
    elseif byte > 10^6
        return @sprintf("%5.0f MB", byte/10^6)
    elseif byte > 10^3
        return @sprintf("%5.0f KB", byte/10^3)
    else
        return @sprintf("%5.0f  B", byte)
    end
end


println("size\t\t", [@sprintf("%9s", "Float$j \t") for j in [32,64,128,256]]...)
for i in union(1,10,100,1000:1000:20000, 30000:10000:100000)
    println("$(i)²\t\t", [KMG(i^2*j/8)*"\t" for j in [32,64,128,256]]...)
end
出力
size		 Float32 	 Float64 	Float128 	Float256 	
1²		    4  B	    8  B	   16  B	   32  B	
10²		  400  B	  800  B	    2 KB	    3 KB	
100²		   40 KB	   80 KB	  160 KB	  320 KB	
1000²		    4 MB	    8 MB	   16 MB	   32 MB	
2000²		   16 MB	   32 MB	   64 MB	  128 MB	
3000²		   36 MB	   72 MB	  144 MB	  288 MB	
4000²		   64 MB	  128 MB	  256 MB	  512 MB	
5000²		  100 MB	  200 MB	  400 MB	  800 MB	
6000²		  144 MB	  288 MB	  576 MB	  1.2 GB	
7000²		  196 MB	  392 MB	  784 MB	  1.6 GB	
8000²		  256 MB	  512 MB	  1.0 GB	  2.0 GB	
9000²		  324 MB	  648 MB	  1.3 GB	  2.6 GB	
10000²		  400 MB	  800 MB	  1.6 GB	  3.2 GB	
11000²		  484 MB	  968 MB	  1.9 GB	  3.9 GB	
12000²		  576 MB	  1.2 GB	  2.3 GB	  4.6 GB	
13000²		  676 MB	  1.4 GB	  2.7 GB	  5.4 GB	
14000²		  784 MB	  1.6 GB	  3.1 GB	  6.3 GB	
15000²		  900 MB	  1.8 GB	  3.6 GB	  7.2 GB	
16000²		  1.0 GB	  2.0 GB	  4.1 GB	  8.2 GB	
17000²		  1.2 GB	  2.3 GB	  4.6 GB	  9.2 GB	
18000²		  1.3 GB	  2.6 GB	  5.2 GB	 10.4 GB	
19000²		  1.4 GB	  2.9 GB	  5.8 GB	 11.6 GB	
20000²		  1.6 GB	  3.2 GB	  6.4 GB	 12.8 GB	
30000²		  3.6 GB	  7.2 GB	 14.4 GB	 28.8 GB	
40000²		  6.4 GB	 12.8 GB	 25.6 GB	 51.2 GB	
50000²		 10.0 GB	 20.0 GB	 40.0 GB	 80.0 GB	
60000²		 14.4 GB	 28.8 GB	 57.6 GB	115.2 GB	
70000²		 19.6 GB	 39.2 GB	 78.4 GB	156.8 GB	
80000²		 25.6 GB	 51.2 GB	102.4 GB	204.8 GB	
90000²		 32.4 GB	 64.8 GB	129.6 GB	259.2 GB	
100000²		 40.0 GB	 80.0 GB	160.0 GB	320.0 GB

参考文献

https://gist.github.com/ohno/239c6452a2633f24a403f2bad21081dc

追伸

Fortran/LAPACKでの実装例を作りました.
https://github.com/ohno/ill-conditioned

Discussion