Improve speed and mean accuracy of pow12_5

`pow12_5` is used in the `RGB{Float64}`-->`XYZ` conversion.
This ensures sufficient accuracy (~2 ULP) in the range of [0.05, 1],
which is typically required for the conversion.

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.485910662571346
 ```
 
 ```julia

src/utilities.jl

@@ -14,17 +14,28 @@ 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"
+function pow12_5(x::Float64)
+    # x^0.4
+    t1 = @evalpoly(@fastmath(min(x, 1.6)),
+        0.17423327702711236,
+        3.035155441155984,  -10.46460992384518,  28.338924866372604, -48.80965034630853,
+        51.970906997398615, -33.061334642907866, 11.494242702033949, -1.678011428260117)
+    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,8 @@ 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(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"