From 4c49e74d9313f1c7360d038b3e19fe3c3e113339 Mon Sep 17 00:00:00 2001
From: searchivairus <leo@boytsov.info>
Date: Wed, 10 Jan 2018 23:19:23 -0500
Subject: [PATCH] No ticket: adding a wrapper for efficient exponentiation
 functions.

---
 similarity_search/include/pow.h    | 60 ++++++++++++++++++++++++--
 similarity_search/include/utils.h  |  8 ++++
 similarity_search/test/test_pow.cc | 68 ++++++++++++++++++++++++++++++
 3 files changed, 133 insertions(+), 3 deletions(-)
 create mode 100644 similarity_search/test/test_pow.cc

diff --git a/similarity_search/include/pow.h b/similarity_search/include/pow.h
index 3dc496b..4d5cfc9 100644
--- a/similarity_search/include/pow.h
+++ b/similarity_search/include/pow.h
@@ -26,6 +26,9 @@
  */
 
 #include <cmath>
+#include <limits>
+
+#include "logging.h"
 
 namespace similarity {
 
@@ -150,9 +153,6 @@ inline T EfficientPow(T Base, unsigned Exp) {
 }
 
 
-
-}
-
 template <class T>
 inline T EfficientFractPowUtil(T Base, uint64_t Exp, uint64_t MaxK) {
     if (Exp == 0)    return 1;     // pow == 0 
@@ -185,5 +185,59 @@ inline T EfficientFractPow(T Base, T FractExp, unsigned NumDig) {
     return EfficientFractPowUtil(Base, Exp, MaxK);
 }
 
+/**
+  * A helper object that does some preprocessing for subsequent
+  * efficient computation of both integer and fractional
+  * powers for exponents x, where x * 2^maxDig is an integer.
+  * In other words, this can be done for exponents that have
+  * zeros beyond maxDig binary digits after the binary point.
+  * When the exponent does not satisfy this property,
+  * we simply default to using the standard std::pow function.
+  */
+template <class T>
+class PowerProxyObject {
+public:
+  /**
+    * Constructor.
+    *
+    * @param  p       is an exponent
+    * @param  maxDig  a maximum number of binary digits to consider (should be <= 31).
+    */
+  PowerProxyObject(const T p, const unsigned maxDig = 18) {
+    CHECK_MSG(p >= 0, "The exponent should be non-negative");
+
+    maxK_ = 1u << maxDig;
+    unsigned pfm = static_cast<unsigned>(std::floor(maxK_ * p));
+
+    p_         = p;
+    isOptim_   = (fabs(maxK_*p_ - pfm) <= 2 * std::numeric_limits<T>::min());
+    intPow_    = pfm >> maxDig;
+    fractPow_  = pfm - (intPow_ << maxDig);
+
+  }
+  
+  /**
+    * Compute pow(base, p_) possibly efficiently. We expect base to be 
+    * non-negative!
+    */
+  T inline pow(const T base) {
+    if (isOptim_) {
+      T mult1  = intPow_ ?    EfficientPow(base, intPow_) : 1;
+      T mult2  = fractPow_ ?  EfficientFractPowUtil(base, fractPow_, maxK_) : 1;
+      return mult1 * mult2;
+    } else {
+      return std::pow(base, p_);
+    }
+  }
+private:
+  T        p_;
+  unsigned maxK_;
+  bool     isOptim_;
+  unsigned intPow_;
+  unsigned fractPow_;
+};
+
+
+}
 
 #endif
diff --git a/similarity_search/include/utils.h b/similarity_search/include/utils.h
index c1aa21c..977950b 100644
--- a/similarity_search/include/utils.h
+++ b/similarity_search/include/utils.h
@@ -319,6 +319,14 @@ inline void ReplaceSomePunct(string &s) {
     if (s[i] == ',' || s[i] == FIELD_DELIMITER) s[i] = ' ';
 }
 
+template <class T>
+T getRelDiff(T val1, T val2) {
+  T diff = std::fabs(val1 - val2);
+  T maxVal = std::max(std::fabs(val1), std::fabs(val2));
+  T diffRel = diff/ std::max(maxVal, numeric_limits<T>::min()) ;
+  return diffRel;
+}
+
 template <class T>
 struct AutoVectDel {
   AutoVectDel(typename std::vector<T *>& Vector) : mVector(Vector) {}
diff --git a/similarity_search/test/test_pow.cc b/similarity_search/test/test_pow.cc
new file mode 100644
index 0000000..2e9131f
--- /dev/null
+++ b/similarity_search/test/test_pow.cc
@@ -0,0 +1,68 @@
+/**
+ * Non-metric Space Library
+ *
+ * Authors: Bilegsaikhan Naidan (https://github.com/bileg), Leonid Boytsov (http://boytsov.info).
+ * With contributions from Lawrence Cayton (http://lcayton.com/) and others.
+ *
+ * For the complete list of contributors and further details see:
+ * https://github.com/searchivarius/NonMetricSpaceLib 
+ * 
+ * Copyright (c) 2018
+ *
+ * This code is released under the * Apache License Version 2.0 http://www.apache.org/licenses/.
+ *
+ */
+
+#include <memory>
+#include <string>
+#include <cmath>
+#include <vector>
+#include <limits>
+
+#include "bunit.h"
+#include "pow.h"
+#include "logging.h"
+#include "utils.h"
+
+namespace similarity {
+
+using namespace std;
+
+const float MAX_REL_DIFF = 1e-6f;
+
+vector<float>   addExps = { 0, 0.125, 0.25, 0.5 };
+vector<float>   vals    = { 0.1, 0.5, 1, 1.5, 2, 4};
+
+template <typename T> bool runTest() {
+  for (unsigned i = 0; i <= 128; ++i) {
+    for (float a : addExps) {
+      T oneExp = a + i; 
+
+      PowerProxyObject<T> obj(oneExp);
+
+      for (float v : vals) {
+        T expRes = pow(T(v), oneExp); 
+        T res = obj.pow(v);
+
+        T absDiff = fabs(res - expRes);
+        T relDiff = getRelDiff(res, expRes);
+        if ( relDiff > MAX_REL_DIFF ) {
+          LOG(LIB_ERROR) << "Mismatch for base=" << v << " exponent=" << oneExp << " expected: " << expRes << " obtained: " << res << " abs diff: " << absDiff << " rel diff: " << relDiff;
+          return false;
+        }
+      }
+    }
+  }
+
+  return true;
+}
+
+TEST(pow_float) {
+  EXPECT_TRUE(runTest<float>());
+}
+
+TEST(pow_double) {
+  EXPECT_TRUE(runTest<double>());
+}
+
+}  // namespace similarity