From bcb46c06c796006bf17a753029042e6ba4aa3d3d Mon Sep 17 00:00:00 2001
From: Mike Reed <reed@google.com>
Date: Mon, 23 Mar 2020 17:51:01 -0400
Subject: [PATCH] Add approx_pow/log2/pow2 to SkVM builder

... in support of programs for colorspacexforms

Change-Id: I72ace09f10511ef8994038a4af3feab8bc1a299e
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/278466
Commit-Queue: Mike Reed <reed@google.com>
Reviewed-by: Mike Klein <mtklein@google.com>
---
 src/core/SkVM.cpp  | 29 ++++++++++++++++++++++
 src/core/SkVM.h    |  4 ++++
 tests/SkVMTest.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 93 insertions(+)

diff --git a/src/core/SkVM.cpp b/src/core/SkVM.cpp
index b59bd21ce3..61d16525e9 100644
--- a/src/core/SkVM.cpp
+++ b/src/core/SkVM.cpp
@@ -926,6 +926,35 @@ namespace skvm {
         return {this->push(Op::sqrt_f32, x.id,NA,NA)};
     }
 
+    // See http://www.machinedlearnings.com/2011/06/fast-approximate-logarithm-exponential.html.
+    F32 Builder::approx_log2(F32 x) {
+        // e - 127 is a fair approximation of log2(x) in its own right...
+        F32 e = mul(to_f32(bit_cast(x)), splat(1.0f / (1<<23)));
+
+        // ... but using the mantissa to refine its error is _much_ better.
+        F32 m = bit_cast(bit_or(bit_and(bit_cast(x), splat(0x007fffff)), splat(0x3f000000)));
+        F32 approx = sub(e,      splat(124.225514990f));
+            approx = sub(approx, mul(splat(1.498030302f), m));
+            approx = sub(approx, div(splat(1.725879990f), add(splat(0.3520887068f), m)));
+
+        return approx;
+    }
+
+    F32 Builder::approx_pow2(F32 x) {
+        F32 f = fract(x);
+        F32 approx = add(x, splat(121.274057500f));
+            approx = sub(approx, mul(splat( 1.490129070f), f));
+            approx = add(approx, div(splat(27.728023300f), sub(splat(4.84252568f), f)));
+
+        return bit_cast(round(mul(splat(1.0f * (1<<23)), approx)));
+    }
+
+    F32 Builder::approx_powf(F32 x, F32 y) {
+        auto is_x = bit_or(eq(x, splat(0.0f)),
+                           eq(x, splat(1.0f)));
+        return select(is_x, x, approx_pow2(mul(approx_log2(x), y)));
+    }
+
     F32 Builder::min(F32 x, F32 y) {
         float X,Y;
         if (this->allImm(x.id,&X, y.id,&Y)) { return this->splat(std::min(X,Y)); }
diff --git a/src/core/SkVM.h b/src/core/SkVM.h
index 1ec4f39b9a..1104bcea15 100644
--- a/src/core/SkVM.h
+++ b/src/core/SkVM.h
@@ -454,6 +454,10 @@ namespace skvm {
         F32 mad(F32 x, F32 y, F32 z) { return this->add(this->mul(x,y), z); }
         F32 sqrt(F32 x);
 
+        F32 approx_log2(F32);
+        F32 approx_pow2(F32);
+        F32 approx_powf(F32 base, F32 exp);
+
         F32 negate(F32 x) {
             return sub(splat(0.0f), x);
         }
diff --git a/tests/SkVMTest.cpp b/tests/SkVMTest.cpp
index 0661ac30a3..3310b176e5 100644
--- a/tests/SkVMTest.cpp
+++ b/tests/SkVMTest.cpp
@@ -1791,3 +1791,63 @@ DEF_TEST(SkVM_Assembler, r) {
         0x20,0x00,0x02,0x4e,
     });
 }
+
+DEF_TEST(SkVM_approx_math, r) {
+    auto eval = [](int N, float values[], auto fn) {
+        skvm::Builder b;
+        skvm::Arg inout  = b.varying<float>();
+
+        b.storeF(inout, fn(&b, b.loadF(inout)));
+
+        b.done().eval(N, values);
+    };
+
+    auto compare = [r](int N, const float values[], const float expected[]) {
+        for (int i = 0; i < N; ++i) {
+            REPORTER_ASSERT(r, SkScalarNearlyEqual(values[i], expected[i], 0.001f));
+        }
+    };
+
+    // log2
+    {
+        float values[] = {0.25f, 0.5f, 1, 2, 4, 8};
+        constexpr int N = SK_ARRAY_COUNT(values);
+        eval(N, values, [](skvm::Builder* b, skvm::F32 v) {
+            return b->approx_log2(v);
+        });
+        const float expected[] = {-2, -1, 0, 1, 2, 3};
+        compare(N, values, expected);
+    }
+
+    // pow2
+    {
+        float values[] = {-2, -1, 0, 1, 2, 3};
+        constexpr int N = SK_ARRAY_COUNT(values);
+        eval(N, values, [](skvm::Builder* b, skvm::F32 v) {
+            return b->approx_pow2(v);
+        });
+        const float expected[] = {0.25f, 0.5f, 1, 2, 4, 8};
+        compare(N, values, expected);
+    }
+
+    // powf -- x^0.5
+    {
+        float bases[] = {0, 1, 4, 9, 16};
+        constexpr int N = SK_ARRAY_COUNT(bases);
+        eval(N, bases, [](skvm::Builder* b, skvm::F32 base) {
+            return b->approx_powf(base, b->splat(0.5f));
+        });
+        const float expected[] = {0, 1, 2, 3, 4};
+        compare(N, bases, expected);
+    }
+    // powf -- 3^x
+    {
+        float exps[] = {-2, -1, 0, 1, 2};
+        constexpr int N = SK_ARRAY_COUNT(exps);
+        eval(N, exps, [](skvm::Builder* b, skvm::F32 exp) {
+            return b->approx_powf(b->splat(3.0f), exp);
+        });
+        const float expected[] = {1/9.0f, 1/3.0f, 1, 3, 9};
+        compare(N, exps, expected);
+    }
+}