diff --git a/Project.toml b/Project.toml
index 7468c32..ee276ed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,8 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
+version = "0.15.8"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.15.7"
-
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -38,7 +37,7 @@ ArrayLayouts = "1.3.1"
 BandedMatrices = "1"
 BlockArrays = "1"
 BlockBandedMatrices = "0.13"
-ContinuumArrays = "0.19"
+ContinuumArrays = "0.20"
 DomainSets = "0.6, 0.7"
 DynamicPolynomials = "0.6"
 FFTW = "1.1"
@@ -52,10 +51,11 @@ IntervalSets = "0.7"
 LazyArrays = "2.8"
 LazyBandedMatrices = "0.11"
 MutableArithmetics = "1"
-QuasiArrays = "0.12"
+QuasiArrays = "0.13"
 RecurrenceRelationshipArrays = "0.1.2"
 RecurrenceRelationships = "0.2"
 SpecialFunctions = "1.0, 2"
+StatsBase = "0.34.6"
 julia = "1.10"
 
 [extras]
@@ -65,7 +65,8 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "StaticArrays", "Random", "DynamicPolynomials"]
+test = ["Base64", "DynamicPolynomials", "ForwardDiff", "Random", "SemiseparableMatrices", "StaticArrays", "StatsBase", "Test"]
diff --git a/examples/detpointprocess.jl b/examples/detpointprocess.jl
index e111538..642b58a 100644
--- a/examples/detpointprocess.jl
+++ b/examples/detpointprocess.jl
@@ -1,5 +1,4 @@
-using ClassicalOrthogonalPolynomials, Plots
-using ClassicalOrthogonalPolynomials: sample
+using ClassicalOrthogonalPolynomials, StatsBase, Plots
 
 x = Inclusion(ChebyshevInterval())
 
diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl
index 6a7db81..8bfbb69 100644
--- a/src/ClassicalOrthogonalPolynomials.jl
+++ b/src/ClassicalOrthogonalPolynomials.jl
@@ -30,7 +30,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu
                     ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                     LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                     _getindex, layout_getindex, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector,
-                    AbstractQuasiFill, equals_layout, QuasiArrayLayout, PolynomialLayout, diff_layout
+                    AbstractQuasiFill, equals_layout, QuasiArrayLayout, PolynomialLayout, diff_layout, findall_layout, minimum_layout, maximum_layout, extrema_layout
 
 import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
 import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!, partialcholesky!, SymTridiagonalConjugation, TridiagonalConjugation
diff --git a/src/roots.jl b/src/roots.jl
index 7d0210c..fd70879 100644
--- a/src/roots.jl
+++ b/src/roots.jl
@@ -18,58 +18,26 @@ function colleaguematrix(P, c)
 end
 
 
-function _findall(::typeof(iszero), ::ExpansionLayout{<:AbstractOPLayout}, f)
+function findall_layout(::ExpansionLayout{<:AbstractOPLayout}, ::typeof(iszero), f)
     C = colleaguematrix(f.args...)
     ax = axes(f,1)
     convert(Vector{eltype(ax)}, filter!(in(ax), eigvals(C)))
 end
-findall(f::Function, v::AbstractQuasiVector) = _findall(f, MemoryLayout(v), v)
-
-# gives a generalization of midpoint for when `a` or `b` is infinite
-function genmidpoint(a::T, b::T) where T
-    if isinf(a) && isinf(b)
-        zero(T)
-    elseif isinf(a)
-        b - 100
-    elseif isinf(b)
-        a + 100
-    else
-        (a+b)/2
-    end
-end
-
-
-function _searchsortedfirst(::ExpansionLayout{<:AbstractOPLayout}, f, x; iterations=47)
-    d = axes(f,1)
-    a,b = first(d), last(d)
-
-    for k=1:iterations  #TODO: decide 47
-        m= genmidpoint(a,b)
-        (f[m] ≤ x) ? (a = m) : (b = m)
-    end
-    (a+b)/2
-end
-searchsortedfirst(f::AbstractQuasiVector, x; kwds...) = _searchsortedfirst(MemoryLayout(f), f, x; kwds...)
-
-function sample(f::AbstractQuasiVector, n...)
-    g = cumsum(f)
-    searchsortedfirst.(Ref(g/last(g)), rand(n...))
-end
 
 ####
 # min/max/extrema
 ####
-function minimum(f::AbstractQuasiVector)
+function minimum_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims)
     r = findall(iszero, diff(f))
     min(first(f), minimum(f[r]), last(f))
 end
 
-function maximum(f::AbstractQuasiVector)
+function maximum_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims)
     r = findall(iszero, diff(f))
     max(first(f), maximum(f[r]), last(f))
 end
 
-function extrema(f::AbstractQuasiVector)
+function extrema_layout(::ExpansionLayout{<:AbstractOPLayout}, f::AbstractQuasiVector, dims...)
     r = findall(iszero, diff(f))
     extrema([first(f); f[r]; last(f)])
 end
\ No newline at end of file
diff --git a/test/test_roots.jl b/test/test_roots.jl
index f316e73..5651a8a 100644
--- a/test/test_roots.jl
+++ b/test/test_roots.jl
@@ -1,5 +1,4 @@
-using ClassicalOrthogonalPolynomials, QuasiArrays, Random, Test
-using ClassicalOrthogonalPolynomials: sample
+using ClassicalOrthogonalPolynomials, QuasiArrays, Random, StatsBase, Test
 
 Random.seed!(5)