CompatHelper: bump compat for "ArrayLayouts" to "0.8" (#100)

* CompatHelper: bump compat for "ArrayLayouts" to "0.8"

* Test on Julia v1.8

* update getindex for orthogonal

* tests pass

* Update Project.toml

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: Sheehan Olver <[email protected]>

Author: github-actions[bot]

.github/workflows/ci.yml

@@ -11,7 +11,8 @@ jobs:
       matrix:
         version:
           - '1.6'
-          - '^1.7.0-0'
+          - '1'
+          - '^1.8.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,6 +1,6 @@
 name = "InfiniteLinearAlgebra"
 uuid = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
-version = "0.6.5"
+version = "0.6.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -17,16 +17,16 @@ MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 
 [compat]
-ArrayLayouts = "0.7.7"
-BandedMatrices = "0.16.9"
+ArrayLayouts = "0.8"
+BandedMatrices = "0.17"
 BlockArrays = "0.16"
 BlockBandedMatrices = "0.11"
 DSP = "0.7"
 FillArrays = "0.12, 0.13"
 InfiniteArrays = "0.12"
 LazyArrays = "0.22"
 LazyBandedMatrices = "0.7.4"
-MatrixFactorizations = "0.8"
+MatrixFactorizations = "0.9"
 SemiseparableMatrices = "0.3"
 julia = "1.6"
 

src/InfiniteLinearAlgebra.jl

@@ -10,7 +10,7 @@ import Base.Broadcast: BroadcastStyle, Broadcasted, broadcasted
 
 import ArrayLayouts: colsupport, rowsupport, triangularlayout, MatLdivVec, triangulardata, TriangularLayout, TridiagonalLayout, 
                         sublayout, _qr, __qr, MatLmulVec, MatLmulMat, AbstractQLayout, materialize!, diagonaldata, subdiagonaldata, supdiagonaldata,
-                        _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul
+                        _bidiag_forwardsub!, mulreduce, RangeCumsum, _factorize, transposelayout, ldiv!, lmul!, mul, CNoPivot
 import BandedMatrices: BandedMatrix, _BandedMatrix, AbstractBandedMatrix, bandeddata, bandwidths, BandedColumns, bandedcolumns,
                         _default_banded_broadcast, banded_similar
 import FillArrays: AbstractFill, getindex_value, axes_print_matrix_row

src/banded/hessenbergq.jl

@@ -141,13 +141,13 @@ end
 # getindex
 ####
 
-function getindex(Q::UpperHessenbergQ, i::Integer, j::Integer)
+function getindex(Q::UpperHessenbergQ, i::Int, j::Int)
     y = zeros(eltype(Q), size(Q, 2))
     y[j] = 1
     lmul!(Q, y)[i]
 end
 
-getindex(Q::LowerHessenbergQ, i::Integer, j::Integer) = (Q')[j,i]'
+getindex(Q::LowerHessenbergQ, i::Int, j::Int) = (Q')[j,i]'
 
 
 ###

src/banded/infqltoeplitz.jl

@@ -4,98 +4,99 @@
 #  0 d e]
 #
 # is the fixed point
-function tail_de(a::AbstractVector{T}; branch=findmax) where T<:Real
+function tail_de(a::AbstractVector{T}; branch=findmax) where {T<:Real}
     m = length(a)
-    C = [view(a,m-1:-1:1) Vcat(-a[end]*Eye(m-2), Zeros{T}(1,m-2))]
+    C = [view(a, m-1:-1:1) Vcat(-a[end] * Eye(m - 2), Zeros{T}(1, m - 2))]
     λ, V = eigen(C)
     n2, j = branch(abs2.(λ))
     isreal(λ[j]) || throw(DomainError(a, "Real-valued QL factorization does not exist. Try ql(complex(A)) to see if a complex-valued QL factorization exists."))
     n2 ≥ a[end]^2 || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
-    c = sqrt((n2 - a[end]^2)/real(V[1,j])^2)
-    c*real(V[end:-1:1,j])
+    c = sqrt((n2 - a[end]^2) / real(V[1, j])^2)
+    c * real(V[end:-1:1, j])
 end
 
-function tail_de(a::AbstractVector{T}; branch=findmax) where T
+function tail_de(a::AbstractVector{T}; branch=findmax) where {T}
     m = length(a)
-    C = [view(a,m-1:-1:1) Vcat(-a[end]*Eye(m-2), Zeros{T}(1,m-2))]
+    C = [view(a, m-1:-1:1) Vcat(-a[end] * Eye(m - 2), Zeros{T}(1, m - 2))]
     λ, V = eigen(C)::Eigen{float(T),float(T),Matrix{float(T)},Vector{float(T)}}
     n2, j = branch(abs2.(λ))
-    n2 ≥ abs2(a[end]) || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
-    c_abs = sqrt((n2 - abs2(a[end]))/abs2(V[1,j]))
-    c_sgn = -sign(λ[j])/sign(V[1,j]*a[end-1] - V[2,j]*a[end])
-    c_sgn*c_abs*V[end:-1:1,j]    
+    n2 ≥ abs2(a[end]) || throw(DomainError(a, "QL factorization does not exist. This could indicate that the operator is not Fredholm or that the dimension of the kernel exceeds that of the co-kernel. Try again with the adjoint."))
+    c_abs = sqrt((n2 - abs2(a[end])) / abs2(V[1, j]))
+    c_sgn = -sign(λ[j]) / sign(V[1, j] * a[end-1] - V[2, j] * a[end])
+    c_sgn * c_abs * V[end:-1:1, j]
 end
 
 
 # this calculates the QL decomposition of X and corrects sign
 function ql_X!(X)
-    s = sign(real(X[2,end])) 
+    s = sign(real(X[2, end]))
     F = ql!(X)
-    if s ≠ sign(real(X[1,end-1])) # we need to normalise the sign if ql! flips it
+    if s ≠ sign(real(X[1, end-1])) # we need to normalise the sign if ql! flips it
         F.τ[1] = 2 - F.τ[1] # 1-F.τ[1] is the sign so this flips it
-        X[1,1:end-1] *= -1
+        X[1, 1:end-1] *= -1
     end
     F
 end
 
 
 
 
-function ql(Op::TriToeplitz{T}; kwds...) where T<:Real
-    Z,A,B = Op.dl.value, Op.d.value, Op.du.value
-    d,e = tail_de([Z,A,B]; kwds...) # fixed point of QL but with two Qs, one that changes sign
+function ql(Op::TriToeplitz{T}; kwds...) where {T<:Real}
+    Z, A, B = Op.dl.value, Op.d.value, Op.du.value
+    d, e = tail_de([Z, A, B]; kwds...) # fixed point of QL but with two Qs, one that changes sign
     X = [Z A B; zero(T) d e]
     F = ql_X!(X)
-    t,ω = F.τ[2],X[1,end]
-    QL(_BandedMatrix(Hcat([zero(T), e, X[2,2], X[2,1]], [ω, X[2,3], X[2,2], X[2,1]] * Ones{T}(1,∞)), ℵ₀, 2, 1), Vcat(F.τ[1],Fill(t,∞)))
+    t, ω = F.τ[2], X[1, end]
+    QL(_BandedMatrix(Hcat([zero(T), e, X[2, 2], X[2, 1]], [ω, X[2, 3], X[2, 2], X[2, 1]] * Ones{T}(1, ∞)), ℵ₀, 2, 1), Vcat(F.τ[1], Fill(t, ∞)))
 end
 
-ql(Op::TriToeplitz{T}) where T = ql(InfToeplitz(Op))
+ql(Op::TriToeplitz{T}) where {T} = ql(InfToeplitz(Op))
 
 # ql for Lower hessenberg InfToeplitz
-function ql_hessenberg(A::InfToeplitz{T}; kwds...) where T
-    l,u = bandwidths(A)
+function ql_hessenberg(A::InfToeplitz{T}; kwds...) where {T}
+    l, u = bandwidths(A)
     @assert u == 1
     a = reverse(A.data.args[1])
     de = tail_de(a; kwds...)
     X = [transpose(a); zero(T) transpose(de)]::Matrix{float(T)}
     F = ql_X!(X) # calculate data for fixed point
-    factors = _BandedMatrix(Hcat([zero(T); X[1,end-1]; X[2,end-1:-1:1]], [0; X[2,end:-1:1]] * Ones{float(T)}(1,∞)), ℵ₀, l+u, 1)
-    QLHessenberg(factors, Fill(F.Q,∞))
+    factors = _BandedMatrix(Hcat([zero(T); X[1, end-1]; X[2, end-1:-1:1]], [0; X[2, end:-1:1]] * Ones{float(T)}(1, ∞)), ℵ₀, l + u, 1)
+    QLHessenberg(factors, Fill(F.Q, ∞))
 end
 
 
 # remove one band of A
 function ql_pruneband(A; kwds...)
-    l,u = bandwidths(A)
-    A_hess = A[:,u:end]
-    Q,L = ql_hessenberg(A_hess; kwds...)
-    p = size(_pertdata(bandeddata(parent(L))),2)+u+1 # pert size
-    dat = (UpperHessenbergQ((Q').q[1:(p+l)])) * A[1:p+l+1,1:p]
-    pert = Array{eltype(dat)}(undef, l+u+1,size(dat,2)-1)
+    l, u = bandwidths(A)
+    A_hess = A[:, u:end]
+    Q, L = ql_hessenberg(A_hess; kwds...)
+    p = size(_pertdata(bandeddata(parent(L))), 2) + u + 1 # pert size
+    dat = (UpperHessenbergQ((Q').q[1:(p+l)])) * A[1:p+l+1, 1:p]
+    pert = Array{eltype(dat)}(undef, l + u + 1, size(dat, 2) - 1)
     for j = 1:u
-        pert[u-j+1:end,j] .= view(dat,1:l+j+1,j)
+        pert[u-j+1:end, j] .= view(dat, 1:l+j+1, j)
     end
-    for j = u+1:size(pert,2)
-        pert[:,j] .= view(dat,j-u+1:j+l+1,j)
+    for j = u+1:size(pert, 2)
+        pert[:, j] .= view(dat, j-u+1:j+l+1, j)
     end
-    H = _BandedMatrix(Hcat(pert, dat[end-l-u:end,end]*Ones{eltype(dat)}(1,∞)), ℵ₀, l+1,u-1)
-    Q,H
+    H = _BandedMatrix(Hcat(pert, dat[end-l-u:end, end] * Ones{eltype(dat)}(1, ∞)), ℵ₀, l + 1, u - 1)
+    Q, H
 end
 
 # represent Q as a product of orthogonal operations
-struct ProductQ{T,QQ<:Tuple} <: AbstractQ{T}
+struct ProductQ{T,QQ<:Tuple} <: LayoutQ{T}
     Qs::QQ
 end
 
 ArrayLayouts.@layoutmatrix ProductQ
 ArrayLayouts.@_layoutlmul ProductQ
 
-ProductQ(Qs::AbstractMatrix...) = ProductQ{mapreduce(eltype,promote_type,Qs),typeof(Qs)}(Qs)
+ProductQ(Qs::AbstractMatrix...) = ProductQ{mapreduce(eltype, promote_type, Qs),typeof(Qs)}(Qs)
 
-adjoint(Q::ProductQ) = ProductQ(reverse(map(adjoint,Q.Qs))...)
+adjoint(Q::ProductQ) = ProductQ(reverse(map(adjoint, Q.Qs))...)
 
 size(Q::ProductQ, dim::Integer) = size(dim == 1 ? Q.Qs[1] : last(Q.Qs), dim == 2 ? 1 : dim)
+axes(Q::ProductQ, dim::Integer) = axes(dim == 1 ? Q.Qs[1] : last(Q.Qs), dim == 2 ? 1 : dim)
 
 function lmul!(Q::ProductQ, v::AbstractVecOrMat)
     for j = length(Q.Qs):-1:1
@@ -104,12 +105,19 @@ function lmul!(Q::ProductQ, v::AbstractVecOrMat)
     v
 end
 
-getindex(Q::ProductQ{<:Any,<:Tuple{Vararg{LowerHessenbergQ}}}, i::Integer, j::Integer) = (Q')[j,i]'
+# Avoid ambiguities
+getindex(Q::ProductQ, i::Int, j::Int) = Q[:, j][i]
 
+function getindex(Q::ProductQ, ::Colon, j::Int)
+    y = zeros(eltype(Q), size(Q, 2))
+    y[j] = 1
+    lmul!(Q, y)
+end
+getindex(Q::ProductQ{<:Any,<:Tuple{Vararg{LowerHessenbergQ}}}, i::Int, j::Int) = (Q')[j, i]'
 
 function _productq_mul(A::ProductQ{T}, x::AbstractVector{S}) where {T,S}
     TS = promote_op(matprod, T, S)
-    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
+    lmul!(A, Base.copymutable(convert(AbstractVector{TS}, x)))
 end
 
 mul(A::ProductQ, x::AbstractVector) = _productq_mul(A, x)
@@ -123,7 +131,7 @@ end
 
 
 
-QLProduct(Qs::Tuple, L::AbstractMatrix{T}) where T = QLProduct{T,typeof(Qs),typeof(L)}(Qs, L)
+QLProduct(Qs::Tuple, L::AbstractMatrix{T}) where {T} = QLProduct{T,typeof(Qs),typeof(L)}(Qs, L)
 QLProduct(F::QLHessenberg) = QLProduct(tuple(F.Q), F.L)
 
 # iteration for destructuring into components
@@ -140,7 +148,8 @@ Matrix(F::QLProduct) = Array(AbstractArray(F))
 Array(F::QLProduct) = Matrix(F)
 
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::QLProduct)
-    summary(io, F); println(io)
+    summary(io, F)
+    println(io)
     println(io, "Q factor:")
     show(io, mime, F.Q)
     println(io, "\nL factor:")
@@ -163,11 +172,11 @@ end
 Base.propertynames(F::QLProduct, private::Bool=false) =
     (:L, :Q, (private ? fieldnames(typeof(F)) : ())...)
 
-function _inf_ql(A::AbstractMatrix{T}; kwds...) where T
-    _,u = bandwidths(A)
+function _inf_ql(A::AbstractMatrix{T}; kwds...) where {T}
+    _, u = bandwidths(A)
     u ≤ 0 && return QLProduct(tuple(Eye{float(T)}(∞)), A)
     u == 1 && return QLProduct(ql_hessenberg(A; kwds...))
-    Q1,H1 = ql_pruneband(A; kwds...)
+    Q1, H1 = ql_pruneband(A; kwds...)
     F̃ = ql(H1; kwds...)
     QLProduct(tuple(Q1, F̃.Qs...), F̃.L)
 end

src/infcholesky.jl

@@ -45,7 +45,7 @@ end
 adaptivecholesky(A) = Cholesky(AdaptiveCholeskyFactors(A), :U, 0)
 
 
-ArrayLayouts._cholesky(::SymmetricLayout{<:AbstractBandedLayout}, ::NTuple{2,OneToInf{Int}}, A) = adaptivecholesky(A)
+ArrayLayouts._cholesky(::SymmetricLayout{<:AbstractBandedLayout}, ::NTuple{2,OneToInf{Int}}, A, ::CNoPivot) = adaptivecholesky(A)
 
 function colsupport(F::AdaptiveCholeskyFactors, j)
     partialcholesky!(F, maximum(j)+bandwidth(F,2))

src/infql.jl

@@ -96,7 +96,7 @@ function ql_hessenberg!(B::InfBandedMatrix{TT}; kwds...) where TT
     QLHessenberg(_BandedMatrix(H, ℵ₀, l, 1), Vcat( LowerHessenbergQ(F.Q).q, F∞.q))
 end
 
-getindex(Q::QLPackedQ{T,<:InfBandedMatrix{T}}, i::Integer, j::Integer) where T =
+getindex(Q::QLPackedQ{T,<:InfBandedMatrix{T}}, i::Int, j::Int) where T =
     (Q'*[Zeros{T}(i-1); one(T); Zeros{T}(∞)])[j]'
 
 getL(Q::QL, ::NTuple{2,InfiniteCardinal{0}}) where T = LowerTriangular(Q.factors)

test/test_hessenbergq.jl

@@ -3,49 +3,49 @@ import InfiniteLinearAlgebra: UpperHessenbergQ, LowerHessenbergQ, tail_de, _Band
 
 @testset "HessenbergQ" begin
     @testset "finite UpperHessenbergQ" begin
-        c,s = cos(0.1), sin(0.1); q = [c s; s -c]; 
+        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
         Q = UpperHessenbergQ(Fill(q,1))
         @test size(Q) == (size(Q,1),size(Q,2)) == (2,2)
         @test Q ≈ q
         Q = UpperHessenbergQ(Fill(q,2))
         @test bandwidths(Q) == (1,2)
         @test Q ≈ [q zeros(2); zeros(1,2) 1] * [1 zeros(1,2); zeros(2) q]
         @test Q' isa LowerHessenbergQ
-        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1] 
+        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1]
 
         q = qr(randn(2,2) + im*randn(2,2)).Q
         Q = UpperHessenbergQ(Fill(q,1))
         @test Q ≈ q
         Q = UpperHessenbergQ(Fill(q,2))
         @test bandwidths(Q) == (1,2)
-        @test Q ≈ [q zeros(2); zeros(1,2) 1] * [1 zeros(1,2); zeros(2) q] 
+        @test Q ≈ [q zeros(2); zeros(1,2) 1] * [1 zeros(1,2); zeros(2) q]
         @test Q' isa LowerHessenbergQ
-        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1] 
+        @test Q' ≈ [1 zeros(1,2); zeros(2) q'] * [q' zeros(2); zeros(1,2) 1]
     end
 
     @testset "finite LowerHessenbergQ" begin
-        c,s = cos(0.1), sin(0.1); q = [c s; s -c]; 
+        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
         Q = LowerHessenbergQ(Fill(q,1))
         @test size(Q) == (size(Q,1),size(Q,2)) == (2,2)
         @test Q ≈ q
         Q = LowerHessenbergQ(Fill(q,2))
         @test bandwidths(Q) == (2,1)
-        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1] 
+        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1]
         @test Q' isa UpperHessenbergQ
-        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q'] 
+        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q']
 
         q = qr(randn(2,2) + im*randn(2,2)).Q
         Q = LowerHessenbergQ(Fill(q,1))
         @test Q ≈ q
         Q = LowerHessenbergQ(Fill(q,2))
         @test bandwidths(Q) == (2,1)
-        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1] 
+        @test Q ≈ [1 zeros(1,2); zeros(2) q] * [q zeros(2); zeros(1,2) 1]
         @test Q' isa UpperHessenbergQ
-        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q'] 
+        @test Q' ≈ [q' zeros(2); zeros(1,2) 1]  * [1 zeros(1,2); zeros(2) q']
     end
 
     @testset "infinite-Q" begin
-        c,s = cos(0.1), sin(0.1); q = [c s; s -c]; 
+        c,s = cos(0.1), sin(0.1); q = [c s; s -c];
         Q = LowerHessenbergQ(Fill(q,∞))
         @test Q[1,1] ≈ 0.9950041652780258
         Q = UpperHessenbergQ(Fill(q,∞))
@@ -58,7 +58,7 @@ import InfiniteLinearAlgebra: UpperHessenbergQ, LowerHessenbergQ, tail_de, _Band
             Q,R = qr(A)
             @test UpperHessenbergQ(Q) ≈ Q
             Q,L = ql(A)
-            @test LowerHessenbergQ(Q) ≈ Q            
+            @test LowerHessenbergQ(Q) ≈ Q
         end
     end
-end 
+end

test/test_infqr.jl

@@ -188,7 +188,7 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
             @test qr(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
             @test factorize(A) isa MatrixFactorizations.QR{Float64,<:InfiniteLinearAlgebra.AdaptiveQRFactors}
 
-            @test (adaptiveqr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (qr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (A \ [ℯ; zeros(∞)])[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test (adaptiveqr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (qr(A) \ [ℯ; zeros(∞)])[1:1000] ≈ (A \ [ℯ; zeros(∞)])[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
         end
 
         @testset "two-bands" begin
@@ -202,22 +202,22 @@ import SemiseparableMatrices: AlmostBandedLayout, VcatAlmostBandedLayout
             x = 0.1
             @test (exp(1 - x)*(-1 + exp(2 + 2x)))/(-1 + exp(4)) ≈ dot(u[1:1000], x.^(0:999))
             u = qr(A) \ Vcat([ℯ,1/ℯ], zeros(∞))
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
             u = qr(A) \ Vcat(ℯ,1/ℯ, zeros(∞))
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
             u = qr(A) \ [ℯ; 1/ℯ; zeros(∞)]
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
 
             A = Vcat(Ones(1,∞), ((-1.0).^(0:∞))', B)
             @test MemoryLayout(typeof(A)) isa VcatAlmostBandedLayout
             u = A \ Vcat(ℯ,1/ℯ, zeros(∞))
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
             u = qr(A) \ [ℯ; 1/ℯ; zeros(∞)]
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
 
             A = Vcat(Ones(1,∞), ((-1).^(0:∞))', B)
             u = A \ [ℯ; 1/ℯ; zeros(∞)]
-            @test u[1:1000] ≈ [1/factorial(1.0k) for k=0:999]
+            @test u[1:1000] ≈ [1/gamma(k+1.0) for k=0:999]
         end
 
         @testset "more bands" begin