Correct reciprocal approximation for power-of-two values.

src/Reactor/Nucleus.cpp

@@ -33,6 +33,7 @@
 #include "Thread.hpp"
 #include "Memory.hpp"
 
+#include <xmmintrin.h>
 #include <fstream>
 
 #if defined(__x86_64__) && defined(_WIN32)
@@ -6069,10 +6070,19 @@ namespace sw
 		return IfThenElse(x < y, x, y);
 	}
 
-	RValue<Float> Rcp_pp(RValue<Float> x)
+	RValue<Float> Rcp_pp(RValue<Float> x, bool exactAtPow2)
+	{
+		if(exactAtPow2)
+		{
+			// rcpss uses a piecewise-linear approximation which minimizes the relative error
+			// but is not exact at power-of-two values. Rectify by multiplying by the inverse.
+			return x86::rcpss(x) * Float(1.0f / _mm_cvtss_f32(_mm_rcp_ss(_mm_set_ps1(1.0f))));
+		}
+		else
 		{
 			return x86::rcpss(x);
 		}
+	}
 
 	RValue<Float> RcpSqrt_pp(RValue<Float> x)
 	{
@@ -6580,10 +6590,19 @@ namespace sw
 		return x86::minps(x, y);
 	}
 
-	RValue<Float4> Rcp_pp(RValue<Float4> x)
+	RValue<Float4> Rcp_pp(RValue<Float4> x, bool exactAtPow2)
+	{
+		if(exactAtPow2)
+		{
+			// rcpps uses a piecewise-linear approximation which minimizes the relative error
+			// but is not exact at power-of-two values. Rectify by multiplying by the inverse.
+			return x86::rcpps(x) * Float4(1.0f / _mm_cvtss_f32(_mm_rcp_ss(_mm_set_ps1(1.0f))));
+		}
+		else
 		{
 			return x86::rcpps(x);
 		}
+	}
 
 	RValue<Float4> RcpSqrt_pp(RValue<Float4> x)
 	{

src/Reactor/Nucleus.hpp

@@ -1899,7 +1899,7 @@ namespace sw
 	RValue<Float> Abs(RValue<Float> x);
 	RValue<Float> Max(RValue<Float> x, RValue<Float> y);
 	RValue<Float> Min(RValue<Float> x, RValue<Float> y);
-	RValue<Float> Rcp_pp(RValue<Float> val);
+	RValue<Float> Rcp_pp(RValue<Float> val, bool exactAtPow2 = false);
 	RValue<Float> RcpSqrt_pp(RValue<Float> val);
 	RValue<Float> Sqrt(RValue<Float> x);
 	RValue<Float> Round(RValue<Float> val);
@@ -2377,7 +2377,7 @@ namespace sw
 	RValue<Float4> Abs(RValue<Float4> x);
 	RValue<Float4> Max(RValue<Float4> x, RValue<Float4> y);
 	RValue<Float4> Min(RValue<Float4> x, RValue<Float4> y);
-	RValue<Float4> Rcp_pp(RValue<Float4> val);
+	RValue<Float4> Rcp_pp(RValue<Float4> val, bool exactAtPow2 = false);
 	RValue<Float4> RcpSqrt_pp(RValue<Float4> val);
 	RValue<Float4> Sqrt(RValue<Float4> x);
 	RValue<Float4> Insert(const Float4 &val, RValue<Float> element, int i);

src/Shader/PixelRoutine.cpp

@@ -140,7 +140,7 @@ namespace sw
 			if(interpolateW())
 			{
 				w = interpolate(xxxx, Dw, rhw, primitive + OFFSET(Primitive,w), false, false);
-				rhw = reciprocal(w);
+				rhw = reciprocal(w, false, false, true);
 
 				if(state.centroid)
 				{

src/Shader/ShaderCore.cpp

@@ -271,7 +271,7 @@ namespace sw
 		return exponential2(log, pp);
 	}
 
-	Float4 reciprocal(RValue<Float4> x, bool pp, bool finite)
+	Float4 reciprocal(RValue<Float4> x, bool pp, bool finite, bool exactAtPow2)
 	{
 		Float4 rcp;
 
@@ -281,7 +281,7 @@ namespace sw
 		}
 		else
 		{
-			rcp = Rcp_pp(x);
+			rcp = Rcp_pp(x, exactAtPow2);
 
 			if(!pp)
 			{

src/Shader/ShaderCore.hpp

@@ -87,7 +87,7 @@ namespace sw
 	Float4 exponential(RValue<Float4> x, bool pp = false);
 	Float4 logarithm(RValue<Float4> x, bool abs, bool pp = false);
 	Float4 power(RValue<Float4> x, RValue<Float4> y, bool pp = false);
-	Float4 reciprocal(RValue<Float4> x, bool pp = false, bool finite = false);
+	Float4 reciprocal(RValue<Float4> x, bool pp = false, bool finite = false, bool exactAtPow2 = false);
 	Float4 reciprocalSquareRoot(RValue<Float4> x, bool abs, bool pp = false);
 	Float4 modulo(RValue<Float4> x, RValue<Float4> y);
 	Float4 sine_pi(RValue<Float4> x, bool pp = false);     // limited to [-pi, pi] range