diff --git a/NEWS.md b/NEWS.md
index 28fff37..be89080 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,12 +2,21 @@
 
 ## Wish list for future developments
 
-* Simplifications that are automatically done by`LazyAlgebra` may change
-  multipliers but must not change the coefficients of the mappings.  Call
-  `simplify(A)` to apply further simplifications that may change the
-  coefficients of the mappings in `A`.  For instance, assuming `a` is an array,
-  `inv(Diag(a))` automatically yields `Inverse(Diag(a))` while
-  `simplify(inv(Diag(a)))` yields `Diag(1 ./ a)`.
+* Traits must be type-stable. That is deciding the value of a specific trait
+  should be done on the sole basis of the the type of the argument(s).
+
+* Optimizations and simplifications of expressions involving mappings belong to
+  two categories:
+  * Automatic simplifications performed by`LazyAlgebra`. These can change the
+    values of multipliers but must not change the coefficients of the mappings.
+    Sadly, it is not possible to have every possible automatic simplifications
+    be type-stable.
+  * The user may explicitly call `simplify(A)` to apply further simplifications
+    that may change the coefficients of the mappings in `A`. For example,
+    assuming `a` is an array, `inv(Diag(a))` automatically yields
+    `Inverse(Diag(a))` while `simplify(inv(Diag(a)))` yields `Diag(1 ./ a)`.
+    These simplifications may not be type-stable. For example the multiplier of
+    a scaled mapping becoming equal to one can be eliminated.
 
 * Calling BLAS should be avoided in some cases, either because BLAS is slower
   than optimized Julia code, or because BLAS may use more than one thread in
@@ -17,8 +26,19 @@
   arguments.  This would be useful to deal with arrays whose elements have
   non-standard numerical types as physical quantities in the `Unitful` package.
 
+## Branch 0.3
+
+### Breaking changes
+
+* Abstract type `Mapping{L}` has a built-in parameter indicating whether the
+  mapping is linbear of not. This breaks compatibility but simplifies a lot
+  many parts of the code and makes the linear trait decidable at compile time.
+
+### Other changes
+
+* Methods `multipler`, `unscaled` can be applied to a mapping type to yield the
+  corresponding type.
 
-## Version 0.2.3
 ## Branch 0.2
 
 ### Version 0.2.5
diff --git a/src/LazyAlgebra.jl b/src/LazyAlgebra.jl
index 0856218..bae8875 100644
--- a/src/LazyAlgebra.jl
+++ b/src/LazyAlgebra.jl
@@ -27,6 +27,7 @@ export
     Jacobian,
     LinearMapping,
     Mapping,
+    NonLinearMapping,
     NonuniformScaling,
     RankOneOperator,
     SingularSystem,
@@ -111,7 +112,8 @@ import Base: Tuple, adjoint, inv, axes,
 
 # Import/using from LinearAlgebra, BLAS and SparseArrays.
 using LinearAlgebra
-import LinearAlgebra: UniformScaling, diag, ⋅, mul!, rmul!
+using LinearAlgebra: UniformScaling
+import LinearAlgebra: diag, ⋅, mul!, rmul!
 using LinearAlgebra.BLAS
 using LinearAlgebra.BLAS: libblas, @blasfunc,
     BlasInt, BlasReal, BlasFloat, BlasComplex
diff --git a/src/fft.jl b/src/fft.jl
index 57e7b5d..9e62c7e 100644
--- a/src/fft.jl
+++ b/src/fft.jl
@@ -426,7 +426,7 @@ MorphismType(::FFTOperator{<:Complex}) = Endomorphism()
 
 ncols(A::FFTOperator) = A.ncols
 ncols(A::Adjoint{<:FFTOperator}) = ncols(unveil(A))
-ncols(A::Inverse{<:FFTOperator}) = ncols(unveil(A))
+ncols(A::Inverse{true,<:FFTOperator}) = ncols(unveil(A))
 ncols(A::InverseAdjoint{<:FFTOperator}) = ncols(unveil(A))
 
 input_size(A::FFTOperator) = A.inpdims # FIXME: input_size(A.forward)
@@ -456,9 +456,9 @@ show(io::IO, A::FFTOperator) = print(io, "FFT")
     (identical(unveil(A), B) ? ncols(A)*Id : compose(A, B))
 *(A::F, B::Adjoint{F}) where {F<:FFTOperator} =
     (identical(A, unveil(B)) ? ncols(A)*Id : compose(A, B))
-*(A::InverseAdjoint{F}, B::Inverse{F}) where {F<:FFTOperator} =
+*(A::InverseAdjoint{F}, B::Inverse{true,F}) where {F<:FFTOperator} =
     (identical(unveil(A), unveil(B)) ? (1//ncols(A))*Id : compose(A, B))
-*(A::Inverse{F}, B::InverseAdjoint{F}) where {F<:FFTOperator} =
+*(A::Inverse{true,F}, B::InverseAdjoint{F}) where {F<:FFTOperator} =
     (identical(unveil(A), unveil(B)) ? (1//ncols(A))*Id : compose(A, B))
 
 function vcreate(P::Type{<:Union{Direct,InverseAdjoint}},
diff --git a/src/mappings.jl b/src/mappings.jl
index e7012b0..1b1e2f2 100644
--- a/src/mappings.jl
+++ b/src/mappings.jl
@@ -46,18 +46,19 @@ end
 #------------------------------------------------------------------------------
 # SYMBOLIC MAPPINGS (FOR TESTS)
 
-struct SymbolicMapping{T} <: Mapping end
-struct SymbolicLinearMapping{T} <: LinearMapping end
-SymbolicMapping(id::AbstractString) = SymbolicMapping(Symbol(id))
-SymbolicMapping(id::Symbol) = SymbolicMapping{Val{id}}()
-SymbolicLinearMapping(id::AbstractString) = SymbolicLinearMapping(Symbol(id))
-SymbolicLinearMapping(x::Symbol) = SymbolicLinearMapping{Val{x}}()
+struct SymbolicMapping{L,S} <: Mapping{L} end
+const SymbolicLinearMapping{S} = SymbolicMapping{true,S}
+SymbolicMapping(id) = SymbolicMapping{false}(id)
+SymbolicMapping{L}(id::AbstractString) where {L} = SymbolicMapping{L}(Symbol(id))
+SymbolicMapping{L}(id::Symbol) where {L} = SymbolicMapping{L,Val{id}}()
+# FIXME: SymbolicLinearMapping(id::AbstractString) = SymbolicLinearMapping(Symbol(id))
+# FIXME: SymbolicLinearMapping(id::Symbol) = SymbolicLinearMapping{Val{id}}()
 
-show(io::IO, A::SymbolicMapping{Val{T}}) where {T} = print(io, T)
-show(io::IO, A::SymbolicLinearMapping{Val{T}}) where {T} = print(io, T)
+show(io::IO, A::SymbolicMapping{L,Val{S}}) where {L,S} = print(io, S)
+# FIXME: show(io::IO, A::SymbolicLinearMapping{L,Val{S}}) where {L,S} = print(io, S)
 
 identical(::T, ::T) where {T<:SymbolicMapping} = true
-identical(::T, ::T) where {T<:SymbolicLinearMapping} = true
+# FIXME: identical(::T, ::T) where {T<:SymbolicLinearMapping} = true
 
 #------------------------------------------------------------------------------
 # NON-UNIFORM SCALING
diff --git a/src/methods.jl b/src/methods.jl
index 7f8aec8..d8de5bc 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -56,7 +56,7 @@ function show(io::IO, A::Scaled)
     show(io, M)
 end
 
-function show(io::IO, A::Scaled{<:Sum})
+function show(io::IO, A::Scaled{<:Any,<:Sum})
     λ, M = multiplier(A), unscaled(A)
     if λ == -1
         print(io, "-(")
@@ -69,36 +69,36 @@ function show(io::IO, A::Scaled{<:Sum})
     end
 end
 
-function show(io::IO, A::Adjoint{<:Mapping})
+function show(io::IO, A::Adjoint)
     show(io, unveil(A))
     print(io, "'")
 end
 
-function show(io::IO, A::Adjoint{T}) where {T<:Union{Scaled,Composition,Sum}}
+function show(io::IO, A::Adjoint{<:Union{Scaled,Composition,Sum}})
     print(io, "(")
     show(io, unveil(A))
     print(io, ")'")
 end
 
-function show(io::IO, A::Inverse{<:Mapping})
+function show(io::IO, A::Inverse)
     print(io, "inv(")
     show(io, unveil(A))
     print(io, ")")
 end
 
-function show(io::IO, A::InverseAdjoint{<:Mapping})
+function show(io::IO, A::InverseAdjoint)
     print(io, "inv(")
     show(io, unveil(A))
     print(io, ")'")
 end
 
-function show(io::IO, A::Jacobian{<:Mapping})
+function show(io::IO, A::Jacobian)
     print(io, "∇(")
     show(io, primitive(A))
     print(io, ",x)")
 end
 
-function show(io::IO, A::Sum{N}) where {N}
+function show(io::IO, A::Sum)
     function show_term(io::IO, A::Sum)
         print(io, "(")
         show(io, A)
@@ -106,7 +106,7 @@ function show(io::IO, A::Sum{N}) where {N}
     end
     show_term(io::IO, A::Mapping) = show(io, A)
 
-    for i in 1:N
+    for i in eachindex(A)
         let B = A[i]
             if isa(B, Scaled)
                 λ, M = multiplier(B), unscaled(B)
@@ -130,8 +130,8 @@ function show(io::IO, A::Sum{N}) where {N}
     end
 end
 
-function show(io::IO, A::Composition{N}) where {N}
-    for i in 1:N
+function show(io::IO, A::Composition)
+    for i in eachindex(A)
         let B = A[i]
             if i > 1
                 print(io, "⋅")
@@ -158,9 +158,21 @@ If `A` is sum or a composition of mappings, `Tuple(A)` yields the same result
 as `terms(A)`.
 
 """
-terms(A::Union{Sum,Composition}) = getfield(A, :ops)
+terms(A::Union{Sum,Composition}) = getfield(A, :terms)
 terms(A::Mapping) = (A,)
 
+"""
+    nterms(A)
+
+yields the number of terms in a sum or composition of mappings `A`. The
+returned value is a *trait*, argument can be a mapping instance or type. If `A`
+is a mapping, `nterms(A)` yields the same result as `length(A)`.
+
+"""
+nterms(x::Union{Sum,Composition}) = nterms(typeof(x))
+nterms(::Type{<:Sum{L,N}}) where {L,N} = N
+nterms(::Type{<:Composition{L,N}}) where {L,N} = N
+
 """
     unveil(A)
 
@@ -172,10 +184,8 @@ As a special case, `A` may be an instance of `LinearAlgebra.UniformScaling` and
 the result is the LazyAlgebra mapping corresponding to `A`.
 
 """
-unveil(A::DecoratedMapping) = getfield(A, :op)
+unveil(A::DecoratedMapping) = getfield(A, :parent)
 unveil(A::Mapping) = A
-unveil(A::UniformScaling) = multiplier(A)*Id
-Mapping(A::UniformScaling) = unveil(A)
 
 """
     unscaled(A)
@@ -185,10 +195,14 @@ otherwise yields `A`. This method also works for intances of
 `LinearAlgebra.UniformScaling`. Call [`multiplier`](@ref) to get the multiplier
 `λ`.
 
+If `A` is a mapping type, yields the corresponding unscaled mapping type.
+
 """
 unscaled(A::Mapping) = A
-unscaled(A::Scaled) = getfield(A, :M)
-unscaled(A::UniformScaling) = Id
+unscaled(A::Scaled) = getfield(A, :mapping)
+
+unscaled(::Type{M}) where {M<:Mapping} = M
+unscaled(::Type{<:Scaled{L,M,S}}) where {L,M,S} = M
 
 """
     multiplier(A)
@@ -198,10 +212,23 @@ yields the multiplier `λ` of the scaled mapping `A = λ*M` (see
 intances of `LinearAlgebra.UniformScaling`. Call [`unscaled`](@ref) to get the
 mapping `M`. `λ`.
 
+If `A` is a mapping type, yields the corresponding multiplier type.
+
 """
-multiplier(A::Scaled) = getfield(A, :λ)
+multiplier(A::Scaled) = getfield(A, :multiplier)
 multiplier(A::Mapping) = 1
+
+multiplier(::Type{M}) where {M<:Mapping} = Int
+multiplier(::Type{<:Scaled{L,M,S}}) where {L,M,S} = S
+
+# Extend a few methods for `LinearAlgebra.UniformScaling` (NOTE: see code in
+# `LinearAlgebra/src/uniformscaling.jl`)
+unveil(A::UniformScaling) = multiplier(A)*Id
+Mapping(A::UniformScaling) = unveil(A)
+unscaled(A::UniformScaling) = Id
+unscaled(::Type{<:UniformScaling}) = typeof(Id)
 multiplier(A::UniformScaling) = getfield(A, :λ)
+multiplier(::Type{<:UniformScaling{T}}) where {T} = T
 
 """
     primitive(J)
@@ -210,7 +237,7 @@ yields the mapping `A` embedded in the Jacobian `J = ∇(A,x)`. Call
 [`variables`](@ref) to get `x` instead.
 
 """
-primitive(J::Jacobian) = getfield(J, :A)
+primitive(J::Jacobian) = getfield(J, :primitive)
 
 """
     variables(J)
@@ -219,7 +246,7 @@ yields the variables `x` embedded in the Jacobian `J = ∇(A,x)`. Call
 [`primitive`](@ref) to get `A` instead.
 
 """
-variables(J::Jacobian) = getfield(J, :x)
+variables(J::Jacobian) = getfield(J, :variables)
 
 """
     identifier(A)
@@ -237,19 +264,19 @@ identifier(A::Mapping) = objectid(unscaled(A))
 Base.isless(A::Mapping, B::Mapping) = isless(identifier(A), identifier(B))
 
 # Extend base methods to simplify the code for reducing expressions.
-first(A::Mapping) = A
-last(A::Mapping) = A
-first(A::Union{Sum,Composition}) = @inbounds A[1]
-last(A::Union{Sum{N},Composition{N}}) where {N} = @inbounds A[N]
-firstindex(A::Union{Sum,Composition}) = 1
-lastindex(A::Union{Sum{N},Composition{N}}) where {N} = N
-length(A::Union{Sum{N},Composition{N}}) where {N} = N
-eltype(::Type{<:Sum{N,T}}) where {N,T} = eltype(T)
-eltype(::Type{<:Composition{N,T}}) where {N,T} = eltype(T)
+# FIXME: Base.first(A::Mapping) = A
+# FIXME: Base.last(A::Mapping) = A
+Base.first(A::Union{Sum,Composition}) = first(terms(A))
+Base.last(A::Union{Sum,Composition}) = last(terms(A))
+Base.firstindex(A::Union{Sum,Composition}) = firstindex(terms(A))
+Base.lastindex(A::Union{Sum,Composition}) = lastindex(terms(A))
+Base.length(A::Union{Sum,Composition}) = length(terms(A))
+Base.keys(A::Union{Sum,Composition}) = keys(terms(A))
+# FIXME: Base.eltype(::Type{<:Sum{L,N,T}}) where {L,N,T} = eltype(T)
+# FIXME: Base.eltype(::Type{<:Composition{L,N,T}}) where {L,N,T} = eltype(T)
 Tuple(A::Union{Sum,Composition}) = terms(A)
 
-@inline @propagate_inbounds getindex(A::Union{Sum,Composition}, i) =
-    getindex(terms(A), i)
+@inline @propagate_inbounds getindex(A::Union{Sum,Composition}, i) = terms(A)[i]
 
 """
     input_type([P=Direct,] A)
@@ -809,8 +836,8 @@ function vcreate(::Type{P}, A::Sum, x,
     vcreate(P, A[1], x, scratch)
 end
 
-function apply!(α::Number, P::Type{<:Union{Direct,Adjoint}}, A::Sum{N},
-                x, scratch::Bool, β::Number, y) where {N}
+function apply!(α::Number, P::Type{<:Union{Direct,Adjoint}}, A::Sum{L,N},
+                x, scratch::Bool, β::Number, y) where {L,N}
     if α == 0
         # Just scale the destination.
         vscale!(y, β)
@@ -860,13 +887,12 @@ throw_unsupported_inverse_of_sum() =
 # To break the type barrier, this is done by a recursion.  The recursion is
 # just done in the other direction for the Adjoint or Inverse operation.
 
-function vcreate(::Type{<:Operations},
-                 A::Composition{N}, x, scratch::Bool) where {N}
-    error("it is not possible to create the output of a composition of mappings")
-end
+vcreate(::Type{<:Operations}, A::Composition, x, scratch::Bool) =
+    error("it is not possible to automatically create the output of a composition of mappings")
 
-function apply!(α::Number, ::Type{P}, A::Composition{N}, x, scratch::Bool,
-                β::Number, y) where {N,P<:Union{Direct,InverseAdjoint}}
+# FIXME: replace `N` by `end` in following methods.
+function apply!(α::Number, ::Type{P}, A::Composition{L,N}, x, scratch::Bool,
+                β::Number, y) where {L,N,P<:Union{Direct,InverseAdjoint}}
     if α == 0
         # Just scale the destination.
         vscale!(y, β)
@@ -879,8 +905,8 @@ function apply!(α::Number, ::Type{P}, A::Composition{N}, x, scratch::Bool,
     return y
 end
 
-function apply(::Type{P}, A::Composition{N}, x,
-               scratch::Bool) where {N,P<:Union{Direct,InverseAdjoint}}
+function apply(::Type{P}, A::Composition, x,
+               scratch::Bool) where {P<:Union{Direct,InverseAdjoint}}
     apply(P, *, terms(A), x, scratch)
 end
 
@@ -892,8 +918,8 @@ function apply(::Type{P}, ::typeof(*), ops::NTuple{N,Mapping}, x,
     apply(P, *, ops[1:N-1], w, scratch)
 end
 
-function apply!(α::Number, ::Type{P}, A::Composition{N}, x, scratch::Bool,
-                β::Number, y) where {N,P<:Union{Adjoint,Inverse}}
+function apply!(α::Number, ::Type{P}, A::Composition{L,N}, x, scratch::Bool,
+                β::Number, y) where {L,N,P<:Union{Adjoint,Inverse}}
     if α == 0
         # Just scale the destination.
         vscale!(y, β)
@@ -906,17 +932,17 @@ function apply!(α::Number, ::Type{P}, A::Composition{N}, x, scratch::Bool,
     return y
 end
 
-function apply(::Type{P}, A::Composition{N}, x,
-               scratch::Bool) where {N,P<:Union{Adjoint,Inverse}}
-    apply(P, *, terms(A), x, scratch)
+function apply(::Type{P}, A::Composition, x,
+               scratch::Bool) where {P<:Union{Adjoint,Inverse}}
+    return apply(P, *, terms(A), x, scratch)
 end
 
 function apply(::Type{P}, ::typeof(*), ops::NTuple{N,Mapping}, x,
-               scratch::Bool) where {N,P<:Union{Adjoint,Inverse}}
+                scratch::Bool) where {N,P<:Union{Adjoint,Inverse}}
     w = apply(P, ops[1], x, scratch)
     N == 1 && return w
     scratch = overwritable(scratch, w, x)
-    apply(P, *, ops[2:N], w, scratch)
+    return apply(P, *, ops[2:N], w, scratch)
 end
 
 # Default rules to apply a Gram operator.  Gram matrices are Hermitian by
diff --git a/src/rules.jl b/src/rules.jl
index b0dc94f..98e2089 100644
--- a/src/rules.jl
+++ b/src/rules.jl
@@ -8,7 +8,7 @@
 # This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
 # released under the MIT "Expat" license.
 #
-# Copyright (c) 2017-2021 Éric Thiébaut.
+# Copyright (c) 2017-2022 Éric Thiébaut.
 #
 
 #------------------------------------------------------------------------------
@@ -37,16 +37,8 @@ isone(::Mapping) = false
 # must hold is that T(A), with T an unqualified type constructor, always yields
 # an instance of T.
 Direct(A::Mapping) = A # provided for completeness
-Adjoint(A::T) where {T<:Mapping} = Adjoint{T}(A)
-Inverse(A::T) where {T<:Mapping} = Inverse{T}(A)
-InverseAdjoint(A::T) where {T<:Mapping} = InverseAdjoint{T}(A)
-Gram(A::T) where {T<:Mapping} = Gram{T}(A)
-Jacobian(A::M, x::T) where {M<:Mapping,T} = Jacobian{M,T}(A, x)
-Scaled(α::S, A::T) where {S<:Number,T<:Mapping} = Scaled{T,S}(α, A)
-Sum(ops::Mapping...) = Sum(ops)
-Sum(ops::T) where {N,T<:NTuple{N,Mapping}} = Sum{N,T}(ops)
-Composition(ops::Mapping...) = Composition(ops)
-Composition(ops::T) where {N,T<:NTuple{N,Mapping}} = Composition{N,T}(ops)
+@inline Sum(terms::Mapping...) = Sum(terms)
+@inline Composition(terms::Mapping...) = Composition(terms)
 
 # Qualified outer constructors to forbid decoration of mappings of specific
 # types when, according to the simplification rules, another more simple
@@ -81,15 +73,13 @@ for (func, blacklist) in ((:Adjoint,        (:Identity,
                                              :Scaled,)),
                           (:Jacobian,       (:Scaled,)),
                           (:Scaled,         (:Scaled,)))
-    for T in blacklist
+    for M in blacklist
         if func === :Scaled
-            @eval $func{T,S}(α::S, A::T) where {S<:Number,T<:$T} =
-                illegal_call_to($func, T)
+            @eval $func(λ::Number, A::$M) = illegal_call_to($func, typeof(A))
         elseif func === :Jacobian
-            @eval $func{M,T}(A::M, x::T) where {M<:$T,T} =
-                illegal_call_to($func, M)
+            @eval $func(A::$M, x) = illegal_call_to($func, typeof(A))
         else
-            @eval $func{T}(A::T) where {T<:$T} = illegal_call_to($func, T)
+            @eval $func(A::$M) = illegal_call_to($func, typeof(A))
         end
     end
 end
@@ -112,7 +102,7 @@ end
 
 @noinline illegal_call_to(::Type{Jacobian}, T::Type) =
     bad_argument("the `Jacobian` constructor cannot be applied to an instance of `",
-                 brief(T), "`, use an expression like `∇(A,x)`")
+                 brief(T), "`, use expressions like `∇(A,x)` or `jacobian(A,x)`")
 
 @noinline illegal_call_to(::Type{Scaled}, T::Type) =
     bad_argument("the `Scaled` constructor cannot be applied to an instance of `",
@@ -127,7 +117,7 @@ brief(::Type{<:Scaled}        ) = "Scaled"
 brief(::Type{<:Sum}           ) = "Sum"
 brief(::Type{<:Composition}   ) = "Composition"
 brief(::Type{<:Identity}      ) = "Identity"
-brief(T::Type) = repr(T)
+brief(::Type{T}) where {T} = repr(T)
 
 #------------------------------------------------------------------------------
 # SCALED TYPE
@@ -144,14 +134,10 @@ brief(T::Type) = repr(T)
 #------------------------------------------------------------------------------
 # ADJOINT TYPE
 
-# Adjoint for non-specific mappings.
-adjoint(A::Mapping) = _adjoint(LinearType(A), A)
-
-_adjoint(::Linear, A::Mapping) = _adjoint(Linear(), SelfAdjointType(A), A)
-_adjoint(::Linear, ::SelfAdjoint, A::Mapping) = A
-_adjoint(::Linear, ::NonSelfAdjoint, A::Mapping) = Adjoint(A)
-_adjoint(::NonLinear, A::Mapping) =
-    throw_forbidden_adjoint_of_non_linear_mapping()
+# Adjoint for non-specific mappings. This method is called to implement he
+# syntax A'. NOTE: Self-Adjoint linear mappings such as `Gram` shall extend
+# this method to behave like the identity.
+adjoint(A::LinearMapping) = Adjoint(A)
 
 # Adjoint for specific mapping types.
 adjoint(A::Identity) = Id
@@ -159,19 +145,20 @@ adjoint(A::Scaled) = conj(multiplier(A))*adjoint(unscaled(A))
 adjoint(A::Adjoint) = unveil(A)
 adjoint(A::Inverse) = inv(adjoint(unveil(A)))
 adjoint(A::InverseAdjoint) = inv(unveil(A))
-adjoint(A::Jacobian) = Jacobian(A)
+adjoint(A::Jacobian) = Jacobian(A,variables(A)) # take the next derivative for the same variables
 adjoint(A::Gram) = A
 adjoint(A::Composition) =
     # It is assumed that the composition has already been simplified, so we
     # just apply the mathematical formula for the adjoint of a composition.
     Composition(reversemap(adjoint, terms(A)))
 
-function adjoint(A::Sum{N}) where {N}
+# FIXME: Not type-stable!
+function adjoint(A::Sum)
     # It is assumed that the sum has already been simplified, so we just apply
     # the mathematical formula for the adjoint of a sum and sort the resulting
     # terms.
-    B = Vector{Mapping}(undef, N)
-    @inbounds for i in 1:N
+    B = Vector{Mapping}(undef, length(A))
+    @inbounds for i in eachindex(A)
         B[i] = adjoint(A[i])
     end
     return Sum(to_tuple(sort!(B)))
@@ -187,8 +174,8 @@ end
     ∇(A, x)
 
 yields a result corresponding to the Jacobian (first partial derivatives) of
-the linear mapping `A` for the variables `x`.  If `A` is a linear mapping,
-`A` is returned whatever `x`.
+the mapping `A` for the variables `x`. If `A` is a linear mapping, `A` is
+returned whatever `x`.
 
 The call
 
@@ -198,9 +185,8 @@ is an alias for `∇(A,x)`.
 
 """
 ∇(A::Mapping, x) = jacobian(A, x)
-jacobian(A::Mapping, x) = _jacobian(LinearType(A), A, x)
-_jacobian(::Linear, A::Mapping, x) = A
-_jacobian(::NonLinear, A::Mapping, x) = Jacobian(A, x)
+jacobian(A::NonLinearMapping, x) = Jacobian(A, x)
+jacobian(A::LinearMapping, x) = A
 jacobian(A::Scaled, x) = multiplier(A)*jacobian(unscaled(A), x)
 
 @doc @doc(∇) jacobian
@@ -209,13 +195,12 @@ jacobian(A::Scaled, x) = multiplier(A)*jacobian(unscaled(A), x)
 # INVERSE TYPE
 
 # Inverse for non-specific mappings (a simple mapping or a sum or mappings).
-inv(A::T) where {T<:Mapping} = Inverse{T}(A)
+inv(A::Mapping) = Inverse(A)
 
 # Inverse for specific mapping types.
 inv(A::Identity) = Id
-inv(A::Scaled) = (is_linear(unscaled(A)) ?
-                  inv(multiplier(A))*inv(unscaled(A)) :
-                  inv(unscaled(A))*(inv(multiplier(A))*Id))
+inv(A::Scaled{true}) = inv(multiplier(A))*inv(unscaled(A))
+inv(A::Scaled) = inv(unscaled(A))*(inv(multiplier(A))*Id)
 inv(A::Inverse) = unveil(A)
 inv(A::InverseAdjoint) = adjoint(unveil(A))
 inv(A::Adjoint) = InverseAdjoint(unveil(A))
@@ -259,15 +244,14 @@ end
     add(A, B)
 
 performs the final stage of simplifying the sum `A + B` of mappings `A` and
-`B`.  This method assumes that any other simplifications than those involving
+`B`. This method assumes that any other simplifications than those involving
 sum of sums have been performed on `A + B` and that `A` and `B` have already
 been simplified (individually).
 
 Trivial simplifications of the composition `A + B` of mappings `A` and `B` must
 be done by specializing the operator `+` and `sum(A,B)` may eventually be
-called to simplify the sum of `A` and `B` when one of these may be a
-sum.  This method just returns `Sum(A,B)` when none of its
-arguments is a composition.
+called to simplify the sum of `A` and `B` when one of these may be a sum. This
+method just returns `Sum(A,B)` when none of its arguments is a composition.
 
 
 yields a simplified sum `A + B` of the mappings `A` and `B`.  This helper
@@ -299,8 +283,8 @@ end
 
 # Add the terms of a sum one-by-one.  Since terms must be re-ordered, there are
 # no obvious better ways to recombine.
-function add!(A::Vector{Mapping}, B::Sum{N}) where {N}
-    @inbounds for i in 1:N
+function add!(A::Vector{Mapping}, B::Sum)
+    @inbounds for i in eachindex(B)
         add!(A, B[i])
     end
     return A
@@ -390,21 +374,21 @@ end
         compose(A, B)
     end
 *(A::Scaled, B::Scaled) =
-    if is_linear(A)
+    if A isa LinearMapping
         (multiplier(A)*multiplier(B))*(unscaled(A)*unscaled(B))
     else
         multiplier(A)*(unscaled(A)*B)
     end
 
 # Simplify compositions involving an inverse mapping.
-*(A::Inverse{T}, B::T) where {T<:Mapping} =
+*(A::Inverse{L,M}, B::M) where {L,M<:Mapping{L}} =
     identical(unveil(A), B) ? Id : compose(A, B)
-*(A::T, B::Inverse{T}) where {T<:Mapping} =
+*(A::M, B::Inverse{L,M}) where {L,M<:Mapping{L}} =
     identical(A, unveil(B)) ? Id : compose(A, B)
 *(A::Inverse, B::Inverse) = compose(A, B)
-*(A::InverseAdjoint{T}, B::Adjoint{T}) where {T<:Mapping} =
+*(A::InverseAdjoint{M}, B::Adjoint{M}) where {M<:LinearMapping} =
     identical(unveil(A), unveil(B)) ? Id : compose(A, B)
-*(A::Adjoint{T}, B::InverseAdjoint{T}) where {T<:Mapping} =
+*(A::Adjoint{M}, B::InverseAdjoint{M}) where {M<:LinearMapping} =
     identical(unveil(A), unveil(B)) ? Id : compose(A, B)
 *(A::InverseAdjoint, B::InverseAdjoint) = compose(A, B)
 
@@ -423,16 +407,14 @@ end
 #
 # In principle, if forming the adjoint has been allowed, it is not needed to
 # check whether operands are linear mappings.
-*(A::Adjoint{T}, B::T) where {T<:Mapping} =
+*(A::Adjoint{M}, B::M) where {M<:LinearMapping} =
     identical(unveil(A), B) ? Gram(B) : compose(A, B)
-*(A::T, B::Adjoint{T}) where {T<:Mapping} =
+*(A::M, B::Adjoint{M}) where {M<:LinearMapping} =
     identical(A, unveil(B)) ? Gram(B) : compose(A, B)
-*(A::Inverse{T}, B::InverseAdjoint{T}) where {T<:Mapping} =
-    identical(unveil(A), unveil(B)) ? Inverse(Gram(unveil(A))) :
-    compose(A, B)
-*(A::InverseAdjoint{T}, B::Inverse{T}) where {T<:Mapping} =
-    identical(unveil(A), unveil(B)) ?
-    Inverse(Gram(Adjoint(unveil(A)))) : compose(A, B)
+*(A::Inverse{true,M}, B::InverseAdjoint{M}) where {M<:LinearMapping} =
+    identical(unveil(A), unveil(B)) ? Inverse(Gram(unveil(A))) : compose(A, B)
+*(A::InverseAdjoint{M}, B::Inverse{true,M}) where {M<:LinearMapping} =
+    identical(unveil(A), unveil(B)) ? Inverse(Gram(Adjoint(unveil(A)))) : compose(A, B)
 
 """
     compose(A,B)
@@ -479,16 +461,16 @@ mapping.
 
 """ compose!
 
-function compose!(A::Vector{Mapping}, B::Composition{N}) where {N}
-    @inbounds for i in 1:N
+function compose!(A::Vector{Mapping}, B::Composition)
+    @inbounds for i in eachindex(B)
         # Build the simplified composition A*B[i].
         compose!(A, B[i])
         if identical(last(A), B[i])
             # The last term of the simplified composition A*B[i] is still B[i],
             # which indicates that composing A with B[i] did not yield any
-            # simplifications.  It is sufficient to append all the other terms
+            # simplifications. It is sufficient to append all the other terms
             # of B to A as no further simplifications are expected.
-            return append_terms!(A, B, (i+1):N)
+            return append_terms!(A, B, (i+1):lastindex(B))
         end
     end
     return A
@@ -497,28 +479,28 @@ end
 function compose!(A::Vector{Mapping}, B::Mapping)
     # Compute the simplified composition of the last term of A with B.  The
     # result is either a simple mapping or a simplified composition.
-    m = length(A); @certify m > 0
-    C = A[m]*B
+    @certify length(A) > 0
+    C = last(A)*B
 
     # Replace the last term of A with C.
-    if C isa Composition && identical(C[1], A[m])
-        # Nothing has changed at the tail of the composition A.  No further
-        # simplifications are expected.  Push all terms of C to A, but the
-        # first term of C which is identical to the last term of A.
-        append_terms!(A, C, 2:length(C))
-    elseif m > 1
-        # Drop the last term of A and compose the remaining terms with C.  This
+    if C isa Composition && identical(last(A), first(C))
+        # Nothing has changed at the tail of the composition A. No further
+        # simplifications are expected. Push all terms of C to A, but the first
+        # term of C which is identical to the last term of A.
+        append_terms!(A, C, firstindex(C)+1:lastindex(C))
+    elseif length(A) > 1
+        # Drop the last term of A and compose the remaining terms with C. This
         # may trigger further simplifications
-        compose!(resize!(A, m - 1), C)
+        compose!(resize!(A, length(A) - 1), C)
     elseif C isa Composition
         # Replace the only term of A by all the terms of the composition C.
-        # This is the same as above but avoids calling `resize!` as a small
+        # This is the same as above but avoids one call to `resize!` as a small
         # optimization.
-        A[1] = C[1]
-        append_terms!(A, C, 2:length(C))
+        A[firstindex(A)] = C[firstindex(C)]
+        append_terms!(A, C, firstindex(C)+1:lastindex(C))
     else
         # Replace the only term of A by the simple mapping C.
-        A[1] = C
+        A[firstindex(A)] = C
     end
     return A
 end
@@ -529,10 +511,10 @@ end
 """
     as_vector(op, A)
 
-yields a vector of mappings with the terms of the mapping `A`.  Argument `op`
-is `+` or `*`.  If `op` is `+` (resp. `*`) and `A` is a sum (resp. a
-composition) of mappings, the terms of `A` are extracted in the returned
-vector; otherwise, the returned vector has just one element which is `A`.
+yields a vector of mappings with the terms of the mapping `A`. Argument `op` is
+`+` or `*`. If `op` is `+` (resp. `*`) and `A` is a sum (resp. a composition)
+of mappings, the terms of `A` are extracted in the returned vector; otherwise,
+the returned vector has just one element which is `A`.
 
 """
 function as_vector(::Union{typeof(+),typeof(*)}, A::Mapping)
@@ -547,14 +529,14 @@ as_vector(::typeof(*), A::Composition) = collect_terms(A)
 """
     collect_terms(A)
 
-collects the terms of the mapping `A` into a vector.  This is similar to
+collects the terms of the mapping `A` into a vector. This is similar to
 `collect(A)` except that the element type of the result is forced to be
 `Mapping`.
 
 """
-function collect_terms(A::Union{Sum{N},Composition{N}}) where {N}
-    V = Vector{Mapping}(undef, N)
-    @inbounds for i in 1:N
+function collect_terms(A::Union{Sum,Composition})
+    V = Vector{Mapping}(undef, length(A))
+    @inbounds for i in eachindex(A)
         V[i] = A[i]
     end
     return V
diff --git a/src/traits.jl b/src/traits.jl
index 1df4a2e..eb228ed 100644
--- a/src/traits.jl
+++ b/src/traits.jl
@@ -8,48 +8,31 @@
 # This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
 # released under the MIT "Expat" license.
 #
-# Copyright (c) 2017-2021 Éric Thiébaut.
+# Copyright (c) 2017-2022 Éric Thiébaut.
 #
 
 """
-
-```julia
-LinearType(A)
-```
+    LinearType(A)
 
 yields the *linear* trait of mapping `A` indicating whether `A` is certainly
-linear.  The returned value is one of the singletons `Linear()` for linear maps
-or `NonLinear()` for other mappings.
+linear. The returned value is one of the singletons `Linear()` for linear
+mappings or `NonLinear()` for other mappings.
 
 See also: [`Trait`](@ref), [`is_linear`](@ref).
 
 """
-LinearType(::Mapping) = NonLinear() # any mapping assumed non-linear by default
-LinearType(::LinearMapping) = Linear()
-LinearType(::Inverse{<:LinearMapping}) = Linear()
-LinearType(::Scaled{<:LinearMapping}) = Linear()
-LinearType(A::Inverse) = LinearType(unveil(A))
-LinearType(A::Scaled) = LinearType(unscaled(A))
-LinearType(A::Union{Sum,Composition}) =
-    (allof(x -> LinearType(x) === Linear(), terms(A)...) ?
-     Linear() : NonLinear())
-LinearType(A::Scaled{T,S}) where {T,S} =
-    # If the multiplier λ of a scaled mapping A = (λ⋅M) is zero, then A behaves
-    # linearly even though M is not a linear mapping.  FIXME: But acknowledging
-    # this as a linear mapping may give rise to troubles later.
-    (multiplier(A) == zero(S) ? Linear() : LinearType(unscaled(A)))
+LinearType(A::Mapping) = LinearType(typeof(A))
+LinearType(::Type{<:LinearMapping}) = Linear()
+LinearType(::Type{<:NonLinearMapping}) = NonLinear()
 
 @doc @doc(LinearType) Linear
 @doc @doc(LinearType) NonLinear
 
 """
-
-```julia
-SelfAdjointType(A)
-```
+    SelfAdjointType(A)
 
 yields the *self-adjoint* trait of mapping `A` indicating whether `A` is
-certainly a self-adjoint linear map.  The returned value is one of the
+certainly a self-adjoint linear map. The returned value is one of the
 singletons `SelfAdjoint()` for self-adjoint linear maps and `NonSelfAdjoint()`
 for other mappings.
 
@@ -68,13 +51,10 @@ SelfAdjointType(A::Gram) = SelfAdjoint()
 @doc @doc(SelfAdjointType) NonSelfAdjoint
 
 """
-
-```julia
-MorphismType(A)
-```
+    MorphismType(A)
 
 yields the *morphism* trait of mapping `A` indicating whether `A` is certainly
-an endomorphism (its input and output spaces are the same).  The returned value
+an endomorphism (its input and output spaces are the same). The returned value
 is one of the singletons `Endomorphism()` for mappings whose input and output
 spaces are the same or `Morphism()` for other mappings.
 
@@ -93,15 +73,12 @@ MorphismType(A::Union{Sum,Composition}) =
 @doc @doc(MorphismType) Endomorphism
 
 """
-
-```julia
-DiagonalType(A)
-```
+    DiagonalType(A)
 
 yields the *diagonal* trait of mapping `A` indicating whether `A` is certainly
-a diagonal linear mapping.  The returned value is one of the singletons
-`DiagonalMapping()` for diagonal linear maps or `NonDiagonalMapping()` for other
-mappings.
+a diagonal linear mapping. The returned value is one of the singletons
+`DiagonalMapping()` for diagonal linear maps or `NonDiagonalMapping()` for
+other mappings.
 
 See also: [`Trait`](@ref), [`is_diagonal`](@ref).
 
@@ -117,37 +94,23 @@ DiagonalType(A::Union{Sum,Composition}) =
 @doc @doc(DiagonalType) DiagonalMapping
 
 """
-```julia
-is_linear(A)
-```
-
-yields whether `A` is certainly a linear mapping.
+    is_linear(A)
 
-!!! note
-    This method is intended to perform certain automatic simplifications or
-    optimizations.  It is guaranted to return `true` when its argument is
-    certainly a linear mapping but it may return `false` even though its
-    argument behaves linearly because it is not always possible to figure out
-    that a complex mapping assemblage has this property.
+yields whether `A` is a linear mapping.
 
 See also: [`LinearType`](@ref).
 
 """
-is_linear(A::LinearMapping) = true
-is_linear(A::Mapping) = _is_linear(LinearType(A))
-_is_linear(::Linear) = true
-_is_linear(::NonLinear) = false
+is_linear(x) = (LinearType(x) === Linear())
 
 """
-```julia
-is_selfadjoint(A)
-```
+    is_selfadjoint(A)
 
 yields whether mapping `A` is certainly a self-adjoint linear mapping.
 
 !!! note
-    This method is intended to perform certain automatic simplifications or
-    optimizations.  It is guaranted to return `true` when its argument is
+    This method is called to perform certain automatic simplifications or
+    optimizations. It is guaranted to return `true` when its argument is
     certainly a self-adjoint linear mapping but it may return `false` even
     though its argument behaves like a self-adjoint linear map because it is
     not always possible to figure out that a complex mapping construction has
@@ -157,20 +120,16 @@ yields whether mapping `A` is certainly a self-adjoint linear mapping.
 See also: [`SelfAdjointType`](@ref).
 
 """
-is_selfadjoint(A::Mapping) = _is_selfadjoint(SelfAdjointType(A))
-_is_selfadjoint(::SelfAdjoint) = true
-_is_selfadjoint(::NonSelfAdjoint) = false
+is_selfadjoint(x) = (SelfAdjointType(x) === SelfAdjoint())
 
 """
-```julia
-is_endomorphism(A)
-```
+    is_endomorphism(A)
 
-yields whether mapping `A` is certainly an endomorphism.
+yields whether `A` is certainly an endomorphism.
 
 !!! note
-    This method is intended to perform certain automatic simplifications or
-    optimizations.  It is guaranted to return `true` when its argument is
+    This method is called to perform certain automatic simplifications or
+    optimizations. It is guaranted to return `true` when its argument is
     certainly an endomorphism but it may return `false` even though its
     argument behaves like an endomorphism because it is not always possible to
     figure out that a complex mapping assemblage has this property.
@@ -178,20 +137,16 @@ yields whether mapping `A` is certainly an endomorphism.
 See also: [`MorphismType`](@ref).
 
 """
-is_endomorphism(A::Mapping) = _is_endomorphism(MorphismType(A))
-_is_endomorphism(::Endomorphism) = true
-_is_endomorphism(::Morphism) = false
+is_endomorphism(x) = (MorphismType(x) === Endomorphism())
 
 """
-```julia
-is_diagonal(A)
-```
+    is_diagonal(A)
 
 yields whether mapping `A` is certainly a diagonal linear map.
 
 !!! note
-    This method is intended to perform certain automatic simplifications or
-    optimizations.  It is guaranted to return `true` when its argument is
+    This method is called to perform certain automatic simplifications or
+    optimizations. It is guaranted to return `true` when its argument is
     certainly a diagonal linear map but it may return `false` even though its
     argument behaves like a diagonal linear map because it is not always
     possible to figure out that a complex mapping assemblage has this property.
@@ -199,6 +154,4 @@ yields whether mapping `A` is certainly a diagonal linear map.
 See also: [`DiagonalType`](@ref).
 
 """
-is_diagonal(A::Mapping) = _is_diagonal(DiagonalType(A))
-_is_diagonal(::DiagonalMapping) = true
-_is_diagonal(::NonDiagonalMapping) = false
+is_diagonal(x) = (DiagonalType(x) === DiagonalMapping())
diff --git a/src/types.jl b/src/types.jl
index db432cd..c9ca98a 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -71,10 +71,15 @@ This definition closely follows the semantic used in the BLAS module that
 const Floats = Union{Reals,Complexes}
 
 """
+    Mapping{L}
 
-A `Mapping` is any function between two variables spaces. Assuming upper case
-Latin letters denote mappings, lower case Latin letters denote variables, and
-Greek letters denote scalars, then:
+is the abstract super-type of mappings in `LazyAlgebra`. Type parameter `L` is
+a boolean indicating whether the mapping is linear. A mapping is any function
+between two variables spaces. Assuming upper case Latin letters denote
+mappings, lower case Latin letters denote variables, and Greek letters denote
+scalars, then:
+
+* `A(x)` yields the result of applying the mapping `A` to `x`;
 
 * `A*x` or `A⋅x` yields the result of applying the mapping `A` to `x`;
 
@@ -137,10 +142,23 @@ See also: [`apply`](@ref), [`apply!`](@ref), [`vcreate`](@ref),
           [`Adjoint`](@ref), [`Inverse`](@ref), [`InverseAdjoint`](@ref).
 
 """
-abstract type Mapping <: Function end
+abstract type Mapping{L} <: Function end
+
+"""
+    LinearMapping
+
+is the abstract super-type of linear mappings in `LazyAlgebra`.
+
+"""
+const LinearMapping = Mapping{true}
 
-abstract type LinearMapping <: Mapping end
-@doc @doc(Mapping) LinearMapping
+"""
+    NonLinearMapping
+
+is the abstract super-type of non-linear mappings in `LazyAlgebra`.
+
+"""
+const NonLinearMapping = Mapping{false}
 
 """
     Identity()
@@ -217,15 +235,9 @@ Call [`unveil(obj)`](@ref) to reveal the linear mapping `A` embedded in `obj`.
 See also [`DecoratedMapping`](@ref).
 
 """
-struct Adjoint{T<:Mapping} <: LinearMapping
-    op::T
-
-    # The outer constructors prevent most illegal calls to `Adjoint(A)` we
-    # just have to check that the argument is a simple linear mapping.
-    function Adjoint{T}(A::T) where {T<:Mapping}
-        is_linear(A) || throw_forbidden_adjoint_of_non_linear_mapping()
-        return new{T}(A)
-    end
+struct Adjoint{M<:LinearMapping} <: LinearMapping
+    parent::M
+    Adjoint(A::M) where {M<:LinearMapping} = new{M}(A)
 end
 
 """
@@ -241,12 +253,9 @@ Call [`unveil(obj)`](@ref) to reveal the mapping `A` embedded in `obj`.
 See also [`DecoratedMapping`](@ref).
 
 """
-struct Inverse{T<:Mapping} <: Mapping
-    op::T
-
-    # The outer constructors prevent all illegal calls to `Inverse(A)` so there
-    # is nothing more to check.
-    Inverse{T}(A::T) where {T<:Mapping} = new{T}(A)
+struct Inverse{L,M<:Mapping{L}} <: Mapping{L}
+    parent::M
+    Inverse(A::M) where {L,M<:Mapping{L}} = new{L,M}(A)
 end
 
 """
@@ -266,44 +275,31 @@ Call [`unveil(obj)`](@ref) to reveal the mapping `A` embedded in `obj`.
 See also [`DecoratedMapping`](@ref).
 
 """
-struct InverseAdjoint{T<:Mapping} <: LinearMapping
-    op::T
-
-    # The outer constructors prevent most illegal calls to `InverseAdjoint(A)`
-    # we just have to check that the argument is a simple linear mapping.
-    function InverseAdjoint{T}(A::T) where {T<:Mapping}
-        is_linear(A) ||
-            bad_argument("taking the inverse adjoint of non-linear mappings is not allowed")
-        return new{T}(A)
-    end
+struct InverseAdjoint{M<:LinearMapping} <: LinearMapping
+    parent::M
+    InverseAdjoint(A::M) where {M<:LinearMapping} = new{M}(A)
 end
 
-const AdjointInverse{T} = InverseAdjoint{T}
+const AdjointInverse{M} = InverseAdjoint{M}
 @doc @doc(InverseAdjoint) AdjointInverse
 
 """
-    Gram(A) -> obj
+    Gram(A) -> B
 
-yields an object instance `obj` representing the composition `A'*A` for the
-linear mapping `A`.
+yields an object `B` representing the composition `A'*A` for the linear mapping
+`A`.
 
 Directly calling this constructor is discouraged, call [`gram(A)`](@ref) or use
 expression `A'*A` instead and benefit from automatic simplification rules.
 
-Call [`unveil(obj)`](@ref) to reveal the linear mapping `A` embedded in `obj`.
+Call [`unveil(B)`](@ref) to reveal the linear mapping `A` embedded in `B`.
 
 See also [`gram`](@ref), [`unveil`](@ref) and [`DecoratedMapping`](@ref).
 
 """
-struct Gram{T<:Mapping} <: LinearMapping
-    op::T
-
-    # The outer constructors prevent most illegal calls to `Gram(A)` we
-    # just have to check that the argument is a simple linear mapping.
-    function Gram{T}(A::T) where {T<:Mapping}
-        is_linear(A) || throw_forbidden_Gram_of_non_linear_mapping()
-        return new{T}(A)
-    end
+struct Gram{M<:LinearMapping} <: LinearMapping
+    parent::M
+    Gram(A::M) where {M<:LinearMapping} = new{M}(A)
 end
 
 """
@@ -328,18 +324,10 @@ Directly calling this constructor is discouraged, call [`jacobian(A,x)`](@ref)
 or [`∇(A,x)`](@ref) instead and benefit from automatic simplification rules.
 
 """
-struct Jacobian{M<:Mapping,T} <: Mapping
-    A::M
-    x::T
-
-    # The outer constructors prevent most illegal calls to `Jacobian(A)` we
-    # just have to check that the argument is not a simple linear mapping.
-    function Jacobian{M,T}(A::M, x::T) where {M<:Mapping,T}
-        is_linear(A) &&
-            bad_argument("the Jacobian of a linear mapping of type `",
-                         M, "` should be the mapping itself")
-        return new{M,T}(A, x)
-    end
+struct Jacobian{M<:NonLinearMapping,V} <: NonLinearMapping
+    primitive::M
+    variables::V
+    Jacobian(A::M, x::V) where {M<:NonLinearMapping,V} = new{M,V}(A, x)
 end
 
 """
@@ -355,70 +343,63 @@ See also: [`apply`](@ref) and [`apply!`](@ref).
 const Operations = Union{Direct,Adjoint,Inverse,InverseAdjoint}
 
 """
-    Scaled(λ, M) -> obj
+    Scaled(λ, A) -> B
 
-yields an object instance `obj` representing `λ*M`, the mapping `M` multiplied
-by a scalar `λ`.
+yields an object `B` representing `λ*A`, that is the mapping `A` multiplied by
+a scalar `λ`.
 
-Directly calling this constructor is discouraged, use expressions like `λ*M`
+Directly calling this constructor is discouraged, use expressions like `λ*A`
 instead and benefit from automatic simplification rules.
 
-Call [`multiplier(obj)`](@ref) and [`unscaled(obj)`](@ref) with a scaled
-mapping `obj = λ*M` to retrieve `λ` and `M` respectively.
+Call [`multiplier(B)`](@ref) and [`unscaled(B)`](@ref) with a scaled
+mapping `B = λ*A` to retrieve `λ` and `A` respectively.
 
 """
-struct Scaled{T<:Mapping,S<:Number} <: Mapping
-    λ::S
-    M::T
-    Scaled{T,S}(λ::S, M::Mapping) where {S<:Number,T<:Mapping} =
-        new{T,S}(λ, M)
+struct Scaled{L,M<:Mapping{L},S<:Number} <: Mapping{L}
+    multiplier::S
+    mapping::M
+    Scaled(λ::S, A::M) where {S<:Number,L,M<:Mapping{L}} = new{L,M,S}(λ, A)
 end
 
 """
-    Sum(A, B, ...) -> obj
+    Sum(A, B...) -> S
 
-yields an object instance `obj` representing the sum `A + B + ...` of the
-mappings `A`, `B`, ...
+yields an object `S` representing the sum `A + B + ...` of the mappings `A`,
+`B...`. The constructor also accepts as argument a tuple of mappings.
 
 Directly calling this constructor is discouraged, use expressions like `A + B +
 ...` instead and benefit from automatic simplification rules.
 
-Call [`terms(obj)`](@ref) retrieve the tuple `(A,B,...)` of the terms of the
-sum stored in `obj`.
+Call [`terms(S)`](@ref) retrieve the tuple `(A,B...)` of the terms of the sum
+stored in `S`.
 
 """
-struct Sum{N,T<:NTuple{N,Mapping}} <: Mapping
-    ops::T
-
-    # The inner constructor ensures that the number of arguments is at least 2.
-    function Sum{N,T}(ops::T) where {N,T<:NTuple{N,Mapping}}
-        N ≥ 2 ||
-            throw(ArgumentError("a sum of mappings has at least 2 components"))
-        new{N,T}(ops)
+struct Sum{L,N,T<:NTuple{N,Mapping}} <: Mapping{L}
+    terms::T
+    function Sum(terms::T) where {N,T<:NTuple{N,Mapping}}
+        N ≥ 2 || throw(ArgumentError("a sum of mappings has at least 2 terms"))
+        return new{T <: NTuple{N,Mapping{true}},N,T}(terms)
     end
 end
 
 """
-    Composition(A, B, ...) -> obj
+    Composition(A, B...) -> C
 
-yields an object instance `obj` representing the composition `A*B*...` of the
-mappings `A`, `B`, ...
+yields an object `C` representing the composition `A*B*...` of the mappings
+`A`, `B...`. The constructor also accepts as argument a tuple of mappings.
 
 Directly calling this constructor is discouraged, use expressions like
 `A*B*...` `A∘B∘...` or `A⋅B⋅...` instead and benefit from automatic
 simplification rules.
 
-Call [`terms(obj)`](@ref) retrieve the tuple `(A,B,...)` of the terms of the
-composition stored in `obj`.
+Call [`terms(C)`](@ref) to retrieve the tuple `(A,B...)` of the terms of the
+composition stored in `C`.
 
 """
-struct Composition{N,T<:NTuple{N,Mapping}} <: Mapping
-    ops::T
-
-    # The inner constructor ensures that the number of arguments is at least 2.
-    function Composition{N,T}(ops::T) where {N,T<:NTuple{N,Mapping}}
-        N ≥ 2 ||
-            throw(ArgumentError("a composition of mappings has at least 2 components"))
-        new{N,T}(ops)
+struct Composition{L,N,T<:NTuple{N,Mapping}} <: Mapping{L}
+    terms::T
+    function Composition(terms::T) where {N,T<:NTuple{N,Mapping}}
+        N ≥ 2 || throw(ArgumentError("a composition of mappings has at least 2 terms"))
+        return new{T <: NTuple{N,Mapping{true}},N,T}(terms)
     end
 end
diff --git a/test/fft-tests.jl b/test/fft-tests.jl
index da4825b..f0193b2 100644
--- a/test/fft-tests.jl
+++ b/test/fft-tests.jl
@@ -70,7 +70,7 @@ end # testset
         show(io, F)
         @test String(take!(io)) == "FFT"
         @test typeof(F') === Adjoint{typeof(F)}
-        @test typeof(inv(F)) === Inverse{typeof(F)}
+        @test typeof(inv(F)) === Inverse{true,typeof(F)}
         @test typeof(inv(F)') === InverseAdjoint{typeof(F)}
         @test typeof(inv(F')) === InverseAdjoint{typeof(F)}
         @test F'*F == length(x)*Identity()
diff --git a/test/rules-tests.jl b/test/rules-tests.jl
index fe84fb8..08067aa 100644
--- a/test/rules-tests.jl
+++ b/test/rules-tests.jl
@@ -111,8 +111,8 @@ end
 
     # Basic methods for sums and compositions of mappings.
     let A1 = A + B + M, A2 = C*B*M*Q
-        @test eltype(A1) <: Mapping
-        @test eltype(A2) <: Mapping
+        # FIXME: @test eltype(A1) <: Mapping
+        # FIXME: @test eltype(A2) <: Mapping
         @test Tuple(A1) === terms(A1)
         @test Tuple(A2) === terms(A2)
         @test length(A1) == 3
@@ -159,10 +159,10 @@ end
     @test ∇(A,x) === A
     @test (3A)' === 3*(A')
     @test (A + 2B)' - A' === 2*B'
-    @test_throws ArgumentError Jacobian(A,x)
-    @test_throws ArgumentError M'
-    @test_throws ArgumentError adjoint(M)
-    @test_throws ArgumentError Adjoint(M)
+    @test_throws MethodError Jacobian(A,x)
+    @test_throws MethodError M'
+    @test_throws MethodError adjoint(M)
+    @test_throws MethodError Adjoint(M)
     @test jacobian(M,x) isa Jacobian
     @test ∇(M,x) === jacobian(M,x)
     @test ∇(3M,x) === 3*∇(M,x)