Skip testing subnormals when flushed to zero

When flush-to-zero (FTZ) mode is enabled, std::nextafter() returns zero
when a subnormal gets passed in. This can cause infinite loops when
trying to iterate over a range.

The inc() function returns the next larger value which does not compare
equal to the input. The only exception is negative zero, which produces
positive zero, since the CPU still discerns between these two when FTZ
has been enabled.

Bug: b/169904252
Change-Id: Idff987a4c14bd4f778e940ec4fc4de26e0a74372
Reviewed-on: https://swiftshader-review.googlesource.com/c/SwiftShader/+/62869
Kokoro-Result: kokoro <noreply+kokoro@google.com>
Tested-by: Nicolas Capens <nicolascapens@google.com>
Reviewed-by: Alexis Hétu <sugoi@google.com>

src/Pipeline/ShaderCore.cpp

diff --git a/src/Pipeline/ShaderCore.cpp b/src/Pipeline/ShaderCore.cpp
index cbc5a6e..f4695cf 100644
--- a/src/Pipeline/ShaderCore.cpp
+++ b/src/Pipeline/ShaderCore.cpp
@@ -166,7 +166,7 @@
 	return (x + x * (x2 * (a2 + x2 * (a4 + x2 * (a6 + x2 * (a8 + x2 * (a10 + x2 * (a12 + x2 * (a14 + x2 * a16)))))))));
 }
 
-// Polynomal approximation of order 5 for sin(x * 2 * pi) in the range [-1/4, 1/4]
+// Polynomial approximation of order 5 for sin(x * 2 * pi) in the range [-1/4, 1/4]
 static Float4 Sin5(Float4 x)
 {
 	// A * x^5 + B * x^3 + C * x

src/System/Math.hpp

diff --git a/src/System/Math.hpp b/src/System/Math.hpp
index 8c6a069..bafa046 100644
--- a/src/System/Math.hpp
+++ b/src/System/Math.hpp
@@ -370,11 +370,36 @@
 }
 
 // Convert floating value v to fixed point with p digits after the decimal point
-constexpr int toFixedPoint(float v, int p)
+inline constexpr int toFixedPoint(float v, int p)
 {
 	return static_cast<int>(v * (1 << p));
 }
 
+// Returns the next floating-point number which is not treated identical to the input.
+// Note that std::nextafter() does not skip representations flushed to zero.
+[[nodiscard]] inline float inc(float x)
+{
+	int x1 = bit_cast<int>(x);
+
+	while(bit_cast<float>(x1) == x)
+	{
+		// Since IEEE 754 uses ones' complement and integers are two's complement,
+		// we need to explicitly hop from negative zero to positive zero.
+		if(x1 == (int)0x80000000)  // -0.0f
+		{
+			// Note that while the comparison -0.0f == +0.0f returns true, this
+			// function returns the next value which can be treated differently.
+			return +0.0f;
+		}
+
+		// Negative ones' complement value are made less negative by subtracting 1
+		// in two's complement representation.
+		x1 += (x1 >= 0) ? 1 : -1;
+	}
+
+	return bit_cast<float>(x1);
+}
+
 }  // namespace sw
 
 #endif  // sw_Math_hpp

tests/MathUnitTests/unittests.cpp

diff --git a/tests/MathUnitTests/unittests.cpp b/tests/MathUnitTests/unittests.cpp
index 0721c6f..5943f1c 100644
--- a/tests/MathUnitTests/unittests.cpp
+++ b/tests/MathUnitTests/unittests.cpp
@@ -13,6 +13,7 @@
 // limitations under the License.
 
 #include "System/Half.hpp"
+#include "System/Math.hpp"
 
 #include <gmock/gmock.h>
 #include <gtest/gtest.h>
@@ -21,7 +22,7 @@
 
 using namespace sw;
 
-// Polynomal approximation of order 5 for sin(x * 2 * pi) in the range [-1/4, 1/4]
+// Polynomial approximation of order 5 for sin(x * 2 * pi) in the range [-1/4, 1/4]
 static float sin5(float x)
 {
 	// A * x^5 + B * x^3 + C * x
@@ -40,7 +41,7 @@
 	const float tolerance = powf(2.0f, -12.0f);  // Vulkan requires absolute error <= 2^−11 inside the range [−pi, pi]
 	const float pi = 3.1415926535f;
 
-	for(float x = -pi; x <= pi; x = nextafterf(x, +INFINITY))
+	for(float x = -pi; x <= pi; x = inc(x))
 	{
 		// Range reduction and mirroring
 		float x_2 = 0.25f - x * (0.5f / pi);
@@ -57,7 +58,7 @@
 	const float tolerance = powf(2.0f, -12.0f);  // Vulkan requires absolute error <= 2^−11 inside the range [−pi, pi]
 	const float pi = 3.1415926535f;
 
-	for(float x = -pi; x <= pi; x = nextafterf(x, +INFINITY))
+	for(float x = -pi; x <= pi; x = inc(x))
 	{
 		// Phase shift, range reduction, and mirroring
 		float x_2 = x * (0.5f / pi);