diff --git a/Project.toml b/Project.toml
index 25c2c05..ceaf35d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,29 +1,27 @@
 name = "ColorVectorSpace"
 uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
-version = "0.8.7"
+version = "0.9.0"
 
 [deps]
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
-Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 FixedPointNumbers = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+TensorCore = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
 
 [compat]
-ColorTypes = "0.8, 0.9, 0.10"
-Colors = "0.9, 0.10, 0.11, 0.12"
-FixedPointNumbers = "0.6, 0.7, 0.8"
-SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
-StatsBase = "0.28, 0.29, 0.30, 0.31, 0.32, 0.33"
+ColorTypes = "0.10"
+FixedPointNumbers = "0.8"
+SpecialFunctions = "0.7, 0.8, 0.9, 0.10"
+TensorCore = "0.1"
 julia = "1"
 
 [extras]
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Statistics", "StatsBase", "LinearAlgebra", "Test"]
+test = ["Colors", "Statistics", "LinearAlgebra", "Test"]
diff --git a/README.md b/README.md
index 13de965..6a047cf 100644
--- a/README.md
+++ b/README.md
@@ -6,18 +6,27 @@
 This package is an add-on to [ColorTypes](https://github.com/JuliaGraphics/ColorTypes.jl), and provides fast
 mathematical operations for objects with types such as `RGB` and
 `Gray`.
+Specifically, with this package both grayscale and `RGB` colors are treated as if they are points
+in a normed vector space.
 
 ## Introduction
 
-Colorspaces such as RGB, unlike XYZ, are technically non-linear; the
-"colorimetrically correct" approach when averaging two RGBs is to
+Colorspaces such as RGB, unlike XYZ, are technically non-linear;
+perhaps the most "colorimetrically correct" approach when averaging two RGBs is to
 first convert each to XYZ, average them, and then convert back to RGB.
+Nor is there a clear definition of computing the sum of two colors.
+As a consequence, Julia's base color package,
+[ColorTypes](https://github.com/JuliaGraphics/ColorTypes.jl),
+does not support mathematical operations on colors.
 
 However, particularly in image processing it is common to ignore this
 concern, and for the sake of performance treat an RGB as if it were a
-3-vector.  This package provides such operations.
+3-vector.  The role of this package is to extend ColorTypes to support such mathematical operations.
+Specifically, it defines `+` and multiplication by a scalar (and by extension, `-` and division by a scalar) for grayscale and `AbstractRGB` colors.
+These are the requirements of a [vector space](https://en.wikipedia.org/wiki/Vector_space).
 
-If you're curious about how much difference it makes, the following
+If you're curious about how much the "colorimetrically correct" and
+"vector space" views differ, the following
 diagram might help. The first 10 `distinguishable_colors` were
 generated, and all pairs were averaged. Each box represents the
 average of the pair of diagonal elements intersected by tracing
@@ -27,16 +36,77 @@ represents the "RGB vector space" version.
 
 ![ColorVectorSpace](images/comparison.png "Comparison")
 
+This package also defines `norm(c)` for RGB and grayscale colors.
+This makes these color spaces [normed vector spaces](https://en.wikipedia.org/wiki/Normed_vector_space).
+Note that `norm` has been designed to satisfy equivalence of grayscale and RGB representations: if
+`x` is a scalar, then `norm(x) == norm(Gray(x)) == norm(RGB(x, x, x))`.
+Effectively, there's a division-by-3 in the `norm(::RGB)` case compared to the Euclidean interpretation of
+the RGB vector space.
+Equivalence is an important principle for the Colors ecosystem, and violations should be reported as likely bugs.
+
 ## Usage
 
 ```julia
 using ColorTypes, ColorVectorSpace
 ```
 
-That's it. Just by loading `ColorVectorSpace`, most basic mathematical
+For the most part, that's it; just by loading `ColorVectorSpace`, most basic mathematical
 operations will "just work" on `AbstractRGB`, `AbstractGray`
 (`Color{T,1}`), `TransparentRGB`, and `TransparentGray` objects.
 (See definitions for the latter inside of `ColorTypes`).
 
-If you discover missing operations, please open an issue, or better
-yet submit a pull request.
+However, there are some additional operations that you may need to distinguish carefully.
+
+### Multiplication
+
+Grayscale values are conceptually similar to scalars, and consequently it seems straightforward to define multiplication of two grayscale values.
+RGB values present more options.
+This package supports three different notions of multiplication: the inner product, the hadamard (elementwise) product, and the tensor product.
+
+```julia
+julia> c1, c2 = RGB(0.2, 0.3, 0.4), RGB(0.5, 0.3, 0.2)
+(RGB{Float64}(0.2,0.3,0.4), RGB{Float64}(0.5,0.3,0.1))
+
+julia> c1⋅c2     # \cdot<TAB> # or dot(c1, c2)
+0.09000000000000001
+
+# This is equivelant to `mapc(*, c1, c2)`
+julia> c1⊙c2     # \odot<TAB> # or hadamard(c1, c2)
+RGB{Float64}(0.1,0.09,0.08000000000000002)
+
+julia> c1⊗c2    # \otimes<TAB> # or tensor(c1, c2)
+RGBRGB{Float64}(
+ 0.1   0.06  0.04000000000000001
+ 0.15  0.09  0.06
+ 0.2   0.12  0.08000000000000002)
+```
+
+Note that `c1⋅c2 = (c1.r*c2.r + c1.g*c2.g + c1.b*c2.b)/3`, where the division by 3 ensures the equivalence `norm(x) == norm(Gray(x)) == norm(RGB(x, x, x))`.
+
+It is designed to not support the ordinary multiplication operation `*` because it is not obvious which one of these should be the default option.
+
+However, `*` is defined for grayscale since all these three multiplication operations (i.e., `⋅`, `⊙` and `⊗`) are equivalent in the 1D vector space.
+
+### Variance
+
+The variance `v = E((c - μ)^2)` (or its bias-corrected version) involves a multiplication,
+and to be consistent with the above you must specify which sense of multiplication you wish to use:
+
+```julia
+julia> cs = [c1, c2]
+2-element Array{RGB{Float64},1} with eltype RGB{Float64}:
+ RGB{Float64}(0.2,0.3,0.4)
+ RGB{Float64}(0.5,0.3,0.2)
+
+julia> varmult(⋅, cs)
+0.021666666666666667
+
+julia> varmult(⊙, cs)
+RGB{Float64}(0.045,0.0,0.020000000000000004)
+
+julia> varmult(⊗, cs)
+RGBRGB{Float64}(
+  0.045  0.0  -0.03
+  0.0    0.0   0.0
+ -0.03   0.0   0.020000000000000004)
+```
diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 0173b5d..0bb46ee 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -1,20 +1,17 @@
 module ColorVectorSpace
 
-using Colors, FixedPointNumbers, SpecialFunctions
+using ColorTypes, FixedPointNumbers, SpecialFunctions
+using TensorCore
+import TensorCore: ⊙, ⊗
+
+using FixedPointNumbers: ShorterThanInt
 
 import Base: ==, +, -, *, /, ^, <, ~
-import Base: abs, abs2, clamp, convert, copy, div, eps, isfinite, isinf,
-    isnan, isless, length, mapreduce, oneunit,
+import Base: abs, clamp, convert, copy, div, eps, float,
+    isfinite, isinf, isnan, isless, length, mapreduce, oneunit,
     promote_op, promote_rule, zero, trunc, floor, round, ceil, bswap,
-    mod, rem, atan, hypot, max, min, real, typemin, typemax
-import LinearAlgebra: norm
-import StatsBase: varm
-import SpecialFunctions: gamma, lgamma, lfact
-import Statistics: middle
-
-export nan
-
-# The unaryOps
+    mod, mod1, rem, atan, hypot, max, min, real, typemin, typemax
+# More unaryOps (mostly math functions)
 import Base:      conj, sin, cos, tan, sinh, cosh, tanh,
                   asin, acos, atan, asinh, acosh, atanh,
                   sec, csc, cot, asec, acsc, acot,
@@ -24,49 +21,30 @@ import Base:      conj, sin, cos, tan, sinh, cosh, tanh,
                   asind, atand, rad2deg, deg2rad,
                   log, log2, log10, log1p, exponent, exp,
                   exp2, exp10, expm1, cbrt, sqrt,
-                  significand,
-                  frexp, modf,
-                  float
-
-export dotc
+                  significand, frexp, modf
+import LinearAlgebra: norm, ⋅, dot, promote_leaf_eltypes  # norm1, norm2, normInf
+import SpecialFunctions: gamma, lgamma, lfact
+using Statistics
+import Statistics: middle, _mean_promote
 
-# TODO: we get rid of these definitions or move them to ColorTypes.jl
-TransparentRGBFloat{C<:AbstractRGB,T<:AbstractFloat} = TransparentColor{C,T,4}
-TransparentGrayFloat{C<:AbstractGray,T<:AbstractFloat} = TransparentColor{C,T,2}
-TransparentRGBNormed{C<:AbstractRGB,T<:Normed} = TransparentColor{C,T,4}
-TransparentGrayNormed{C<:AbstractGray,T<:Normed} = TransparentColor{C,T,2}
+export RGBRGB, nan, dotc, dot, ⋅, hadamard, ⊙, tensor, ⊗, norm, varmult
 
 MathTypes{T,C} = Union{AbstractRGB{T},TransparentRGB{C,T},AbstractGray{T},TransparentGray{C,T}}
 
-# convert(RGB{Float32}, NaN) doesn't and shouldn't work, so we need to reintroduce nan
-nan(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
-nan(::Type{C}) where {C<:MathTypes} = _nan(eltype(C), C)
-_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:AbstractGray} = (x = convert(T, NaN); C(x))
-_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:TransparentGray} = (x = convert(T, NaN); C(x,x))
-_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:AbstractRGB} = (x = convert(T, NaN); C(x,x,x))
-_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:TransparentRGB} = (x = convert(T, NaN); C(x,x,x,x))
-
-## Generic algorithms
-mapreduce(f, op::Union{typeof(&), typeof(|)}, a::MathTypes) = f(a)  # ambiguity
-mapreduce(f, op, a::MathTypes) = f(a)
-Base.add_sum(c1::MathTypes,c2::MathTypes) = mapc(Base.add_sum, c1, c2)
-Base.reduce_first(::typeof(Base.add_sum), c::MathTypes) = mapc(x->Base.reduce_first(Base.add_sum, x), c)
-function Base.reduce_empty(::typeof(Base.add_sum), ::Type{T}) where {T<:MathTypes}
-    z = Base.reduce_empty(Base.add_sum, eltype(T))
-    return zero(base_colorant_type(T){typeof(z)})
-end
-
-for f in (:trunc, :floor, :round, :ceil, :eps, :bswap)
-    @eval $f(g::Gray{T}) where {T} = Gray{T}($f(gray(g)))
-end
-eps(::Type{Gray{T}}) where {T} = Gray(eps(T))
+## Version compatibility with ColorTypes
 
-for f in (:trunc, :floor, :round, :ceil)
-    @eval $f(::Type{T}, g::Gray) where {T<:Integer} = Gray{T}($f(T, gray(g)))
+if !hasmethod(zero, (Type{TransparentGray},))
+    zero(::Type{C}) where {C<:TransparentGray} = C(0,0)
+    zero(::Type{C}) where {C<:AbstractRGB}     = C(0,0,0)
+    zero(::Type{C}) where {C<:TransparentRGB}  = C(0,0,0,0)
+    zero(p::Colorant) = zero(typeof(p))
 end
 
-for f in (:mod, :rem, :mod1)
-    @eval $f(x::Gray, m::Gray) = Gray($f(gray(x), gray(m)))
+if !hasmethod(one, (Type{TransparentGray},))
+    Base.one(::Type{C}) where {C<:TransparentGray} = C(1,1)
+    Base.one(::Type{C}) where {C<:AbstractRGB}     = C(1,1,1)
+    Base.one(::Type{C}) where {C<:TransparentRGB}  = C(1,1,1,1)
+    Base.one(p::Colorant) = one(typeof(p))
 end
 
 # Real values are treated like grays
@@ -74,24 +52,24 @@ if !hasmethod(gray, (Number,))
     ColorTypes.gray(x::Real) = x
 end
 
-dotc(x::T, y::T) where {T<:Real} = acc(x)*acc(y)
-dotc(x::Real, y::Real) = dotc(promote(x, y)...)
+## Traits and key utilities
 
 # Return types for arithmetic operations
 multype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)*zero(B)))
 sumtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)+zero(B)))
 divtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)/oneunit(B)))
 powtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)^zero(B)))
-multype(a::Colorant, b::Colorant) = multype(eltype(a),eltype(b))
-sumtype(a::Colorant, b::Colorant) = sumtype(eltype(a),eltype(b))
-divtype(a::Colorant, b::Colorant) = divtype(eltype(a),eltype(b))
-powtype(a::Colorant, b::Colorant) = powtype(eltype(a),eltype(b))
+sumtype(a::Colorant, b::Colorant) = coltype(sumtype(eltype(a),eltype(b)))
 
 coltype(::Type{T}) where {T<:Fractional} = T
-coltype(::Type{T}) where {T}             = Float64
+coltype(::Type{T}) where {T<:Number}     = floattype(T)
+
+acctype(::Type{T}) where {T<:FixedPoint}        = floattype(T)
+acctype(::Type{T}) where {T<:ShorterThanInt}    = Int
+acctype(::Type{Rational{T}}) where {T<:Integer} = typeof(zero(T)/oneunit(T))
+acctype(::Type{T}) where {T<:Real}              = T
 
-acctype(::Type{T}) where {T<:FixedPoint} = FixedPointNumbers.floattype(T)
-acctype(::Type{T}) where {T<:Number} = T
+acctype(::Type{T1}, ::Type{T2}) where {T1,T2} = acctype(promote_type(T1, T2))
 
 acc(x::Number) = convert(acctype(typeof(x)), x)
 
@@ -101,7 +79,6 @@ color_rettype(::Type{A}, ::Type{B}) where {A<:AbstractRGB,B<:AbstractRGB} = _col
 color_rettype(::Type{A}, ::Type{B}) where {A<:AbstractGray,B<:AbstractGray} = _color_rettype(base_colorant_type(A), base_colorant_type(B))
 color_rettype(::Type{A}, ::Type{B}) where {A<:TransparentRGB,B<:TransparentRGB} = _color_rettype(base_colorant_type(A), base_colorant_type(B))
 color_rettype(::Type{A}, ::Type{B}) where {A<:TransparentGray,B<:TransparentGray} = _color_rettype(base_colorant_type(A), base_colorant_type(B))
-_color_rettype(::Type{A}, ::Type{B}) where {A<:Colorant,B<:Colorant} = error("binary operation with $A and $B, return type is ambiguous")
 _color_rettype(::Type{C}, ::Type{C}) where {C<:Colorant} = C
 
 color_rettype(c1::Colorant, c2::Colorant) = color_rettype(typeof(c1), typeof(c2))
@@ -119,6 +96,47 @@ parametric(::Type{RGB24}, ::Type{N0f8}) = RGB24
 parametric(::Type{AGray32}, ::Type{N0f8}) = AGray32
 parametric(::Type{ARGB32}, ::Type{N0f8}) = ARGB32
 
+# Useful for leveraging iterator algorithms. Don't use this externally, as the implementation may change.
+channels(c::AbstractGray)    = (gray(c),)
+channels(c::TransparentGray) = (gray(c), alpha(c))
+channels(c::AbstractRGB)     = (red(c), green(c), blue(c))
+channels(c::TransparentRGB)  = (red(c), green(c), blue(c), alpha(c))
+
+nan(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
+nan(::Type{C}) where {C<:MathTypes} = _nan(eltype(C), C)
+_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:AbstractGray} = (x = convert(T, NaN); C(x))
+_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:TransparentGray} = (x = convert(T, NaN); C(x,x))
+_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:AbstractRGB} = (x = convert(T, NaN); C(x,x,x))
+_nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:TransparentRGB} = (x = convert(T, NaN); C(x,x,x,x))
+
+
+## Generic algorithms
+Base.add_sum(c1::MathTypes,c2::MathTypes) = mapc(Base.add_sum, c1, c2)
+Base.reduce_first(::typeof(Base.add_sum), c::MathTypes) = mapc(x->Base.reduce_first(Base.add_sum, x), c)
+function Base.reduce_empty(::typeof(Base.add_sum), ::Type{T}) where {T<:MathTypes}
+    z = Base.reduce_empty(Base.add_sum, eltype(T))
+    return zero(base_colorant_type(T){typeof(z)})
+end
+
+
+## Rounding & mod
+for f in (:trunc, :floor, :round, :ceil, :eps, :bswap)
+    @eval $f(g::Gray{T}) where {T} = Gray{T}($f(gray(g)))
+end
+eps(::Type{Gray{T}}) where {T} = Gray(eps(T))
+
+for f in (:trunc, :floor, :round, :ceil)
+    @eval $f(::Type{T}, g::Gray) where {T<:Integer} = Gray{T}($f(T, gray(g)))
+end
+
+for f in (:mod, :rem, :mod1)
+    @eval $f(x::Gray, m::Gray) = Gray($f(gray(x), gray(m)))
+end
+
+dotc(x::T, y::T) where {T<:Real} = acc(x)*acc(y)
+dotc(x::Real, y::Real) = dotc(promote(x, y)...)
+
+
 ## Math on Colors. These implementations encourage inlining and,
 ## for the case of Normed types, nearly halve the number of multiplications (for RGB)
 
@@ -159,32 +177,23 @@ end
 (/)(c::AbstractRGB, f::Integer) = (one(eltype(c))/f)*c
 (/)(c::TransparentRGB, f::Integer) = (one(eltype(c))/f)*c
 
+
+# New multiplication operators
+(⋅)(x::AbstractRGB, y::AbstractRGB)  = (T = acctype(eltype(x), eltype(y)); T(red(x))*T(red(y)) + T(green(x))*T(green(y)) + T(blue(x))*T(blue(y)))/3
+(⊙)(x::C, y::C) where C<:AbstractRGB = base_color_type(C)(red(x)*red(y), green(x)*green(y), blue(x)*blue(y))
+(⊙)(x::AbstractRGB, y::AbstractRGB)  = ⊙(promote(x, y)...)
+# ⊗ defined below
+
 isfinite(c::Colorant{T}) where {T<:Normed} = true
 isfinite(c::Colorant) = mapreducec(isfinite, &, true, c)
 isnan(c::Colorant{T}) where {T<:Normed} = false
 isnan(c::Colorant) = mapreducec(isnan, |, false, c)
 isinf(c::Colorant{T}) where {T<:Normed} = false
 isinf(c::Colorant) = mapreducec(isinf, |, false, c)
-abs(c::AbstractRGB) = abs(red(c))+abs(green(c))+abs(blue(c)) # should this have a different name?
-abs(c::AbstractRGB{T}) where {T<:Normed} = Float32(red(c))+Float32(green(c))+Float32(blue(c)) # should this have a different name?
-abs(c::TransparentRGB) = abs(red(c))+abs(green(c))+abs(blue(c))+abs(alpha(c)) # should this have a different name?
-abs(c::TransparentRGB{T}) where {T<:Normed} = Float32(red(c))+Float32(green(c))+Float32(blue(c))+Float32(alpha(c)) # should this have a different name?
-abs2(c::AbstractRGB) = red(c)^2+green(c)^2+blue(c)^2
-abs2(c::AbstractRGB{T}) where {T<:Normed} = Float32(red(c))^2+Float32(green(c))^2+Float32(blue(c))^2
-abs2(c::TransparentRGB) = (ret = abs2(color(c)); ret + convert(typeof(ret), alpha(c))^2)
-norm(c::AbstractRGB) = sqrt(abs2(c))
-norm(c::TransparentRGB) = sqrt(abs2(c))
-
-oneunit(::Type{C}) where {C<:AbstractRGB}     = C(1,1,1)
-oneunit(::Type{C}) where {C<:TransparentRGB}  = C(1,1,1,1)
-
-zero(::Type{C}) where {C<:AbstractRGB}    = C(0,0,0)
-zero(::Type{C}) where {C<:TransparentRGB} = C(0,0,0,0)
-zero(::Type{C}) where {C<:YCbCr} = C(0,0,0)
-zero(::Type{C}) where {C<:HSV} = C(0,0,0)
-oneunit(p::Colorant) = oneunit(typeof(p))
-Base.one(c::Colorant) = Base.one(typeof(c))
-zero(p::Colorant) = zero(typeof(p))
+abs(c::MathTypes) = mapc(abs, c)
+norm(c::MathTypes, p::Real=2) = (cc = channels(c); norm(cc, p)/(p == 0 ? length(cc) : length(cc)^(1/p)))
+
+promote_leaf_eltypes(x::Union{AbstractArray{T},Tuple{T,Vararg{T}}}) where {T<:MathTypes} = eltype(T)
 
 # These constants come from squaring the conversion to grayscale
 # (rec601 luma), and normalizing
@@ -212,18 +221,20 @@ const unaryOps = (:~, :conj, :abs,
                   :(SpecialFunctions.besselj0), :(SpecialFunctions.besselj1), :(SpecialFunctions.bessely0), :(SpecialFunctions.bessely1),
                   :(SpecialFunctions.eta), :(SpecialFunctions.zeta), :(SpecialFunctions.digamma))
 for op in unaryOps
-    @eval ($op)(c::AbstractGray) = $op(gray(c))
+    @eval ($op)(c::AbstractGray) = Gray($op(gray(c)))
 end
 
 middle(c::AbstractGray) = arith_colorant_type(c)(middle(gray(c)))
 middle(x::C, y::C) where {C<:AbstractGray} = arith_colorant_type(C)(middle(gray(x), gray(y)))
 
+_mean_promote(x::MathTypes, y::MathTypes) = mapc(FixedPointNumbers.Treduce, y)
+
 (*)(f::Real, c::AbstractGray) = arith_colorant_type(c){multype(typeof(f),eltype(c))}(f*gray(c))
 (*)(f::Real, c::TransparentGray) = arith_colorant_type(c){multype(typeof(f),eltype(c))}(f*gray(c), f*alpha(c))
 (*)(c::AbstractGray, f::Real) = (*)(f, c)
 (*)(c::TransparentGray, f::Real) = (*)(f, c)
 (/)(c::AbstractGray, f::Real) = (one(f)/f)*c
-(/)(n::Number, c::AbstractGray) = n/gray(c)
+(/)(n::Number, c::AbstractGray) = base_color_type(c)(n/gray(c))
 (/)(c::TransparentGray, f::Real) = (one(f)/f)*c
 (/)(c::AbstractGray, f::Integer) = (one(eltype(c))/f)*c
 (/)(c::TransparentGray, f::Integer) = (one(eltype(c))/f)*c
@@ -238,12 +249,18 @@ middle(x::C, y::C) where {C<:AbstractGray} = arith_colorant_type(C)(middle(gray(
 (+)(c::TransparentGray) = c
 (-)(c::AbstractGray) = typeof(c)(-gray(c))
 (-)(c::TransparentGray) = typeof(c)(-gray(c),-alpha(c))
-(/)(a::AbstractGray, b::AbstractGray) = gray(a)/gray(b)
-div(a::AbstractGray, b::AbstractGray) = div(gray(a), gray(b))
-(+)(a::AbstractGray, b::Number) = gray(a)+b
-(-)(a::AbstractGray, b::Number) = gray(a)-b
-(+)(a::Number, b::AbstractGray) = a+gray(b)
-(-)(a::Number, b::AbstractGray) = a-gray(b)
+(/)(a::C, b::C) where C<:AbstractGray = base_color_type(C)(gray(a)/gray(b))
+(/)(a::AbstractGray, b::AbstractGray) = /(promote(a, b)...)
+(+)(a::AbstractGray, b::Number) = base_color_type(a)(gray(a)+b)
+(+)(a::Number, b::AbstractGray) = b+a
+(-)(a::AbstractGray, b::Number) = base_color_type(a)(gray(a)-b)
+(-)(a::Number, b::AbstractGray) = base_color_type(b)(a-gray(b))
+
+(⋅)(x::AbstractGray, y::AbstractGray) = gray(x)*gray(y)
+(⊙)(x::C, y::C) where C<:AbstractGray = base_color_type(C)(gray(x)*gray(y))
+(⊙)(x::AbstractGray, y::AbstractGray) = ⊙(promote(x, y)...)
+(⊗)(x::AbstractGray, y::AbstractGray) = ⊙(x, y)
+
 max(a::T, b::T) where {T<:AbstractGray} = T(max(gray(a),gray(b)))
 max(a::AbstractGray, b::AbstractGray) = max(promote(a,b)...)
 max(a::Number, b::AbstractGray) = max(promote(a,b)...)
@@ -253,58 +270,26 @@ min(a::AbstractGray, b::AbstractGray) = min(promote(a,b)...)
 min(a::Number, b::AbstractGray) = min(promote(a,b)...)
 min(a::AbstractGray, b::Number) = min(promote(a,b)...)
 
-norm(c::AbstractGray) = abs(gray(c))
-abs(c::TransparentGray) = abs(gray(c))+abs(alpha(c)) # should this have a different name?
-abs(c::TransparentGrayNormed) = Float32(gray(c)) + Float32(alpha(c)) # should this have a different name?
-abs2(c::AbstractGray) = gray(c)^2
-abs2(c::AbstractGray{T}) where {T<:Normed} = Float32(gray(c))^2
-abs2(c::TransparentGray) = gray(c)^2+alpha(c)^2
-abs2(c::TransparentGrayNormed) = Float32(gray(c))^2 + Float32(alpha(c))^2
-atan(x::Gray, y::Gray) = atan(convert(Real, x), convert(Real, y))
-hypot(x::Gray, y::Gray) = hypot(convert(Real, x), convert(Real, y))
-norm(c::TransparentGray) = sqrt(abs2(c))
-
-(<)(g1::AbstractGray, g2::AbstractGray) = gray(g1) < gray(g2)
-(<)(c::AbstractGray, r::Real) = gray(c) < r
-(<)(r::Real, c::AbstractGray) = r < gray(c)
-if !hasmethod(isless, Tuple{AbstractGray,AbstractGray})  # this was moved to ColorTypes 0.10
-    isless(g1::AbstractGray, g2::AbstractGray) = isless(gray(g1), gray(g2))
+atan(x::AbstractGray, y::AbstractGray)  = atan(gray(x), gray(y))
+hypot(x::AbstractGray, y::AbstractGray) = hypot(gray(x), gray(y))
+
+if which(<, Tuple{AbstractGray,AbstractGray}).module === Base  # planned for ColorTypes 0.11
+    (<)(g1::AbstractGray, g2::AbstractGray) = gray(g1) < gray(g2)
+    (<)(c::AbstractGray, r::Real) = gray(c) < r
+    (<)(r::Real, c::AbstractGray) = r < gray(c)
 end
-isless(c::AbstractGray, r::Real) = isless(gray(c), r)
-isless(r::Real, c::AbstractGray) = isless(r, gray(c))
-
-function Base.isapprox(x::AbstractArray{Cx},
-                       y::AbstractArray{Cy};
-                       rtol::Real=Base.rtoldefault(eltype(Cx),eltype(Cy),0),
-                       atol::Real=0,
-                       norm::Function=norm) where {Cx<:MathTypes,Cy<:MathTypes}
-    d = norm(x - y)
-    if isfinite(d)
-        return d <= atol + rtol*max(norm(x), norm(y))
-    else
-        # Fall back to a component-wise approximate comparison
-        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol), zip(x, y))
-    end
+if !hasmethod(isless, Tuple{AbstractGray,Real})  # planned for ColorTypes 0.11
+    isless(c::AbstractGray, r::Real) = isless(gray(c), r)
+    isless(r::Real, c::AbstractGray) = isless(r, gray(c))
 end
 
-zero(::Type{C}) where {C<:TransparentGray} = C(0,0)
-oneunit(::Type{C}) where {C<:TransparentGray} = C(1,1)
-
 dotc(x::T, y::T) where {T<:AbstractGray} = acc(gray(x))*acc(gray(y))
 dotc(x::AbstractGray, y::AbstractGray) = dotc(promote(x, y)...)
 
-float(::Type{T}) where {T<:Gray} = typeof(float(zero(T)))
-
 # Mixed types
 (+)(a::MathTypes, b::MathTypes) = (+)(promote(a, b)...)
 (-)(a::MathTypes, b::MathTypes) = (-)(promote(a, b)...)
 
-# Arrays---necessary methods
-+(A::AbstractArray{C}) where {C<:MathTypes} = A
-+(A::Array{C}) where {C<:MathTypes} = A
-
-varm(v::AbstractArray{C}, s::AbstractGray; corrected::Bool=true) where {C<:AbstractGray} =
-        varm(map(gray,v),gray(s); corrected=corrected)
 real(::Type{C}) where {C<:AbstractGray} = real(eltype(C))
 
 # To help type inference
@@ -316,46 +301,124 @@ typemax(::Type{T}) where {T<:ColorTypes.AbstractGray} = T(typemax(eltype(T)))
 typemin(::T) where {T<:ColorTypes.AbstractGray} = T(typemin(eltype(T)))
 typemax(::T) where {T<:ColorTypes.AbstractGray} = T(typemax(eltype(T)))
 
-include("precompile.jl")
-_precompile_()
-
-## Deprecations
-
-@deprecate (+)(A::AbstractArray{CV}, b::AbstractRGB) where {CV<:AbstractRGB} (.+)(A, b)
-@deprecate (+)(b::AbstractRGB, A::AbstractArray{CV}) where {CV<:AbstractRGB} (.+)(b, A)
-@deprecate (-)(A::AbstractArray{CV}, b::AbstractRGB) where {CV<:AbstractRGB} (.-)(A, b)
-@deprecate (-)(b::AbstractRGB, A::AbstractArray{CV}) where {CV<:AbstractRGB} (.-)(b, A)
-@deprecate (*)(A::AbstractArray{T}, b::AbstractRGB) where {T<:Number} A.*b
-@deprecate (*)(b::AbstractRGB, A::AbstractArray{T}) where {T<:Number} A.*b
-
-@deprecate (+)(A::AbstractArray{CV}, b::TransparentRGB) where {CV<:TransparentRGB} (.+)(A, b)
-@deprecate (+)(b::TransparentRGB, A::AbstractArray{CV}) where {CV<:TransparentRGB} (.+)(b, A)
-@deprecate (-)(A::AbstractArray{CV}, b::TransparentRGB) where {CV<:TransparentRGB} (.-)(A, b)
-@deprecate (-)(b::TransparentRGB, A::AbstractArray{CV}) where {CV<:TransparentRGB} (.-)(b, A)
-@deprecate (*)(A::AbstractArray{T}, b::TransparentRGB) where {T<:Number} A.*b
-@deprecate (*)(b::TransparentRGB, A::AbstractArray{T}) where {T<:Number} A.*b
-
-@deprecate (+)(A::AbstractArray{CV}, b::AbstractGray) where {CV<:AbstractGray} (.+)(A, b)
-@deprecate (+)(b::AbstractGray, A::AbstractArray{CV}) where {CV<:AbstractGray} (.+)(b, A)
-@deprecate (-)(A::AbstractArray{CV}, b::AbstractGray) where {CV<:AbstractGray} (.-)(A, b)
-@deprecate (-)(b::AbstractGray, A::AbstractArray{CV}) where {CV<:AbstractGray} (.-)(b, A)
-@deprecate (*)(A::AbstractArray{T}, b::AbstractGray) where {T<:Number} A.*b
-@deprecate (*)(b::AbstractGray, A::AbstractArray{T}) where {T<:Number} A.*b
-@deprecate (/)(A::AbstractArray{C}, b::AbstractGray) where {C<:AbstractGray} A./b
-
-@deprecate (+)(A::AbstractArray{CV}, b::TransparentGray) where {CV<:TransparentGray} (.+)(A, b)
-@deprecate (+)(b::TransparentGray, A::AbstractArray{CV}) where {CV<:TransparentGray} (.+)(b, A)
-@deprecate (-)(A::AbstractArray{CV}, b::TransparentGray) where {CV<:TransparentGray} (.-)(A, b)
-@deprecate (-)(b::TransparentGray, A::AbstractArray{CV}) where {CV<:TransparentGray} (.-)(b, A)
-@deprecate (*)(A::AbstractArray{T}, b::TransparentGray) where {T<:Number} A.*b
-@deprecate (*)(b::TransparentGray, A::AbstractArray{T}) where {T<:Number} A.*b
-
-## Deprecations
-
-## From 2020-Sept-9
-# Since ImageContrastAdjustment is now doing its own binning, I think this can be safely deprecated and we can eliminate
-# the dependency on StatsBase.
-import StatsBase: histrange
-@deprecate histrange(v::AbstractArray{Gray{T}}, n::Integer) where {T} histrange(convert(Array{Float32}, map(gray, v)), n, :right)
+## RGB tensor products
+
+"""
+    RGBRGB(rr, gr, br, rg, gg, bg, rb, gb, bb)
+
+Represent the [tensor product](https://en.wikipedia.org/wiki/Tensor_product) of two RGB values.
+
+# Example
+
+```jldoctest
+julia> a, b = RGB(0.2f0, 0.3f0, 0.5f0), RGB(0.77f0, 0.11f0, 0.22f0)
+(RGB{Float32}(0.2f0,0.3f0,0.5f0), RGB{Float32}(0.77f0,0.11f0,0.22f0))
+
+julia> a ⊗ b
+RGBRGB{Float32}(
+ 0.154f0  0.022f0  0.044f0
+ 0.231f0  0.033f0  0.066f0
+ 0.385f0  0.055f0  0.11f0 )
+"""
+struct RGBRGB{T}
+    rr::T
+    gr::T
+    br::T
+    rg::T
+    gg::T
+    bg::T
+    rb::T
+    gb::T
+    bb::T
+end
+Base.eltype(::Type{RGBRGB{T}}) where T = T
+Base.Matrix{T}(p::RGBRGB) where T = T[p.rr p.rg p.rb;
+                                      p.gr p.gg p.gb;
+                                      p.br p.bg p.bb]
+Base.Matrix(p::RGBRGB{T}) where T = Matrix{T}(p)
+
+function Base.show(io::IO, p::RGBRGB)
+    print(io, "RGBRGB{", eltype(p), "}(\n")
+    Base.print_matrix(io, Matrix(p))
+    print(io, ')')
+end
+
+Base.zero(::Type{RGBRGB{T}}) where T = (z = zero(T); RGBRGB(z, z, z, z, z, z, z, z, z))
+Base.zero(a::RGBRGB) = zero(typeof(a))
+
++(a::RGBRGB) = a
+-(a::RGBRGB) = RGBRGB(-a.rr, -a.gr, -a.br, -a.rg, -a.gg, -a.bg, -a.rb, -a.gb, -a.bb)
++(a::RGBRGB, b::RGBRGB) = RGBRGB(a.rr + b.rr, a.gr + b.gr, a.br + b.br,
+                                 a.rg + b.rg, a.gg + b.gg, a.bg + b.bg,
+                                 a.rb + b.rb, a.gb + b.gb, a.bb + b.bb)
+-(a::RGBRGB, b::RGBRGB) = +(a, -b)
+*(α::Real, a::RGBRGB) = RGBRGB(α*a.rr, α*a.gr, α*a.br, α*a.rg, α*a.gg, α*a.bg, α*a.rb, α*a.gb, α*a.bb)
+*(a::RGBRGB, α::Real) = α*a
+/(a::RGBRGB, α::Real) = (1/α)*a
+
+function ⊗(a::AbstractRGB, b::AbstractRGB)
+    ar, ag, ab = red(a), green(a), blue(a)
+    br, bg, bb = red(b), green(b), blue(b)
+    agbr, abbg, arbb, abbr, arbg, agbb = ag*br, ab*bg, ar*bb, ab*br, ar*bg, ag*bb
+    return RGBRGB(ar*br, agbr, abbr, arbg, ag*bg, abbg, arbb, agbb, ab*bb)
+end
+
+"""
+    varmult(op, itr; corrected::Bool=true, mean=Statistics.mean(itr), dims=:)
+
+Compute the variance of elements of `itr`, using `op` as the multiplication operator.
+The keyword arguments behave identically to those of `Statistics.var`.
+
+# Example
+
+```julia
+julia> cs = [RGB(0.2, 0.3, 0.4), RGB(0.5, 0.3, 0.2)]
+2-element Array{RGB{Float64},1} with eltype RGB{Float64}:
+ RGB{Float64}(0.2,0.3,0.4)
+ RGB{Float64}(0.5,0.3,0.2)
+
+julia> varmult(⋅, cs)
+0.021666666666666667
+
+julia> varmult(⊙, cs)
+RGB{Float64}(0.045,0.0,0.020000000000000004)
+
+julia> varmult(⊗, cs)
+RGBRGB{Float64}(
+  0.045  0.0  -0.03
+  0.0    0.0   0.0
+ -0.03   0.0   0.020000000000000004)
+```
+"""
+function varmult(op, itr; corrected::Bool=true, dims=:, mean=Statistics.mean(itr; dims=dims))
+    if dims === (:)
+        v = mapreduce(c->(Δc = c-mean; op(Δc, Δc)), +, itr; dims=dims)
+        n = length(itr)
+    else
+        # TODO: avoid temporary creation
+        v = mapreduce(Δc->op(Δc, Δc), +, itr .- mean; dims=dims)
+        n = length(itr) // length(v)
+    end
+    return v / (corrected ? max(1, n-1) : max(1, n))
+end
+
+function __init__()
+    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :register_error_hint)
+        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+            if exc.f === _color_rettype && length(argtypes) >= 2
+                # Color is not necessary, this is just to show it's possible.
+                A, B = argtypes
+                A !== B && print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
+            end
+        end
+    end
+end
+
+## Precompilation
+
+if Base.VERSION >= v"1.4.2"
+    include("precompile.jl")
+    _precompile_()
+end
 
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index a75b7ac..c707eac 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,18 +1,14 @@
-module ColorVectorSpaceTests
-
-using LinearAlgebra, Statistics
-using ColorVectorSpace, Colors, FixedPointNumbers, StatsBase
+using LinearAlgebra, Statistics, SpecialFunctions
+using ColorVectorSpace, Colors, FixedPointNumbers
 
 using Test
 
-const var = Statistics.var
+n8sum(x,y) = Float64(N0f8(x)) + Float64(N0f8(y))
 
 macro test_colortype_approx_eq(a, b)
     :(test_colortype_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
 end
 
-n8sum(x,y) = Float64(N0f8(x)) + Float64(N0f8(y))
-
 function test_colortype_approx_eq(a::Colorant, b::Colorant, astr, bstr)
     @test typeof(a) == typeof(b)
     n = length(fieldnames(typeof(a)))
@@ -21,8 +17,29 @@ function test_colortype_approx_eq(a::Colorant, b::Colorant, astr, bstr)
     end
 end
 
+struct RatRGB <: AbstractRGB{Rational{Int}}
+    r::Rational{Int}
+    g::Rational{Int}
+    b::Rational{Int}
+end
+ColorTypes.red(c::RatRGB)   = c.r
+ColorTypes.green(c::RatRGB) = c.g
+ColorTypes.blue(c::RatRGB)  = c.b
+
 @testset "Colortypes" begin
 
+    @testset "convert" begin
+        for x in (0.5, 0.5f0, NaN, NaN32, N0f8(0.5))
+            @test @inferred(convert(Gray{typeof(x)}, x))  === @inferred(convert(Gray, x))  === Gray(x)
+            @test @inferred(convert(RGB{typeof(x)}, x))   === @inferred(convert(RGB, x))   === RGB(x, x, x)
+            # These should be fixed by a future release of ColorTypes
+            @test_broken @inferred(convert(AGray{typeof(x)}, x)) === @inferred(convert(AGray, x)) === AGray(x, 1)
+            @test_broken @inferred(convert(ARGB{typeof(x)}, x))  === @inferred(convert(ARGB, x))  === ARGB(x, x, x, 1)
+            @test_broken @inferred(convert(GrayA{typeof(x)}, x)) === @inferred(convert(GrayA, x)) === GrayA(x, 1)
+            @test_broken @inferred(convert(RGBA{typeof(x)}, x))  === @inferred(convert(RGBA, x))  === RGBA(x, x, x, 1)
+        end
+    end
+
     @testset "nan" begin
         function make_checked_nan(::Type{T}) where T
             x = nan(T)
@@ -39,34 +56,47 @@ end
         end
     end
 
+    @testset "traits" begin
+        @test floattype(Gray{N0f8}) === Gray{float(N0f8)}
+    end
+
     @testset "Arithmetic with Gray" begin
         cf = Gray{Float32}(0.1)
-        @test +cf == cf
-        @test -cf == Gray(-0.1f0)
+        @test @inferred(+cf) === cf
+        @test @inferred(-cf) === Gray(-0.1f0)
+        @test @inferred(one(cf)*cf) === cf
+        @test oneunit(cf) === Gray(1.0f0)
         ccmp = Gray{Float32}(0.2)
-        @test 2*cf == ccmp
-        @test cf*2 == ccmp
-        @test ccmp/2 == cf
-        @test 2.0f0*cf == ccmp
-        @test_colortype_approx_eq cf*cf Gray{Float32}(0.01)
-        @test_colortype_approx_eq cf^2 Gray{Float32}(0.01)
-        @test_colortype_approx_eq cf^3.0f0 Gray{Float32}(0.001)
-        @test eltype(2.0*cf) == Float64
-        @test abs2(ccmp) == 0.2f0^2
-        @test norm(cf) == 0.1f0
-        @test sum(abs2, ccmp) == 0.2f0^2
+        @test @inferred(2*cf) === cf*2 === 2.0f0*cf === cf*2.0f0 === ccmp
+        @test @inferred(ccmp/2) === cf
+        @test @inferred(cf*cf) === Gray{Float32}(0.1f0*0.1f0)
+        @test @inferred(Gray{N0f32}(0.5)*Gray(0.5f0)) === Gray(Float64(N0f32(0.5)) * 0.5)
+        @test @inferred(cf^2 ) === Gray{Float32}(0.1f0*0.1f0)
+        @test @inferred(cf^3.0f0) === Gray{Float32}(0.1f0^3.0f0)
+        @test @inferred(2.0*cf) === cf*2.0 === Gray(2.0*0.1f0)
+        cf64 = Gray(0.2)
+        @test cf / cf64 === Gray(0.1f0/0.2)
+        @test_throws MethodError cf ÷ cf
+        @test cf + 0.1     === 0.1 + cf        === Gray(Float64(0.1f0) + 0.1)
+        @test cf64 - 0.1f0 === -(0.1f0 - cf64) === Gray( 0.2 - Float64(0.1f0))
+        @test_throws MethodError abs2(ccmp)
+        @test norm(cf) == norm(cf, 2) == norm(gray(cf))
+        @test norm(cf, 1)   == norm(gray(cf), 1)
+        @test norm(cf, Inf) == norm(gray(cf), Inf)
+        @test @inferred(abs(cf)) === Gray(0.1f0)
         cu = Gray{N0f8}(0.1)
-        @test 2*cu == Gray(2*cu.val)
-        @test 2.0f0*cu == Gray(2.0f0*cu.val)
+        @test @inferred(2*cu) === cu*2 === Gray(2*gray(cu))
+        @test @inferred(2.0f0*cu) === cu*2.0f0 === Gray(2.0f0*gray(cu))
         f = N0f8(0.5)
-        @test (f*cu).val ≈ f*cu.val
-        @test cf/2.0f0 == Gray{Float32}(0.05)
-        @test cu/2 == Gray(cu.val/2)
-        @test cu/0.5f0 == Gray(cu.val/0.5f0)
-        @test cf+cf == ccmp
+        @test @inferred(gray(f*cu)) === gray(cu*f) ===f*gray(cu)
+        @test @inferred(cf/2.0f0) === Gray{Float32}(0.05)
+        @test @inferred(cu/2) === Gray(cu.val/2)
+        @test @inferred(cu/0.5f0) === Gray(cu.val/0.5f0)
+        @test @inferred(cf+cf) === ccmp
         @test isfinite(cf)
         @test isfinite(Gray(true))
         @test !isinf(cf)
+        @test !isinf(Gray(f))
         @test !isnan(cf)
         @test !isfinite(Gray(NaN))
         @test !isinf(Gray(NaN))
@@ -74,13 +104,20 @@ end
         @test !isfinite(Gray(Inf))
         @test Gray(Inf) |> isinf
         @test !isnan(Gray(Inf))
-        @test abs(Gray(0.1)) ≈ 0.1
-        @test eps(Gray{N0f8}) == Gray(eps(N0f8))  # #282
+        @test abs(Gray(0.1)) === Gray(0.1)
+        @test eps(Gray{N0f8}) === Gray(eps(N0f8))  # #282
         @test atan(Gray(0.1), Gray(0.2)) == atan(0.1, 0.2)
+        @test hypot(Gray(0.2), Gray(0.3)) === hypot(0.2, 0.3)
+        # Multiplication
+        @test cf ⋅ cf   === gray(cf)^2
+        @test cf ⋅ cf64 === gray(cf) * gray(cf64)
+        @test cf ⊙ cf   === Gray(gray(cf)^2)
+        @test cf ⊙ cf64 === Gray(gray(cf) * gray(cf64))
+        @test cf ⊗ cf   === Gray(gray(cf)^2)
+        @test cf ⊗ cf64 === Gray(gray(cf) * gray(cf64))
 
         acu = Gray{N0f8}[cu]
         acf = Gray{Float32}[cf]
-        @test +acu === acu
         @test @inferred(acu./trues(1)) == acu
         @test typeof(acu./trues(1)) == Vector{typeof(cu/true)}
         @test @inferred(ones(Int, 1)./acu) == [1/cu]
@@ -102,22 +139,23 @@ end
         @test (acf./Gray{Float32}(2))[1] ≈ 0.05f0
         @test (acu/2)[1] == Gray(gray(acu[1])/2)
         @test (acf/2)[1] ≈ Gray{Float32}(0.05f0)
-        @test sum(abs2, [cf, ccmp]) ≈ 0.05f0
 
-        @test gray(0.8) == 0.8
+        @test gray(0.8) === 0.8
 
         a = Gray{N0f8}[0.8,0.7]
-        @test sum(a) == Gray(n8sum(0.8,0.7))
-        @test abs( var(a) - (a[1]-a[2])^2 / 2 ) <= 0.001
+        @test a == a
+        @test a === a
         @test isapprox(a, a)
+        @test sum(a) == Gray(n8sum(0.8,0.7))
+        @test sum(a[1:1]) == a[1]
+        @test abs( varmult(*, a) - (a[1]-a[2])^2 / 2 ) <= 0.001
+
         @test real(Gray{Float32}) <: Real
         @test zero(ColorTypes.Gray) == 0
         @test oneunit(ColorTypes.Gray) == 1
-        a = Gray{N0f8}[0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
-        @test StatsBase.histrange(a,10) == 0.1f0:0.1f0:1f0
 
         @test typeof(float(Gray{N0f16}(0.5))) <: AbstractFloat
-        @test quantile( Gray{N0f16}[0.0,0.5,1.0], 0.1) ≈ 0.10000152590218968
+        @test quantile( Gray{N0f16}[0.0,0.5,1.0], 0.1) ≈ 0.1 atol=eps(N0f16)
         @test middle(Gray(0.2)) === Gray(0.2)
         @test middle(Gray(0.2), Gray(0.4)) === Gray((0.2+0.4)/2)
 
@@ -144,13 +182,23 @@ end
         @test !(isless(0.5, g1))
         @test g1 < 0.5
         @test !(0.5 < g1)
-        @test (@inferred(max(g1, g2)) ) == g2
-        @test max(g1, 0.1) == 0.2
+        @test @inferred(max(g1, g2)) === g2
+        @test @inferred(max(g1, Gray(0.3))) === Gray(0.3)
+        @test max(g1, 0.1) === max(0.1, g1) === Float64(gray(g1))
         @test (@inferred(min(g1, g2)) ) == g1
-        @test min(g1, 0.1) == 0.1
+        @test min(g1, 0.1) === min(0.1, g1) === 0.1
         a = Gray{Float64}(0.9999999999999999)
         b = Gray{Float64}(1.0)
 
+        @test (Gray(0.3) < Gray(NaN)) == (0.3 < NaN)
+        @test (Gray(NaN) < Gray(0.3)) == (NaN < 0.3)
+        @test isless(Gray(0.3), Gray(NaN)) == isless(0.3, NaN)
+        @test isless(Gray(NaN), Gray(0.3)) == isless(NaN, 0.3)
+        @test isless(Gray(0.3), NaN) == isless(0.3, NaN)
+        @test isless(Gray(NaN), 0.3) == isless(NaN, 0.3)
+        @test isless(0.3, Gray(NaN)) == isless(0.3, NaN)
+        @test isless(NaN, Gray(0.3)) == isless(NaN, 0.3)
+
         @test isapprox(a, b)
         a = Gray{Float64}(0.99)
         @test !(isapprox(a, b, rtol = 0.01))
@@ -158,18 +206,38 @@ end
     end
 
     @testset "Unary operations with Gray" begin
+        ntested = 0
         for g in (Gray(0.4), Gray{N0f8}(0.4))
             @test @inferred(zero(g)) === typeof(g)(0)
             @test @inferred(oneunit(g)) === typeof(g)(1)
-            for op in ColorVectorSpace.unaryOps
-                try
-                    v = @eval $op(gray(g))  # if this fails, don't bother
-                    @show op
-                    @test op(g) == v
-                catch
+            for opgroup in (ColorVectorSpace.unaryOps, (:trunc, :floor, :round, :ceil, :eps, :bswap))
+                for op in opgroup
+                    op ∈ (:frexp, :exponent, :modf, :lfact) && continue
+                    op === :~ && eltype(g) === Float64 && continue
+                    op === :significand && eltype(g) === N0f8 && continue
+                    try
+                        v = @eval $op(gray($g))  # if this fails, don't bother with the next test
+                        @test @eval($op($g)) === Gray(v)
+                        ntested += 1
+                    catch ex
+                        @test ex isa Union{DomainError,ArgumentError}
+                    end
                 end
             end
         end
+        @test ntested > 130
+        for g in (Gray{N0f8}(0.4), Gray{N0f8}(0.6))
+            for op in (:trunc, :floor, :round, :ceil)
+                v = @eval $op(Bool, gray($g))
+                @test @eval($op(Bool, $g)) === Gray(v)
+            end
+        end
+        for (g1, g2) in ((Gray(0.4), Gray(0.3)), (Gray(N0f8(0.4)), Gray(N0f8(0.3))))
+            for op in (:mod, :rem, :mod1)
+                v = @eval $op(gray($g1), gray($g2))
+                @test @eval($op($g1, $g2)) === Gray(v)
+            end
+        end
         u = N0f8(0.4)
         @test ~Gray(u) == Gray(~u)
         @test -Gray(u) == Gray(-u)
@@ -248,9 +316,12 @@ end
         @test !isfinite(RGB(1, Inf, 0.5))
         @test isinf(RGB(1, Inf, 0.5))
         @test !isnan(RGB(1, Inf, 0.5))
-        @test abs(RGB(0.1,0.2,0.3)) ≈ 0.6
-        @test sum(abs2, RGB(0.1,0.2,0.3)) ≈ 0.14
-        @test norm(RGB(0.1,0.2,0.3)) ≈ sqrt(0.14)
+        @test abs(RGB(0.1,0.2,0.3)) == RGB(0.1,0.2,0.3)
+        @test_throws MethodError abs2(RGB(0.1,0.2,0.3))
+        @test_throws MethodError sum(abs2, RGB(0.1,0.2,0.3))
+        @test norm(RGB(0.1,0.2,0.3)) ≈ sqrt(0.14)/sqrt(3)
+
+        @test_throws MethodError RGBX(0, 0, 1) + XRGB(1, 0, 0)
 
         acu = RGB{N0f8}[cu]
         acf = RGB{Float32}[cf]
@@ -269,6 +340,7 @@ end
 
         a = RGB{N0f8}[RGB(1,0,0), RGB(1,0.8,0)]
         @test sum(a) == RGB(2.0,0.8,0)
+        @test sum(typeof(a)()) == RGB(0.0,0.0,0)
         @test isapprox(a, a)
         a = RGB{Float64}(1.0, 1.0, 0.9999999999999999)
         b = RGB{Float64}(1.0, 1.0, 1.0)
@@ -284,6 +356,30 @@ end
         @test RGB24(1,0,0)/2.0 === RGB(0.5,0,0)
         # issue #133
         @test RGB24(1, 0, 0) + RGB24(0, 0, 1) === RGB24(1, 0, 1)
+
+        # Multiplication
+        @test_throws MethodError cf*cf
+        cf64 = mapc(Float64, cf)
+        @test cf ⋅ cf   === (red(cf)^2 + green(cf)^2 + blue(cf)^2)/3
+        @test cf ⋅ cf64 === (red(cf)*red(cf64) + green(cf)*green(cf64) + blue(cf)*blue(cf64))/3
+        @test cf ⊙ cf   === RGB(red(cf)^2, green(cf)^2, blue(cf)^2)
+        @test cf ⊙ cf64 === RGB(red(cf)*red(cf64), green(cf)*green(cf64), blue(cf)*blue(cf64))
+        c2 = rand(RGB{Float64})
+        rr = cf ⊗ c2
+        @test Matrix(rr) == [red(cf)*red(c2)   red(cf)*green(c2)   red(cf)*blue(c2);
+                             green(cf)*red(c2) green(cf)*green(c2) green(cf)*blue(c2);
+                             blue(cf)*red(c2)  blue(cf)*green(c2)  blue(cf)*blue(c2)]
+        @test +rr === rr
+        @test -rr === RGBRGB(-rr.rr, -rr.gr, -rr.br, -rr.rg, -rr.gg, -rr.bg, -rr.rb, -rr.gb, -rr.bb)
+        @test rr + rr == 2*rr == rr*2
+        @test rr - rr == zero(rr)
+        io = IOBuffer()
+        print(io, N0f8)
+        Tstr = String(take!(io))
+        cfn = RGB{N0f8}(0.1, 0.2, 0.3)
+        show(io, cfn ⊗ cfn)
+        spstr = Base.VERSION >= v"1.5" ? "" : " "
+        @test String(take!(io)) == "RGBRGB{$Tstr}(\n 0.012N0f8  0.02N0f8   0.031N0f8\n 0.02N0f8   0.039N0f8  0.059N0f8\n 0.031N0f8  0.059N0f8  0.09N0f8$spstr)"
     end
 
     @testset "Arithemtic with RGBA" begin
@@ -324,7 +420,7 @@ end
         @test !isfinite(RGBA(0.2, 1, 0.5, Inf))
         @test RGBA(0.2, 1, 0.5, Inf) |> isinf
         @test !isnan(RGBA(0.2, 1, 0.5, Inf))
-        @test abs(RGBA(0.1,0.2,0.3,0.2)) ≈ 0.8
+        @test abs(RGBA(0.1,0.2,0.3,0.2)) === RGBA(0.1,0.2,0.3,0.2)
 
         acu = RGBA{N0f8}[cu]
         acf = RGBA{Float32}[cf]
@@ -372,8 +468,14 @@ end
         @test ARGB32(0.4, 0, 0.2, 0.5) + AGray32(0.4, 0.2) === ARGB32(0.8, 0.4, 0.6, 0.5N0f8+0.2N0f8)
     end
 
+    @testset "Custom RGB arithmetic" begin
+        cf = RatRGB(1//10, 2//10, 3//10)
+        @test cf ⋅ cf   === (Float64(red(cf))^2 + Float64(green(cf))^2 + Float64(blue(cf))^2)/3
+    end
+
     @testset "dotc" begin
         @test dotc(0.2, 0.2) == 0.2^2
+        @test dotc(Int8(3), Int16(6)) === 18
         @test dotc(0.2, 0.3f0) == 0.2*0.3f0
         @test dotc(N0f8(0.2), N0f8(0.3)) == Float32(N0f8(0.2))*Float32(N0f8(0.3))
         @test dotc(Gray{N0f8}(0.2), Gray24(0.3)) == Float32(N0f8(0.2))*Float32(N0f8(0.3))
@@ -396,12 +498,42 @@ end
 
     @testset "Colors issue #326" begin
         A = rand(RGB{N0f8}, 2, 2)
-        if VERSION >= v"1.5"
-            @test_broken @inferred mean(A) == mean(map(c->mapc(FixedPointNumbers.Treduce, c), A))
-        else
-            @test @inferred mean(A) == mean(map(c->mapc(FixedPointNumbers.Treduce, c), A))
+        @test @inferred(mean(A)) == mean(map(c->mapc(FixedPointNumbers.Treduce, c), A))
+    end
+
+    @testset "Equivalence" begin
+        x = 0.4
+        g = Gray(x)
+        c = RGB(g)
+        for p in (0, 1, 2, Inf)
+            @test norm(x, p) == norm(g, p) ≈ norm(c, p)
         end
+        @test dot(x, x) == dot(g, g) ≈ dot(c, c)
+        @test_throws MethodError mapreduce(x->x^2, +, c)   # this risks breaking equivalence & noniterability
+    end
+
+    @testset "varmult" begin
+        cs = [RGB(0.2, 0.3, 0.4), RGB(0.5, 0.3, 0.2)]
+        @test varmult(⋅, cs) ≈ 2*(0.15^2 + 0.1^2)/3    # the /3 is for the 3 color channels, i.e., equivalence
+        @test varmult(⋅, cs; corrected=false) ≈ (0.15^2 + 0.1^2)/3
+        @test varmult(⋅, cs; mean=RGB(0, 0, 0)) ≈ (0.2^2+0.3^2+0.4^2 + 0.5^2+0.3^2+0.2^2)/3
+        @test varmult(⊙, cs) ≈ 2*RGB(0.15^2, 0, 0.1^2)
+        @test Matrix(varmult(⊗, cs)) ≈ 2*[0.15^2 0 -0.1*0.15; 0 0 0; -0.1*0.15 0 0.1^2]
+
+        cs = [RGB(0.1, 0.2,  0.3)  RGB(0.3, 0.5, 0.3);
+              RGB(0.2, 0.21, 0.33) RGB(0.4, 0.51, 0.33);
+              RGB(0.3, 0.22, 0.36) RGB(0.5, 0.52, 0.36)]
+        v1 = RGB(0.1^2, 0.15^2, 0)
+        @test varmult(⊙, cs, dims=2) ≈ 2*[v1, v1, v1]
+        v2 = RGB(0.1^2, 0.01^2, 0.03^2)
+        @test varmult(⊙, cs, dims=1) ≈ [v2 v2]
+    end
+
+    @testset "copy" begin
+        g = Gray{N0f8}(0.2)
+        @test copy(g) === g
+        c = RGB(0.1, 0.2, 0.3)
+        @test copy(c) === c
     end
-end
 
 end