From 1b84ef2b50958bd22a2706f0cfbdf2d623bb59ac Mon Sep 17 00:00:00 2001 From: Mike Reed Date: Mon, 13 Apr 2020 17:56:24 -0400 Subject: [PATCH] [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 Commit-Queue: Mike Reed --- src/core/SkVM.cpp | 48 ++++++++++++++++++++++++++++++++++++++-------- src/core/SkVM.h | 2 ++ tests/SkVMTest.cpp | 31 ++++++++++++++++++++++++++++-- 3 files changed, 71 insertions(+), 10 deletions(-) diff --git a/src/core/SkVM.cpp b/src/core/SkVM.cpp index 7467e98cc3..cadc4a53be 100644 --- a/src/core/SkVM.cpp +++ b/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)}; diff --git a/src/core/SkVM.h b/src/core/SkVM.h index 728d543b8a..5b730fb427 100644 --- a/src/core/SkVM.h +++ b/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); } diff --git a/tests/SkVMTest.cpp b/tests/SkVMTest.cpp index 41d8cebdd1..c89542c9d7 100644 --- a/tests/SkVMTest.cpp +++ b/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(); + skvm::Arg in1 = b.varying(); + skvm::Arg out = b.varying(); + 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); } } }