Adds a generalized shift operation to the Fixed16 class. (#508)

* Adds a generalized shift operation to the Fixed16 class.

* Simplifies the code a bit.

src/utils/Fixed16.cxxtest

@@ -213,3 +213,39 @@ TEST(Fixed16Test, SignedZero)
     // Adding and subtracting does not preserve the sign.
     EXPECT_TRUE(v1.is_positive());
 }
+
+/// Helper function to test mulpow2 operation. Computes base mulpow2 shift with
+/// fixed16 and with double, and verifies that the results is within the given
+/// relative precision form each other.
+/// @param base left operand of mulpow2
+/// @param shift right operand of mulpow2
+/// @param precision relative precisoin of how close the FIxed16 computation
+/// should be to the real result. Example 0.01 for 1% precision.
+void mulpow2_test(double base, double shift, double precision)
+{
+    double result = base * pow(2.0, shift);
+    string expl = StringPrintf("while computing %lg mulpow2 %lg", base, shift);
+    SCOPED_TRACE(expl);
+    Fixed16 fbase {Fixed16::FROM_DOUBLE, base};
+    Fixed16 fshift {Fixed16::FROM_DOUBLE, shift};
+    fbase.mulpow2(fshift);
+    EXPECT_NEAR(result, fbase.to_float(), result * precision);
+}
+
+TEST(Fixed16Test, MulPow2)
+{
+    // General case is within 1% precision.
+    mulpow2_test(13, 1.385, 0.01);
+    mulpow2_test(13, -2.442, 0.01);
+
+    // Integer shifts should work very well.
+    mulpow2_test(1, 4, 0.00001);
+    mulpow2_test(1.77429, 4, 0.00001);
+    mulpow2_test(1, -4, 0.00001);
+    mulpow2_test(1.77429, -4, 0.0001);
+
+    // When the fractional part does not have a lot of binary digits, the
+    // result should be pretty good too.
+    mulpow2_test(30481, -12.5, 0.00001);
+    mulpow2_test(4377, 2.375, 0.0001);
+}

src/utils/Fixed16.hxx

@@ -132,6 +132,56 @@ public:
         return round();
     }
 
+    /// Multiplies *this with pow(2, o). This is effectively a generalized
+    /// shift operation that works on fractional numbers too. The precision is
+    /// limited.
+    ///
+    /// Modifies *this.
+    /// @param o number of "bits" to shift with. May be positive or negative.
+    /// @return *this = *this * pow(2, o);
+    Fixed16 &mulpow2(Fixed16 o)
+    {
+        const Fixed16* coeffs;
+        uint16_t f;
+        if (o.is_positive())
+        {
+            // multiplying
+            value_ <<= o.trunc();
+            f = o.frac();
+            static constexpr const Fixed16 pown[6] = {
+                {FROM_DOUBLE, 1.4142135623730951},  // 2^(1/2)
+                {FROM_DOUBLE, 1.1892071150027210},  // 2^(1/4)
+                {FROM_DOUBLE, 1.0905077326652577},  // 2^(1/8)
+                {FROM_DOUBLE, 1.0442737824274138},  // 2^(1/16)
+                {FROM_DOUBLE, 1.0218971486541166},  // 2^(1/32)
+                {FROM_DOUBLE, 1.0108892860517005}}; // 2^(1/64)
+            coeffs = pown;
+        }
+        else
+        {
+            // dividing
+            o.negate();
+            value_ >>= o.trunc();
+            f = o.frac();
+            static constexpr const Fixed16 pown[6] = {
+                {FROM_DOUBLE, 0.7071067811865476},  // 2^(-1/2)
+                {FROM_DOUBLE, 0.8408964152537145},  // 2^(-1/4)
+                {FROM_DOUBLE, 0.9170040432046712},  // 2^(-1/8)
+                {FROM_DOUBLE, 0.9576032806985737},  // 2^(-1/16)
+                {FROM_DOUBLE, 0.9785720620877001},  // 2^(-1/32)
+                {FROM_DOUBLE, 0.9892280131939755}}; // 2^(-1/64)
+            coeffs = pown;
+        }
+        for (unsigned idx = 0, bit = 0x8000; idx < 6; ++idx, bit >>= 1)
+        {
+            if (f & bit)
+            {
+                *this *= coeffs[idx];
+            }
+        }
+        return *this;
+    }
+
     /// @return the value rounded to nearest integer.
     int16_t round() const {
         int16_t b = (value_ + 0x8000) >> 16;