Improve speed and mean accuracy of pow12_5

docs/src/colordifferences.md

@@ -10,7 +10,7 @@ julia> colordiff(colorant"red", colorant"blue")
 52.88136782250768
 
 julia> colordiff(HSV(0, 0.75, 0.5), HSL(0, 0.75, 0.5))
-19.48591066257135
+19.48591066257134
 ```
 
 ```julia

src/conversions.jl

@@ -278,9 +278,11 @@ cnvt(::Type{HSI{T}}, c::Color) where {T} = cnvt(HSI{T}, convert(RGB{T}, c))
 # Everything to XYZ
 # -----------------
 
-function invert_srgb_compand(v::Fractional)
+@inline function invert_srgb_compand(v)
+    F = typeof(0.5 * v)
+    vf = F(v)
     # `pow12_5` is an optimized function to get `x^2.4`
-    v <= 0.04045 ? 1/12.92 * v : pow12_5(1/1.055 * (v + 0.055))
+    vf > F(0.04045) ? pow12_5(muladd(F(1000/1055), vf, F(55/1055))) : F(100/1292) * vf
 end
 
 # lookup table for `N0f8` (the extra two elements are for `Float32` splines)
@@ -306,7 +308,9 @@ function invert_srgb_compand(v::Float32)
 end
 
 function cnvt(::Type{XYZ{T}}, c::AbstractRGB) where T
-    r, g, b = invert_srgb_compand(red(c)), invert_srgb_compand(green(c)), invert_srgb_compand(blue(c))
+    r = invert_srgb_compand(red(c))
+    g = invert_srgb_compand(green(c))
+    b = invert_srgb_compand(blue(c))
     XYZ{T}(0.4124564*r + 0.3575761*g + 0.1804375*b,
            0.2126729*r + 0.7151522*g + 0.0721750*b,
            0.0193339*r + 0.1191920*g + 0.9503041*b)

src/utilities.jl

@@ -14,17 +14,27 @@ end
 # mod6 supports the input `x` in [-2^28, 2^29]
 mod6(::Type{T}, x::Int32) where T = unsafe_trunc(T, x - 6 * ((widemul(x, 0x2aaaaaaa) + Int64(0x20000000)) >> 0x20))
 
-# TODO: move `pow7` from "src/differences.jl" to here
-
 pow3_4(x) = (y = @fastmath(sqrt(x)); y*@fastmath(sqrt(y))) # x^(3/4)
 
 # `pow5_12` is called from `srgb_compand`.
 # `cbrt` generates a function call, so there is little benefit of `@fastmath`.
 pow5_12(x) = pow3_4(x) / cbrt(x) # 5/12 == 1/2 + 1/4 - 1/3 == 3/4 - 1/3
 
 # `pow12_5` is called from `invert_srgb_compand`.
-# x^y ≈ exp(y*log(x)) ≈ exp2(y*log2(y)); the middle form is faster
-@noinline pow12_5(x) = x^2 * exp(0.4 * log(x)) # 12/5 == 2.4 == 2 + 0.4
+pow12_5(x) = pow12_5(Float64(x))
+pow12_5(x::BigFloat) = x^big"2.4"
+@inline function pow12_5(x::Float64)
+    # x^0.4
+    t1 = @evalpoly(@fastmath(min(x, 1.75)), 0.24295462640373672,
+        1.7489099720303518, -1.9919942887850166, 1.3197188815160004, -0.3257258790067756)
+    t2 = muladd(2/5, muladd(x / t1^2, @fastmath(sqrt(t1)), -t1), t1) # Newton's method
+    t3 = muladd(2/5, muladd(x / t2^2, @fastmath(sqrt(t2)), -t2), t2)
+    t4 = muladd(2/5, muladd(x / t3^2, @fastmath(sqrt(t3)), -t3), t3)
+    # x^0.4 * x^2
+    rx = reinterpret(Float64, reinterpret(UInt64, x) & 0xffffffff_f8000000) # hi
+    e = x - rx # lo
+    muladd(t4, rx^2, t4 * (rx + rx + e) * e)
+end
 
 pow7(x) = (y = x*x*x; y*y*x)
 pow7(x::Integer) = pow7(Float64(x)) # avoid overflow

test/utilities.jl

@@ -19,6 +19,9 @@ using InteractiveUtils # for `subtypes`
         @warn "Optimization technique in `pow5_12` may have the opposite effect."
     end
 
+    @test Colors.pow12_5(0.6) ≈ Colors.pow12_5(big"0.6") atol=1e-6
+    @test Colors.pow12_5(0.6N0f16) ≈ Colors.pow12_5(big"0.6") atol=1e-6
+
     @testset "hex" begin
         @test hex(RGB(1,0.5,0)) == "FF8000"
         @test hex(RGBA(1,0.5,0,0.25)) == "FF800040"