From 08dbd9e385b4bcfd49df1ebe53a7025e307fb9dc Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Tue, 17 Mar 2020 07:48:41 -0500
Subject: [PATCH 1/9] Redesign color arithmetic

This notably defines 3 multiplication operators for RGB colors.
It also un-defines `abs2`, because how that should work is a bit ambiguous.
Finally, it defines a new `varmult` function, which allows one to
compute variance using a specific multiplication operator.

There are some compatibility definitions for current releases of ColorTypes.

Co-authored-by: Johnny Chen <johnnychen94@hotmail.com>
---
 Project.toml            |  11 +-
 README.md               |  81 +++++++-
 src/ColorVectorSpace.jl | 411 ++++++++++++++++++++++++----------------
 test/runtests.jl        | 133 ++++++++-----
 4 files changed, 419 insertions(+), 217 deletions(-)

diff --git a/Project.toml b/Project.toml
index 25c2c05..75fab88 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ColorVectorSpace"
 uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
-version = "0.8.7"
+version = "1.0.0"
 
 [deps]
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
@@ -9,21 +9,20 @@ 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"
+SpecialFunctions = "0.7, 0.8, 0.9, 0.10"
+TensorCore = "0.1"
 julia = "1"
 
 [extras]
 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 = ["Statistics", "LinearAlgebra", "Test"]
diff --git a/README.md b/README.md
index 13de965..8480ef4 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,74 @@ 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))`.
+
 ## 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..f519bb0 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -1,20 +1,17 @@
 module ColorVectorSpace
 
 using Colors, 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, abs2, 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
+# 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,38 @@ 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))
+## Version compatibility with ColorTypes
 
-## 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)))
+if !hasmethod(zero, (Type{AGray{N0f8}},))
+    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
-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)))
+if !hasmethod(one, (Gray{N0f8},))
+    Base.one(::Type{C}) where {C<:AbstractGray}    = C(1)
+    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
 
-for f in (:mod, :rem, :mod1)
-    @eval $f(x::Gray, m::Gray) = Gray($f(gray(x), gray(m)))
+if !hasmethod(oneunit, (Type{AGray{N0f8}},))
+    oneunit(::Type{C}) where {C<:TransparentGray} = C(1,1)
+    oneunit(::Type{C}) where {C<:AbstractRGB}     = C(1,1,1)
+    oneunit(::Type{C}) where {C<:TransparentRGB}  = C(1,1,1,1)
+    oneunit(p::Colorant) = oneunit(typeof(p))
 end
 
 # Real values are treated like grays
@@ -74,24 +60,27 @@ 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))
+multype(a::Colorant, b::Colorant) = coltype(multype(eltype(a),eltype(b)))
+sumtype(a::Colorant, b::Colorant) = coltype(sumtype(eltype(a),eltype(b)))
+divtype(a::Colorant, b::Colorant) = coltype(divtype(eltype(a),eltype(b)))
+powtype(a::Colorant, b::Colorant) = coltype(powtype(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)
 
@@ -119,6 +108,49 @@ 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))
+
+## Math on colors and color types
+
+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))
+
+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 +191,33 @@ 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)       # if we make abs2 return a scalar, this could be confusing
+# abs2(c::MathTypes) = mapc(abs2, c)   # or mapreducec(abs2, +, float(zero(eltype(c))), c); or is it better to make this undefined?
+norm(c::MathTypes, p::Real=2) = (cc = channels(c); norm(cc, p)/(p == 0 ? length(cc) : length(cc)^(1/p)))
+# norm1(c::MathTypes)   =      mapreducec(abs∘float,  +,   float(zero(eltype(c))), c)
+# norm2(c::MathTypes)   = sqrt(mapreducec(abs2∘float, +,   float(zero(eltype(c))), c))
+# normInf(c::MathTypes) =      mapreducec(abs,        max, zero(eltype(c)),        c)
+
+# function Base.rtoldefault(::Union{C1,Type{C1}}, ::Union{C2,Type{C2}}, atol::Real) where {C1<:MathTypes,C2<:MathTypes}
+#     T1, T2 = eltype(C1), eltype(C2)
+#     @show T1, T2
+#     return Base.rtoldefault(eltype(C1), eltype(C2), atol)
+# end
+
+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 +245,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 +273,20 @@ 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)...)
+div(a::C, b::C) where C<:AbstractGray = base_color_type(C)(div(gray(a), gray(b)))
+div(a::AbstractGray, b::AbstractGray) = div(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,42 +296,35 @@ 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)
+atan(x::AbstractGray, y::AbstractGray)  = atan(gray(x), gray(y))
+hypot(x::AbstractGray, y::AbstractGray) = hypot(gray(x), gray(y))
+
+if !hasmethod(<, Tuple{AbstractGray,AbstractGray})  # 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
 if !hasmethod(isless, Tuple{AbstractGray,AbstractGray})  # this was moved to ColorTypes 0.10
     isless(g1::AbstractGray, g2::AbstractGray) = isless(gray(g1), gray(g2))
 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)
+# 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
+# end
 
 dotc(x::T, y::T) where {T<:AbstractGray} = acc(gray(x))*acc(gray(y))
 dotc(x::AbstractGray, y::AbstractGray) = dotc(promote(x, y)...)
@@ -299,12 +335,6 @@ float(::Type{T}) where {T<:Gray} = typeof(float(zero(T)))
 (+)(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 +346,109 @@ 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
+
++(a::RGBRGB) = a
+-(a::RGBRGB) = RGB(-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
+
+## 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..ff7622d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,18 +1,14 @@
-module ColorVectorSpaceTests
-
 using LinearAlgebra, Statistics
-using ColorVectorSpace, Colors, FixedPointNumbers, StatsBase
+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)))
@@ -23,6 +19,18 @@ end
 
 @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)
@@ -41,29 +49,30 @@ 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)
         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(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)
+        # @test @inferred(abs2(ccmp)) === Gray(0.2f0^2)
+        @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)
+        # @test @inferred(sum(abs2, ccmp)) == Gray(0.2f0^2)
         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)
@@ -74,13 +83,12 @@ 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)
 
         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 +110,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 sum(abs2, [cf, ccmp]) ≈ Gray(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 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)
 
@@ -248,9 +257,9 @@ 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 sum(abs2, RGB(0.1,0.2,0.3)) == RGB(0.1^2,0.2^2,0.3^2)
+        @test norm(RGB(0.1,0.2,0.3)) ≈ sqrt(0.14)/sqrt(3)
 
         acu = RGB{N0f8}[cu]
         acf = RGB{Float32}[cf]
@@ -284,6 +293,18 @@ 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})
+        @test Matrix(cf ⊗ c2) == [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)]
     end
 
     @testset "Arithemtic with RGBA" begin
@@ -324,7 +345,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]
@@ -396,12 +417,34 @@ 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)
+    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
-end
 
 end

From f2dfaaddc46e595a426e2d6282d717acee3bedee Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 06:48:58 -0600
Subject: [PATCH 2/9] Clean cruft, dep only on ColorTypes

A single test depends on Colors, otherwise this package is now
independent of Colors.
---
 Project.toml            |  9 ++++-----
 src/ColorVectorSpace.jl | 10 +++-------
 test/runtests.jl        |  7 +++----
 3 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/Project.toml b/Project.toml
index 75fab88..d6aaeea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ version = "1.0.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"
@@ -12,17 +11,17 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 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"
+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"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Statistics", "LinearAlgebra", "Test"]
+test = ["Colors", "Statistics", "LinearAlgebra", "Test"]
diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index f519bb0..6ed0e59 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -1,13 +1,13 @@
 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, float,
+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
@@ -204,12 +204,8 @@ 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::MathTypes) = mapc(abs, c)       # if we make abs2 return a scalar, this could be confusing
-# abs2(c::MathTypes) = mapc(abs2, c)   # or mapreducec(abs2, +, float(zero(eltype(c))), c); or is it better to make this undefined?
+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)))
-# norm1(c::MathTypes)   =      mapreducec(abs∘float,  +,   float(zero(eltype(c))), c)
-# norm2(c::MathTypes)   = sqrt(mapreducec(abs2∘float, +,   float(zero(eltype(c))), c))
-# normInf(c::MathTypes) =      mapreducec(abs,        max, zero(eltype(c)),        c)
 
 # function Base.rtoldefault(::Union{C1,Type{C1}}, ::Union{C2,Type{C2}}, atol::Real) where {C1<:MathTypes,C2<:MathTypes}
 #     T1, T2 = eltype(C1), eltype(C2)
diff --git a/test/runtests.jl b/test/runtests.jl
index ff7622d..b50c419 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -58,12 +58,11 @@ end
         @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)
-        # @test @inferred(abs2(ccmp)) === Gray(0.2f0^2)
+        @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)
-        # @test @inferred(sum(abs2, ccmp)) == Gray(0.2f0^2)
         cu = Gray{N0f8}(0.1)
         @test @inferred(2*cu) === cu*2 === Gray(2*gray(cu))
         @test @inferred(2.0f0*cu) === cu*2.0f0 === Gray(2.0f0*gray(cu))
@@ -110,7 +109,6 @@ 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]) ≈ Gray(0.05f0)
 
         @test gray(0.8) === 0.8
 
@@ -258,7 +256,8 @@ end
         @test isinf(RGB(1, Inf, 0.5))
         @test !isnan(RGB(1, Inf, 0.5))
         @test abs(RGB(0.1,0.2,0.3)) == RGB(0.1,0.2,0.3)
-        # @test sum(abs2, RGB(0.1,0.2,0.3)) == RGB(0.1^2,0.2^2,0.3^2)
+        @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)
 
         acu = RGB{N0f8}[cu]

From 04685f876f061d66548d131a4173ebcbf0f56ddc Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 09:02:39 -0600
Subject: [PATCH 3/9] explain

---
 README.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/README.md b/README.md
index 8480ef4..6a047cf 100644
--- a/README.md
+++ b/README.md
@@ -40,6 +40,9 @@ 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
 

From 0cf0121ea75a1c560bf8abe12322d5bd041678a0 Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 09:06:39 -0600
Subject: [PATCH 4/9] Improve test coverage

---
 src/ColorVectorSpace.jl | 17 +++-----
 test/runtests.jl        | 89 ++++++++++++++++++++++++++++++++++++-----
 2 files changed, 85 insertions(+), 21 deletions(-)

diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 6ed0e59..8d0ba82 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -10,7 +10,7 @@ import Base: ==, +, -, *, /, ^, <, ~
 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
+    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,
@@ -70,7 +70,6 @@ powtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)^zero(B)))
 multype(a::Colorant, b::Colorant) = coltype(multype(eltype(a),eltype(b)))
 sumtype(a::Colorant, b::Colorant) = coltype(sumtype(eltype(a),eltype(b)))
 divtype(a::Colorant, b::Colorant) = coltype(divtype(eltype(a),eltype(b)))
-powtype(a::Colorant, b::Colorant) = coltype(powtype(eltype(a),eltype(b)))
 
 coltype(::Type{T}) where {T<:Fractional} = T
 coltype(::Type{T}) where {T<:Number}     = floattype(T)
@@ -126,8 +125,6 @@ _nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:TransparentRGB} = (x = con
 
 
 ## 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}
@@ -271,8 +268,6 @@ _mean_promote(x::MathTypes, y::MathTypes) = mapc(FixedPointNumbers.Treduce, y)
 (-)(c::TransparentGray) = typeof(c)(-gray(c),-alpha(c))
 (/)(a::C, b::C) where C<:AbstractGray = base_color_type(C)(gray(a)/gray(b))
 (/)(a::AbstractGray, b::AbstractGray) = /(promote(a, b)...)
-div(a::C, b::C) where C<:AbstractGray = base_color_type(C)(div(gray(a), gray(b)))
-div(a::AbstractGray, b::AbstractGray) = div(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)
@@ -295,14 +290,11 @@ min(a::AbstractGray, b::Number) = min(promote(a,b)...)
 atan(x::AbstractGray, y::AbstractGray)  = atan(gray(x), gray(y))
 hypot(x::AbstractGray, y::AbstractGray) = hypot(gray(x), gray(y))
 
-if !hasmethod(<, Tuple{AbstractGray,AbstractGray})  # planned for ColorTypes 0.11
+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
-if !hasmethod(isless, Tuple{AbstractGray,AbstractGray})  # this was moved to ColorTypes 0.10
-    isless(g1::AbstractGray, g2::AbstractGray) = isless(gray(g1), gray(g2))
-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))
@@ -384,8 +376,11 @@ function Base.show(io::IO, p::RGBRGB)
     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) = RGB(-a.rr, -a.gr, -a.br, -a.rg, -a.gg, -a.bg, -a.rb, -a.gb, -a.bb)
+-(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)
diff --git a/test/runtests.jl b/test/runtests.jl
index b50c419..d9e54da 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-using LinearAlgebra, Statistics
+using LinearAlgebra, Statistics, SpecialFunctions
 using ColorVectorSpace, Colors, FixedPointNumbers
 
 using Test
@@ -17,6 +17,15 @@ 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
@@ -51,13 +60,21 @@ end
         cf = Gray{Float32}(0.1)
         @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 @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)
@@ -85,6 +102,14 @@ end
         @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]
@@ -158,6 +183,15 @@ end
         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))
@@ -165,18 +199,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)
@@ -301,9 +355,18 @@ end
         @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})
-        @test Matrix(cf ⊗ c2) == [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)]
+        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()
+        cfn = RGB{N0f8}(0.1, 0.2, 0.3)
+        show(io, cfn ⊗ cfn)
+        @test String(take!(io)) == "RGBRGB{N0f8}(\n 0.012N0f8  0.02N0f8   0.031N0f8\n 0.02N0f8   0.039N0f8  0.059N0f8\n 0.031N0f8  0.059N0f8  0.09N0f8)"
     end
 
     @testset "Arithemtic with RGBA" begin
@@ -392,6 +455,11 @@ 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(0.2, 0.3f0) == 0.2*0.3f0
@@ -427,6 +495,7 @@ end
             @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

From 9ff7e7e3076f6242a2a19d61250c4d950683cae9 Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 10:19:55 -0600
Subject: [PATCH 5/9] improve coverage and register err hint

---
 src/ColorVectorSpace.jl | 27 ++++++++++++---------------
 test/runtests.jl        | 24 +++++++++++++++++++++---
 2 files changed, 33 insertions(+), 18 deletions(-)

diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 8d0ba82..37c6e63 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -33,28 +33,20 @@ MathTypes{T,C} = Union{AbstractRGB{T},TransparentRGB{C,T},AbstractGray{T},Transp
 
 ## Version compatibility with ColorTypes
 
-if !hasmethod(zero, (Type{AGray{N0f8}},))
+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
 
-if !hasmethod(one, (Gray{N0f8},))
-    Base.one(::Type{C}) where {C<:AbstractGray}    = C(1)
+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
 
-if !hasmethod(oneunit, (Type{AGray{N0f8}},))
-    oneunit(::Type{C}) where {C<:TransparentGray} = C(1,1)
-    oneunit(::Type{C}) where {C<:AbstractRGB}     = C(1,1,1)
-    oneunit(::Type{C}) where {C<:TransparentRGB}  = C(1,1,1,1)
-    oneunit(p::Colorant) = oneunit(typeof(p))
-end
-
 # Real values are treated like grays
 if !hasmethod(gray, (Number,))
     ColorTypes.gray(x::Real) = x
@@ -67,9 +59,7 @@ 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) = coltype(multype(eltype(a),eltype(b)))
 sumtype(a::Colorant, b::Colorant) = coltype(sumtype(eltype(a),eltype(b)))
-divtype(a::Colorant, b::Colorant) = coltype(divtype(eltype(a),eltype(b)))
 
 coltype(::Type{T}) where {T<:Fractional} = T
 coltype(::Type{T}) where {T<:Number}     = floattype(T)
@@ -89,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))
@@ -317,8 +306,6 @@ end
 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)...)
@@ -435,6 +422,16 @@ function varmult(op, itr; corrected::Bool=true, dims=:, mean=Statistics.mean(itr
     return v / (corrected ? max(1, n-1) : max(1, n))
 end
 
+function __init__()
+    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
+             print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
+        end
+    end
+end
+
 ## Precompilation
 
 if Base.VERSION >= v"1.4.2"
diff --git a/test/runtests.jl b/test/runtests.jl
index d9e54da..f3c84a1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -56,6 +56,10 @@ ColorTypes.blue(c::RatRGB)  = c.b
         end
     end
 
+    @testset "traits" begin
+        @test floattype(Gray{N0f8}) === Gray{float(N0f8)}
+    end
+
     @testset "Arithmetic with Gray" begin
         cf = Gray{Float32}(0.1)
         @test @inferred(+cf) === cf
@@ -92,6 +96,7 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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))
@@ -142,6 +147,7 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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
@@ -176,10 +182,11 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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)
 
@@ -314,6 +321,8 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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]
         @test typeof(acu+acf) == Vector{RGB{Float32}}
@@ -331,6 +340,7 @@ ColorTypes.blue(c::RatRGB)  = c.b
 
         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)
@@ -462,6 +472,7 @@ ColorTypes.blue(c::RatRGB)  = c.b
 
     @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))
@@ -515,4 +526,11 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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

From ec4457b4e83035031ba373e996eb48ad69d94aa4 Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 10:23:47 -0600
Subject: [PATCH 6/9] Protect err hint registration

---
 src/ColorVectorSpace.jl | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 37c6e63..2a94607 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -423,11 +423,13 @@ function varmult(op, itr; corrected::Bool=true, dims=:, mean=Statistics.mean(itr
 end
 
 function __init__()
-    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
-             print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
+    if 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
+                print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
+            end
         end
     end
 end

From 610f1482e44b09728e5ea495d14b91b8be292b4d Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 10:34:29 -0600
Subject: [PATCH 7/9] Fix Julia 1.0

---
 src/ColorVectorSpace.jl | 2 +-
 test/runtests.jl        | 5 ++++-
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 2a94607..58feb8e 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -423,7 +423,7 @@ function varmult(op, itr; corrected::Bool=true, dims=:, mean=Statistics.mean(itr
 end
 
 function __init__()
-    if isdefined(Base.Experimental, :register_error_hint)
+    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.
diff --git a/test/runtests.jl b/test/runtests.jl
index f3c84a1..c707eac 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -374,9 +374,12 @@ ColorTypes.blue(c::RatRGB)  = c.b
         @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)
-        @test String(take!(io)) == "RGBRGB{N0f8}(\n 0.012N0f8  0.02N0f8   0.031N0f8\n 0.02N0f8   0.039N0f8  0.059N0f8\n 0.031N0f8  0.059N0f8  0.09N0f8)"
+        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

From 24cdcb51016da7049dccd7577cff8a4253a95e4c Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 10:40:58 -0600
Subject: [PATCH 8/9] Call this 0.9 "in an abundance of caution"

I *hate* that phrase...
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d6aaeea..ceaf35d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ColorVectorSpace"
 uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
-version = "1.0.0"
+version = "0.9.0"
 
 [deps]
 ColorTypes = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"

From d61e3638171b9e660de325dc9fea197bb71f24be Mon Sep 17 00:00:00 2001
From: Tim Holy <tim.holy@gmail.com>
Date: Sat, 16 Jan 2021 13:12:01 -0600
Subject: [PATCH 9/9] cruft

---
 src/ColorVectorSpace.jl | 28 ++++------------------------
 1 file changed, 4 insertions(+), 24 deletions(-)

diff --git a/src/ColorVectorSpace.jl b/src/ColorVectorSpace.jl
index 58feb8e..0bb46ee 100644
--- a/src/ColorVectorSpace.jl
+++ b/src/ColorVectorSpace.jl
@@ -102,8 +102,6 @@ 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))
 
-## Math on colors and color types
-
 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))
@@ -112,7 +110,6 @@ _nan(::Type{T}, ::Type{C}) where {T<:AbstractFloat,C<:AbstractRGB} = (x = conver
 _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)
@@ -121,6 +118,8 @@ function Base.reduce_empty(::typeof(Base.add_sum), ::Type{T}) where {T<:MathType
     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
@@ -137,6 +136,7 @@ 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)
 
@@ -193,12 +193,6 @@ isinf(c::Colorant) = mapreducec(isinf, |, false, c)
 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)))
 
-# function Base.rtoldefault(::Union{C1,Type{C1}}, ::Union{C2,Type{C2}}, atol::Real) where {C1<:MathTypes,C2<:MathTypes}
-#     T1, T2 = eltype(C1), eltype(C2)
-#     @show T1, T2
-#     return Base.rtoldefault(eltype(C1), eltype(C2), atol)
-# end
-
 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
@@ -289,20 +283,6 @@ if !hasmethod(isless, Tuple{AbstractGray,Real})  # planned for ColorTypes 0.11
     isless(r::Real, c::AbstractGray) = isless(r, gray(c))
 end
 
-# 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
-# end
-
 dotc(x::T, y::T) where {T<:AbstractGray} = acc(gray(x))*acc(gray(y))
 dotc(x::AbstractGray, y::AbstractGray) = dotc(promote(x, y)...)
 
@@ -428,7 +408,7 @@ function __init__()
             if exc.f === _color_rettype && length(argtypes) >= 2
                 # Color is not necessary, this is just to show it's possible.
                 A, B = argtypes
-                print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
+                A !== B && print(io, "\nIn binary operation with $A and $B, the return type is ambiguous")
             end
         end
     end