Use function barriers for MulAddAccumulate

src/SemiclassicalOrthogonalPolynomials.jl

@@ -7,7 +7,7 @@ import Base: getindex, axes, size, \, /, *, +, -, summary, show, ==, copy, sum,
 
 import ArrayLayouts: MemoryLayout, ldiv, diagonaldata, subdiagonaldata, supdiagonaldata
 import BandedMatrices: bandwidths, AbstractBandedMatrix, BandedLayout, _BandedMatrix
-import LazyArrays: resizedata!, paddeddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, ApplyArray,
+import LazyArrays: resizedata!, paddeddata, maybe_cacheddata, CachedVector, CachedMatrix, CachedAbstractVector, LazyMatrix, LazyVector, arguments, ApplyLayout, colsupport, AbstractCachedVector, ApplyArray,
                     AccumulateAbstractVector, LazyVector, AbstractCachedMatrix, BroadcastLayout, simplifiable
 import ClassicalOrthogonalPolynomials: OrthogonalPolynomial, recurrencecoefficients, jacobimatrix, normalize, _p0, UnitInterval, orthogonalityweight, NormalizedOPLayout, MappedOPLayout,
                                         Bidiagonal, Tridiagonal, SymTridiagonal, symtridiagonalize, normalizationconstant, LanczosPolynomial,

src/derivatives.jl

@@ -1,7 +1,7 @@
 """
-MulAddAccumulate(μ, A, B)
+    MulAddAccumulate(μ, A, B)
 
-represents the vector satisfying v[k+1] == A[k]*v[k]+B[k] with v[1] == μ
+represents the vector satisfying v[k+1] == A[k+1]*v[k]+B[k+1] with v[1] == μ
 """
 mutable struct MulAddAccumulate{T} <: AbstractCachedVector{T}
     data::Vector{T}
@@ -21,11 +21,17 @@ MulAddAccumulate(data::Vector{T}, A::AbstractVector, B::AbstractVector, datasize
 MulAddAccumulate(μ, A, B) = MulAddAccumulate([μ], A, B, (1,))
 MulAddAccumulate(A, B) = MulAddAccumulate(A[1]+B[1], A, B)
 
-function LazyArrays.cache_filldata!(K::MulAddAccumulate, inds)
-    A,B = K.A,K.B
+LazyArrays.cache_filldata!(K::MulAddAccumulate, inds) = _muladdaccumulate_filldata!(K.data, K.A, K.B, inds)
+function _muladdaccumulate_filldata!(K, A, B, inds)
+    n = maximum(inds)
+    resizedata!(A, n)
+    resizedata!(B, n)
+    __muladdaccumulate_filldata!(K, maybe_cacheddata(A), maybe_cacheddata(B), inds)
+end
+function __muladdaccumulate_filldata!(data, A, B, inds)
     @inbounds for k in inds
-        K.data[k] = muladd(A[k], K.data[k-1], B[k])
-    end
+        data[k] = muladd(A[k], data[k-1], B[k])
+    end 
 end
 
 """

test/test_derivative.jl

@@ -1,7 +1,15 @@
 using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, LazyArrays, Test
 import ClassicalOrthogonalPolynomials: recurrencecoefficients, _BandedMatrix, _p0, Weighted
 import LazyArrays: Accumulate, AccumulateAbstractVector
-import SemiclassicalOrthogonalPolynomials: MulAddAccumulate, HalfWeighted, toclassical
+import SemiclassicalOrthogonalPolynomials: MulAddAccumulate, HalfWeighted, toclassical, MulAddAccumulate
+
+@testset "MulAddAccumulate" begin 
+    A = Vcat(1, cache(2:3), cache(4:∞) ./ cache(3:∞))
+    B = Vcat(cache(1:5), [1, 2], cache(3:∞))
+    v = MulAddAccumulate(A, B)
+    @test v[1] ≈ 2
+    @test all(k -> v[k+1] ≈ A[k+1] * v[k] + B[k+1], 2:1000)
+end
 
 @testset "Derivative" begin
     @testset "basics" begin