diff --git a/Project.toml b/Project.toml
index fd7d007..62afe87 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BasisFunctions"
 uuid = "4343a256-5453-507d-8aad-01a9d7189916"
 authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
-version = "0.6.6"
+version = "0.7"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -40,7 +40,7 @@ BasisFunctionsPGFPlotsXExt = "PGFPlotsX"
 BlockArrays = "0.16"
 Calculus = "0.5"
 CompositeTypes = "0.1.3"
-DomainIntegrals = "0.4.5"
+DomainIntegrals = "0.5"
 DomainSets = "0.7"
 DoubleFloats = "1"
 FFTW = "1.6"
@@ -63,7 +63,7 @@ SpecialFunctions = "1.0, 2.0"
 StaticArrays = "1.4"
 Test = "1"
 ToeplitzMatrices = "0.7"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
diff --git a/src/BasisFunctions.jl b/src/BasisFunctions.jl
index 21e58b3..9d2adac 100644
--- a/src/BasisFunctions.jl
+++ b/src/BasisFunctions.jl
@@ -56,7 +56,6 @@ import LinearAlgebra:
     norm, pinv, normalize, cross, ×, dot, adjoint, mul!, rank,
     diag, isdiag, eigvals, issymmetric, svdvals
 import SparseArrays: sparse
-@deprecate sparse_matrix sparse
 
 
 import BlockArrays.BlockVector
@@ -158,6 +157,7 @@ export GenericFunctionSpace, space, MeasureSpace, FourierSpace, ChebyshevTSpace,
 # from bases/dictionary/dictionary.jl
 export Dictionary, Dictionary1d, Dictionary2d, Dictionary3d,
     interpolation_grid, left, right, domain, codomain,
+    domaintype, codomaintype, coefficienttype,
     measure, hasmeasure,
     eval_expansion, eval_element, eval_element_derivative, eval_gradient,
     resize,
@@ -356,8 +356,13 @@ include("bases/dictionary/expansions.jl")
 
 include("util/products.jl")
 
-include("computations/generic_operators.jl")
+include("computations/transform.jl")
+include("computations/evaluation.jl")
+include("computations/approximation.jl")
+include("computations/differentiation.jl")
 include("computations/conversion.jl")
+include("computations/generic_operators.jl")
+include("computations/fun.jl")
 
 ################################################################
 # Generic dictionaries
@@ -418,4 +423,7 @@ include("util/plot.jl")
 include("test/Test.jl")
 
 include("examples/pwconstants.jl")
+
+include("deprecated.jl")
+
 end # module
diff --git a/src/bases/dictionary/dict_evaluation.jl b/src/bases/dictionary/dict_evaluation.jl
index eb5e975..c930be1 100644
--- a/src/bases/dictionary/dict_evaluation.jl
+++ b/src/bases/dictionary/dict_evaluation.jl
@@ -217,8 +217,6 @@ default_eval_expansion(dict::Dictionary, coefficients, x) =
 default_unsafe_eval_expansion(dict::Dictionary, coefficients, x) =
     sum(coefficients[idx]*val for (idx,val) in unsafe_pointvalues(dict, x))
 
-@deprecate eval_expansion(dict::Dictionary, coefficients, grid::AbstractGrid; options...) eval_expansion.(Ref(dict), Ref(coefficients), grid; options...)
-
 function eval_expansion_derivative(dict::Dictionary, coefficients, x, order; options...)
     @assert size(coefficients) == size(dict)
     in_support(dict, x) ? unsafe_eval_expansion_derivative(dict, coefficients, x, order) : zero(span_codomaintype(dict, coefficients))
diff --git a/src/bases/dictionary/dictionary.jl b/src/bases/dictionary/dictionary.jl
index d94c37b..14b0440 100644
--- a/src/bases/dictionary/dictionary.jl
+++ b/src/bases/dictionary/dictionary.jl
@@ -1,10 +1,4 @@
 
-export domaintype,
-    codomaintype,
-    coefficienttype,
-    prectype
-
-
 """
 A `Dictionary{S,T}` is an ordered family of functions, in which each function
 maps a variable of type `S` to a variable of type `T`. The dictionary can be
@@ -101,6 +95,11 @@ lastindex(d::Dictionary) = last(eachindex(d))
 
 ≈(d1::Dictionary, d2::Dictionary) = d1==d2
 
+isequaldict(d1, d2) = isequaldict1(d1, d2)
+isequaldict1(d1, d2) = isequaldict2(d1, d2)
+isequaldict2(d1, d2) = default_isequaldict(d1, d2)
+default_isequaldict(d1, d2) = d1 == d2
+
 #############################
 # Domain and codomain type
 #############################
diff --git a/src/bases/dictionary/expansions.jl b/src/bases/dictionary/expansions.jl
index 18daf4d..e751fad 100644
--- a/src/bases/dictionary/expansions.jl
+++ b/src/bases/dictionary/expansions.jl
@@ -54,6 +54,7 @@ function expansion(dict::Dictionary, coefficients)
     T = promote_type(coefficienttype(dict), eltype(coef))
     Expansion(ensure_coefficienttype(T, dict), coef)
 end
+expansion(e::Expansion) = e
 
 similar(e::Expansion, coefficients) = Expansion(dictionary(e), coefficients)
 
@@ -64,7 +65,6 @@ coefficients(e::Expansion) = e.coefficients
 mapped_expansion(e::Expansion, m) =
     Expansion(mapped_dict(dictionary(e), m), coefficients(e))
 
-export random_expansion
 random_expansion(dict::Dictionary) = Expansion(dict, rand(dict))
 
 # For expansions of composite types, return a Expansion of a subdict
@@ -119,17 +119,6 @@ ei(dim,i, coefficients) = tuple((coefficients*Matrix{Int}(I, dim, dim)[:,i])...)
 differentiate(f::Expansion, var, order) = differentiate(f, ei(dimension(f), var, order))
 antidifferentiate(f::Expansion, var, order) = antidifferentiate(f, ei(dimension(f), var, order))
 
-function Base.:^(f::Expansion, i::Int)
-    @assert i >= 0
-    if i == 1
-        f
-    elseif i == 2
-        f * f
-    else
-        f^(i-1) * f
-    end
-end
-
 # To be implemented: Laplacian (needs multiplying functions)
 ## Δ(f::Expansion)
 
@@ -140,8 +129,6 @@ adjoint(f::Expansion) = differentiate(f)
 
 roots(f::Expansion) = expansion_roots(dictionary(f), coefficients(f))
 
-Base.@deprecate roots(dict::Dictionary, coef::AbstractVector) expansion_roots(dict, coef)
-
 show(io::IO, mime::MIME"text/plain", fun::Expansion) = composite_show(io, mime, fun)
 Display.displaystencil(fun::Expansion) = ["Expansion(", dictionary(fun), ", ", coefficients(fun), ")"]
 
@@ -199,56 +186,12 @@ function to_compatible_expansions(e1::Expansion, e2::Expansion)
     end
 end
 
-function (+)(f::Expansion, g::Expansion)
-    dict, coef = expansion_sum(dictionary(f),dictionary(g),coefficients(f),coefficients(g))
-    Expansion(dict, coef)
-end
-(-)(e1::Expansion, e2::Expansion) = e1 + similar(e2, -coefficients(e2))
-
-(-)(e::Expansion) = similar(e, -coefficients(e))
-
-@deprecate dict_multiply expansion_multiply
-expansion_multiply(dict1, dict2, coef1, coef2) =
-    expansion_multiply1(dict1, dict2, coef1, coef2)
-expansion_multiply1(dict1, dict2, coef1, coef2) =
-    expansion_multiply2(dict1, dict2, coef1, coef2)
-expansion_multiply2(dict1, dict2, coef1, coef2) =
-    error("Multiplication of expansions in these dictionaries is not implemented.")
-
-expansion_sum(dict1, dict2, coef1, coef2) =
-    expansion_sum1(dict1, dict2, coef1, coef2)
-expansion_sum1(dict1, dict2, coef1, coef2) =
-    expansion_sum2(dict1, dict2, coef1, coef2)
-function expansion_sum2(dict1, dict2, coef1, coef2)
-    if iscompatible(dict1, dict2)
-        default_expansion_sum(dict1, dict2, coef1, coef2)
-    else
-        error("Summation of expansions in these dictionaries is not implemented.")
-    end
-end
-
-function default_expansion_sum(dict1, dict2, coef1, coef2)
-    @assert iscompatible(dict1, dict2)
-    e1, e2 = to_compatible_expansions(Expansion(dict1, coef1), Expansion(dict2, coef2))
-    dictionary(e1), coefficients(e1)+coefficients(e2)
-end
-
-function (*)(f::Expansion, g::Expansion)
-    dict, coef = expansion_multiply(dictionary(f),dictionary(g),coefficients(f),coefficients(g))
-    Expansion(dict, coef)
-end
-
 (*)(op::DictionaryOperator, e::Expansion) = apply(op, e)
 
 Base.complex(e::Expansion) = Expansion(complex(dictionary(e)), complex(coefficients(e)))
-function Base.real(e::Expansion)
-    dict, coef = expansion_real(dictionary(e), coefficients(e))
-    Expansion(dict, coef)
-end
-function Base.imag(e::Expansion)
-    dict, coef = expansion_imag(dictionary(e), coefficients(e))
-    Expansion(dict, coef)
-end
+Base.real(e::Expansion) = Expansion(expansion_real(dictionary(e), coefficients(e))...)
+Base.imag(e::Expansion) = Expansion(expansion_imag(dictionary(e), coefficients(e))...)
+
 function expansion_real(dict::Dictionary, coefficients)
     @assert isreal(dict)
     real(dict), real(coefficients)
@@ -263,14 +206,6 @@ Base.one(e::Expansion) = similar(e, coefficients_of_one(dictionary(e)))
 expansion_of_one(dict::Dictionary) = Expansion(dict, coefficients_of_one(dict))
 expansion_of_x(dict::Dictionary) = Expansion(dict, coefficients_of_x(dict))
 
-for op in (:+, :-)
-    @eval Base.$op(a::Number, e::Expansion) = $op(a*one(e), e)
-    @eval Base.$op(e::Expansion, a::Number) = $op(e, a*one(e))
-end
-Base.:*(e::Expansion, a::Number) = expansion(dictionary(e), coefficients(e)*a)
-Base.:*(a::Number, e::Expansion) = expansion(dictionary(e), a*coefficients(e))
-Base.:/(e::Expansion, a::Number) = expansion(dictionary(e), coefficients(e)/a)
-
 apply(op::DictionaryOperator, e::Expansion) = Expansion(dest(op), op * coefficients(e))
 
 # TODO: we may not want to identify an Expansion with its array of coefficients
@@ -281,22 +216,3 @@ Base.collect(e::Expansion) = coefficients(e)
 
 Base.BroadcastStyle(e::Expansion) = Base.Broadcast.DefaultArrayStyle{dimension(e)}()
 
-# Interoperate with basis functions.
-# TODO: implement cleaner using promotion mechanism
-
-Base.:*(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) * expansion(φ2)
-Base.:*(f::Expansion, φ2::AbstractBasisFunction) = f * expansion(φ2)
-Base.:*(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) * f
-
-Base.:*(A::Number, φ::AbstractBasisFunction) = A * expansion(φ)
-Base.:*(φ::AbstractBasisFunction, A::Number) = A * expansion(φ)
-
-Base.:+(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) + expansion(φ2)
-Base.:+(f::Expansion, φ2::AbstractBasisFunction) = f + expansion(φ2)
-Base.:+(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) + f
-
-Base.:-(φ1::AbstractBasisFunction, φ2::AbstractBasisFunction) = expansion(φ1) - expansion(φ2)
-Base.:-(f::Expansion, φ2::AbstractBasisFunction) = f - expansion(φ2)
-Base.:-(φ1::AbstractBasisFunction, f::Expansion) = expansion(φ1) - f
-
-
diff --git a/src/bases/dictionary/gram.jl b/src/bases/dictionary/gram.jl
index efda1af..727f278 100644
--- a/src/bases/dictionary/gram.jl
+++ b/src/bases/dictionary/gram.jl
@@ -8,12 +8,6 @@ hasmeasure(Φ::Dictionary) = false
 # Shortcut: Dictionaries of the same type have just one measure
 defaultmeasure(Φ1::D, Φ2::D) where {D <: Dictionary} = measure(Φ1)
 
-@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) dict_innerproduct(Φ1, i, Φ2, j; options...)
-@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct_native(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct_native(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct1(Φ1, i, Φ2, j, measure; options...)
-@deprecate innerproduct2(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct2(Φ1, i, Φ2, j, measure; options...)
-
 dict_innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) =
     dict_innerproduct(Φ1, i, Φ2, j, defaultmeasure(Φ1, Φ2); options...)
 
diff --git a/src/bases/dictionary/promotion.jl b/src/bases/dictionary/promotion.jl
index 7886627..1aad01f 100644
--- a/src/bases/dictionary/promotion.jl
+++ b/src/bases/dictionary/promotion.jl
@@ -90,3 +90,20 @@ widen_prectype(::Type{Tuple{A,B}}, ::Type{U}) where {A,B,U} =
     Tuple{widen_prectype(A,U),widen_prectype(B,U)}
 widen_prectype(::Type{Tuple{A,B,C}}, ::Type{U}) where {A,B,C,U} =
     Tuple{widen_prectype(A,U),widen_prectype(B,U),widen_prectype(C,U)}
+
+
+## Promotion with conversion
+
+"""
+Promote two dictionaries to a common family, in such a way that expansions in both given
+dictionaries can be converted to expansions in the promoted dictionaries.
+"""
+promote_convertible(d1, d2) = promote_convertible1(promote_domaintype(d1, d2)...)
+promote_convertible1(d1, d2) = promote_convertible2(d1, d2)
+promote_convertible2(d1, d2) = default_promote_convertible(d1, d2)
+
+nopromotion(d1,d2) = error("Don't know how to promote and convert between dictionaries $(d1) and $(d2).")
+
+default_promote_convertible(d1::D, d2::D) where D = d1, d2
+default_promote_convertible(d1, d2) =
+    isequaldict(d1,d2) ? d1 : nopromotion(d1,d2)
diff --git a/src/bases/fourier/cosineseries.jl b/src/bases/fourier/cosineseries.jl
index 8c7b264..5dcea12 100644
--- a/src/bases/fourier/cosineseries.jl
+++ b/src/bases/fourier/cosineseries.jl
@@ -17,6 +17,7 @@ show(io::IO, b::CosineSeries{Float64}) = print(io, "CosineSeries($(length(b)))")
 show(io::IO, b::CosineSeries{T}) where T = print(io, "CosineSeries{$(T)}($(length(b)))")
 
 similar(b::CosineSeries, ::Type{T}, n::Int) where {T} = CosineSeries{T}(n)
+isequaldict(b1::CosineSeries, b2::CosineSeries) = length(b1)==length(b2)
 
 isbasis(b::CosineSeries) = true
 isorthogonal(b::CosineSeries, ::FourierWeight) = true
diff --git a/src/bases/fourier/fourier.jl b/src/bases/fourier/fourier.jl
index 7f1df5d..e206575 100644
--- a/src/bases/fourier/fourier.jl
+++ b/src/bases/fourier/fourier.jl
@@ -93,6 +93,7 @@ tofourier(d::FourierLike{T}) where {T} = convert(Fourier{T}, d)
 size(b::Fourier) = (b.n,)
 
 similar(b::Fourier, ::Type{T}, n::Int) where {T} = Fourier{T}(n)
+isequaldict(b1::Fourier, b2::Fourier) = length(b1)==length(b2)
 
 show(io::IO, b::Fourier{Float64}) = print(io, "Fourier($(length(b)))")
 show(io::IO, b::Fourier{T}) where T = print(io, "Fourier{$(T)}($(length(b)))")
diff --git a/src/bases/fourier/sinc.jl b/src/bases/fourier/sinc.jl
index 2d4914c..b765470 100644
--- a/src/bases/fourier/sinc.jl
+++ b/src/bases/fourier/sinc.jl
@@ -16,6 +16,7 @@ show(io::IO, b::PeriodicSincFunctions{T}) where T =
     print(io, "PeriodicSincFunctions{$(T)}($(length(b)))")
 
 similar(b::PeriodicSincFunctions, ::Type{T}, n::Int) where {T} = PeriodicSincFunctions{T}(n)
+isequaldict(b1::PeriodicSincFunctions, b2::PeriodicSincFunctions) = length(b1)==length(b2)
 
 size(b::PeriodicSincFunctions) = (b.n,)
 
diff --git a/src/bases/fourier/sineseries.jl b/src/bases/fourier/sineseries.jl
index d5ddf61..074cfeb 100644
--- a/src/bases/fourier/sineseries.jl
+++ b/src/bases/fourier/sineseries.jl
@@ -14,6 +14,7 @@ show(io::IO, b::SineSeries{T}) where T = print(io, "SineSeries{$(T)}($(length(b)
 SineSeries(n::Int) = SineSeries{Float64}(n)
 
 similar(b::SineSeries, ::Type{T}, n::Int) where {T} = SineSeries{T}(n)
+isequaldict(b1::SineSeries, b2::SineSeries) = length(b1)==length(b2)
 
 isbasis(b::SineSeries) = true
 isorthogonal(b::SineSeries, ::FourierWeight) = true
diff --git a/src/bases/fourier/trig.jl b/src/bases/fourier/trig.jl
index b51aeb9..771551f 100644
--- a/src/bases/fourier/trig.jl
+++ b/src/bases/fourier/trig.jl
@@ -21,6 +21,7 @@ coshalf(n::Int) = 1+nhalf(n)
 sinhalf(n::Int) = n-coshalf(n)
 
 similar(d::TrigSeries, ::Type{T}, n::Int) where {T} = TrigSeries{T}(n)
+isequaldict(b1::TrigSeries, b2::TrigSeries) = length(b1)==length(b2)
 
 show(io::IO, d::TrigSeries{Float64}) = print(io, "TrigSeries($(length(d)))")
 show(io::IO, d::TrigSeries{T}) where T = print(io, "TrigSeries{$(T)}($(length(d)))")
diff --git a/src/bases/modifiers/complexified_dict.jl b/src/bases/modifiers/complexified_dict.jl
index 4e7eeb2..eab0d75 100644
--- a/src/bases/modifiers/complexified_dict.jl
+++ b/src/bases/modifiers/complexified_dict.jl
@@ -22,6 +22,9 @@ Base.real(dict::ComplexifiedDict) = superdict(dict)
 
 similardictionary(s::ComplexifiedDict, s2::Dictionary) = ComplexifiedDict(s2)
 
+isequaldict1(d1::ComplexifiedDict, d2) = isequaldict(superdict(d1), d2)
+isequaldict2(d1, d2::ComplexifiedDict) = isequaldict(d1, superdict(d2))
+
 components(dict::ComplexifiedDict) = map(ComplexifiedDict, components(superdict(dict)))
 component(dict::ComplexifiedDict, i) = ComplexifiedDict(component(superdict(dict), i))
 
@@ -58,6 +61,25 @@ span_issubset2(d1, d2::ComplexifiedDict) = span_issubset(d1, superdict(d2))
 iscompatible(d1::Dictionary, d2::ComplexifiedDict) = iscompatible(d1, superdict(d2))
 iscompatible(d1::ComplexifiedDict, d2::Dictionary) = iscompatible(superdict(d1), d2)
 iscompatible(d1::ComplexifiedDict, d2::ComplexifiedDict) = iscompatible(superdict(d1), superdict(d2))
+
+function promote_convertible1(d1::ComplexifiedDict, d2::ComplexifiedDict)
+    reald1, reald2 = promote_convertible(superdict(d1), superdict(d2))
+    complex(reald1), complex(reald2)
+end
+function promote_convertible1(d1::ComplexifiedDict, d2)
+    reald1, reald2 = promote_convertible(superdict(d1), d2)
+    complex(reald1), complex(reald2)
+end
+function promote_convertible2(d1, d2::ComplexifiedDict)
+    reald1, reald2 = promote_convertible(d1, superdict(d2))
+    complex(reald1), complex(reald2)
+end
+
+function conversion1(T, d1::ComplexifiedDict, d2; options...)
+    @assert !isreal(T)
+    op = conversion(T, superdict(d1), d2; options...)
+    wrap_operator(d1, d2, op)
+end
 function conversion2(T, d1, d2::ComplexifiedDict; options...)
     @assert !isreal(T)
     op = conversion(T, d1, superdict(d2); options...)
diff --git a/src/bases/modifiers/mapped_dict.jl b/src/bases/modifiers/mapped_dict.jl
index bbb354f..38b8259 100644
--- a/src/bases/modifiers/mapped_dict.jl
+++ b/src/bases/modifiers/mapped_dict.jl
@@ -25,9 +25,6 @@ mapped_dict2(dict, map) = MappedDict(dict, map)
 mapped_dict1(dict::MappedDict, map) = mapped_dict(superdict(dict), map ∘ forward_map(dict))
 mapped_dict2(dict, map::IdentityMap) = dict
 
-# Convenience function, similar to apply_map for grids etcetera
-@deprecate apply_map(dict::Dictionary, map) mapped_dict(dict, map)
-
 components(dict::MappedDict) = map(d->mapped_dict(d, forward_map(dict)), components(superdict(dict)))
 
 # In the constructor we check the domain and codomain types.
diff --git a/src/bases/modifiers/operated_dict.jl b/src/bases/modifiers/operated_dict.jl
index 5d039ee..e51d546 100644
--- a/src/bases/modifiers/operated_dict.jl
+++ b/src/bases/modifiers/operated_dict.jl
@@ -283,6 +283,9 @@ function mixedgram2(T, dict1, dict2::OperatedDict, measure; options...)
 	wrap_operator(dict2, dict1, G *  operator(dict2))
 end
 
+promote_convertible1(d1::OperatedDict, d2) = promote_convertible(superdict(d1), d2)
+promote_convertible2(d1, d2::OperatedDict) = promote_convertible(d1, superdict(d2))
+
 function conversion1(::Type{T}, src::OperatedDict, dest; options...) where {T}
 	op = conversion(T, superdict(src), dest; options...) * operator(src)
 	wrap_operator(src, dest, op)
diff --git a/src/bases/modifiers/subdicts.jl b/src/bases/modifiers/subdicts.jl
index b58b64e..bc0d252 100644
--- a/src/bases/modifiers/subdicts.jl
+++ b/src/bases/modifiers/subdicts.jl
@@ -99,6 +99,20 @@ dict_innerproduct1(d1::Subdictionary, i, d2, j, measure; options...) =
 dict_innerproduct2(d1, i, d2::Subdictionary, j, measure; options...) =
     dict_innerproduct(d1, i, superdict(d2), superindices(d2, j), measure; options...)
 
+promote_convertible1(d1::Subdictionary, d2) =
+    d1 == d2 ? (d1, d2) : promote_convertible(superdict(d1), d2)
+promote_convertible2(d1, d2::Subdictionary) =
+    d1 == d2 ? (d1, d2) : promote_convertible(d1, superdict(d2))
+
+function conversion1(::Type{T}, d1::Subdictionary, d2) where T
+    if d1 == d2
+        IdentityOperator{T}(d1, d1)
+    else
+        E = extension(T, d1, superdict(d1))
+        op = conversion(superdict(d1), d2)
+        op * E
+    end
+end
 
 
 #####################
diff --git a/src/bases/poly/chebyshev/ChebyshevT.jl b/src/bases/poly/chebyshev/ChebyshevT.jl
index d73c241..1f43e21 100644
--- a/src/bases/poly/chebyshev/ChebyshevT.jl
+++ b/src/bases/poly/chebyshev/ChebyshevT.jl
@@ -15,8 +15,9 @@ ChebyshevT(d::PolynomialDegree) = ChebyshevT(value(d)+1)
 ChebyshevT{T}(d::PolynomialDegree) where T = ChebyshevT{T}(value(d)+1)
 
 similar(b::ChebyshevT, ::Type{T}, n::Int) where {T} = ChebyshevT{T}(n)
+isequaldict(b1::ChebyshevT, b2::ChebyshevT) = length(b1)==length(b2)
 
-to_chebyshev_dict(dict::IntervalOPS{T}) where T = ChebyshevT{T}(length(dict))
+to_chebyshev_dict(dict::Dictionary{T,T}) where T = ChebyshevT{T}(length(dict))
 to_chebyshev(f) = to_chebyshev(expansion(f))
 to_chebyshev(f::Expansion) = to_chebyshev(dictionary(f), coefficients(f))
 to_chebyshev(dict::Dictionary, coef) =
diff --git a/src/bases/poly/chebyshev/ChebyshevU.jl b/src/bases/poly/chebyshev/ChebyshevU.jl
index 6e9ed54..da05189 100644
--- a/src/bases/poly/chebyshev/ChebyshevU.jl
+++ b/src/bases/poly/chebyshev/ChebyshevU.jl
@@ -10,6 +10,7 @@ ChebyshevU(d::PolynomialDegree) = ChebyshevU(value(d)+1)
 ChebyshevU{T}(d::PolynomialDegree) where T = ChebyshevU{T}(value(d)+1)
 
 similar(b::ChebyshevU, ::Type{T}, n::Int) where {T} = ChebyshevU{T}(n)
+isequaldict(b1::ChebyshevU, b2::ChebyshevU) = length(b1)==length(b2)
 
 show(io::IO, d::ChebyshevU{Float64}) = print(io, "ChebyshevU($(length(d)))")
 show(io::IO, d::ChebyshevU{T}) where T = print(io, "ChebyshevU{$(T)}($(length(d)))")
diff --git a/src/bases/poly/monomials.jl b/src/bases/poly/monomials.jl
index c7eb1ca..49e347b 100644
--- a/src/bases/poly/monomials.jl
+++ b/src/bases/poly/monomials.jl
@@ -65,6 +65,7 @@ function differentiation(::Type{T}, src::Monomials, dest::Monomials, order::Int;
 end
 
 similar(b::Monomials, ::Type{T}, n::Int) where {T} = Monomials{T}(n)
+isequaldict(b1::Monomials, b2::Monomials) = length(b1)==length(b2)
 
 complex(d::Monomials{T}) where {T} = similar(d, complex(T))
 
diff --git a/src/bases/poly/ops/hermite.jl b/src/bases/poly/ops/hermite.jl
index b76d2e9..23498ee 100644
--- a/src/bases/poly/ops/hermite.jl
+++ b/src/bases/poly/ops/hermite.jl
@@ -14,6 +14,7 @@ Hermite(d::PolynomialDegree) = Hermite(value(d)+1)
 Hermite{T}(d::PolynomialDegree) where T = Hermite{T}(value(d)+1)
 
 similar(b::Hermite, ::Type{T}, n::Int) where {T} = Hermite{T}(n)
+isequaldict(b1::Hermite, b2::Hermite) = length(b1)==length(b2)
 
 support(b::Hermite{T}) where {T} = DomainSets.FullSpace(T)
 
diff --git a/src/bases/poly/ops/jacobi.jl b/src/bases/poly/ops/jacobi.jl
index 1d0ef70..f684879 100644
--- a/src/bases/poly/ops/jacobi.jl
+++ b/src/bases/poly/ops/jacobi.jl
@@ -28,6 +28,9 @@ const JacobiExpansion{T,C} = Expansion{T,T,Jacobi{T},C}
 
 similar(b::Jacobi, ::Type{T}, n::Int) where {T} = Jacobi{T}(n, b.α, b.β)
 
+iscompatible(d1::Jacobi, d2::Jacobi) = d1.α == d2.α && d1.β == d2.β
+isequaldict(b1::Jacobi, b2::Jacobi) = length(b1)==length(b2) && iscompatible(b1,b2)
+
 first_moment(b::Jacobi{T}) where {T} = (b.α+b.β+1≈0) ?
     T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1) :
     T(2).^(b.α+b.β+1)*gamma(b.α+1)*gamma(b.β+1)/(b.α+b.β+1)/gamma(b.α+b.β+1)
@@ -38,7 +41,6 @@ jacobi_β(b::Jacobi) = b.β
 
 measure(b::Jacobi) = JacobiWeight(b.α, b.β)
 
-iscompatible(d1::Jacobi, d2::Jacobi) = d1.α == d2.α && d1.β == d2.β
 
 isorthogonal(dict::Jacobi, measure::JacobiWeight) =
 	dict.α == measure.α && dict.β == measure.β
@@ -88,21 +90,21 @@ function dict_innerproduct_native(d1::Jacobi, i::PolynomialDegree, d2::Jacobi, j
 	end
 end
 
-function expansion_sum(src1::Jacobi, src2::Jacobi, coef1, coef2)
-	if iscompatible(src1, src2)
-		default_expansion_sum(src1, src2, coef1, coef2)
-	else
-		to_chebyshev_expansion_sum(src1, src2, coef1, coef2)
-	end
-end
-
-function expansion_multiply(src1::Jacobi, src2::Jacobi, coef1, coef2)
-	if iscompatible(src1, src2)
-		default_poly_expansion_multiply(src1, src2, coef1, coef2)
-	else
-		to_chebyshev_expansion_multiply(src1, src2, coef1, coef2)
-	end
-end
+# function expansion_sum(src1::Jacobi, src2::Jacobi, coef1, coef2)
+# 	if iscompatible(src1, src2)
+# 		default_expansion_sum(src1, src2, coef1, coef2)
+# 	else
+# 		to_chebyshev_expansion_sum(src1, src2, coef1, coef2)
+# 	end
+# end
+
+# function expansion_multiply(src1::Jacobi, src2::Jacobi, coef1, coef2)
+# 	if iscompatible(src1, src2)
+# 		default_poly_expansion_multiply(src1, src2, coef1, coef2)
+# 	else
+# 		to_chebyshev_expansion_multiply(src1, src2, coef1, coef2)
+# 	end
+# end
 
 ## Printing
 function show(io::IO, b::Jacobi{Float64})
diff --git a/src/bases/poly/ops/laguerre.jl b/src/bases/poly/ops/laguerre.jl
index 29b2d0b..01421ce 100644
--- a/src/bases/poly/ops/laguerre.jl
+++ b/src/bases/poly/ops/laguerre.jl
@@ -21,6 +21,8 @@ Laguerre{T}(d::PolynomialDegree, args...; options...) where T =
 	Laguerre{T}(value(d)+1, args...; options...)
 
 similar(b::Laguerre, ::Type{T}, n::Int) where {T} = Laguerre{T}(n, b.α)
+isequaldict(b1::Laguerre, b2::Laguerre) = length(b1)==length(b2) &&
+	b1.α == b2.α
 
 support(b::Laguerre{T}) where {T} = HalfLine{T}()
 
diff --git a/src/bases/poly/ops/legendre.jl b/src/bases/poly/ops/legendre.jl
index 35252e9..1c202b9 100644
--- a/src/bases/poly/ops/legendre.jl
+++ b/src/bases/poly/ops/legendre.jl
@@ -13,8 +13,9 @@ Legendre(d::PolynomialDegree) = Legendre(value(d)+1)
 Legendre{T}(d::PolynomialDegree) where T = Legendre{T}(value(d)+1)
 
 similar(b::Legendre, ::Type{T}, n::Int) where {T} = Legendre{T}(n)
+isequaldict(b1::Legendre, b2::Legendre) = length(b1)==length(b2)
 
-to_legendre_dict(dict::IntervalOPS{T}) where T = Legendre{T}(length(dict))
+to_legendre_dict(dict::Dictionary{T,T}) where T = Legendre{T}(length(dict))
 to_legendre(f) = to_legendre(expansion(f))
 to_legendre(f::Expansion) = to_legendre(dictionary(f), coefficients(f))
 to_legendre(dict::Dictionary, coef) =
diff --git a/src/bases/poly/ops/orthopoly.jl b/src/bases/poly/ops/orthopoly.jl
index eed258f..5512240 100644
--- a/src/bases/poly/ops/orthopoly.jl
+++ b/src/bases/poly/ops/orthopoly.jl
@@ -73,6 +73,10 @@ function unsafe_eval_element_derivative(b::OPS, idx::PolynomialDegree, x, order)
 end
 
 ## Some routines for OPS on [-1,1] are implemented for Chebyshev polynomials
+
+promote_convertible(d1::IntervalOPS, d2::IntervalOPS) =
+    iscompatible(d1,d2) ? (d1, d2) : (to_chebyshev_dict(d1), to_chebyshev_dict(d2))
+
 expansion_roots(dict::IntervalOPS, coef) = roots(to_chebyshev(dict, coef))
 
 expansion_multiply(src1::IntervalOPS, src2::IntervalOPS, coef1, coef2) =
@@ -87,15 +91,6 @@ function expansion_multiply(src1::I, src2::I, coef1, coef2) where {I <: Interval
     dictionary(result2), coefficients(result2)
 end
 
-expansion_sum(src1::IntervalOPS, src2::IntervalOPS, coef1, coef2) =
-    to_chebyshev_expansion_sum(src1, src2, coef1, coef2)
-function to_chebyshev_expansion_sum(src1, src2, coef1, coef2)
-    result = to_chebyshev(src1, coef1) + to_chebyshev(src2, coef2)
-    dictionary(result), coefficients(result)
-end
-expansion_sum(src1::I, src2::I, coef1, coef2) where {I <: IntervalOPS} =
-    default_expansion_sum(src1, src2, coef1, coef2)
-
 
 hasmeasure(dict::OPS) = true
 
@@ -134,8 +129,6 @@ function symmetric_jacobi_matrix(b::OPS)
     SymTridiagonal(J)
 end
 
-Base.@deprecate roots(b::OPS) ops_roots(b)
-
 "Return the roots of the (N+1)st orthogonal polynomial."
 function ops_roots(b::OPS{T}) where {T<:Number}
     J = symmetric_jacobi_matrix(b)
diff --git a/src/computations/conversion.jl b/src/computations/conversion.jl
index 5bb1e89..57c19ee 100644
--- a/src/computations/conversion.jl
+++ b/src/computations/conversion.jl
@@ -12,15 +12,13 @@ export conversion,
 
 """
 ```
-conversion(::Type{T}, src::Dictionary, dest::Dictionary; options...)
+conversion([::Type{T}, ]src::Dictionary, dest::Dictionary; options...)
 ```
 
-Construct an operator with element type `T` that converts coefficients from the
-source dictionary to coefficients of the destination dictionary.
+Construct an operator that converts coefficients from the source dictionary to
+coefficients of the destination dictionary.
 
-The conversion is often exact, but this can not always be achieved with
-certainty. For example, it is hard to know in advance whether converting to
-evaluations in a grid is invertible.
+The conversion is exact, up to side-effects of finite-precision calculations.
 """
 conversion
 
@@ -158,20 +156,25 @@ end
 change_basis(F1::Expansion; dest::Dictionary) =
     conversion(dictionary(F1), dest) * F1
 
-conversion(::Type{T}, src, dest; options...) where {T} =
+noconversion(src, dest) = error("No known exact conversion from $(src) to $(dest).")
+
+conversion(::Type{T}, src, dest; options...) where T =
     conversion1(T, src, dest; options...)
-conversion1(T, src::Dictionary, dest; options...) =
-    conversion2(T, src, dest; options...)
-function conversion2(T, src, dest::Dictionary; verbose=false, options...)
-    if issubset(Span(src), Span(dest))
+conversion1(T, src, dest; options...) = conversion2(T, src, dest; options...)
+conversion2(T, src, dest; options...) = default_conversion(T, src, dest; options...)
+
+function default_conversion(::Type{T}, src, dest; verbose=false, options...) where T
+    if isequaldict(src, dest)
+        IdentityOperator{T}(src, dest)
+    elseif issubset(Span(src), Span(dest))
         verbose && println("WARN: using default conversion from $(src) to $(dest)")
-        default_conversion(T, src, dest; verbose, options...)
+        explicit_conversion(T, src, dest; verbose, options...)
     else
-        error("No known exact conversion from $(src) to $(dest)")
+        noconversion(src, dest)
     end
 end
 
-function default_conversion(T, src, dest; options...)
+function explicit_conversion(::Type{T}, src, dest; options...) where T
     if hasmeasure(dest)
         projection(T, src, dest, measure(dest); options...)
     elseif hasinterpolationgrid(dest)
@@ -185,7 +188,7 @@ function default_conversion(T, src, dest; options...)
         A = evaluation(T, src, pts; options...)
         I*A
     else
-        error("Don't know how to reliably convert from $(src) to $(dest).")
+        noconversion(src, dest)
     end
 end
 
@@ -193,6 +196,7 @@ end
 # To grid: we invoke `evaluation`
 # From grid: we invoke `interpolation`.
 # Problem is that we don't know in general whether these would be exact.
+# TODO: these should be removed (breaking change)
 conversion(::Type{T}, src::Dictionary, dest::GridBasis; options...) where {T} =
     evaluation(T, src, dest; options...)
 
diff --git a/src/computations/evaluation.jl b/src/computations/evaluation.jl
index ae3e308..334bc6a 100644
--- a/src/computations/evaluation.jl
+++ b/src/computations/evaluation.jl
@@ -19,8 +19,6 @@ function evaluation_matrix!(A::AbstractMatrix, Φ::Dictionary, pts)
     A
 end
 
-@deprecate dense_evaluation default_evaluation
-
 function default_evaluation(::Type{T}, Φ::Dictionary, gb::GridBasis; options...) where {T}
     A = evaluation_matrix(Φ, grid(gb), T)
     ArrayOperator(A, Φ, gb)
diff --git a/src/computations/fun.jl b/src/computations/fun.jl
new file mode 100644
index 0000000..776e00c
--- /dev/null
+++ b/src/computations/fun.jl
@@ -0,0 +1,85 @@
+# Code to compute with functions
+
+const FUN = Union{<:TypedFunction,<:Expansion}
+
+"""
+Promote the given two functions such that they have compatible type, suitable for further
+operations such as addition and multiplication.
+
+In most cases the results are two expansions with a dictionary of the same family.
+
+See also: [`funpromote_samelength`](@ref)
+"""
+funpromote(f, g) = funpromote(expansion(f), expansion(g))
+
+funpromote(f::Expansion, g::Expansion) =
+    funpromote(dictionary(f), dictionary(g), coefficients(f), coefficients(g))
+
+funpromote(d1, d2, c1, c2) = _funpromote(d1, d2, c1, c2, promote_convertible(d1, d2)...)
+_funpromote(d1::D1, d2::D2, c1, c2, d3::D1, d4::D2) where {D1,D2} =
+    expansion(d1, c1), expansion(d2, c2)    # no conversion happened
+function _funpromote(d1, d2, c1, c2, d3, d4)
+    C1 = conversion(d1, d3)
+    C2 = conversion(d2, d4)
+    C1 * expansion(d1, c1), C2 * expansion(d2, c2)
+end
+
+
+Base.:+(f::FUN, g::FUN) = expansion_sum(funpromote(f, g)...)
+Base.:-(f::FUN) = -expansion(f)
+Base.:-(f::Expansion) = expansion(dictionary(f), -coefficients(f))
+Base.:-(f::FUN, g::FUN) = f + (-g)
+
+expansion_sum(f::Expansion, g::Expansion) = expansion_sum(dictionary(f), dictionary(g), coefficients(f), coefficients(g))
+
+function expansion_sum(d1, d2, c1, c2)
+    # we assume that d1 and d2 are compatible, just looking for a possible size mismatch
+    if size(d1) == size(d2)
+        expansion(d1, c1+c2)
+    else
+        if length(d1) < length(d2)
+            e12 = extension(d1, d2)
+            expansion(d2, e12*c1 + c2)
+        elseif length(d1) > length(d2)
+            e21 = extension(d2, d1)
+            expansion(d1, c1 + e21*c2)
+        else
+            error("Sizes of dictionaries don't match, but lengths do.")
+        end
+    end
+end
+
+
+Base.:*(f::FUN, g::FUN) = expansion_multiply(funpromote(f, g)...)
+
+function expansion_multiply(f::Expansion, g::Expansion)
+    dict, coef = expansion_multiply(dictionary(f), dictionary(g), coefficients(f), coefficients(g))
+    expansion(dict, coef)
+end
+
+# Multiplication by scalars
+
+Base.:*(f::FUN, a::Number) = a*f
+Base.:*(a::Number, f::FUN) = a*expansion(f)
+Base.:*(a::Number, f::Expansion) = expansion(dictionary(f), a*coefficients(f))
+
+Base.:/(f::FUN, a::Number) = expansion(f) / a
+Base.:/(f::Expansion, a::Number) = expansion(dictionary(f), coefficients(f)/a)
+
+for op in (:+, :-)
+    @eval Base.$op(a::Number, f::FUN) = $op(a*one(expansion(f)), f)
+    @eval Base.$op(f::FUN, a::Number) = $op(f, a*one(expansion(f)))
+end
+
+# exponentiation
+
+function Base.:^(f::FUN, i::Int)
+    @assert i >= 0
+    if i == 1
+        f
+    elseif i == 2
+        f * f
+    else
+        f^(i-1) * f
+    end
+end
diff --git a/src/computations/generic_operators.jl b/src/computations/generic_operators.jl
index 6dbcab6..acdbee2 100644
--- a/src/computations/generic_operators.jl
+++ b/src/computations/generic_operators.jl
@@ -9,14 +9,6 @@ function (*)(op::FunOperator, args...)
     Expansion(dest(op.op), coef)
 end
 
-# In this file we define the interface for a number of generic functions:
-# See the individual files for details on the interfaces.
-
-include("transform.jl")
-include("evaluation.jl")
-include("approximation.jl")
-include("differentiation.jl")
-
 
 
 #####################################
diff --git a/src/deprecated.jl b/src/deprecated.jl
new file mode 100644
index 0000000..8ca2568
--- /dev/null
+++ b/src/deprecated.jl
@@ -0,0 +1,18 @@
+
+@deprecate sparse_matrix sparse
+
+@deprecate eval_expansion(dict::Dictionary, coefficients, grid::AbstractGrid; options...) eval_expansion.(Ref(dict), Ref(coefficients), grid; options...)
+
+@deprecate roots(dict::Dictionary, coef::AbstractVector) expansion_roots(dict, coef)
+
+@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j; options...) dict_innerproduct(Φ1, i, Φ2, j; options...)
+@deprecate innerproduct(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct(Φ1, i, Φ2, j, measure; options...)
+@deprecate innerproduct_native(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct_native(Φ1, i, Φ2, j, measure; options...)
+@deprecate innerproduct1(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct1(Φ1, i, Φ2, j, measure; options...)
+@deprecate innerproduct2(Φ1::Dictionary, i, Φ2::Dictionary, j, measure; options...) dict_innerproduct2(Φ1, i, Φ2, j, measure; options...)
+
+@deprecate apply_map(dict::Dictionary, map) mapped_dict(dict, map)
+
+@deprecate roots(b::OPS) ops_roots(b)
+
+@deprecate dense_evaluation default_evaluation