diff --git a/ext/DynamicQuantitiesUnitfulExt.jl b/ext/DynamicQuantitiesUnitfulExt.jl
index 045c027e..16f09420 100644
--- a/ext/DynamicQuantitiesUnitfulExt.jl
+++ b/ext/DynamicQuantitiesUnitfulExt.jl
@@ -1,6 +1,7 @@
 module DynamicQuantitiesUnitfulExt
 
-import DynamicQuantities
+using DynamicQuantities: DynamicQuantities, ABSTRACT_QUANTITY_TYPES
+
 import Unitful
 import Unitful: @u_str
 
@@ -23,30 +24,33 @@ function unitful_equivalences()
     return NamedTuple((k => si_units[k] for k in keys(si_units)))
 end
 
-Base.convert(::Type{Unitful.Quantity}, x::DynamicQuantities.Quantity) =
-    let
-        validate_upreferred()
-        cumulator = DynamicQuantities.ustrip(x)
-        dims = DynamicQuantities.dimension(x)
-        if dims isa DynamicQuantities.SymbolicDimensions
-            throw(ArgumentError("Conversion of a `DynamicQuantities.Quantity` to a `Unitful.Quantity` is not defined with dimensions of type `SymbolicDimensions`. Instead, you can first use the `uexpand` function to convert the dimensions to their base SI form of type `Dimensions`, then convert this quantity to a `Unitful.Quantity`."))
+for (_, _, Q) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        function Base.convert(::Type{Unitful.Quantity}, x::$Q)
+            validate_upreferred()
+            cumulator = DynamicQuantities.ustrip(x)
+            dims = DynamicQuantities.dimension(x)
+            if dims isa DynamicQuantities.SymbolicDimensions
+                throw(ArgumentError("Conversion of a `DynamicQuantities." * string($Q) * "` to a `Unitful.Quantity` is not defined with dimensions of type `SymbolicDimensions`. Instead, you can first use the `uexpand` function to convert the dimensions to their base SI form of type `Dimensions`, then convert this quantity to a `Unitful.Quantity`."))
+            end
+            equiv = unitful_equivalences()
+            for dim in keys(dims)
+                value = dims[dim]
+                iszero(value) && continue
+                cumulator *= equiv[dim]^value
+            end
+            cumulator
         end
-        equiv = unitful_equivalences()
-        for dim in keys(dims)
-            value = dims[dim]
-            iszero(value) && continue
-            cumulator *= equiv[dim]^value
+        function Base.convert(::Type{$Q}, x::Unitful.Quantity{T}) where {T}
+            return convert($Q{T,DynamicQuantities.DEFAULT_DIM_TYPE}, x)
+        end
+        function Base.convert(::Type{$Q{T,D}}, x::Unitful.Quantity) where {T,R,D<:DynamicQuantities.AbstractDimensions{R}}
+            value = Unitful.ustrip(Unitful.upreferred(x))
+            dimension = convert(D, Unitful.dimension(x))
+            return $Q(convert(T, value), dimension)
         end
-        cumulator
-    end
-
-Base.convert(::Type{DynamicQuantities.Quantity}, x::Unitful.Quantity{T}) where {T} = convert(DynamicQuantities.Quantity{T,DynamicQuantities.DEFAULT_DIM_TYPE}, x)
-Base.convert(::Type{DynamicQuantities.Quantity{T,D}}, x::Unitful.Quantity) where {T,R,D<:DynamicQuantities.AbstractDimensions{R}} =
-    let
-        value = Unitful.ustrip(Unitful.upreferred(x))
-        dimension = convert(D, Unitful.dimension(x))
-        return DynamicQuantities.Quantity(convert(T, value), dimension)
     end
+end
 
 Base.convert(::Type{DynamicQuantities.Dimensions}, d::Unitful.Dimensions) = convert(DynamicQuantities.DEFAULT_DIM_TYPE, d)
 Base.convert(::Type{DynamicQuantities.Dimensions{R}}, d::Unitful.Dimensions{D}) where {R,D} =
diff --git a/src/DynamicQuantities.jl b/src/DynamicQuantities.jl
index 08d06121..bd0e8d8b 100644
--- a/src/DynamicQuantities.jl
+++ b/src/DynamicQuantities.jl
@@ -1,8 +1,8 @@
 module DynamicQuantities
 
 export Units, Constants
-export AbstractDimensions, AbstractQuantity, AbstractGenericQuantity, UnionAbstractQuantity
-export Quantity, GenericQuantity, Dimensions, SymbolicDimensions, QuantityArray, DimensionError
+export AbstractDimensions, AbstractQuantity, AbstractGenericQuantity, AbstractRealQuantity, UnionAbstractQuantity
+export Quantity, GenericQuantity, RealQuantity, Dimensions, SymbolicDimensions, QuantityArray, DimensionError
 export ustrip, dimension
 export ulength, umass, utime, ucurrent, utemperature, uluminosity, uamount
 export uparse, @u_str, sym_uparse, @us_str, uexpand, uconvert
diff --git a/src/arrays.jl b/src/arrays.jl
index 775b7ce7..b01d8ced 100644
--- a/src/arrays.jl
+++ b/src/arrays.jl
@@ -52,12 +52,13 @@ struct QuantityArray{T,N,D<:AbstractDimensions,Q<:UnionAbstractQuantity{T,D},V<:
     end
 end
 
-# Construct with a Quantity (easier, as you can use the units):
 QuantityArray(v::AbstractArray; kws...) = QuantityArray(v, DEFAULT_DIM_TYPE(; kws...))
 for (type, base_type, default_type) in ABSTRACT_QUANTITY_TYPES
-    @eval begin
-        QuantityArray(v::AbstractArray{<:$base_type}, q::$type) = QuantityArray(v .* ustrip(q), dimension(q), typeof(q))
-        QuantityArray(v::AbstractArray{<:$base_type}, d::AbstractDimensions) = QuantityArray(v, d, $default_type)
+    @eval QuantityArray(v::AbstractArray{<:$base_type}, q::$type) = QuantityArray(v .* ustrip(q), dimension(q), typeof(q))
+
+    # Only define defaults for Quantity and GenericQuantity. Other types, the user needs to declare explicitly.
+    if type in (AbstractQuantity, AbstractGenericQuantity)
+        @eval QuantityArray(v::AbstractArray{<:$base_type}, d::AbstractDimensions) = QuantityArray(v, d, $default_type)
     end
 end
 QuantityArray(v::QA) where {Q<:UnionAbstractQuantity,QA<:AbstractArray{Q}} =
diff --git a/src/constants.jl b/src/constants.jl
index 736a71b4..6a39ca2f 100644
--- a/src/constants.jl
+++ b/src/constants.jl
@@ -1,7 +1,6 @@
 module Constants
 
 import ..DEFAULT_QUANTITY_TYPE
-import ..Quantity
 import ..Units as U
 import ..Units: _add_prefixes
 
diff --git a/src/disambiguities.jl b/src/disambiguities.jl
index a0571754..49a3cf2b 100644
--- a/src/disambiguities.jl
+++ b/src/disambiguities.jl
@@ -1,14 +1,108 @@
-Base.isless(::AbstractQuantity, ::Missing) = missing
-Base.isless(::Missing, ::AbstractQuantity) = missing
-Base.:(==)(::AbstractQuantity, ::Missing) = missing
-Base.:(==)(::Missing, ::AbstractQuantity) = missing
-Base.isapprox(::AbstractQuantity, ::Missing; kws...) = missing
-Base.isapprox(::Missing, ::AbstractQuantity; kws...) = missing
+for op in (:isless, :(==), :isequal, :(<)), (type, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        Base.$(op)(::$type, ::Missing) = missing
+        Base.$(op)(::Missing, ::$type) = missing
+    end
+end
+for op in (:isapprox,), (type, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        Base.$(op)(::$type, ::Missing; kws...) = missing
+        Base.$(op)(::Missing, ::$type; kws...) = missing
+    end
+end
 
-Base.:(==)(::AbstractQuantity, ::WeakRef) = error("Cannot compare a quantity to a weakref")
-Base.:(==)(::WeakRef, ::AbstractQuantity) = error("Cannot compare a weakref to a quantity")
+for (type, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        Base.:(==)(::$type, ::WeakRef) = error("Cannot compare a quantity to a weakref")
+        Base.:(==)(::WeakRef, ::$type) = error("Cannot compare a weakref to a quantity")
+    end
+end
 
 Base.:*(l::AbstractDimensions, r::Number) = error("Please use an `UnionAbstractQuantity` for multiplication. You used multiplication on types: $(typeof(l)) and $(typeof(r)).")
 Base.:*(l::Number, r::AbstractDimensions) = error("Please use an `UnionAbstractQuantity` for multiplication. You used multiplication on types: $(typeof(l)) and $(typeof(r)).")
 Base.:/(l::AbstractDimensions, r::Number) = error("Please use an `UnionAbstractQuantity` for division. You used division on types: $(typeof(l)) and $(typeof(r)).")
 Base.:/(l::Number, r::AbstractDimensions) = error("Please use an `UnionAbstractQuantity` for division. You used division on types: $(typeof(l)) and $(typeof(r)).")
+
+# Promotion ambiguities
+function Base.promote_rule(::Type{F}, ::Type{Bool}) where {F<:FixedRational}
+    return F
+end
+function Base.promote_rule(::Type{Bool}, ::Type{F}) where {F<:FixedRational}
+    return F
+end
+function Base.promote_rule(::Type{F}, ::Type{BigFloat}) where {F<:FixedRational}
+    return promote_type(Rational{eltype(F)}, BigFloat)
+end
+function Base.promote_rule(::Type{BigFloat}, ::Type{F}) where {F<:FixedRational}
+    return promote_type(Rational{eltype(F)}, BigFloat)
+end
+function Base.promote_rule(::Type{F}, ::Type{T}) where {F<:FixedRational,T<:AbstractIrrational}
+    return promote_type(Rational{eltype(F)}, T)
+end
+function Base.promote_rule(::Type{T}, ::Type{F}) where {F<:FixedRational,T<:AbstractIrrational}
+    return promote_type(Rational{eltype(F)}, T)
+end
+
+################################################################################
+# Assorted calls found by Aqua: ################################################
+################################################################################
+
+for type in (Signed, Float64, Float32, Rational), op in (:flipsign, :copysign)
+    @eval function Base.$(op)(x::$type, y::AbstractRealQuantity)
+        return $(op)(x, ustrip(y))
+    end
+end
+for type in (:(Complex), :(Complex{Bool}))
+    @eval begin
+        function Base.:*(l::$type, r::AbstractRealQuantity)
+            new_quantity(typeof(r), l * ustrip(r), dimension(r))
+        end
+        function Base.:*(l::AbstractRealQuantity, r::$type)
+            new_quantity(typeof(l), ustrip(l) * r, dimension(l))
+        end
+    end
+end
+function Complex{T}(q::AbstractRealQuantity) where {T<:Real}
+    @assert iszero(dimension(q)) "$(typeof(q)): $(q) has dimensions! Use `ustrip` instead."
+    return Complex{T}(ustrip(q))
+end
+for type in (:Bool, :Complex)
+    @eval function $type(q::AbstractRealQuantity)
+        @assert iszero(dimension(q)) "$(typeof(q)): $(q) has dimensions! Use `ustrip` instead."
+        return $type(ustrip(q))
+    end
+end
+function Base.:/(l::Complex, r::AbstractRealQuantity)
+    new_quantity(typeof(r), l / ustrip(r), inv(dimension(r)))
+end
+function Base.:/(l::AbstractRealQuantity, r::Complex)
+    new_quantity(typeof(l), ustrip(l) / r, dimension(l))
+end
+for op in (:(==), :isequal), base_type in (AbstractIrrational, AbstractFloat)
+    @eval begin
+        function Base.$(op)(l::AbstractRealQuantity, r::$base_type)
+            return $(op)(ustrip(l), r) && iszero(dimension(l))
+        end
+        function Base.$(op)(l::$base_type, r::AbstractRealQuantity)
+            return $(op)(l, ustrip(r)) && iszero(dimension(r))
+        end
+    end
+end
+function Base.isless(l::AbstractRealQuantity, r::AbstractFloat)
+    iszero(dimension(l)) || throw(DimensionError(l, r))
+    return isless(ustrip(l), r)
+end
+function Base.isless(l::AbstractFloat, r::AbstractRealQuantity)
+    iszero(dimension(r)) || throw(DimensionError(l, r))
+    return isless(l, ustrip(r))
+end
+for (type, _, _) in ABSTRACT_QUANTITY_TYPES, numeric_type in (Bool, BigFloat)
+    @eval begin
+        function Base.promote_rule(::Type{Q}, ::Type{$numeric_type}) where {T,D,Q<:$type{T,D}}
+            return with_type_parameters(promote_quantity_on_value(Q, $numeric_type), promote_type(T, $numeric_type), D)
+        end
+        function Base.promote_rule(::Type{$numeric_type}, ::Type{Q}) where {T,D,Q<:$type{T,D}}
+            return with_type_parameters(promote_quantity_on_value(Q, $numeric_type), promote_type(T, $numeric_type), D)
+        end
+    end
+end
diff --git a/src/fixed_rational.jl b/src/fixed_rational.jl
index a6a3313e..7b50f616 100644
--- a/src/fixed_rational.jl
+++ b/src/fixed_rational.jl
@@ -14,6 +14,7 @@ struct FixedRational{T<:Integer,den} <: Real
     num::T
     global unsafe_fixed_rational(num::Integer, ::Type{T}, ::Val{den}) where {T,den} = new{T,den}(num)
 end
+@inline _denom(::Type{F}) where {T,den,F<:FixedRational{T,den}} = den
 
 """
     denom(F::FixedRational)
@@ -21,14 +22,14 @@ end
 Since `den` can be a different type than `T`, this function
 is used to get the denominator as a `T`.
 """
-denom(::Type{F}) where {T,den,F<:FixedRational{T,den}} = convert(T, den)
+denom(::Type{<:F}) where {T,F<:FixedRational{T}} = convert(T, _denom(F))
 denom(x::FixedRational) = denom(typeof(x))
 
 # But, for Val(den), we need to use the same type as at init.
 # Otherwise, we would have type instability.
-val_denom(::Type{F}) where {T,den,F<:FixedRational{T,den}} = Val(den)
+val_denom(::Type{<:F}) where {F<:FixedRational} = Val(_denom(F))
 
-Base.eltype(::Type{F}) where {T,F<:FixedRational{T}} = T
+Base.eltype(::Type{<:FixedRational{T}}) where {T} = T
 
 const DEFAULT_NUMERATOR_TYPE = Int32
 const DEFAULT_DENOM = DEFAULT_NUMERATOR_TYPE(2^4 * 3^2 * 5^2 * 7)
@@ -73,23 +74,34 @@ Rational(x::F) where {F<:FixedRational} = Rational{eltype(F)}(x)
 Base.round(::Type{T}, x::F, r::RoundingMode=RoundNearest) where {T,F<:FixedRational} = div(convert(T, x.num), convert(T, denom(F)), r)
 Base.decompose(x::F) where {T,F<:FixedRational{T}} = (x.num, zero(T), denom(F))
 
-# Promotion rules:
-function Base.promote_rule(::Type{<:FixedRational{T1,den1}}, ::Type{<:FixedRational{T2,den2}}) where {T1,T2,den1,den2}
-    return error("Refusing to promote `FixedRational` types with mixed denominators. Use `Rational` instead.")
+# Promotion with self or rational-like
+function Base.promote_rule(::Type{F1}, ::Type{F2}) where {F1<:FixedRational,F2<:FixedRational}
+    _denom(F1) == _denom(F2) ||
+        error("Refusing to promote `FixedRational` types with mixed denominators. Use `Rational` instead.")
+    return FixedRational{promote_type(eltype(F1), eltype(F2)), _denom(F1)}
 end
-function Base.promote_rule(::Type{<:FixedRational{T1,den}}, ::Type{<:FixedRational{T2,den}}) where {T1,T2,den}
-    return FixedRational{promote_type(T1,T2),den}
+function Base.promote_rule(::Type{F}, ::Type{Rational{T2}}) where {F<:FixedRational,T2}
+    return Rational{promote_type(eltype(F),T2)}
 end
-function Base.promote_rule(::Type{<:FixedRational{T1}}, ::Type{Rational{T2}}) where {T1,T2}
-    return Rational{promote_type(T1,T2)}
-end
-function Base.promote_rule(::Type{<:FixedRational{T1}}, ::Type{T2}) where {T1,T2<:Real}
-    return promote_type(Rational{T1}, T2)
+function Base.promote_rule(::Type{Rational{T2}}, ::Type{F}) where {F<:FixedRational,T2}
+    return Rational{promote_type(eltype(F),T2)}
 end
+
+# We want to consume integers
 function Base.promote_rule(::Type{F}, ::Type{<:Integer}) where {F<:FixedRational}
-    # Want to consume integers:
     return F
 end
+function Base.promote_rule(::Type{<:Integer}, ::Type{F}) where {F<:FixedRational}
+    return F
+end
+
+# Promotion with general types promotes like a rational
+function Base.promote_rule(::Type{T}, ::Type{T2}) where {T2<:Real,T<:FixedRational}
+    return promote_type(Rational{eltype(T)}, T2)
+end
+function Base.promote_rule(::Type{T2}, ::Type{T}) where {T2<:Real,T<:FixedRational}
+    return promote_type(Rational{eltype(T)}, T2)
+end
 
 Base.string(x::FixedRational) =
     let
diff --git a/src/math.jl b/src/math.jl
index 3fb8e03b..67608330 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -1,4 +1,6 @@
 for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
+    # For div, we don't want to go more generic than `Number`
+    div_base_type = base_type <: Number ? base_type : Number
     @eval begin
         function Base.:*(l::$type, r::$type)
             l, r = promote_except_value(l, r)
@@ -20,7 +22,7 @@ for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
         function Base.:/(l::$type, r::$base_type)
             new_quantity(typeof(l), ustrip(l) / r, dimension(l))
         end
-        function Base.div(x::$type, y::Number, r::RoundingMode=RoundToZero)
+        function Base.div(x::$type, y::$div_base_type, r::RoundingMode=RoundToZero)
             new_quantity(typeof(x), div(ustrip(x), y, r), dimension(x))
         end
 
@@ -30,7 +32,7 @@ for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
         function Base.:/(l::$base_type, r::$type)
             new_quantity(typeof(r), l / ustrip(r), inv(dimension(r)))
         end
-        function Base.div(x::Number, y::$type, r::RoundingMode=RoundToZero)
+        function Base.div(x::$div_base_type, y::$type, r::RoundingMode=RoundToZero)
             new_quantity(typeof(y), div(x, ustrip(y), r), inv(dimension(y)))
         end
 
@@ -53,8 +55,10 @@ end
 Base.:*(l::AbstractDimensions, r::AbstractDimensions) = map_dimensions(+, l, r)
 Base.:/(l::AbstractDimensions, r::AbstractDimensions) = map_dimensions(-, l, r)
 
-# Defines + and -
-for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES, op in (:+, :-)
+# Defines +, -, and mod
+for (type, true_base_type, _) in ABSTRACT_QUANTITY_TYPES, op in (:+, :-, :mod)
+    # Only define `mod` on `Number` types:
+    base_type = (op == :mod && !(true_base_type <: Number)) ? Number : true_base_type
     @eval begin
         function Base.$op(l::$type, r::$type)
             l, r = promote_except_value(l, r)
@@ -75,13 +79,17 @@ end
 Base.:-(l::UnionAbstractQuantity) = new_quantity(typeof(l), -ustrip(l), dimension(l))
 
 # Combining different abstract types
-for op in (:*, :/, :+, :-, :div, :atan, :atand, :copysign, :flipsign, :mod),
+for op in (:*, :/, :+, :-, :atan, :atand, :copysign, :flipsign, :div, :mod),
     (t1, _, _) in ABSTRACT_QUANTITY_TYPES,
     (t2, _, _) in ABSTRACT_QUANTITY_TYPES
 
     t1 == t2 && continue
 
-    @eval Base.$op(l::$t1, r::$t2) = $op(promote_except_value(l, r)...)
+    if op == :div
+        @eval Base.$op(x::$t1, y::$t2, r::RoundingMode=RoundToZero) = $op(promote_except_value(x, y)..., r)
+    else
+        @eval Base.$op(l::$t1, r::$t2) = $op(promote_except_value(l, r)...)
+    end
 end
 
 # We don't promote on the dimension types:
@@ -118,6 +126,10 @@ for (type, _, _) in ABSTRACT_QUANTITY_TYPES
         Base.:^(l::$type, r::Integer) = _pow_int(l, r)
         Base.:^(l::$type, r::Number) = _pow(l, r)
         Base.:^(l::$type, r::Rational) = _pow(l, r)
+        function Base.:^(l::$type, r::Complex)
+            iszero(dimension(l)) || throw(DimensionError(l))
+            return new_quantity(typeof(l), ustrip(l)^r, dimension(l))
+        end
     end
 end
 @inline Base.literal_pow(::typeof(^), l::AbstractDimensions, ::Val{p}) where {p} = map_dimensions(Base.Fix1(*, p), l)
@@ -175,13 +187,13 @@ end
 ############################## Same dimension as input ##################################
 for f in (
     :float, :abs, :real, :imag, :conj, :adjoint, :unsigned,
-    :nextfloat, :prevfloat, :identity, :transpose, :significand
+    :nextfloat, :prevfloat, :transpose, :significand
 )
     @eval function Base.$f(q::UnionAbstractQuantity)
         return new_quantity(typeof(q), $f(ustrip(q)), dimension(q))
     end
 end
-for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES, f in (:copysign, :flipsign, :mod)
+for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES, f in (:copysign, :flipsign,)
     # These treat the x as the magnitude, so we take the dimensions from there,
     # and ignore any dimensions on y, since those will cancel out.
     @eval begin
@@ -197,6 +209,33 @@ for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES, f in (:copysign, :flipsign,
         end
     end
 end
+# Define :rem (unfortunately we have to create a method for each rounding mode to avoid ambiguity)
+for (type, true_base_type, _) in ABSTRACT_QUANTITY_TYPES, rounding_mode in (RoundingMode, RoundingMode{:ToZero}, RoundingMode{:Down}, RoundingMode{:Up}, RoundingMode{:FromZero})
+
+    # We don't want to go more generic than `Number` for mod and rem
+    base_type = true_base_type <: Number ? true_base_type : Number
+    # Add extra args:
+    param = rounding_mode === RoundingMode ? (()) : (:(::$rounding_mode),)
+    extra_f_args = rounding_mode === RoundingMode ? (:RoundToZero,) : (:($rounding_mode()),)
+
+    for (type2, _, _) in ABSTRACT_QUANTITY_TYPES
+        @eval function Base.rem(x::$type, y::$type2, $(param...))
+            x, y = promote_except_value(x, y)
+            dimension(x) == dimension(y) || throw(DimensionError(x, y))
+            return new_quantity(typeof(x), rem(ustrip(x), ustrip(y), $(extra_f_args...)), dimension(x))
+        end
+    end
+    @eval begin
+        function Base.rem(x::$type, y::$base_type, $(param...))
+            iszero(dimension(x)) || throw(DimensionError(x))
+            return new_quantity(typeof(x), rem(ustrip(x), y, $(extra_f_args...)), dimension(x))
+        end
+        function Base.rem(x::$base_type, y::$type, $(param...))
+            iszero(dimension(y)) || throw(DimensionError(y))
+            return new_quantity(typeof(y), rem(x, ustrip(y), $(extra_f_args...)), dimension(y))
+        end
+    end
+end
 function Base.ldexp(x::UnionAbstractQuantity, n::Integer)
     return new_quantity(typeof(x), ldexp(ustrip(x), n), dimension(x))
 end
diff --git a/src/symbolic_dimensions.jl b/src/symbolic_dimensions.jl
index 8638f705..46435978 100644
--- a/src/symbolic_dimensions.jl
+++ b/src/symbolic_dimensions.jl
@@ -104,9 +104,9 @@ end
     uexpand(q::UnionAbstractQuantity{<:Any,<:SymbolicDimensions})
 
 Expand the symbolic units in a quantity to their base SI form.
-In other words, this converts a `Quantity` with `SymbolicDimensions`
+In other words, this converts a quantity with `SymbolicDimensions`
 to one with `Dimensions`. The opposite of this function is `uconvert`,
-for converting to specific symbolic units, or `convert(Quantity{<:Any,<:SymbolicDimensions}, q)`,
+for converting to specific symbolic units, or, e.g., `convert(Quantity{<:Any,<:SymbolicDimensions}, q)`,
 for assuming SI units as the output symbols.
 """
 function uexpand(q::Q) where {T,R,D<:SymbolicDimensions{R},Q<:UnionAbstractQuantity{T,D}}
@@ -264,6 +264,8 @@ function map_dimensions(op::O, l::SymbolicDimensions{L}, r::SymbolicDimensions{R
     return SymbolicDimensions(I, V)
 end
 
+const DEFAULT_SYMBOLIC_QUANTITY_TYPE = with_type_parameters(DEFAULT_QUANTITY_TYPE, DEFAULT_VALUE_TYPE, SymbolicDimensions{DEFAULT_DIM_BASE_TYPE})
+
 """
     SymbolicUnitsParse
 
@@ -277,7 +279,8 @@ module SymbolicUnitsParse
     import ..SYMBOL_CONFLICTS
     import ..SymbolicDimensions
 
-    import ...Quantity
+    import ...constructorof
+    import ...DEFAULT_SYMBOLIC_QUANTITY_TYPE
     import ...DEFAULT_VALUE_TYPE
     import ...DEFAULT_DIM_BASE_TYPE
 
@@ -287,7 +290,8 @@ module SymbolicUnitsParse
         import ..SYMBOL_CONFLICTS
         import ..SymbolicDimensions
 
-        import ..Quantity
+        import ..constructorof
+        import ..DEFAULT_SYMBOLIC_QUANTITY_TYPE
         import ..DEFAULT_VALUE_TYPE
         import ..DEFAULT_DIM_BASE_TYPE
 
@@ -299,11 +303,11 @@ module SymbolicUnitsParse
             CONSTANT_SYMBOLS_EXIST[] || lock(CONSTANT_SYMBOLS_LOCK) do
                 CONSTANT_SYMBOLS_EXIST[] && return nothing
                 for unit in setdiff(CONSTANT_SYMBOLS, SYMBOL_CONFLICTS)
-                    @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+                    @eval const $unit = constructorof(DEFAULT_SYMBOLIC_QUANTITY_TYPE)(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
                 end
                 # Evaluate conflicting symbols to non-symbolic form:
                 for unit in SYMBOL_CONFLICTS
-                    @eval const $unit = convert(Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions}, EagerConstants.$unit)
+                    @eval const $unit = convert(DEFAULT_SYMBOLIC_QUANTITY_TYPE, EagerConstants.$unit)
                 end
                 CONSTANT_SYMBOLS_EXIST[] = true
             end
@@ -318,7 +322,7 @@ module SymbolicUnitsParse
         UNIT_SYMBOLS_EXIST[] || lock(UNIT_SYMBOLS_LOCK) do
             UNIT_SYMBOLS_EXIST[] && return nothing
             for unit in UNIT_SYMBOLS
-                @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+                @eval const $unit = constructorof(DEFAULT_SYMBOLIC_QUANTITY_TYPE)(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
             end
             UNIT_SYMBOLS_EXIST[] = true
         end
@@ -345,11 +349,11 @@ module SymbolicUnitsParse
         _generate_unit_symbols()
         Constants._generate_unit_symbols()
         raw_result = eval(Meta.parse(raw_string))
-        return copy(as_quantity(raw_result))::Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}}
+        return copy(as_quantity(raw_result))::DEFAULT_SYMBOLIC_QUANTITY_TYPE
     end
 
-    as_quantity(q::Quantity) = q
-    as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE})
+    as_quantity(q::DEFAULT_SYMBOLIC_QUANTITY_TYPE) = q
+    as_quantity(x::Number) = convert(DEFAULT_SYMBOLIC_QUANTITY_TYPE, x)
     as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
 end
 
@@ -380,3 +384,6 @@ end
 function Base.promote_rule(::Type{SymbolicDimensions{R1}}, ::Type{Dimensions{R2}}) where {R1,R2}
     return Dimensions{promote_type(R1,R2)}
 end
+function Base.promote_rule(::Type{Dimensions{R2}}, ::Type{SymbolicDimensions{R1}}) where {R1,R2}
+    return Dimensions{promote_type(R1,R2)}
+end
diff --git a/src/types.jl b/src/types.jl
index 116dd883..68f519b0 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -56,6 +56,14 @@ _as well as any other future abstract quantity types_,
 """
 abstract type AbstractGenericQuantity{T,D} end
 
+"""
+    AbstractRealQuantity{T,D} <: Real
+
+This has the same behavior as `AbstractQuantity` but is subtyped to `Real` rather
+than `Number`.
+"""
+abstract type AbstractRealQuantity{T,D} <: Real end
+
 """
     UnionAbstractQuantity{T,D}
 
@@ -64,7 +72,7 @@ It is used throughout the library to declare methods which can take both types.
 You should generally specialize on this type, rather than its constituents,
 as it will also include future abstract quantity types.
 """
-const UnionAbstractQuantity{T,D} = Union{AbstractQuantity{T,D},AbstractGenericQuantity{T,D}}
+const UnionAbstractQuantity{T,D} = Union{AbstractQuantity{T,D},AbstractGenericQuantity{T,D},AbstractRealQuantity{T,D}}
 
 """
     Dimensions{R<:Real} <: AbstractDimensions{R}
@@ -164,6 +172,18 @@ struct GenericQuantity{T,D<:AbstractDimensions} <: AbstractGenericQuantity{T,D}
     GenericQuantity(x::_T, dimensions::_D) where {_T,_D<:AbstractDimensions} = new{_T,_D}(x, dimensions)
 end
 
+"""
+    RealQuantity{T<:Real,D<:AbstractDimensions} <: AbstractRealQuantity{T,D} <: Real
+
+This has the same behavior as `Quantity` but is subtyped to `AbstractRealQuantity <: Real`.
+"""
+struct RealQuantity{T<:Real,D<:AbstractDimensions} <: AbstractRealQuantity{T,D}
+    value::T
+    dimensions::D
+
+    RealQuantity(x::_T, dimensions::_D) where {_T,_D<:AbstractDimensions} = new{_T,_D}(x, dimensions)
+end
+
 """
     ABSTRACT_QUANTITY_TYPES
 
@@ -171,7 +191,12 @@ A constant tuple of the existing abstract quantity types,
 each as a tuple with (1) the abstract type,
 (2) the base type, and (3) the default exported concrete type.
 """
-const ABSTRACT_QUANTITY_TYPES = ((AbstractQuantity, Number, Quantity), (AbstractGenericQuantity, Any, GenericQuantity))
+const ABSTRACT_QUANTITY_TYPES = (
+    (AbstractQuantity, Number, Quantity),
+    (AbstractGenericQuantity, Any, GenericQuantity),
+    (AbstractRealQuantity, Real, RealQuantity)
+)
+
 
 for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
     @eval begin
@@ -191,8 +216,13 @@ end
 
 const DEFAULT_QUANTITY_TYPE = Quantity{DEFAULT_VALUE_TYPE, DEFAULT_DIM_TYPE}
 
-new_dimensions(::Type{D}, dims...) where {D<:AbstractDimensions} = constructorof(D)(dims...)
-new_quantity(::Type{Q}, l, r) where {Q<:UnionAbstractQuantity} = constructorof(Q)(l, r)
+@inline function new_dimensions(::Type{D}, dims...) where {D<:AbstractDimensions}
+    return constructorof(D)(dims...)
+end
+@inline function new_quantity(::Type{Q}, val, dims) where {Q<:UnionAbstractQuantity}
+    Qout = promote_quantity_on_value(Q, typeof(val))
+    return constructorof(Qout)(val, dims)
+end
 
 dim_type(::Type{Q}) where {T,D<:AbstractDimensions,Q<:UnionAbstractQuantity{T,D}} = D
 dim_type(::Type{<:UnionAbstractQuantity}) = DEFAULT_DIM_TYPE
@@ -208,6 +238,7 @@ if you need custom behavior.
 constructorof(::Type{<:Dimensions}) = Dimensions
 constructorof(::Type{<:Quantity}) = Quantity
 constructorof(::Type{<:GenericQuantity}) = GenericQuantity
+constructorof(::Type{<:RealQuantity}) = RealQuantity
 
 """
     with_type_parameters(::Type{<:AbstractDimensions}, ::Type{R})
@@ -226,6 +257,9 @@ end
 function with_type_parameters(::Type{<:GenericQuantity}, ::Type{T}, ::Type{D}) where {T,D}
     return GenericQuantity{T,D}
 end
+function with_type_parameters(::Type{<:RealQuantity}, ::Type{T}, ::Type{D}) where {T,D}
+    return RealQuantity{T,D}
+end
 
 # The following functions should be overloaded for special types
 function constructorof(::Type{T}) where {T<:Union{UnionAbstractQuantity,AbstractDimensions}}
diff --git a/src/units.jl b/src/units.jl
index 4b8b831e..81d0d3a5 100644
--- a/src/units.jl
+++ b/src/units.jl
@@ -3,7 +3,6 @@ module Units
 import ..DEFAULT_DIM_TYPE
 import ..DEFAULT_VALUE_TYPE
 import ..DEFAULT_QUANTITY_TYPE
-import ..Quantity
 
 @assert DEFAULT_VALUE_TYPE == Float64 "`units.jl` must be updated to support a different default value type."
 
diff --git a/src/uparse.jl b/src/uparse.jl
index 68ee5b52..e913d044 100644
--- a/src/uparse.jl
+++ b/src/uparse.jl
@@ -1,6 +1,7 @@
 module UnitsParse
 
-import ..Quantity
+import ..constructorof
+import ..DEFAULT_QUANTITY_TYPE
 import ..DEFAULT_DIM_TYPE
 import ..DEFAULT_VALUE_TYPE
 import ..Units: UNIT_SYMBOLS
@@ -33,11 +34,11 @@ the quantity corresponding to the speed of light multiplied by Hertz,
 squared.
 """
 function uparse(s::AbstractString)
-    return as_quantity(eval(Meta.parse(s)))::Quantity{DEFAULT_VALUE_TYPE,DEFAULT_DIM_TYPE}
+    return as_quantity(eval(Meta.parse(s)))::DEFAULT_QUANTITY_TYPE
 end
 
-as_quantity(q::Quantity) = q
-as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), DEFAULT_DIM_TYPE)
+as_quantity(q::DEFAULT_QUANTITY_TYPE) = q
+as_quantity(x::Number) = convert(DEFAULT_QUANTITY_TYPE,  x)
 as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
 
 """
diff --git a/src/utils.jl b/src/utils.jl
index d31791a4..2cb32ac1 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -25,24 +25,55 @@ end
     return output
 end
 
-Base.convert(::Type{Number}, q::AbstractQuantity) = q
-function Base.convert(::Type{T}, q::UnionAbstractQuantity) where {T<:Number}
-    @assert iszero(dimension(q)) "$(typeof(q)): $(q) has dimensions! Use `ustrip` instead."
-    return convert(T, ustrip(q))
-end
 function Base.promote_rule(::Type{Dimensions{R1}}, ::Type{Dimensions{R2}}) where {R1,R2}
     return Dimensions{promote_type(R1,R2)}
 end
-function Base.promote_rule(::Type{<:GenericQuantity{T1,D1}}, ::Type{<:GenericQuantity{T2,D2}}) where {T1,T2,D1,D2}
-    return GenericQuantity{promote_type(T1,T2),promote_type(D1,D2)}
-end
-function Base.promote_rule(::Type{<:Quantity{T1,D1}}, ::Type{<:GenericQuantity{T2,D2}}) where {T1,T2,D1,D2}
-    return GenericQuantity{promote_type(T1,T2),promote_type(D1,D2)}
-end
-function Base.promote_rule(::Type{<:Quantity{T1,D1}}, ::Type{<:Quantity{T2,D2}}) where {T1,T2,D1,D2}
-    return Quantity{promote_type(T1,T2),promote_type(D1,D2)}
+
+# Define all the quantity x quantity promotion rules
+"""
+    promote_quantity_on_value(Q::Type, T::Type)
+
+Find the next quantity type in the hierarchy that can accommodate the type `T`.
+If the current quantity type can already accommodate `T`, then the current type is returned.
+For example, `promote_quantity_on_value(Quantity, Float64)` would return `Quantity`, and
+`promote_quantity_on_value(RealQuantity, String)` would return `GenericQuantity`.
+The user should overload this function to define a custom type hierarchy.
+
+Also see `promote_quantity_on_quantity`.
+"""
+@inline promote_quantity_on_value(::Type{<:Union{GenericQuantity,Quantity,RealQuantity}}, ::Type{<:Any}) = GenericQuantity
+@inline promote_quantity_on_value(::Type{<:Union{Quantity,RealQuantity}}, ::Type{<:Number}) = Quantity
+@inline promote_quantity_on_value(::Type{<:RealQuantity}, ::Type{<:Real}) = RealQuantity
+@inline promote_quantity_on_value(T, _) = T
+
+"""
+    promote_quantity_on_quantity(Q1, Q2)
+
+Defines the type hierarchy for quantities, returning the most specific type
+that is compatible with both input quantity types. For example,
+`promote_quantity_on_quantity(Quantity, GenericQuantity)` would return `GenericQuantity`,
+as it can store both `Quantity` and `GenericQuantity` values.
+Similarly, `promote_quantity_on_quantity(RealQuantity, RealQuantity)` would return `RealQuantity`,
+as that is the most specific type.
+
+Also see `promote_quantity_on_value`.
+"""
+@inline promote_quantity_on_quantity(::Type{<:Union{GenericQuantity,Quantity,RealQuantity}}, ::Type{<:Union{GenericQuantity,Quantity,RealQuantity}}) = GenericQuantity
+@inline promote_quantity_on_quantity(::Type{<:Union{Quantity,RealQuantity}}, ::Type{<:Union{Quantity,RealQuantity}}) = Quantity
+@inline promote_quantity_on_quantity(::Type{<:RealQuantity}, ::Type{<:RealQuantity}) = RealQuantity
+@inline promote_quantity_on_quantity(::Type{Q}, ::Type{Q}) where {Q<:UnionAbstractQuantity} = Q
+
+for (type1, _, _) in ABSTRACT_QUANTITY_TYPES, (type2, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval function Base.promote_rule(::Type{Q1}, ::Type{Q2}) where {T1,T2,D1,D2,Q1<:$type1{T1,D1},Q2<:$type2{T2,D2}}
+        return with_type_parameters(
+            promote_quantity_on_quantity(Q1, Q2),
+            promote_type(T1, T2),
+            promote_type(D1, D2),
+        )
+    end
 end
 
+
 # Define promotion rules for all basic numeric types, individually.
 # We don't want to define an opinionated promotion on <:Number,
 # or even <:AbstractFloat, as it could conflict with other
@@ -57,16 +88,31 @@ const BASE_NUMERIC_TYPES = Union{
     Rational{Int64}, Rational{UInt64}, Rational{Int128}, Rational{UInt128},
     Rational{BigInt},
 }
-for (type, _, _) in ABSTRACT_QUANTITY_TYPES
-    @eval function Base.promote_rule(::Type{Q}, ::Type{T2}) where {T,D,Q<:$type{T,D},T2<:BASE_NUMERIC_TYPES}
-        return with_type_parameters(Q, promote_type(T, T2), D)
+
+for (type, base_type, _) in ABSTRACT_QUANTITY_TYPES
+    !(base_type <: Number) && @eval begin
+        function Base.convert(::Type{Q}, x::BASE_NUMERIC_TYPES) where {T,D,Q<:$type{T,D}}
+            return new_quantity(Q, convert(T, x), D())
+        end
     end
-    @eval function Base.convert(::Type{Q}, x::BASE_NUMERIC_TYPES) where {T,D,Q<:$type{T,D}}
-        return new_quantity(Q, convert(T, x), D())
+    @eval begin
+        function Base.promote_rule(::Type{Q}, ::Type{T2}) where {T,D,Q<:$type{T,D},T2<:BASE_NUMERIC_TYPES}
+            return with_type_parameters(promote_quantity_on_value(Q, T2), promote_type(T, T2), D)
+        end
+        function Base.promote_rule(::Type{T2}, ::Type{Q}) where {T,D,Q<:$type{T,D},T2<:BASE_NUMERIC_TYPES}
+            return with_type_parameters(promote_quantity_on_value(Q, T2), promote_type(T, T2), D)
+        end
     end
 end
-function Base.promote_rule(::Type{<:AbstractQuantity}, ::Type{<:Number})
-    return Number
+
+for (type, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        function (::Type{T})(q::$type) where {T<:Number}
+            q isa T && return q
+            @assert iszero(dimension(q)) "$(typeof(q)): $(q) has dimensions! Use `ustrip` instead."
+            return convert(T, ustrip(q))
+        end
+    end
 end
 
 """
@@ -118,47 +164,89 @@ Base.keys(q::UnionAbstractQuantity) = keys(ustrip(q))
 
 
 # Numeric checks
-function Base.isapprox(l::UnionAbstractQuantity, r::UnionAbstractQuantity; kws...)
-    l, r = promote_except_value(l, r)
-    return isapprox(ustrip(l), ustrip(r); kws...) && dimension(l) == dimension(r)
-end
-function Base.isapprox(l::Number, r::UnionAbstractQuantity; kws...)
-    iszero(dimension(r)) || throw(DimensionError(l, r))
-    return isapprox(l, ustrip(r); kws...)
-end
-function Base.isapprox(l::UnionAbstractQuantity, r::Number; kws...)
-    iszero(dimension(l)) || throw(DimensionError(l, r))
-    return isapprox(ustrip(l), r; kws...)
-end
-Base.iszero(d::AbstractDimensions) = all_dimensions(iszero, d)
-function Base.:(==)(l::UnionAbstractQuantity, r::UnionAbstractQuantity)
-    l, r = promote_except_value(l, r)
-    ustrip(l) == ustrip(r) && dimension(l) == dimension(r)
+for op in (:(<=), :(<), :(>=), :(>), :isless), (type, true_base_type, _) in ABSTRACT_QUANTITY_TYPES
+    # Avoid creating overly generic operations on these:
+    base_type = true_base_type <: Number ? true_base_type : Number
+    @eval begin
+        function Base.$(op)(l::$type, r::$type)
+            l, r = promote_except_value(l, r)
+            dimension(l) == dimension(r) || throw(DimensionError(l, r))
+            return $(op)(ustrip(l), ustrip(r))
+        end
+        function Base.$(op)(l::$type, r::$base_type)
+            iszero(dimension(l)) || throw(DimensionError(l, r))
+            return $(op)(ustrip(l), r)
+        end
+        function Base.$(op)(l::$base_type, r::$type)
+            iszero(dimension(r)) || throw(DimensionError(l, r))
+            return $(op)(l, ustrip(r))
+        end
+    end
 end
-Base.:(==)(l::Number, r::UnionAbstractQuantity) = ustrip(l) == ustrip(r) && iszero(dimension(r))
-Base.:(==)(l::UnionAbstractQuantity, r::Number) = ustrip(l) == ustrip(r) && iszero(dimension(l))
-Base.:(==)(l::AbstractDimensions, r::AbstractDimensions) = all_dimensions(==, l, r)
-function Base.isless(l::UnionAbstractQuantity, r::UnionAbstractQuantity)
-    l, r = promote_except_value(l, r)
-    dimension(l) == dimension(r) || throw(DimensionError(l, r))
-    return isless(ustrip(l), ustrip(r))
+for op in (:isequal, :(==)), (type, true_base_type, _) in ABSTRACT_QUANTITY_TYPES
+    # Avoid creating overly generic operations on these:
+    base_type = true_base_type <: Number ? true_base_type : Number
+    @eval begin
+        function Base.$(op)(l::$type, r::$type)
+            l, r = promote_except_value(l, r)
+            return $(op)(ustrip(l), ustrip(r)) && dimension(l) == dimension(r)
+        end
+        function Base.$(op)(l::$type, r::$base_type)
+            return $(op)(ustrip(l), r) && iszero(dimension(l))
+        end
+        function Base.$(op)(l::$base_type, r::$type)
+            return $(op)(l, ustrip(r)) && iszero(dimension(r))
+        end
+    end
 end
-function Base.isless(l::UnionAbstractQuantity, r::Number)
-    iszero(dimension(l)) || throw(DimensionError(l, r))
-    return isless(ustrip(l), r)
+for op in (:(<=), :(<), :(>=), :(>), :isless, :isgreater, :isequal, :(==)),
+    (t1, _, _) in ABSTRACT_QUANTITY_TYPES,
+    (t2, _, _) in ABSTRACT_QUANTITY_TYPES
+
+    t1 == t2 && continue
+
+    @eval function Base.$(op)(l::$t1, r::$t2)
+        return $(op)(promote_except_value(l, r)...)
+    end
 end
-function Base.isless(l::Number, r::UnionAbstractQuantity)
-    iszero(dimension(r)) || throw(DimensionError(l, r))
-    return isless(l, ustrip(r))
+# Define isapprox:
+for (type, true_base_type, _) in ABSTRACT_QUANTITY_TYPES
+    # Avoid creating overly generic operations on these:
+    base_type = true_base_type <: Number ? true_base_type : Number
+    @eval begin
+        function Base.isapprox(l::$type, r::$type; kws...)
+            dimension(l) == dimension(r) || throw(DimensionError(l, r))
+            return isapprox(ustrip(l), ustrip(r); kws...)
+        end
+        function Base.isapprox(l::$base_type, r::$type; kws...)
+            iszero(dimension(r)) || throw(DimensionError(l, r))
+            return isapprox(l, ustrip(r); kws...)
+        end
+        function Base.isapprox(l::$type, r::$base_type; kws...)
+            iszero(dimension(l)) || throw(DimensionError(l, r))
+            return isapprox(ustrip(l), r; kws...)
+        end
+    end
+    for (type2, _, _) in ABSTRACT_QUANTITY_TYPES
+
+        type == type2 && continue
+
+        @eval function Base.isapprox(l::$type, r::$type2; kws...)
+            return isapprox(promote_except_value(l, r)...; kws...)
+        end
+    end
 end
 
+
 # Simple flags:
 for f in (
-    :iszero, :isfinite, :isinf, :isnan, :isreal, :signbit,
+    :isone, :iszero, :isfinite, :isinf, :isnan, :isreal, :signbit,
     :isempty, :iseven, :isodd, :isinteger, :ispow2
 )
     @eval Base.$f(q::UnionAbstractQuantity) = $f(ustrip(q))
 end
+Base.iszero(d::AbstractDimensions) = all_dimensions(iszero, d)
+Base.:(==)(l::AbstractDimensions, r::AbstractDimensions) = all_dimensions(==, l, r)
 
 
 # Base.one, typemin, typemax
@@ -177,7 +265,7 @@ end
 Base.one(::Type{D}) where {D<:AbstractDimensions} = D()
 Base.one(::D) where {D<:AbstractDimensions} = one(D)
 
-# Additive identities (zero)
+# Additive identities (zero). We have to invalidate these due to different behavior with conversion
 Base.zero(q::Q) where {Q<:UnionAbstractQuantity} = new_quantity(Q, zero(ustrip(q)), dimension(q))
 Base.zero(::AbstractDimensions) = error("There is no such thing as an additive identity for a `AbstractDimensions` object, as + is only defined for `UnionAbstractQuantity`.")
 Base.zero(::Type{<:UnionAbstractQuantity}) = error("Cannot create an additive identity for a `UnionAbstractQuantity` type, as the dimensions are unknown. Please use `zero(::UnionAbstractQuantity)` instead.")
@@ -219,9 +307,17 @@ tryrationalize(::Type{R}, x) where {R} = isinteger(x) ? convert(R, round(Int, x)
 Base.showerror(io::IO, e::DimensionError) = print(io, "DimensionError: ", e.q1, " and ", e.q2, " have incompatible dimensions")
 Base.showerror(io::IO, e::DimensionError{<:Any,Nothing}) = print(io, "DimensionError: ", e.q1, " is not dimensionless")
 
-Base.convert(::Type{Q}, q::UnionAbstractQuantity) where {Q<:UnionAbstractQuantity} = q
-Base.convert(::Type{Q}, q::UnionAbstractQuantity) where {T,Q<:UnionAbstractQuantity{T}} = new_quantity(Q, convert(T, ustrip(q)), dimension(q))
-Base.convert(::Type{Q}, q::UnionAbstractQuantity) where {T,D,Q<:UnionAbstractQuantity{T,D}} = new_quantity(Q, convert(T, ustrip(q)), convert(D, dimension(q)))
+for (type, _, _) in ABSTRACT_QUANTITY_TYPES, (type2, _, _) in ABSTRACT_QUANTITY_TYPES
+    @eval begin
+        Base.convert(::Type{Q}, q::$type) where {Q<:$type2} = q
+        Base.convert(::Type{Q}, q::$type) where {T,Q<:$type2{T}} = new_quantity(Q, convert(T, ustrip(q)), dimension(q))
+        Base.convert(::Type{Q}, q::$type) where {T,D,Q<:$type2{T,D}} = new_quantity(Q, convert(T, ustrip(q)), convert(D, dimension(q)))
+    end
+    # TODO: This invalidates some methods. But we have to, because
+    # the conversion in `number.jl` has a type assertion step, whereas
+    # we want to allow things like `convert(Quantity{Float64}, 1.0u"m")`,
+    # with the type for the dimensions being inferred.
+end
 
 Base.convert(::Type{D}, d::AbstractDimensions) where {D<:AbstractDimensions} = d
 Base.convert(::Type{D}, d::AbstractDimensions) where {R,D<:AbstractDimensions{R}} = D(d)
diff --git a/test/test_measurements.jl b/test/test_measurements.jl
index a02ab642..b76bc809 100644
--- a/test/test_measurements.jl
+++ b/test/test_measurements.jl
@@ -2,7 +2,7 @@ using DynamicQuantities
 using Measurements
 using Measurements: value, uncertainty
 
-for Q in (Quantity, GenericQuantity)
+for Q in (RealQuantity, Quantity, GenericQuantity)
     x = Q(1.0u"m/s") ± Q(0.1u"m/s")
 
     @test ustrip(x^2) == ustrip(x)^2
@@ -11,7 +11,9 @@ for Q in (Quantity, GenericQuantity)
     @test dimension(x)^2 == dimension(x^2)
     @test_throws DimensionError 0.5u"m" ± 0.1u"s"
 
-    # Mixed types:
-    y = Q{Float16}(0.1u"m/s") ± Q{Float32}(0.1u"m/s")
-    @test typeof(y) <: Q{Measurement{Float32}}
+    if Q in (Quantity, GenericQuantity)
+        # Mixed types:
+        y = Q{Float16}(0.1u"m/s") ± Q{Float32}(0.1u"m/s")
+        @test typeof(y) <: Q{Measurement{Float32}}
+    end
 end
diff --git a/test/test_scitypes.jl b/test/test_scitypes.jl
index 5f1897a7..51afe69d 100644
--- a/test/test_scitypes.jl
+++ b/test/test_scitypes.jl
@@ -1,4 +1,5 @@
 using DynamicQuantities
+using DynamicQuantities: DEFAULT_QUANTITY_TYPE, constructorof
 using ScientificTypes
 import ScientificTypes as ST
 
@@ -18,6 +19,6 @@ sch = schema(X)
 
 @test first(sch.names) == :x
 @test first(sch.scitypes) == Continuous
-@test first(sch.types) <: Quantity{Float64}
+@test first(sch.types) <: constructorof(DEFAULT_QUANTITY_TYPE){Float64}
 
 @test first(schema((; x=rand(1:10, 5) .* Quantity{Int}(u"m/s"))).scitypes) == Count
diff --git a/test/test_unitful.jl b/test/test_unitful.jl
index 8da1945e..d0dfddcd 100644
--- a/test/test_unitful.jl
+++ b/test/test_unitful.jl
@@ -6,10 +6,10 @@ import Ratios: SimpleRatio
 import SaferIntegers: SafeInt16
 using Test
 
-risapprox(x::Unitful.Quantity, y::Unitful.Quantity; kws...) =
-    let (xfloat, yfloat) = (Unitful.ustrip ∘ Unitful.upreferred).((x, y))
-        return isapprox(xfloat, yfloat; kws...)
-    end
+function risapprox(x::Unitful.Quantity, y::Unitful.Quantity; kws...)
+    (xfloat, yfloat) = (Unitful.ustrip ∘ Unitful.upreferred).((x, y))
+    return isapprox(xfloat, yfloat; kws...)
+end
 
 for T in [DEFAULT_VALUE_TYPE, Float16, Float32, Float64], R in [DEFAULT_DIM_BASE_TYPE, Rational{Int16}, Rational{Int32}, SimpleRatio{Int}, SimpleRatio{SafeInt16}]
     D = DynamicQuantities.Dimensions{R}
diff --git a/test/unittests.jl b/test/unittests.jl
index f7636d19..e952ac3a 100644
--- a/test/unittests.jl
+++ b/test/unittests.jl
@@ -1,8 +1,9 @@
 using DynamicQuantities
 using DynamicQuantities: FixedRational
-using DynamicQuantities: DEFAULT_DIM_BASE_TYPE, DEFAULT_DIM_TYPE, DEFAULT_VALUE_TYPE
+using DynamicQuantities: DEFAULT_QUANTITY_TYPE, DEFAULT_DIM_BASE_TYPE, DEFAULT_DIM_TYPE, DEFAULT_VALUE_TYPE
 using DynamicQuantities: array_type, value_type, dim_type, quantity_type
-using DynamicQuantities: GenericQuantity
+using DynamicQuantities: GenericQuantity, with_type_parameters, constructorof
+using DynamicQuantities: promote_quantity_on_quantity, promote_quantity_on_value
 using Ratios: SimpleRatio
 using SaferIntegers: SafeInt16
 using StaticArrays: SArray, MArray
@@ -14,10 +15,13 @@ function record_show(s, f=show)
     f(io, s)
     return String(take!(io))
 end
+function unsafe_isapprox(x, y; kwargs...)
+    return isapprox(ustrip(x), ustrip(y); kwargs...) && dimension(x) == dimension(y)
+end
 
 @testset "Basic utilities" begin
 
-    for Q in [Quantity, GenericQuantity], T in [DEFAULT_VALUE_TYPE, Float16, Float32, Float64], R in [DEFAULT_DIM_BASE_TYPE, Rational{Int16}, Rational{Int32}, SimpleRatio{Int}, SimpleRatio{SafeInt16}]
+    for Q in [Quantity, GenericQuantity, RealQuantity], T in [DEFAULT_VALUE_TYPE, Float16, Float32, Float64], R in [DEFAULT_DIM_BASE_TYPE, Rational{Int16}, Rational{Int32}, SimpleRatio{Int}, SimpleRatio{SafeInt16}]
         D = Dimensions{R}
         x = Q(T(0.2), D, length=1, mass=2.5)
 
@@ -94,7 +98,7 @@ end
 
         y = Q(T(2 // 10), D, length=1, mass=6 // 2)
 
-        @test !(y ≈ x)
+        @test !(unsafe_isapprox(y, x))
 
         y = x * Inf32
 
@@ -166,6 +170,8 @@ end
     @test iseven(Quantity(3, length=1)) == false
     @test isodd(Quantity(2, length=1)) == false
     @test isodd(Quantity(3, length=1)) == true
+    @test isone(Quantity(1, length=1)) == true
+    @test isone(Quantity(2, length=1)) == false
     @test isinteger(Quantity(2, length=1)) == true
     @test isinteger(Quantity(2.1, length=1)) == false
     @test ispow2(Quantity(2, length=1)) == true
@@ -183,6 +189,12 @@ end
 
 end
 
+@testset "Ranges" begin
+    x = [xi for xi in 0.0u"km/s":0.1u"km/s":1.0u"km/s"]
+    @test x[2] == 0.1u"km/s"
+    @test x[end] == 1.0u"km/s"
+end
+
 @testset "Complex numbers" begin
     x = (0.5 + 0.6im) * u"km/s"
     @test string(x) == "(500.0 + 600.0im) m s⁻¹"
@@ -197,6 +209,19 @@ end
     @test conj(x) == (0.5 - 0.6im) * u"km/s"
     @test angle(x) == angle(ustrip(x))
     @test adjoint(ustrip(x^2)) ≈ adjoint(x^2) / u"m/s"^2
+
+    # Can create by division as well:
+    x = RealQuantity(1.0u"km/s") / (1.0 + 0.5im)
+    @test typeof(x) == Quantity{Complex{Float64}, DEFAULT_DIM_TYPE}
+    @test ustrip(x) ≈ 1000.0 / (1.0 + 0.5im)
+    @test ulength(x) == 1.0
+    @test utime(x) == -1.0
+
+    x = (1.0 + 0.5im) / RealQuantity(1.0u"km/s")
+    @test typeof(x) == Quantity{Complex{Float64}, DEFAULT_DIM_TYPE}
+    @test ustrip(x) ≈ (1.0 + 0.5im) / 1000.0
+    @test ulength(x) == -1.0
+    @test utime(x) == 1.0
 end
 
 @testset "Fallbacks" begin
@@ -387,18 +412,18 @@ end
     @test utime(x) == -2
 
     y = 0.9u"sqrt(mΩ)"
-    @test typeof(y) == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(y) == with_type_parameters(DEFAULT_QUANTITY_TYPE, Float64, DEFAULT_DIM_TYPE)
     @test ustrip(y) ≈ 0.02846049894151541
     @test ucurrent(y) == -1
     @test ulength(y) == 1
 
     y = BigFloat(0.3) * u"mΩ"
-    @test typeof(y) == Quantity{BigFloat,DEFAULT_DIM_TYPE}
+    @test typeof(y) == with_type_parameters(DEFAULT_QUANTITY_TYPE, BigFloat, DEFAULT_DIM_TYPE)
     @test ustrip(y) ≈ 0.0003
     @test ulength(y) == 2
 
-    y32 = convert(Quantity{Float32,Dimensions{Rational{Int16}}}, y)
-    @test typeof(y32) == Quantity{Float32,Dimensions{Rational{Int16}}}
+    y32 = convert(with_type_parameters(DEFAULT_QUANTITY_TYPE, Float32, Dimensions{Rational{Int16}}), y)
+    @test typeof(y32) == with_type_parameters(DEFAULT_QUANTITY_TYPE, Float32, Dimensions{Rational{Int16}})
     @test ustrip(y32) ≈ 0.0003
 
     z = u"yr"
@@ -406,13 +431,13 @@ end
     @test ustrip(z) ≈ 60 * 60 * 24 * 365.25
 
     # Test type stability of extreme range of units
-    @test typeof(u"1") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"1f0") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"s"^2) == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"Ω") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"Gyr") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"fm") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"fm"^2) == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(u"1") == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"1f0") == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"s"^2) == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"Ω") == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"Gyr") == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"fm") == DEFAULT_QUANTITY_TYPE
+    @test typeof(u"fm"^2) == DEFAULT_QUANTITY_TYPE
 
     @test_throws LoadError eval(:(u":x"))
 end
@@ -480,13 +505,13 @@ end
     @eval struct MyNumber <: Real
         x::Float64
     end
-    a = 0.5u"km/s"
+    a = RealQuantity(0.5u"km/s")
     b = MyNumber(0.5)
     ar = [a, b]
-    @test ar isa Vector{Number}
+    @test ar isa Vector{Real}
     @test a === ar[1]
     @test b === ar[2]
-    @test promote_type(MyNumber, typeof(a)) == Number
+    @test promote_type(MyNumber, typeof(a)) == Real
 
     # Explicit conversion so coverage can see it:
     D = DEFAULT_DIM_TYPE
@@ -496,6 +521,8 @@ end
     @test promote_type(GenericQuantity{Float32,D}, GenericQuantity{Float64,D}) == GenericQuantity{Float64,D}
     @test promote_type(SymbolicDimensions{Rational{Int}}, SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}) == SymbolicDimensions{Rational{Int}}
     @test promote_type(Dimensions{Rational{Int}}, SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}) == Dimensions{Rational{Int}}
+
+    @test promote_quantity_on_quantity(RealQuantity, RealQuantity) == RealQuantity
 end
 
 struct MyDimensions{R} <: AbstractDimensions{R}
@@ -558,6 +585,9 @@ end
 
     # But, we always need to use a quantity when mixing with mathematical operations:
     @test_throws ErrorException MyQuantity(0.1) + 0.1 * MyDimensions()
+
+    # Explicitly test that `promote_quantity_on_quantity` has a reasonable default
+    @test promote_quantity_on_quantity(typeof(MyQuantity(0.1)), typeof(MyQuantity(0.1))) == MyQuantity{Float64,DEFAULT_DIM_TYPE}
 end
 
 @testset "Symbolic dimensions" begin
@@ -606,7 +636,7 @@ end
 
     q = 1.5us"km/s"
     @test q == 1.5 * us"km" / us"s"
-    @test typeof(q) <: Quantity{Float64,<:SymbolicDimensions}
+    @test typeof(q) <: with_type_parameters(DEFAULT_QUANTITY_TYPE, Float64, SymbolicDimensions{DEFAULT_DIM_BASE_TYPE})
     @test string(dimension(q)) == "s⁻¹ km"
     @test uexpand(q) == 1.5u"km/s"
     @test string(dimension(us"Constants.au^1.5")) == "au³ᐟ²"
@@ -680,13 +710,13 @@ end
     @test_throws DimensionError uconvert(us"nm * J", 5e-9u"m")
 
     # Types:
-    @test typeof(uconvert(us"nm", 5e-9u"m")) <: Quantity{Float64,<:SymbolicDimensions}
+    @test typeof(uconvert(us"nm", 5e-9u"m")) <: constructorof(DEFAULT_QUANTITY_TYPE){Float64,<:SymbolicDimensions}
     @test typeof(uconvert(us"nm", GenericQuantity(5e-9u"m"))) <: GenericQuantity{Float64,<:SymbolicDimensions}
     @test uconvert(GenericQuantity(us"nm"), GenericQuantity(5e-9u"m")) ≈ 5us"nm"
     @test uconvert(GenericQuantity(us"nm"), GenericQuantity(5e-9u"m")) ≈ GenericQuantity(5us"nm")
 
     # We only want to convert the dimensions, and ignore the quantity type:
-    @test typeof(uconvert(GenericQuantity(us"nm"), 5e-9u"m")) <: Quantity{Float64,<:SymbolicDimensions}
+    @test typeof(uconvert(GenericQuantity(us"nm"), 5e-9u"m")) <: constructorof(DEFAULT_QUANTITY_TYPE){Float64,<:SymbolicDimensions}
 
     q = 1.5u"Constants.M_sun"
     qs = uconvert(us"Constants.M_sun", 5.0 * q)
@@ -701,7 +731,7 @@ end
     VERSION >= v"1.8" &&
         @test_throws "You passed a quantity" uconvert(1.2us"m", 1.0u"m")
 
-    for Q in (Quantity, GenericQuantity)
+    for Q in (RealQuantity, Quantity, GenericQuantity)
         # Different types require converting both arguments:
         q = convert(Q{Float16}, 1.5u"g")
         qs = uconvert(convert(Q{Float16}, us"g"), 5 * q)
@@ -746,30 +776,102 @@ end
 
 
 @testset "Test ambiguities" begin
-    R = DEFAULT_DIM_BASE_TYPE
-    x = convert(R, 10)
-    y = convert(R, 5)
-    @test promote(x, y) == (x, y)
-    @test_throws ErrorException promote(x, convert(FixedRational{Int32,100}, 10))
-    @test promote_type(typeof(u"km/s"), typeof(convert(Quantity{Float32}, u"km/s"))) <: Quantity{Float64}
+    @testset "FixedRational" begin
+        R = DEFAULT_DIM_BASE_TYPE
+        x = convert(R, 10)
+        y = convert(R, 5)
+        @test promote(x, y) == (x, y)
+        @test_throws ErrorException promote(x, convert(FixedRational{Int32,100}, 10))
+        @test promote_type(typeof(u"km/s"), typeof(convert(Quantity{Float32}, u"km/s"))) <: Quantity{Float64}
+
+        x = FixedRational{Int32,100}(1)
+        # Need explicit `promote_rule` calls here so coverage picks it up
+        @test promote_rule(typeof(x), typeof(true)) == typeof(x)
+        @test promote_rule(typeof(true), typeof(x)) == typeof(x)
+        @test promote_rule(typeof(x), typeof(BigFloat(1))) == promote_type(Rational{Int32}, BigFloat)
+        @test promote_rule(typeof(BigFloat(1)), typeof(x)) == promote_type(Rational{Int32}, BigFloat)
+        @test promote_rule(typeof(x), typeof(π)) == promote_type(Rational{Int32}, typeof(π))
+        @test promote_rule(typeof(π), typeof(x)) == promote_type(Rational{Int32}, typeof(π))
+    end
 
-    x = 1.0u"m"
-    s = "test"
-    y = WeakRef(s)
-    @test_throws ErrorException x == y
-    @test_throws ErrorException y == x
+    @testset "Weakref" begin
+        x = 1.0u"m"
+        s = "test"
+        y = WeakRef(s)
+        @test_throws ErrorException x == y
+        @test_throws ErrorException y == x
+    end
+
+    @testset "Arrays" begin
+        qarr1 = QuantityArray(randn(3), u"km/s")
+        qarr2 = qarr1
+        @test convert(typeof(qarr2), qarr2) === qarr1
+    end
 
-    qarr1 = QuantityArray(randn(3), u"km/s")
-    qarr2 = qarr1
-    @test convert(typeof(qarr2), qarr2) === qarr1
+    @testset "Rational power law" begin
+        x = RealQuantity(1.0u"m")
+        y = x ^ (3//2)
+        @test y == Quantity(1.0, length=3//2)
+        @test typeof(y) == RealQuantity{Float64,DEFAULT_DIM_TYPE}
+    end
 
-    x = 1.0u"m"
-    y = x ^ (3//2)
-    @test y == Quantity(1.0, length=3//2)
-    @test typeof(y) == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @testset "Numeric promotion rules" begin
+        for Q in (RealQuantity, Quantity, GenericQuantity)
+            x = Q(1.0u"m")
+            @test promote_type(typeof(x), Bool) == typeof(x)
+            @test promote_type(Bool, typeof(x)) == typeof(x)
+            @test promote_type(typeof(x), BigFloat) == with_type_parameters(Q, BigFloat, DEFAULT_DIM_TYPE)
+            @test promote_type(BigFloat, typeof(x)) == with_type_parameters(Q, BigFloat, DEFAULT_DIM_TYPE)
+        end
+    end
+
+    @testset "Complex numbers" begin
+        for Q in (RealQuantity, Quantity, GenericQuantity)
+            x = 1.0im
+            y = Q(0.5u"m")
+            @test typeof(x * y) == with_type_parameters(promote_quantity_on_value(Q, ComplexF64), Complex{Float64}, DEFAULT_DIM_TYPE)
+            @test typeof(y * x) == with_type_parameters(promote_quantity_on_value(Q, ComplexF64), Complex{Float64}, DEFAULT_DIM_TYPE)
+            @test ustrip(x * y) == 0.5im
+            @test ustrip(y * x) == 0.5im
+
+            # Bool version
+            x = true * im
+            y = Q(0.5u"m")
+            @test typeof(x * y) == with_type_parameters(promote_quantity_on_value(Q, ComplexF64), Complex{Float64}, DEFAULT_DIM_TYPE)
+            @test typeof(y * x) == with_type_parameters(promote_quantity_on_value(Q, ComplexF64), Complex{Float64}, DEFAULT_DIM_TYPE)
+            @test ustrip(x * y) == 0.5im
+            @test ustrip(y * x) == 0.5im
+
+            # Complex powers
+            x = Q(0.5u"1")
+            out = x ^ (1 + 2im)
+            @test typeof(out) == with_type_parameters(promote_quantity_on_value(Q, ComplexF64), Complex{Float64}, DEFAULT_DIM_TYPE)
+            @test ustrip(out) ≈ 0.5 ^ (1 + 2im)
+
+            for CT in (Complex, Complex{Bool})
+                x = Q(1.0)
+                @test CT(x) == CT(1.0)
+                @test typeof(CT(x)) <: CT
+                x = Q(1.0, length=1)
+                @test_throws AssertionError CT(x)
+            end
+        end
+    end
+
+    @testset "Bool" begin
+        for Q in (RealQuantity, Quantity, GenericQuantity)
+            x = Q(1.0u"1")
+            @test Bool(x) == true
+            @test Bool(ustrip(x)) == true
+            @test Bool(Q(0.0u"1")) == false
+            @test Bool(ustrip(Q(0.0u"1"))) == false
+            x = Q(1.0u"m")
+            @test_throws AssertionError Bool(x)
+        end
+    end
 end
 
-for Q in (Quantity, GenericQuantity)
+for Q in (RealQuantity, Quantity, GenericQuantity)
     @testset "Arrays" begin
         @testset "Basics" begin
             x = QuantityArray(randn(32), Q(u"km/s"))
@@ -795,7 +897,7 @@ for Q in (Quantity, GenericQuantity)
 
             # Test default constructors:
             @test QuantityArray(ones(3), u"m/s") == QuantityArray(ones(3), length=1, time=-1)
-            @test typeof(QuantityArray(ones(3), u"m/s")) <: QuantityArray{Float64,1,<:Dimensions,<:Quantity,<:Array}
+            @test typeof(QuantityArray(ones(3), u"m/s")) <: QuantityArray{Float64,1,<:Dimensions,<:constructorof(DEFAULT_QUANTITY_TYPE),<:Array}
 
             # We can create quantity arrays with generic quantity
             @test typeof(QuantityArray([[1.0], [2.0, 3.0]], dimension(u"m/s"))) <: QuantityArray{<:Any,1,<:Dimensions,<:GenericQuantity,<:Array}
@@ -961,14 +1063,14 @@ for Q in (Quantity, GenericQuantity)
             @test ustrip(x .* y) == ustrip(x) .* ustrip(y)
         end
 
-        Q == Quantity && @testset "Broadcast different arrays" begin
+        Q in (Quantity, RealQuantity) && @testset "Broadcast different arrays" begin
             f(x, y, z, w) = x * y + z * w
             g(x, y, z, w) = f.(x, y, z, w)
 
             x = randn(32)
             y = QuantityArray(randn(32), u"km/s")
             z = rand(1:10, 32)
-            w = Quantity{Float32}(u"m/s")
+            w = Q{Float32}(u"m/s")
             @test typeof(g(x, y, z, w)) <: QuantityArray{Float64}
 
             y32 = QuantityArray(ustrip(y), dimension(y))
@@ -987,7 +1089,7 @@ for Q in (Quantity, GenericQuantity)
             @test typeof(b .* y) <: QuantityArray{Float64}
         end
 
-        Q == Quantity && @testset "Broadcast scalars" begin
+        Q in (RealQuantity, Quantity) && @testset "Broadcast scalars" begin
             for (x, qx) in ((0.5, 0.5u"s"), ([0.5, 0.2], GenericQuantity([0.5, 0.2], time=1)))
                 @test size(qx) == size(x)
                 @test length(qx) == length(x)
@@ -1149,21 +1251,21 @@ end
         :log, :log2, :log10, :log1p, :exp, :exp2, :exp10, :expm1, :frexp, :exponent,
         :atan, :atand
     )
-    for Q in (Quantity, GenericQuantity), D in (Dimensions, SymbolicDimensions), f in functions
+    for Q in (RealQuantity, Quantity, GenericQuantity), D in (Dimensions, SymbolicDimensions), f in functions
         # Only test on valid domain
         valid_inputs = filter(
             x -> is_input_valid(eval(f), x),
             5rand(100) .- 2.5
         )
         for x in valid_inputs[1:3]
-            qx_dimensionless = Quantity(x, D)
-            qx_dimensions = Quantity(x, convert(D, dimension(u"m/s")))
+            qx_dimensionless = Q(x, D)
+            qx_dimensions = Q(x, convert(D, dimension(u"m/s")))
             @eval @test $f($qx_dimensionless) == $f($x)
             @eval @test_throws DimensionError $f($qx_dimensions)
             if f in (:atan, :atand)
                 for y in valid_inputs[end-3:end]
-                    qy_dimensionless = Quantity(y, D)
-                    qy_dimensions = Quantity(y, convert(D, dimension(u"m/s")))
+                    qy_dimensionless = Q(y, D)
+                    qy_dimensions = Q(y, convert(D, dimension(u"m/s")))
                     @eval @test $f($y, $qx_dimensionless) == $f($y, $x)
                     @eval @test $f($qy_dimensionless, $x) == $f($y, $x)
                     @eval @test $f($qy_dimensionless, $qx_dimensionless) == $f($y, $x)
@@ -1182,12 +1284,13 @@ end
     functions = (
         :float, :abs, :real, :imag, :conj, :adjoint, :unsigned,
         :nextfloat, :prevfloat, :identity, :transpose,
-        :copysign, :flipsign, :mod, :modf,
+        :copysign, :flipsign, :modf,
         :floor, :trunc, :ceil, :significand,
-        :ldexp, :round,
+        :ldexp, :round, :mod, :rem
     )
-    for Q in (Quantity, GenericQuantity), D in (Dimensions, SymbolicDimensions), f in functions
+    for Q in (RealQuantity, Quantity, GenericQuantity), D in (Dimensions, SymbolicDimensions), f in functions
         T = f in (:abs, :real, :imag, :conj) ? ComplexF64 : Float64
+        T <: Complex && Q == RealQuantity && continue
         if f == :modf  # Functions that return multiple outputs
             for x in 5rand(T, 3) .- 2.5
                 dim = convert(D, dimension(u"m/s"))
@@ -1204,10 +1307,24 @@ end
                     qx_dimensions = Q(x, dim)
                     qy_dimensions = Q(y, dim)
                     @eval @test $f($qx_dimensions, $qy_dimensions) == $Q($f($x, $y), $dim)
-                    if f in (:copysign, :flipsign, :mod)
+                    if f in (:copysign, :flipsign)
                         # Also do test without dimensions
                         @eval @test $f($x, $qy_dimensions) == $f($x, $y)
                         @eval @test $f($qx_dimensions, $y) == $Q($f($x, $y), $dim)
+                    elseif f in (:rem, :mod)
+                        # Also do test without dimensions (need dimensionless)
+                        qx_dimensionless = Q(x, D)
+                        qy_dimensionless = Q(y, D)
+                        @eval @test $f($x, $qy_dimensionless) ≈ $Q($f($x, $y), $D)
+                        @eval @test $f($qx_dimensionless, $y) ≈ $Q($f($x, $y), $D)
+                        @eval @test_throws DimensionError $f($qx_dimensions, $y)
+                        @eval @test_throws DimensionError $f($x, $qy_dimensions)
+                        if f == :rem && VERSION >= v"1.9"
+                            # Can also do other rounding modes
+                            for r in (:RoundFromZero, :RoundNearest, :RoundUp, :RoundDown)
+                                @eval @test $f($qx_dimensions, $qy_dimensions, $r) ≈ $Q($f($x, $y, $r), $dim)
+                            end
+                        end
                     end
                 end
             end
@@ -1247,8 +1364,42 @@ end
     end
 end
 
+@testset "Assorted comparison functions" begin
+    functions = (
+        :(<=), :(<), :(>=), :(>), :isless, :isequal, :(==),
+    )
+    x = 5randn(10) .- 2.5
+    y = 5randn(10) .- 2.5
+    for Q in (RealQuantity, Quantity, GenericQuantity), D in (Dimensions, SymbolicDimensions), f in functions
+        ground_truth = @eval $f.($x, $y)
+        dim = convert(D, dimension(u"m/s"))
+        qx_dimensions = [Q(xi, dim) for xi in x]
+        qy_dimensions = [Q(yi, dim) for yi in y]
+        @eval @test all($f.($qx_dimensions, $qy_dimensions) .== $ground_truth)
+        if f in (:isequal, :(==))
+            # These include a dimension check in the result, rather than
+            # throwing an error
+            @eval @test !any($f.($qx_dimensions, $y))
+            @eval @test !any($f.($x, $qy_dimensions))
+        else
+            @eval @test_throws DimensionError $f($qx_dimensions[1], $y[1])
+            @eval @test_throws DimensionError $f($x[1], $qy_dimensions[1])
+        end
+        qx_dimensionless = [Q(xi, D) for xi in x]
+        qy_dimensionless = [Q(yi, D) for yi in y]
+        @eval @test all($f.($qx_dimensionless, $y) .== $ground_truth)
+        @eval @test all($f.($x, $qy_dimensionless) .== $ground_truth)
+
+        qx_real_dimensions = [RealQuantity(xi, dim) for xi in x]
+        qy_real_dimensions = [RealQuantity(yi, dim) for yi in y]
+        # Mixed quantity input
+        @eval @test all($f.($qx_real_dimensions, $qy_dimensions) .== $ground_truth)
+        @eval @test all($f.($qx_dimensions, $qy_real_dimensions) .== $ground_truth)
+    end
+end
+
 @testset "Test div" begin
-    for Q in (Quantity, GenericQuantity)
+    for Q in (RealQuantity, Quantity, GenericQuantity)
         x = Q{Int}(10, length=1)
         y = Q{Int}(3, mass=-1)
         @test div(x, y) == Q{Int}(3, length=1, mass=1)
@@ -1260,4 +1411,12 @@ end
             @test div(10, y, RoundFromZero) == Q{Int}(4, mass=1)
         end
     end
+    # Also test mixed quantities:
+    x = RealQuantity{Int}(10, length=1)
+    y = Quantity{Int}(3, mass=-1)
+    @test div(x, y) == Quantity{Int}(3, length=1, mass=1)
+    @test typeof(div(x, y)) <: Quantity{Int}
+    if VERSION >= v"1.9"
+        @test div(x, y, RoundFromZero) == Quantity{Int}(4, length=1, mass=1)
+    end
 end