diff --git a/Project.toml b/Project.toml
index ff6bdda..1636f13 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,20 +17,28 @@ Aqua = "0.8"
 ArrayLayouts = "0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1"
 BandedMatrices = "0.15, 0.16, 0.17, 1"
 CircularArrays = "1"
+ClassicalOrthogonalPolynomials = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.10, 0.11, 0.12"
 Documenter = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 1"
+InfiniteArrays = "0.10, 0.11, 0.12, 0.13"
 Infinities = "0.1"
 LazyArrays = "0.19, 0.20, 0.21, 0.22, 1"
 LinearAlgebra = "0, 1"
 QuadGK = "1, 2"
+SpecialFunctions = "1, 2"
 StaticArrays = "0.8, 0.9, 0.10, 0.11, 0.12, 1"
 Test = "1"
 julia = "1.5"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Aqua", "Documenter", "QuadGK"]
+test = ["Test", "Aqua", "BandedMatrices", "ClassicalOrthogonalPolynomials", "Documenter", "InfiniteArrays", "LinearAlgebra", "QuadGK", "SpecialFunctions"]
diff --git a/src/BandedSylvesters.jl b/src/BandedSylvesters.jl
index a654a2e..79365be 100644
--- a/src/BandedSylvesters.jl
+++ b/src/BandedSylvesters.jl
@@ -2,6 +2,8 @@ import ArrayLayouts: colsupport, rowsupport
 import BandedMatrices: bandwidth
 import Base: size
 
+export BandedSylvesterRecurrencePlan, BandedSylvesterRecurrence
+
 """
     BandedSylvesterRecurrence{T, TA<:AbstractMatrix{T}, TB<:AbstractMatrix{T}, TC<:AbstractMatrix{T}, TX<:AbstractMatrix{T}} <: AbstractLinearRecurrence{slicetype(TX)}
 
@@ -56,28 +58,31 @@ size(P::BandedSylvesterRecurrencePlan) = P.size
 
 function init(P::BandedSylvesterRecurrencePlan; T=Float64, init=:default)
     if init == :default
-        buffer = tocircular(P.init(T, size(P)...))
+        buffer = tocircular(method_to_precision(P.init, (Int,))(T, size(P)[1]))
     elseif init == :rand
-        buffer = CircularArray(rand(T, front(size(P))..., P.firstind))
+        buffer = CircularArray(rand(T, size(P)[1], +(P.bandB...)))
     end
-    A = P.fA(T)
-    B = P.fB(T)
-    BandedSylvesterRecurrence(A, B, buffer, P.firstind, P.size[end])
+    A = method_to_precision(P.fA,())(T)
+    B = method_to_precision(P.fB,())(T)
+    C = method_to_precision(P.fC,())(T)
+    BandedSylvesterRecurrence(A, B, C, buffer, P.firstind, P.size[end]), eachslice(view(buffer.data, fill(:, ndims(buffer)-1)..., axes(buffer)[end][1:P.firstind-1]), dims=ndims(buffer))
 end
 
 function step!(R::BandedSylvesterRecurrence)
+    # A[m,:]X[:,n] + X[m,:]B[:,n] + C[m,n] = 0
     if R.sliceind > R.lastind
         return nothing
     end
-    A, B, X = R.A, R.B, R.buffer
+    A, B, C, X = R.A, R.B, R.C, R.buffer
     nx = R.sliceind
     n = nx - bandwidth(B,1)
     ind = axes(X, 1)
     Bs = colsupport(B, n)
-    Bv = B[(Base.OneTo(nx-1)) ∩ Bs, n]
+    I2 = (Base.OneTo(nx-1)) ∩ Bs
+    Bv = B[I2, n]
     for m in ind
         I1 = rowsupport(A, m) ∩ ind
-        X[m, nx] = (_dotu(view(A,m,I1), view(X,I1,n)) + _dotu(view(X,m,I2), Bv) + C[nx,n]) / B[nx,n]
+        X[m, nx] = -(_dotu(view(A,m,I1), view(X,I1,n)) + _dotu(view(X,m,I2), Bv) + C[m,n]) / B[nx,n]
     end
     R.sliceind += 1
     return view(X, :, nx)
diff --git a/src/EltypeExtensions.jl b/src/EltypeExtensions.jl
index d0d4fe0..19dbd8a 100644
--- a/src/EltypeExtensions.jl
+++ b/src/EltypeExtensions.jl
@@ -123,6 +123,7 @@ julia> precisionconvert(Float16, [[m/n for n in 1:3] for m in 1:3])
 """
 precisionconvert(T,A) = precisionconvert(T,A,precision(T))
 precisionconvert(::Type{T}, A::S, prec) where {T,S} = convert(_to_precisiontype(T,S), A)
+precisionconvert(::Type{BigFloat}, A::S) where {S} = convert(_to_precisiontype(BigFloat,S), A)
 precisionconvert(::Type{BigFloat}, x::Real, prec) = BigFloat(x, precision=prec)
 precisionconvert(::Type{BigFloat}, x::Complex, prec) = Complex(BigFloat(real(x), precision=prec), BigFloat(imag(x), precision=prec))
 precisionconvert(::Type{BigFloat}, A, prec) = precisionconvert.(BigFloat, A, precision=prec)
\ No newline at end of file
diff --git a/test/fractional_integral_operator.jl b/test/fractional_integral_operator.jl
new file mode 100644
index 0000000..bf01dd8
--- /dev/null
+++ b/test/fractional_integral_operator.jl
@@ -0,0 +1,47 @@
+using ClassicalOrthogonalPolynomials
+using LinearAlgebra
+using BandedMatrices
+using SpecialFunctions
+using InfiniteArrays
+
+OpR(γ,δ,α,β) = Jacobi(γ,δ) \ Jacobi(α,β); # R_{(α,β)}^{(γ,δ)}
+OpL(α,β,k,j) = Jacobi(α,β) \ (JacobiWeight(k,j).*Jacobi(α+k,β+j)); # L_{(α+k,β+j)}^{(α,β)}
+OpW(β) = Diagonal(β+1:∞); # W^{(β)}
+OpInvW(β) = Diagonal(inv.(β+1:∞)); # W^{(β)}^{-1}
+OpΛ(r,μ,k) = BandedMatrix(-k=>(gamma.((r+1):μ:∞)./gamma.((r+k*μ+1):μ:∞)))
+
+OpI(α,β,b,p) = p//2^(p-1)*OpL(α,β,0,p)*OpR(α,β+p,α-1,b+p)*OpInvW(b+p-1)*OpR(α,b+p-1,α,β); # I^{(α,β)}_{b,p}
+function OpI(α,β,b,p,N) 
+    L = p>1 ? *(map(x->x[1:N,1:N],OpL(α,β,0,p).args)...) : OpL(α,β,0,p)[1:N,1:N]
+    R1 = b!=β ? *(map(x->x[1:N,1:N],OpR(α,β+p,α-1,b+p).args)...) : OpR(α,β+p,α-1,b+p)[1:N,1:N]
+    R2 = b+p-1-β>1 ? *(map(x->x[1:N,1:N],OpR(α,b+p-1,α,β).args)...) : OpR(α,b+p-1,α,β)[1:N,1:N]
+    p/2^(p-1)*L*R1*Diagonal(OpInvW(b+p-1).diag[1:N])*R2
+end
+
+fracpochhammer(a,b,n) = n==0 ? one(promote_type(typeof(a),typeof(b))) : prod(x/y for (x,y) in zip(range(a,length=n),range(b,length=n)))
+fracpochhammer(a,b,stepa,stepb,n) = n==0 ? one(promote_type(typeof(a),typeof(b))) : prod(x/y for (x,y) in zip(range(a,step=stepa,length=n),range(b,step=stepb,length=n)))
+OpCsingle(α,β,k,n) = (-1)^(n-k)*fracpochhammer(k+β+1,one(β),n-k)*fracpochhammer(n+α+β+1,2*one(β),1,2,k); # [k,n]
+function OpCcolumn(α,β,n) # [:,n]
+    (α,β)=promote(α,β)
+    ret=zeros(typeof(α),n+1)
+    ret[1]=(-1)^n*fracpochhammer(β+1,1,n);
+    for k=1:n
+        ret[k+1]=ret[k]*(n+α+β+k)*(n-k+1)/(-2*k*(k+β))
+    end
+    ret
+end
+function OpC(α,β,N)
+    (α,β)=promote(α,β)
+    ret=zeros(typeof(α),N+1,N+1)
+    for n=0:N
+        ret[1:n+1,n+1] = OpCcolumn(α,β,n)
+    end
+    ret
+end
+
+function OpI11(α,β,b,p,μ,N)
+    k=Int(round(p*μ))
+    (α,β,b,p,μ) = promote(α,β,b,p,μ)
+    A = OpC(α,β,N+k)
+    return A\(2^(μ-k)*OpΛ(b/p,1/p,k)[1:N+k+1,1:N+1]*A[1:N+1,1:N+1])
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 615d18a..ba4747f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -29,6 +29,23 @@ using QuadGK
         ret = hcat(stable_recurrence(P)...)
         @test ret[92:end, 92:end] ≈ [quadgk(x->(cos(x)/(1+sin(x)))^m*exp(im*n*x), 0, π)[1] for m in 91:100, n in 91:100]
     end
+    @testset "Fractional integral operator" begin
+        include("fractional_integral_operator.jl")
+        α, β, b, p, μ = 0.0, 0.0, 0.0, 2, 1//2
+        N = 100
+        k = Int(μ*p)
+        fA(T) = OpI(T(α), T(β), T(b), T(p), N+k+p+1)
+        fB(T) = -fA(T)
+        fC(T) = Zeros{T}(∞,∞)
+        function init(T, N)
+            ret = zeros(T, N, 2*p+1)
+            ret[1:k+2, 1:k+1] = OpI11(T(α), T(β), T(b), T(p), T(μ), k)
+            ret
+        end
+        P = BandedSylvesterRecurrencePlan(fA, fB, fC, init, (N+k+1,N+1), (p, p))
+        FIO = hcat(stable_recurrence(P)...)
+        @test FIO[1:N,:] * FIO[1:N+1,1:N] ≈ OpI(α,β,b,p,N)
+    end
 end
 
 @testset "Extensions" begin