diff --git a/Project.toml b/Project.toml
index 2959561e..d45d668b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.17"
-
+version = "0.18.0"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -19,6 +18,13 @@ RecurrenceRelationships = "807425ed-42ea-44d6-a357-6771516d7b2c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
+[weakdeps]
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+
+[extensions]
+FastTransformsInfiniteArraysExt = "InfiniteArrays"
+
+
 [compat]
 AbstractFFTs = "1.0"
 ArrayLayouts = "1.10"
@@ -28,15 +34,17 @@ FastGaussQuadrature = "0.4, 0.5, 1"
 FastTransforms_jll = "0.6.2"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
 GenericFFT = "0.1"
+InfiniteArrays = "0.15"
 LazyArrays = "2.2"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "0.10, 1, 2"
 ToeplitzMatrices = "0.7.1, 0.8"
-julia = "1.7"
+julia = "1.10"
 
 [extras]
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random"]
+test = ["InfiniteArrays", "Random", "Test"]
diff --git a/src/GramMatrix.jl b/src/GramMatrix.jl
index 8bf6dc42..98d66378 100644
--- a/src/GramMatrix.jl
+++ b/src/GramMatrix.jl
@@ -59,35 +59,13 @@ In the standard (classical) normalization, ``p_0(x) = 1``, so that the moments
 The recurrence is built from ``X^\\top W = WX``.
 """
 GramMatrix(μ::AbstractVector{T}, X::XT) where {T, XT <: AbstractMatrix{T}} = GramMatrix(μ, X, one(T))
-function GramMatrix(μ::AbstractVector{T}, X::XT, p0::T) where {T, XT <: AbstractMatrix{T}}
-    N = length(μ)
-    n = (N+1)÷2
-    @assert N == size(X, 1) == size(X, 2)
-    @assert bandwidths(X) == (1, 1)
-    W = LowerTriangular(Matrix{T}(undef, N, N))
-    if n > 0
-        @inbounds for m in 1:N
-            W[m, 1] = p0*μ[m]
-        end
-    end
-    if n > 1
-        @inbounds for m in 2:N-1
-            W[m, 2] = (X[m-1, m]*W[m-1, 1] + (X[m, m]-X[1, 1])*W[m, 1] + X[m+1, m]*W[m+1, 1])/X[2, 1]
-        end
-    end
-    @inbounds @simd for n in 3:n
-        for m in n:N-n+1
-            W[m, n] = (X[m-1, m]*W[m-1, n-1] + (X[m, m]-X[n-1, n-1])*W[m, n-1] + X[m+1, m]*W[m+1, n-1] - X[n-2, n-1]*W[m, n-2])/X[n, n-1]
-        end
-    end
-    return GramMatrix(Symmetric(W[1:n, 1:n], :L), eval(XT.name.name)(view(X, 1:n, 1:n)))
-end
 
-function GramMatrix(μ::PaddedVector{T}, X::XT, p0::T) where {T, XT <: AbstractMatrix{T}}
-    N = length(μ)
-    b = length(μ.args[2])-1
+
+function GramMatrix(μ::AbstractVector{T}, X::XT, p0::T) where {T, XT <: AbstractMatrix{T}}
+    N = size(X, 1)
+    b = length(μ)-1
     n = (N+1)÷2
-    @assert N == size(X, 1) == size(X, 2)
+    @assert N == size(X, 2)
     @assert bandwidths(X) == (1, 1)
     W = BandedMatrix{T}(undef, (N, N), (b, 0))
     if n > 0
@@ -100,12 +78,12 @@ function GramMatrix(μ::PaddedVector{T}, X::XT, p0::T) where {T, XT <: AbstractM
             W[m, 2] = (X[m-1, m]*W[m-1, 1] + (X[m, m]-X[1, 1])*W[m, 1] + X[m+1, m]*W[m+1, 1])/X[2, 1]
         end
     end
-    @inbounds @simd for n in 3:n
-        for m in n:min(N-n+1, b+n)
-            W[m, n] = (X[m-1, m]*W[m-1, n-1] + (X[m, m]-X[n-1, n-1])*W[m, n-1] + X[m+1, m]*W[m+1, n-1] - X[n-2, n-1]*W[m, n-2])/X[n, n-1]
+    @inbounds @simd for j in 3:N
+        for m in j:min(N, b+j)
+            W[m, j] = (X[m-1, m]*W[m-1, j-1] + (X[m, m]-X[j-1, j-1])*W[m, j-1] + (m < N ? X[m+1, m]*W[m+1, j-1] : zero(T)) - X[j-2, j-1]*W[m, j-2])/X[j, j-1]
         end
     end
-    return GramMatrix(Symmetric(W[1:n, 1:n], :L), eval(XT.name.name)(view(X, 1:n, 1:n)))
+    return GramMatrix(Symmetric(W, :L), X)
 end
 
 #
diff --git a/test/grammatrixtests.jl b/test/grammatrixtests.jl
index 3ac9d626..a017ddd0 100644
--- a/test/grammatrixtests.jl
+++ b/test/grammatrixtests.jl
@@ -23,12 +23,14 @@ using FastTransforms, BandedMatrices, LazyArrays, LinearAlgebra, Test
         X = BandedMatrix(SymTridiagonal(zeros(T, n+b), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n+b-1])) # normalized Legendre X
         W = I+X^2+X^4
         W = Symmetric(W[1:n, 1:n])
-        X = BandedMatrix(SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1])) # normalized Legendre X
-        G = GramMatrix(W, X)
+        X̃ = BandedMatrix(SymTridiagonal(zeros(T, n), [sqrt(T(n)^2/(4*n^2-1)) for n in 1:n-1])) # normalized Legendre X
+        G = GramMatrix(W, X̃)
         @test bandwidths(G) == (b, b)
         F = cholesky(G)
         @test F.L*F.L' ≈ W
 
+        @test G ≈ GramMatrix(W[1:5, 1], X̃)
+
         X = BandedMatrix(SymTridiagonal(T[2n-1 for n in 1:n+b], T[-n for n in 1:n+b-1])) # Laguerre X, tests nonzero diagonal
         W = I+X^2+X^4
         W = Symmetric(W[1:n, 1:n])
diff --git a/test/inffasttransformstests.jl b/test/inffasttransformstests.jl
new file mode 100644
index 00000000..e69de29b