-
Notifications
You must be signed in to change notification settings - Fork 11
/
hilbert.jl
127 lines (105 loc) · 3.59 KB
/
hilbert.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
#The m-by-n Hilbert matrix has matrix elements
# H_{ij} = 1/(i+j-1)
#
export Hilbert, InverseHilbert
"""
Hilbert(m [,n])
Construct `m × m` or `m × n`
[`Hilbert` matrix](http://en.wikipedia.org/wiki/Hilbert_matrix)
from its specified dimensions,
where element `i,j` equal to `1 / (i+j-1)`.
```jldoctest hilbert1
julia> A = Hilbert(5)
5×5 Hilbert{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> Matrix(A)
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
```
Inverses are also integer matrices:
```jldoctest hilbert1
julia> inv(A)
5×5 InverseHilbert{Rational{Int64}}:
25//1 -300//1 1050//1 -1400//1 630//1
-300//1 4800//1 -18900//1 26880//1 -12600//1
1050//1 -18900//1 79380//1 -117600//1 56700//1
-1400//1 26880//1 -117600//1 179200//1 -88200//1
630//1 -12600//1 56700//1 -88200//1 44100//1
```
"""
struct Hilbert{T} <: AbstractMatrix{T}
m :: Int
n :: Int
end
Hilbert(::Type{T}, m::Integer, n::Integer) where {T <: Number} = Hilbert{T}(m, n)
Hilbert(::Type{T}, m::Integer, n::Integer) where {T <: Integer} = Hilbert(Rational{T}, m, n)
Hilbert(m::Integer, n::Integer) = Hilbert(Int, m, n)
Hilbert(n::Integer) = Hilbert(n, n)
Hilbert(::Type{T}, n::Integer) where {T <: Number} = Hilbert(T, n, n)
# Define its size
size(H::Hilbert, dim::Integer) = dim==1 ? H.m : dim==2 ? H.n : 1
size(H::Hilbert) = size(H,1), size(H,2)
# Index into a Hilbert matrix
function getindex(H::Hilbert{T}, i::Integer, j::Integer) where {T}
return one(T)/(i+j-1)
end
# Dense version (provided in Core)
#if VERSION < v"1.6"
# Matrix(H::Hilbert) = [H[i,j] for i in 1:size(H,1), j in 1:size(H,2)]
#end
# Some properties
ishermitian(H::Hilbert) = H.m==H.n
isposdef(H::Hilbert) = H.m==H.n
det(H::Hilbert) = inv(det(inv(H)))
# Inverse of a Hilbert matrix
struct InverseHilbert{T} <: AbstractMatrix{T}
n :: Int
end
InverseHilbert(::Type{T}, n::Integer) where {T <: Number} = InverseHilbert{T}(n)
InverseHilbert(n::Integer) = InverseHilbert(Int, n)
# Define its size
size(A::InverseHilbert, dim::Integer) = dim==1 || dim==2 ? A.n : 1
size(A::InverseHilbert) = size(A,1), size(A,2)
# Properties
ishermitian(::InverseHilbert) = true
isposdef(::InverseHilbert) = true
# Index into a inverse Hilbert matrix
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T}
out = (-1)^(i+j) * (i+j-1) * binomial(A.n+i-1, A.n-j) *
binomial(A.n+j-1, A.n-i) * binomial(i+j-2, i-1)^2
return T(out)
end
# Use bigger numbers in case of bigger types
function getindex(A::InverseHilbert{T}, i::Integer, j::Integer) where {T<:Union{BigInt,BigFloat}}
N = big(A.n)
out = (-1)^(i+j) * (i+j-1) * binomial(N+i-1, N-j) *
binomial(N+j-1, N-i) * binomial(big(i+j-2), i-1)^2
return T(out)
end
"""
det(A::InverseHilbert)
Explicit formula for the determinant.
Caution: this function overflows easily.
"""
det(A::InverseHilbert{T}) where {T} = prod(T, (2k+1) * binomial(2k,k)^2 for k in 1:A.n-1)
# Dense version (in Core)
#if VERSION < v"1.6"
# Matrix(A::InverseHilbert) = [A[i,j] for i in 1:size(A,1), j in 1:size(A,2)]
#end
# Define Base.inv
function inv(H::Hilbert{T}) where {T}
H.m == H.n || throw(ArgumentError("Works only for square Hilbert matrices."))
return InverseHilbert(T, H.n)
end
function inv(A::InverseHilbert{T}) where {T}
HT = promote_type(T,typeof(Rational(one(T))))
return Hilbert(HT, A.n)
end