From 0b44e0131998806baf01fd675454072d8820dfcf Mon Sep 17 00:00:00 2001 From: "Steinar H. Gunderson" Date: Sat, 21 Jul 2018 10:58:59 +0200 Subject: [PATCH] Pack the two b values in the equations with shared exponent. --- equations.frag | 40 +++++++++++++++++++++++++++++++++++++++- sor.frag | 19 ++++++++++++++++++- 2 files changed, 57 insertions(+), 2 deletions(-) diff --git a/equations.frag b/equations.frag index 190be60..3dc21ec 100644 --- a/equations.frag +++ b/equations.frag @@ -11,6 +11,44 @@ uniform sampler2D smoothness_x_tex, smoothness_y_tex; // TODO: Consider a specialized version for the case where we know that du = dv = 0, // since we run so few iterations. +// Similar to packHalf2x16, but the two values share exponent, and are stored +// as 12-bit fixed point numbers multiplied by that exponent (the leading one +// can't be implicit in this kind of format). This allows us to store a much +// greater range of numbers (8-bit, ie., full fp32 range), and also gives us an +// extra mantissa bit. (Well, ostensibly two, but because the numbers have to +// be stored denormalized, we only really gain one.) +// +// The price we pay is that if the numbers are of very different magnitudes, +// the smaller number gets less precision. +uint pack_floats_shared(float a, float b) +{ + float greatest = max(abs(a), abs(b)); + + // Find the exponent, increase it by one, and negate it. + // E.g., if the nonbiased exponent is 3, the number is between + // 2^3 and 2^4, so our normalization factor to get within -1..1 + // is going to be 2^-4. + // + // exponent -= 127; + // exponent = -(exponent + 1); + // exponent += 127; + // + // is the same as + // + // exponent = 252 - exponent; + uint e = floatBitsToUint(greatest) & 0x7f800000u; + float normalizer = uintBitsToFloat((252 << 23) - e); + + // The exponent is the same range as fp32, so just copy it + // verbatim, shifted up to where the sign bit used to be. + e <<= 1; + + // Quantize to 12 bits. + uint qa = uint(int(round(a * (normalizer * 2047.0)))); + uint qb = uint(int(round(b * (normalizer * 2047.0)))); + + return (qa & 0xfffu) | ((qb & 0xfffu) << 12) | e; +} void main() { @@ -106,5 +144,5 @@ void main() equation.x = floatBitsToUint(1.0 / A11); equation.y = floatBitsToUint(A12); equation.z = floatBitsToUint(1.0 / A22); - equation.w = packHalf2x16(vec2(b1, b2)); + equation.w = pack_floats_shared(b1, b2); } diff --git a/sor.frag b/sor.frag index dc4a5b8..25f0ce5 100644 --- a/sor.frag +++ b/sor.frag @@ -6,13 +6,30 @@ out vec2 diff_flow; uniform sampler2D diff_flow_tex, smoothness_x_tex, smoothness_y_tex; uniform usampler2D equation_tex; +// See pack_floats_shared() in equations.frag. +vec2 unpack_floats_shared(uint c) +{ + // Recover the exponent, and multiply it in. Add one because + // we have denormalized mantissas, then another one because we + // already reduced the exponent by one. Then subtract 20, because + // we are going to shift up the number by 20 below to recover the sign bits. + float normalizer = uintBitsToFloat(((c >> 1) & 0x7f800000u) - (18 << 23)); + normalizer *= (1.0 / 2047.0); + + // Shift the values up so that we recover the sign bit, then normalize. + float a = int(uint(c & 0x000fffu) << 20) * normalizer; + float b = int(uint(c & 0xfff000u) << 8) * normalizer; + + return vec2(a, b); +} + void main() { uvec4 equation = texture(equation_tex, tc); float inv_A11 = uintBitsToFloat(equation.x); float A12 = uintBitsToFloat(equation.y); float inv_A22 = uintBitsToFloat(equation.z); - vec2 b = unpackHalf2x16(equation.w); + vec2 b = unpack_floats_shared(equation.w); // Subtract the missing terms from the right-hand side // (it couldn't be done earlier, because we didn't know -- 2.39.2