Copy optimized transcendental intrinsics to ShaderCore
SwiftShader currently uses Reactor's implementations of transcendental
functions for SPIR-V GLSL.std.450 extended instructions. This puts
graphics-specific intrinsics in Reactor, and we have to specify the
desired precision through a parameter which only informally reflects
the Vulkan precision requirements.
These implementations belong in Vulkan-specific code, so as a first step
this change copies them over. The GLSLstd450 code explicitly uses the
rr namespace functions to not alter any behavior at this point.
Bug: b/169755552
Change-Id: Ia2c9e71a2688c077d358d37641917ab40f46d4be
Reviewed-on: https://swiftshader-review.googlesource.com/c/SwiftShader/+/62231
Kokoro-Result: kokoro <noreply+kokoro@google.com>
Tested-by: Nicolas Capens <nicolascapens@google.com>
Reviewed-by: Alexis Hétu <sugoi@google.com>
diff --git a/src/Pipeline/ShaderCore.cpp b/src/Pipeline/ShaderCore.cpp
index 81a4251..b796868 100644
--- a/src/Pipeline/ShaderCore.cpp
+++ b/src/Pipeline/ShaderCore.cpp
@@ -150,6 +150,282 @@
return x;
}
+static Float4 Reciprocal(RValue<Float4> x, bool pp = false, bool finite = false, bool exactAtPow2 = false)
+{
+ Float4 rcp = Rcp_pp(x, exactAtPow2);
+
+ if(!pp)
+ {
+ rcp = (rcp + rcp) - (x * rcp * rcp);
+ }
+
+ return rcp;
+}
+
+static Float4 SinOrCos(RValue<Float4> x, bool sin)
+{
+ // Reduce to [-0.5, 0.5] range
+ Float4 y = x * Float4(1.59154943e-1f); // 1/2pi
+ y = y - Round(y);
+
+ // From the paper: "A Fast, Vectorizable Algorithm for Producing Single-Precision Sine-Cosine Pairs"
+ // This implementation passes OpenGL ES 3.0 precision requirements, at the cost of more operations:
+ // !pp : 17 mul, 7 add, 1 sub, 1 reciprocal
+ // pp : 4 mul, 2 add, 2 abs
+
+ Float4 y2 = y * y;
+ Float4 c1 = y2 * (y2 * (y2 * Float4(-0.0204391631f) + Float4(0.2536086171f)) + Float4(-1.2336977925f)) + Float4(1.0f);
+ Float4 s1 = y * (y2 * (y2 * (y2 * Float4(-0.0046075748f) + Float4(0.0796819754f)) + Float4(-0.645963615f)) + Float4(1.5707963235f));
+ Float4 c2 = (c1 * c1) - (s1 * s1);
+ Float4 s2 = Float4(2.0f) * s1 * c1;
+ Float4 r = Reciprocal(s2 * s2 + c2 * c2);
+
+ if(sin)
+ {
+ return Float4(2.0f) * s2 * c2 * r;
+ }
+ else
+ {
+ return ((c2 * c2) - (s2 * s2)) * r;
+ }
+}
+
+// Approximation of atan in [0..1]
+static Float4 Atan_01(Float4 x)
+{
+ // From 4.4.49, page 81 of the Handbook of Mathematical Functions, by Milton Abramowitz and Irene Stegun
+ const Float4 a2(-0.3333314528f);
+ const Float4 a4(0.1999355085f);
+ const Float4 a6(-0.1420889944f);
+ const Float4 a8(0.1065626393f);
+ const Float4 a10(-0.0752896400f);
+ const Float4 a12(0.0429096138f);
+ const Float4 a14(-0.0161657367f);
+ const Float4 a16(0.0028662257f);
+ Float4 x2 = x * x;
+ return (x + x * (x2 * (a2 + x2 * (a4 + x2 * (a6 + x2 * (a8 + x2 * (a10 + x2 * (a12 + x2 * (a14 + x2 * a16)))))))));
+}
+
+Float4 Sin(RValue<Float4> x)
+{
+ return SinOrCos(x, true);
+}
+
+Float4 Cos(RValue<Float4> x)
+{
+ return SinOrCos(x, false);
+}
+
+Float4 Tan(RValue<Float4> x)
+{
+ return SinOrCos(x, true) / SinOrCos(x, false);
+}
+
+static Float4 Asin_4_terms(RValue<Float4> x)
+{
+ // From 4.4.45, page 81 of the Handbook of Mathematical Functions, by Milton Abramowitz and Irene Stegun
+ // |e(x)| <= 5e-8
+ const Float4 half_pi(1.57079632f);
+ const Float4 a0(1.5707288f);
+ const Float4 a1(-0.2121144f);
+ const Float4 a2(0.0742610f);
+ const Float4 a3(-0.0187293f);
+ Float4 absx = Abs(x);
+ return As<Float4>(As<Int4>(half_pi - Sqrt(Float4(1.0f) - absx) * (a0 + absx * (a1 + absx * (a2 + absx * a3)))) ^
+ (As<Int4>(x) & Int4(0x80000000)));
+}
+
+static Float4 Asin_8_terms(RValue<Float4> x)
+{
+ // From 4.4.46, page 81 of the Handbook of Mathematical Functions, by Milton Abramowitz and Irene Stegun
+ // |e(x)| <= 0e-8
+ const Float4 half_pi(1.5707963268f);
+ const Float4 a0(1.5707963050f);
+ const Float4 a1(-0.2145988016f);
+ const Float4 a2(0.0889789874f);
+ const Float4 a3(-0.0501743046f);
+ const Float4 a4(0.0308918810f);
+ const Float4 a5(-0.0170881256f);
+ const Float4 a6(0.006700901f);
+ const Float4 a7(-0.0012624911f);
+ Float4 absx = Abs(x);
+ return As<Float4>(As<Int4>(half_pi - Sqrt(Float4(1.0f) - absx) * (a0 + absx * (a1 + absx * (a2 + absx * (a3 + absx * (a4 + absx * (a5 + absx * (a6 + absx * a7)))))))) ^
+ (As<Int4>(x) & Int4(0x80000000)));
+}
+
+RValue<Float4> Asin(RValue<Float4> x, Precision p)
+{
+ // TODO(b/169755566): Surprisingly, deqp-vk's precision.acos.highp/mediump tests pass when using the 4-term polynomial
+ // approximation version of acos, unlike for Asin, which requires higher precision algorithms.
+
+ if(p == Precision::Full)
+ {
+ return rr::Asin(x, p);
+ }
+
+ return Asin_8_terms(x);
+}
+
+RValue<Float4> Acos(RValue<Float4> x, Precision p)
+{
+ // pi/2 - arcsin(x)
+ return Float4(1.57079632e+0f) - Asin_4_terms(x);
+}
+
+Float4 Atan(RValue<Float4> x)
+{
+ Float4 absx = Abs(x);
+ Int4 O = CmpNLT(absx, Float4(1.0f));
+ Float4 y = As<Float4>((O & As<Int4>(Float4(1.0f) / absx)) | (~O & As<Int4>(absx))); // FIXME: Vector select
+
+ const Float4 half_pi(1.57079632f);
+ Float4 theta = Atan_01(y);
+ return As<Float4>(((O & As<Int4>(half_pi - theta)) | (~O & As<Int4>(theta))) ^ // FIXME: Vector select
+ (As<Int4>(x) & Int4(0x80000000)));
+}
+
+Float4 Atan2(RValue<Float4> y, RValue<Float4> x)
+{
+ const Float4 pi(3.14159265f); // pi
+ const Float4 minus_pi(-3.14159265f); // -pi
+ const Float4 half_pi(1.57079632f); // pi/2
+ const Float4 quarter_pi(7.85398163e-1f); // pi/4
+
+ // Rotate to upper semicircle when in lower semicircle
+ Int4 S = CmpLT(y, Float4(0.0f));
+ Float4 theta = As<Float4>(S & As<Int4>(minus_pi));
+ Float4 x0 = As<Float4>((As<Int4>(y) & Int4(0x80000000)) ^ As<Int4>(x));
+ Float4 y0 = Abs(y);
+
+ // Rotate to right quadrant when in left quadrant
+ Int4 Q = CmpLT(x0, Float4(0.0f));
+ theta += As<Float4>(Q & As<Int4>(half_pi));
+ Float4 x1 = As<Float4>((Q & As<Int4>(y0)) | (~Q & As<Int4>(x0))); // FIXME: Vector select
+ Float4 y1 = As<Float4>((Q & As<Int4>(-x0)) | (~Q & As<Int4>(y0))); // FIXME: Vector select
+
+ // Mirror to first octant when in second octant
+ Int4 O = CmpNLT(y1, x1);
+ Float4 x2 = As<Float4>((O & As<Int4>(y1)) | (~O & As<Int4>(x1))); // FIXME: Vector select
+ Float4 y2 = As<Float4>((O & As<Int4>(x1)) | (~O & As<Int4>(y1))); // FIXME: Vector select
+
+ // Approximation of atan in [0..1]
+ Int4 zero_x = CmpEQ(x2, Float4(0.0f));
+ Int4 inf_y = IsInf(y2); // Since x2 >= y2, this means x2 == y2 == inf, so we use 45 degrees or pi/4
+ Float4 atan2_theta = Atan_01(y2 / x2);
+ theta += As<Float4>((~zero_x & ~inf_y & ((O & As<Int4>(half_pi - atan2_theta)) | (~O & (As<Int4>(atan2_theta))))) | // FIXME: Vector select
+ (inf_y & As<Int4>(quarter_pi)));
+
+ // Recover loss of precision for tiny theta angles
+ // This combination results in (-pi + half_pi + half_pi - atan2_theta) which is equivalent to -atan2_theta
+ Int4 precision_loss = S & Q & O & ~inf_y;
+
+ return As<Float4>((precision_loss & As<Int4>(-atan2_theta)) | (~precision_loss & As<Int4>(theta))); // FIXME: Vector select
+}
+
+Float4 Exp2(RValue<Float4> x)
+{
+ // This implementation is based on 2^(i + f) = 2^i * 2^f,
+ // where i is the integer part of x and f is the fraction.
+
+ // For 2^i we can put the integer part directly in the exponent of
+ // the IEEE-754 floating-point number. Clamp to prevent overflow
+ // past the representation of infinity.
+ Float4 x0 = x;
+ x0 = Min(x0, As<Float4>(Int4(0x43010000))); // 129.00000e+0f
+ x0 = Max(x0, As<Float4>(Int4(0xC2FDFFFF))); // -126.99999e+0f
+
+ Int4 i = RoundInt(x0 - Float4(0.5f));
+ Float4 ii = As<Float4>((i + Int4(127)) << 23); // Add single-precision bias, and shift into exponent.
+
+ // For the fractional part use a polynomial
+ // which approximates 2^f in the 0 to 1 range.
+ Float4 f = x0 - Float4(i);
+ Float4 ff = As<Float4>(Int4(0x3AF61905)); // 1.8775767e-3f
+ ff = ff * f + As<Float4>(Int4(0x3C134806)); // 8.9893397e-3f
+ ff = ff * f + As<Float4>(Int4(0x3D64AA23)); // 5.5826318e-2f
+ ff = ff * f + As<Float4>(Int4(0x3E75EAD4)); // 2.4015361e-1f
+ ff = ff * f + As<Float4>(Int4(0x3F31727B)); // 6.9315308e-1f
+ ff = ff * f + Float4(1.0f);
+
+ return ii * ff;
+}
+
+Float4 Log2(RValue<Float4> x)
+{
+ Float4 x0;
+ Float4 x1;
+ Float4 x2;
+ Float4 x3;
+
+ x0 = x;
+
+ x1 = As<Float4>(As<Int4>(x0) & Int4(0x7F800000));
+ x1 = As<Float4>(As<UInt4>(x1) >> 8);
+ x1 = As<Float4>(As<Int4>(x1) | As<Int4>(Float4(1.0f)));
+ x1 = (x1 - Float4(1.4960938f)) * Float4(256.0f); // FIXME: (x1 - 1.4960938f) * 256.0f;
+ x0 = As<Float4>((As<Int4>(x0) & Int4(0x007FFFFF)) | As<Int4>(Float4(1.0f)));
+
+ x2 = (Float4(9.5428179e-2f) * x0 + Float4(4.7779095e-1f)) * x0 + Float4(1.9782813e-1f);
+ x3 = ((Float4(1.6618466e-2f) * x0 + Float4(2.0350508e-1f)) * x0 + Float4(2.7382900e-1f)) * x0 + Float4(4.0496687e-2f);
+ x2 /= x3;
+
+ x1 += (x0 - Float4(1.0f)) * x2;
+
+ Int4 pos_inf_x = CmpEQ(As<Int4>(x), Int4(0x7F800000));
+ return As<Float4>((pos_inf_x & As<Int4>(x)) | (~pos_inf_x & As<Int4>(x1)));
+}
+
+Float4 Exp(RValue<Float4> x)
+{
+ // TODO: Propagate the constant
+ return sw::Exp2(Float4(1.44269504f) * x); // 1/ln(2)
+}
+
+Float4 Log(RValue<Float4> x)
+{
+ // TODO: Propagate the constant
+ return Float4(6.93147181e-1f) * sw::Log2(x); // ln(2)
+}
+
+Float4 Pow(RValue<Float4> x, RValue<Float4> y)
+{
+ Float4 log = sw::Log2(x);
+ log *= y;
+ return sw::Exp2(log);
+}
+
+Float4 Sinh(RValue<Float4> x)
+{
+ return (sw::Exp(x) - sw::Exp(-x)) * Float4(0.5f);
+}
+
+Float4 Cosh(RValue<Float4> x)
+{
+ return (sw::Exp(x) + sw::Exp(-x)) * Float4(0.5f);
+}
+
+Float4 Tanh(RValue<Float4> x)
+{
+ Float4 e_x = sw::Exp(x);
+ Float4 e_minus_x = sw::Exp(-x);
+ return (e_x - e_minus_x) / (e_x + e_minus_x);
+}
+
+Float4 Asinh(RValue<Float4> x)
+{
+ return sw::Log(x + Sqrt(x * x + Float4(1.0f)));
+}
+
+Float4 Acosh(RValue<Float4> x)
+{
+ return sw::Log(x + Sqrt(x + Float4(1.0f)) * Sqrt(x - Float4(1.0f)));
+}
+
+Float4 Atanh(RValue<Float4> x)
+{
+ return sw::Log((Float4(1.0f) + x) / (Float4(1.0f) - x)) * Float4(0.5f);
+}
+
Float4 exponential2(RValue<Float4> x, bool pp)
{
// This implementation is based on 2^(i + f) = 2^i * 2^f,
diff --git a/src/Pipeline/ShaderCore.hpp b/src/Pipeline/ShaderCore.hpp
index 8817d06..ad7a4bc 100644
--- a/src/Pipeline/ShaderCore.hpp
+++ b/src/Pipeline/ShaderCore.hpp
@@ -183,6 +183,28 @@
} // namespace SIMD
+// Vulkan 'SPIR-V Extended Instructions for GLSL' (GLSL.std.450) compliant transcendental functions
+Float4 Sin(RValue<Float4> x);
+Float4 Cos(RValue<Float4> x);
+Float4 Tan(RValue<Float4> x);
+RValue<Float4> Asin(RValue<Float4> x, rr::Precision p); // TODO(b/169755552): Remove rr::Precision
+RValue<Float4> Acos(RValue<Float4> x, rr::Precision p); // TODO(b/169755552): Remove rr::Precision
+Float4 Atan(RValue<Float4> x);
+Float4 Atan2(RValue<Float4> y, RValue<Float4> x);
+Float4 Exp2(RValue<Float4> x);
+Float4 Log2(RValue<Float4> x);
+Float4 Exp(RValue<Float4> x);
+Float4 Log(RValue<Float4> x);
+Float4 Pow(RValue<Float4> x, RValue<Float4> y);
+Float4 Sinh(RValue<Float4> x);
+Float4 Cosh(RValue<Float4> x);
+Float4 Tanh(RValue<Float4> x);
+Float4 Asinh(RValue<Float4> x);
+Float4 Acosh(RValue<Float4> x);
+Float4 Atanh(RValue<Float4> x);
+
+// Legacy transcendental functions
+// TODO(b/169755552): Consolidate with the functions above
Float4 exponential2(RValue<Float4> x, bool pp = false);
Float4 logarithm2(RValue<Float4> x, bool pp = false);
Float4 exponential(RValue<Float4> x, bool pp = false);
diff --git a/src/Pipeline/SpirvShaderGLSLstd450.cpp b/src/Pipeline/SpirvShaderGLSLstd450.cpp
index 7672201..355dc19 100644
--- a/src/Pipeline/SpirvShaderGLSLstd450.cpp
+++ b/src/Pipeline/SpirvShaderGLSLstd450.cpp
@@ -595,7 +595,7 @@
auto radians = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Sin(radians.Float(i)));
+ dst.move(i, rr::Sin(radians.Float(i)));
}
}
break;
@@ -604,7 +604,7 @@
auto radians = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Cos(radians.Float(i)));
+ dst.move(i, rr::Cos(radians.Float(i)));
}
}
break;
@@ -613,7 +613,7 @@
auto radians = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Tan(radians.Float(i)));
+ dst.move(i, rr::Tan(radians.Float(i)));
}
}
break;
@@ -624,7 +624,7 @@
ApplyDecorationsForId(&d, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Asin(val.Float(i), d.RelaxedPrecision ? Precision::Relaxed : Precision::Full));
+ dst.move(i, rr::Asin(val.Float(i), d.RelaxedPrecision ? Precision::Relaxed : Precision::Full));
}
}
break;
@@ -635,7 +635,7 @@
ApplyDecorationsForId(&d, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Acos(val.Float(i), d.RelaxedPrecision ? Precision::Relaxed : Precision::Full));
+ dst.move(i, rr::Acos(val.Float(i), d.RelaxedPrecision ? Precision::Relaxed : Precision::Full));
}
}
break;
@@ -644,7 +644,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Atan(val.Float(i)));
+ dst.move(i, rr::Atan(val.Float(i)));
}
}
break;
@@ -653,7 +653,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Sinh(val.Float(i)));
+ dst.move(i, rr::Sinh(val.Float(i)));
}
}
break;
@@ -662,7 +662,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Cosh(val.Float(i)));
+ dst.move(i, rr::Cosh(val.Float(i)));
}
}
break;
@@ -671,7 +671,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Tanh(val.Float(i)));
+ dst.move(i, rr::Tanh(val.Float(i)));
}
}
break;
@@ -680,7 +680,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Asinh(val.Float(i)));
+ dst.move(i, rr::Asinh(val.Float(i)));
}
}
break;
@@ -689,7 +689,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Acosh(val.Float(i)));
+ dst.move(i, rr::Acosh(val.Float(i)));
}
}
break;
@@ -698,7 +698,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Atanh(val.Float(i)));
+ dst.move(i, rr::Atanh(val.Float(i)));
}
}
break;
@@ -708,7 +708,7 @@
auto y = Operand(this, state, insn.word(6));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Atan2(x.Float(i), y.Float(i)));
+ dst.move(i, rr::Atan2(x.Float(i), y.Float(i)));
}
}
break;
@@ -718,7 +718,7 @@
auto y = Operand(this, state, insn.word(6));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Pow(x.Float(i), y.Float(i)));
+ dst.move(i, rr::Pow(x.Float(i), y.Float(i)));
}
}
break;
@@ -727,7 +727,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Exp(val.Float(i)));
+ dst.move(i, rr::Exp(val.Float(i)));
}
}
break;
@@ -736,7 +736,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Log(val.Float(i)));
+ dst.move(i, rr::Log(val.Float(i)));
}
}
break;
@@ -745,7 +745,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Exp2(val.Float(i)));
+ dst.move(i, rr::Exp2(val.Float(i)));
}
}
break;
@@ -754,7 +754,7 @@
auto val = Operand(this, state, insn.word(5));
for(auto i = 0u; i < type.componentCount; i++)
{
- dst.move(i, Log2(val.Float(i)));
+ dst.move(i, rr::Log2(val.Float(i)));
}
}
break;
diff --git a/src/Reactor/Reactor.hpp b/src/Reactor/Reactor.hpp
index d4ca19a..f99c545 100644
--- a/src/Reactor/Reactor.hpp
+++ b/src/Reactor/Reactor.hpp
@@ -2411,7 +2411,6 @@
RValue<Float4> Ceil(RValue<Float4> x);
// Trigonometric functions
-// TODO: Currently unimplemented for Subzero.
RValue<Float4> Sin(RValue<Float4> x);
RValue<Float4> Cos(RValue<Float4> x);
RValue<Float4> Tan(RValue<Float4> x);
@@ -2427,7 +2426,6 @@
RValue<Float4> Atan2(RValue<Float4> x, RValue<Float4> y);
// Exponential functions
-// TODO: Currently unimplemented for Subzero.
RValue<Float4> Pow(RValue<Float4> x, RValue<Float4> y);
RValue<Float4> Exp(RValue<Float4> x);
RValue<Float4> Log(RValue<Float4> x);
