[skvm] approx_atan2

I propose landing this, but then pause on extending math functions until
we develop a clearer migration story for sksl.

Change-Id: Id42ec37071da058e6e7809abe1ed0570d48df8e7
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/283229
Reviewed-by: Mike Klein <mtklein@google.com>
Commit-Queue: Mike Reed <reed@google.com>

src/core/SkVM.cpp

@@ -956,26 +956,58 @@ namespace skvm {
          return x;
      }
 
-    /*  Use symmetry (atan is odd)
-     *  Use identity atan(x) = pi/2 - atan(1/x) for x > 1
-     *  Use 4th order polynomial approximation from https://arachnoid.com/polysolve/
+    /*  Use 4th order polynomial approximation from https://arachnoid.com/polysolve/
      *      with 129 values of x,atan(x) for x:[0...1]
+     *  This only works for 0 <= x <= 1
+     */
+    static F32 approx_atan_unit(F32 x) {
+        // for now we might be given NaN... let that through
+        x->assert_true((x != x) | ((x >= 0) & (x <= 1)));
+        return poly(x, 0.14130025741326729f,
+                      -0.34312835980675116f,
+                      -0.016172900528248768f,
+                       1.0037696976200385f,
+                      -0.00014758242182738969f);
+    }
+
+    /*  Use identity atan(x) = pi/2 - atan(1/x) for x > 1
      */
     F32 Builder::approx_atan(F32 x) {
         I32 neg = (x < 0.0f);
         x = select(neg, -x, x);
         I32 flip = (x > 1.0f);
         x = select(flip, 1/x, x);
-        x = poly(x, 0.14130025741326729f,
-                   -0.34312835980675116f,
-                   -0.016172900528248768f,
-                    1.0037696976200385f,
-                   -0.00014758242182738969f);
+        x = approx_atan_unit(x);
         x = select(flip, SK_ScalarPI/2 - x, x);
         x = select(neg, -x, x);
         return x;
     }
 
+    /*  Use identity atan(x) = pi/2 - atan(1/x) for x > 1
+     *  By swapping y,x to ensure the ratio is <= 1, we can safely call atan_unit()
+     *  which avoids a 2nd divide instruction if we had instead called atan().
+     */
+    F32 Builder::approx_atan(F32 y0, F32 x0) {
+
+        I32 flip = (abs(y0) > abs(x0));
+        F32 y = select(flip, x0, y0);
+        F32 x = select(flip, y0, x0);
+        F32 arg = y/x;
+
+        I32 neg = (arg < 0.0f);
+        arg = select(neg, -arg, arg);
+
+        F32 r = approx_atan_unit(arg);
+        r = select(flip, SK_ScalarPI/2 - r, r);
+        r = select(neg, -r, r);
+
+        // handle quadrant distinctions
+        r = select((y0 >= 0) & (x0  < 0), r + SK_ScalarPI, r);
+        r = select((y0  < 0) & (x0 <= 0), r - SK_ScalarPI, r);
+        // Note: we don't try to handle 0,0 or infinities (yet)
+        return r;
+    }
+
     F32 Builder::min(F32 x, F32 y) {
         if (float X,Y; this->allImm(x.id,&X, y.id,&Y)) { return splat(std::min(X,Y)); }
         return {this, this->push(Op::min_f32, x.id, y.id)};

src/core/SkVM.h

@@ -572,6 +572,7 @@ namespace skvm {
         F32 approx_asin(F32 x);
         F32 approx_acos(F32 x) { return sub(SK_ScalarPI/2, approx_asin(x)); }
         F32 approx_atan(F32 x);
+        F32 approx_atan(F32 y, F32 x);
 
         F32 lerp(F32  lo, F32  hi, F32  t) { return mad(sub(hi, lo), t, lo); }
         F32 lerp(F32a lo, F32a hi, F32a t) { return lerp(_(lo), _(hi), _(t)); }
@@ -918,6 +919,7 @@ namespace skvm {
     static inline F32 approx_asin(F32 x) { return x->approx_asin(x); }
     static inline F32 approx_acos(F32 x) { return x->approx_acos(x); }
     static inline F32 approx_atan(F32 x) { return x->approx_atan(x); }
+    static inline F32 approx_atan(F32 y, F32 x) { return x->approx_atan(y, x); }
 
     static inline F32 clamp01(F32 x) { return x->clamp01(x); }
     static inline F32     abs(F32 x) { return x->    abs(x); }

tests/SkVMTest.cpp

@@ -1787,8 +1787,26 @@ DEF_TEST(SkVM_approx_math, r) {
 
         if (err > tolerance) {
     //        SkDebugf("arg %g, expected %g, actual %g\n", arg, expected, actual);
+            REPORTER_ASSERT(r, true);
+        }
+        return err;
+    };
+
+    auto test2 = [r](float arg0, float arg1, float expected, float tolerance, auto prog) {
+        skvm::Builder b;
+        skvm::Arg in0  = b.varying<float>();
+        skvm::Arg in1  = b.varying<float>();
+        skvm::Arg out  = b.varying<float>();
+        b.storeF(out, prog(b.loadF(in0), b.loadF(in1)));
+        float actual;
+        b.done().eval(1, &arg0, &arg1, &actual);
+
+        float err = std::abs(actual - expected);
+
+        if (err > tolerance) {
+    //        SkDebugf("[%g, %g]: expected %g, actual %g\n", arg0, arg1, expected, actual);
+            REPORTER_ASSERT(r, true);
         }
-        REPORTER_ASSERT(r, err <= tolerance);
         return err;
     };
 
@@ -1839,12 +1857,21 @@ DEF_TEST(SkVM_approx_math, r) {
         if (0) { SkDebugf("asin error %g\n", err); }
 
         err = 0;
-        for (float x = -10; x <= 10; x += 0.1f) {
+        for (float x = -10; x <= 10; x += 1.0f/16) {
             err += test(x, atan(x), tol, [](skvm::F32 x) {
                 return approx_atan(x);
             });
         }
         if (0) { SkDebugf("atan error %g\n", err); }
+
+        for (float y = -3; y <= 3; y += 1) {
+            for (float x = -3; x <= 3; x += 1) {
+                err += test2(y, x, atan2(y,x), tol, [](skvm::F32 y, skvm::F32 x) {
+                    return approx_atan(y,x);
+                });
+            }
+        }
+        if (0) { SkDebugf("atan2 error %g\n", err); }
     }
 }