Finish repairing the deferred KxQ rebuild, and pin the estimator constants

LICENSE

@@ -296,3 +296,17 @@ APPENDIX B: Additional licenses relevant to this product.
     Code Locations:
       * common/include/ceiling_power_of_2.hpp
     that is adapted from the above.
+
+    =============================================================
+    FDLIBM
+    =============================================================
+    Original source code:
+      Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+      Developed at SunSoft, a Sun Microsystems, Inc. business.
+
+      Permission to use, copy, modify, and distribute this software is freely
+      granted, provided that this notice is preserved.
+
+    Code Locations:
+      common/include/fdlibm_log.hpp
+    that is adapted from __ieee754_log in FDLIBM 5.3.

common/CMakeLists.txt

@@ -40,6 +40,7 @@ install(FILES
       include/conditional_back_inserter.hpp
       include/conditional_forward.hpp
       include/count_zeros.hpp
+      include/fdlibm_log.hpp
       include/inv_pow2_table.hpp
       include/kolmogorov_smirnov_impl.hpp
       include/kolmogorov_smirnov.hpp

common/include/fdlibm_log.hpp

@@ -0,0 +1,101 @@
+// fdlibm __ieee754_log, used by Java's StrictMath.log (and by Math.log on most JVMs).
+// Derived from FDLIBM 5.3, Copyright (C) 1993 by Sun Microsystems, Inc.
+// "Permission to use, copy, modify, and distribute this software is freely granted,
+//  provided that this notice is preserved."
+#ifndef FDLIBM_LOG_HPP
+#define FDLIBM_LOG_HPP
+#include <cstdint>
+#include <cstring>
+#include <limits>
+
+namespace fdlibm {
+
+inline int32_t hi_word(double x){ uint64_t u; std::memcpy(&u,&x,8); return (int32_t)(uint32_t)(u>>32); }
+inline uint32_t lo_word(double x){ uint64_t u; std::memcpy(&u,&x,8); return (uint32_t)u; }
+inline void set_hi_word(double& x, uint32_t hi){ uint64_t u; std::memcpy(&u,&x,8);
+  u = (u & 0x00000000ffffffffULL) | ((uint64_t)hi<<32); std::memcpy(&x,&u,8); }
+
+// Forces a value to be rounded to a double before it is used again. fdlibm needs strict
+// IEEE-754 evaluation: a fused multiply-add anywhere in the polynomial below changes the
+// result. Pragmas are overridden by an explicit -ffp-contract=fast, so use a volatile
+// round-trip, which the standard requires the compiler to honour.
+inline double rnd(double v) { volatile double t = v; return t; }
+
+inline double log(double x) {
+// fdlibm depends on strict IEEE-754 evaluation: a fused multiply-add would change the result
+// of the polynomial evaluation below, so contraction must be off for this function.
+#if defined(__clang__)
+#pragma clang fp contract(off)
+#endif
+  static const double
+    ln2_hi = 6.93147180369123816490e-01,  /* 3fe62e42 fee00000 */
+    ln2_lo = 1.90821492927058770002e-10,  /* 3dea39ef 35793c76 */
+    two54  = 1.80143985094819840000e+16,  /* 43500000 00000000 */
+    Lg1 = 6.666666666666735130e-01,       /* 3FE55555 55555593 */
+    Lg2 = 3.999999999940941908e-01,       /* 3FD99999 9997FA04 */
+    Lg3 = 2.857142874366239149e-01,       /* 3FD24924 94229359 */
+    Lg4 = 2.222219843214978396e-01,       /* 3FCC71C5 1D8E78AF */
+    Lg5 = 1.818357216161805012e-01,       /* 3FC74664 96CB03DE */
+    Lg6 = 1.531383769920937332e-01,       /* 3FC39A09 D078C69F */
+    Lg7 = 1.479819860511658591e-01,       /* 3FC2F112 DF3E5244 */
+    zero = 0.0;
+
+  double hfsq,f,s,z,R,w,t1,t2,dk;
+  int32_t k,hx,i,j;
+  uint32_t lx;
+
+  hx = hi_word(x);
+  lx = lo_word(x);
+
+  k = 0;
+  if (hx < 0x00100000) {                       /* x < 2**-1022 */
+    // fdlibm writes these as -two54/zero and (x-x)/zero, which also raise the divide-by-zero
+    // and invalid flags. MSVC rejects a compile-time division by a zero constant (C2124), so
+    // return the same values directly. The estimators never call log() with these inputs.
+    if (((hx & 0x7fffffff) | lx) == 0) {                          /* log(+-0) = -inf */
+      return -std::numeric_limits<double>::infinity();
+    }
+    if (hx < 0) {                                                 /* log(-#) = NaN */
+      return std::numeric_limits<double>::quiet_NaN();
+    }
+    k -= 54; x *= two54;                       /* subnormal: scale up */
+    hx = hi_word(x);
+  }
+  if (hx >= 0x7ff00000) { return x+x; }
+  k += (hx>>20) - 1023;
+  hx &= 0x000fffff;
+  i = (hx + 0x95f64) & 0x100000;
+  set_hi_word(x, (uint32_t)(hx | (i ^ 0x3ff00000)));   /* normalize x or x/2 */
+  k += (i>>20);
+  f = x - 1.0;
+  if ((0x000fffff & (2+hx)) < 3) {             /* |f| < 2**-20 */
+    if (f == zero) {
+      if (k == 0) { return zero; }
+      dk = (double)k; return rnd(dk*ln2_hi) + rnd(dk*ln2_lo);
+    }
+    R = rnd(rnd(f*f)*rnd(0.5 - rnd(0.33333333333333333*f)));
+    if (k == 0) { return f-R; }
+    dk = (double)k; return rnd(dk*ln2_hi) - (rnd(R - rnd(dk*ln2_lo)) - f);
+  }
+  s = f/(2.0+f);
+  dk = (double)k;
+  z = s*s;
+  i = hx - 0x6147a;
+  w = z*z;
+  j = 0x6b851 - hx;
+  t1 = rnd(w*rnd(Lg2 + rnd(w*rnd(Lg4 + rnd(w*Lg6)))));
+  t2 = rnd(z*rnd(Lg1 + rnd(w*rnd(Lg3 + rnd(w*rnd(Lg5 + rnd(w*Lg7)))))));
+  i |= j;
+  R = t2 + t1;
+  if (i > 0) {
+    hfsq = rnd(0.5*f)*f;
+    if (k == 0) { return f - rnd(hfsq - rnd(s*(hfsq+R))); }
+    return rnd(dk*ln2_hi) - (rnd(hfsq - rnd(rnd(s*(hfsq+R)) + rnd(dk*ln2_lo))) - f);
+  } else {
+    if (k == 0) { return f - rnd(s*(f-R)); }
+    return rnd(dk*ln2_hi) - (rnd(rnd(s*(f-R)) - rnd(dk*ln2_lo)) - f);
+  }
+}
+
+} // namespace fdlibm
+#endif

hll/include/CouponList-internal.hpp

@@ -33,15 +33,15 @@ namespace datasketches {
 
 template<typename A>
 CouponList<A>::CouponList(uint8_t lgConfigK, target_hll_type tgtHllType, hll_mode mode, const A& allocator):
-HllSketchImpl<A>(lgConfigK, tgtHllType, mode, false),
+HllSketchImpl<A>(lgConfigK, tgtHllType, mode),
 couponCount_(0),
 oooFlag_(false),
 coupons_(1ULL << (mode == hll_mode::LIST ? hll_constants::LG_INIT_LIST_SIZE : hll_constants::LG_INIT_SET_SIZE), 0, allocator)
 {}
 
 template<typename A>
 CouponList<A>::CouponList(const CouponList& that, const target_hll_type tgtHllType):
-HllSketchImpl<A>(that.lgConfigK_, tgtHllType, that.mode_, false),
+HllSketchImpl<A>(that.lgConfigK_, tgtHllType, that.mode_),
 couponCount_(that.couponCount_),
 oooFlag_(that.oooFlag_),
 coupons_(that.coupons_)

hll/include/HarmonicNumbers-internal.hpp

@@ -21,6 +21,7 @@
 #define _HARMONICNUMBERS_INTERNAL_HPP_
 
 #include "HarmonicNumbers.hpp"
+#include "fdlibm_log.hpp"
 
 #include <cmath>
 
@@ -70,7 +71,7 @@ double HarmonicNumbers<A>::harmonicNumber(const uint64_t x_i) {
   } else {
     double x = static_cast<double>(x_i);
     double invSq = 1.0 / (x * x);
-    double sum = log(x) + EULER_MASCHERONI_CONSTANT + (1.0 / (2.0 * x));
+    double sum = fdlibm::log(x) + EULER_MASCHERONI_CONSTANT + (1.0 / (2.0 * x));
     /* note: the number of terms included from this series expansion is appropriate
        for the size of the exact table (25) and the precision of doubles */
     double pow = invSq; // now n^-2

hll/include/Hll4Array-internal.hpp

@@ -30,8 +30,8 @@
 namespace datasketches {
 
 template<typename A>
-Hll4Array<A>::Hll4Array(uint8_t lgConfigK, bool startFullSize, const A& allocator):
-HllArray<A>(lgConfigK, target_hll_type::HLL_4, startFullSize, allocator),
+Hll4Array<A>::Hll4Array(uint8_t lgConfigK, const A& allocator):
+HllArray<A>(lgConfigK, target_hll_type::HLL_4, allocator),
 auxHashMap_(nullptr)
 {
   const uint32_t numBytes = this->hll4ArrBytes(lgConfigK);
@@ -53,7 +53,7 @@ Hll4Array<A>::Hll4Array(const Hll4Array<A>& that) :
 
 template<typename A>
 Hll4Array<A>::Hll4Array(const HllArray<A>& other) :
-  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_4, other.isStartFullSize(), other.getAllocator()),
+  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_4, other.getAllocator()),
   auxHashMap_(nullptr)
 {
   const int numBytes = this->hll4ArrBytes(this->lgConfigK_);

hll/include/Hll4Array.hpp

@@ -28,7 +28,7 @@ namespace datasketches {
 template<typename A>
 class Hll4Array final : public HllArray<A> {
   public:
-    explicit Hll4Array(uint8_t lgConfigK, bool startFullSize, const A& allocator);
+    explicit Hll4Array(uint8_t lgConfigK, const A& allocator);
     explicit Hll4Array(const Hll4Array<A>& that);
     explicit Hll4Array(const HllArray<A>& that);
 

hll/include/Hll6Array-internal.hpp

@@ -27,16 +27,16 @@
 namespace datasketches {
 
 template<typename A>
-Hll6Array<A>::Hll6Array(uint8_t lgConfigK, bool startFullSize, const A& allocator):
-HllArray<A>(lgConfigK, target_hll_type::HLL_6, startFullSize, allocator)
+Hll6Array<A>::Hll6Array(uint8_t lgConfigK, const A& allocator):
+HllArray<A>(lgConfigK, target_hll_type::HLL_6, allocator)
 {
   const int numBytes = this->hll6ArrBytes(lgConfigK);
   this->hllByteArr_.resize(numBytes, 0);
 }
 
 template<typename A>
 Hll6Array<A>::Hll6Array(const HllArray<A>& other) :
-  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_6, other.isStartFullSize(), other.getAllocator())
+  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_6, other.getAllocator())
 {
   const int numBytes = this->hll6ArrBytes(this->lgConfigK_);
   this->hllByteArr_.resize(numBytes, 0);

hll/include/Hll6Array.hpp

@@ -30,7 +30,7 @@ class Hll6Iterator;
 template<typename A>
 class Hll6Array final : public HllArray<A> {
   public:
-    Hll6Array(uint8_t lgConfigK, bool startFullSize, const A& allocator);
+    Hll6Array(uint8_t lgConfigK, const A& allocator);
     explicit Hll6Array(const HllArray<A>& that);
 
     virtual ~Hll6Array() = default;

hll/include/Hll8Array-internal.hpp

@@ -25,16 +25,16 @@
 namespace datasketches {
 
 template<typename A>
-Hll8Array<A>::Hll8Array(uint8_t lgConfigK, bool startFullSize, const A& allocator):
-HllArray<A>(lgConfigK, target_hll_type::HLL_8, startFullSize, allocator)
+Hll8Array<A>::Hll8Array(uint8_t lgConfigK, const A& allocator):
+HllArray<A>(lgConfigK, target_hll_type::HLL_8, allocator)
 {
   const int numBytes = this->hll8ArrBytes(lgConfigK);
   this->hllByteArr_.resize(numBytes, 0);
 }
 
 template<typename A>
 Hll8Array<A>::Hll8Array(const HllArray<A>& other):
-  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_8, other.isStartFullSize(), other.getAllocator())
+  HllArray<A>(other.getLgConfigK(), target_hll_type::HLL_8, other.getAllocator())
 {
   const int numBytes = this->hll8ArrBytes(this->lgConfigK_);
   this->hllByteArr_.resize(numBytes, 0);

hll/include/Hll8Array.hpp

@@ -30,7 +30,7 @@ class Hll8Iterator;
 template<typename A>
 class Hll8Array final : public HllArray<A> {
   public:
-    Hll8Array(uint8_t lgConfigK, bool startFullSize, const A& allocator);
+    Hll8Array(uint8_t lgConfigK, const A& allocator);
     explicit Hll8Array(const HllArray<A>& that);
 
     virtual ~Hll8Array() = default;

hll/include/HllArray-internal.hpp

@@ -21,6 +21,7 @@
 #define _HLLARRAY_INTERNAL_HPP_
 
 #include "HllArray.hpp"
+#include "fdlibm_log.hpp"
 #include "HllUtil.hpp"
 #include "HarmonicNumbers.hpp"
 #include "CubicInterpolation.hpp"
@@ -35,8 +36,8 @@
 namespace datasketches {
 
 template<typename A>
-HllArray<A>::HllArray(uint8_t lgConfigK, target_hll_type tgtHllType, bool startFullSize, const A& allocator):
-HllSketchImpl<A>(lgConfigK, tgtHllType, hll_mode::HLL, startFullSize),
+HllArray<A>::HllArray(uint8_t lgConfigK, target_hll_type tgtHllType, const A& allocator):
+HllSketchImpl<A>(lgConfigK, tgtHllType, hll_mode::HLL),
 hipAccum_(0.0),
 kxq0_(1 << lgConfigK),
 kxq1_(0.0),
@@ -49,7 +50,7 @@ rebuild_kxq_curmin_(false)
 
 template<typename A>
 HllArray<A>::HllArray(const HllArray& other, target_hll_type tgtHllType) :
-  HllSketchImpl<A>(other.getLgConfigK(), tgtHllType, hll_mode::HLL, other.isStartFullSize()),
+  HllSketchImpl<A>(other.getLgConfigK(), tgtHllType, hll_mode::HLL),
   // remaining fields are initialized to empty sketch defaults
   // and left to subclass constructor to populate
   hipAccum_(0.0),
@@ -109,7 +110,6 @@ HllArray<A>* HllArray<A>::newHll(const void* bytes, size_t len, const A& allocat
   const target_hll_type tgtHllType = HllSketchImpl<A>::extractTgtHllType(data[hll_constants::MODE_BYTE]);
   const bool oooFlag = ((data[hll_constants::FLAGS_BYTE] & hll_constants::OUT_OF_ORDER_FLAG_MASK) ? true : false);
   const bool comapctFlag = ((data[hll_constants::FLAGS_BYTE] & hll_constants::COMPACT_FLAG_MASK) ? true : false);
-  const bool startFullSizeFlag = ((data[hll_constants::FLAGS_BYTE] & hll_constants::FULL_SIZE_FLAG_MASK) ? true : false);
 
   const uint8_t lgK = data[hll_constants::LG_K_BYTE];
   const uint8_t curMin = data[hll_constants::HLL_CUR_MIN_BYTE];
@@ -139,7 +139,7 @@ HllArray<A>* HllArray<A>::newHll(const void* bytes, size_t len, const A& allocat
     aux_ptr = aux_hash_map_ptr(auxHashMap, auxHashMap->make_deleter());
   }
 
-  HllArray<A>* sketch = HllSketchImplFactory<A>::newHll(lgK, tgtHllType, startFullSizeFlag, allocator);
+  HllArray<A>* sketch = HllSketchImplFactory<A>::newHll(lgK, tgtHllType, allocator);
   sketch->putCurMin(curMin);
   sketch->putOutOfOrderFlag(oooFlag);
   if (!oooFlag) { sketch->putHipAccum(hip); }
@@ -180,12 +180,11 @@ HllArray<A>* HllArray<A>::newHll(std::istream& is, const A& allocator) {
   const target_hll_type tgtHllType = HllSketchImpl<A>::extractTgtHllType(listHeader[hll_constants::MODE_BYTE]);
   const bool oooFlag = ((listHeader[hll_constants::FLAGS_BYTE] & hll_constants::OUT_OF_ORDER_FLAG_MASK) ? true : false);
   const bool comapctFlag = ((listHeader[hll_constants::FLAGS_BYTE] & hll_constants::COMPACT_FLAG_MASK) ? true : false);
-  const bool startFullSizeFlag = ((listHeader[hll_constants::FLAGS_BYTE] & hll_constants::FULL_SIZE_FLAG_MASK) ? true : false);
 
   const uint8_t lgK = listHeader[hll_constants::LG_K_BYTE];
   const uint8_t curMin = listHeader[hll_constants::HLL_CUR_MIN_BYTE];
 
-  HllArray* sketch = HllSketchImplFactory<A>::newHll(lgK, tgtHllType, startFullSizeFlag, allocator);
+  HllArray* sketch = HllSketchImplFactory<A>::newHll(lgK, tgtHllType, allocator);
   typedef std::unique_ptr<HllArray<A>, std::function<void(HllSketchImpl<A>*)>> hll_array_ptr;
   hll_array_ptr sketch_ptr(sketch, sketch->get_deleter());
   sketch->putCurMin(curMin);
@@ -570,7 +569,7 @@ double HllArray<A>::getHllBitMapEstimate() const {
 
   //This will eventually go away.
   if (numUnhitBuckets == 0) {
-    return configK * log(configK / 0.5);
+    return configK * fdlibm::log(configK / 0.5);
   }
 
   const uint32_t numHitBuckets = configK - numUnhitBuckets;
@@ -603,34 +602,37 @@ bool HllArray<A>::isRebuildKxqCurminFlag() const {
 template<typename A>
 void HllArray<A>::check_rebuild_kxq_cur_min() {
   if (!rebuild_kxq_curmin_) { return; }
+  // the deferred rebuild is only ever set on an HLL_8 union gadget. Guarding here also keeps
+  // this from ever rewriting curMin_ on an HLL_4 array, whose nibbles are stored relative to it.
+  if (this->getCurMode() != hll_mode::HLL || this->getTgtHllType() != target_hll_type::HLL_8) {
+    rebuild_kxq_curmin_ = false;
+    return;
+  }
 
-  uint8_t cur_min = 64;
-  uint32_t num_at_cur_min = 0;
+  uint32_t num_zeros = 0;
   double kxq0 = 1 << this->lgConfigK_;
   double kxq1 = 0;
 
-  auto it = this->begin(true); // want all points to adjust cur_min
+  auto it = this->begin(true); // want all slots, including the empty ones
   const auto end = this->end();
   while (it != end) {
     uint8_t v = HllUtil<A>::getValue(*it);
     if (v > 0) {
       if (v < 32) { kxq0 += INVERSE_POWERS_OF_2[v] - 1.0; }
       else        { kxq1 += INVERSE_POWERS_OF_2[v] - 1.0; }
-    }
-    if (v > cur_min) { ++it; continue; }
-    if (v < cur_min) {
-      cur_min = v;
-      num_at_cur_min = 1;
     } else {
-      ++num_at_cur_min;
-    }    
+      ++num_zeros;
+    }
     ++it;
   }
 
   kxq0_ = kxq0;
   kxq1_ = kxq1;
-  curMin_ = cur_min;
-  numAtCurMin_ = num_at_cur_min;
+  // HLL_8 convention: curMin is always 0 and numAtCurMin is the number of zero registers.
+  // That is the representation the incremental update path maintains, so the rebuilt state is
+  // indistinguishable from it and the timing of this rebuild is not observable in the image.
+  curMin_ = 0;
+  numAtCurMin_ = num_zeros;
   rebuild_kxq_curmin_ = false;
   // HipAccum is not affected
 

hll/include/HllArray.hpp

@@ -33,7 +33,7 @@ class HllArray : public HllSketchImpl<A> {
   public:
     using vector_bytes = std::vector<uint8_t, typename std::allocator_traits<A>::template rebind_alloc<uint8_t>>;
 
-    HllArray(uint8_t lgConfigK, target_hll_type tgtHllType, bool startFullSize, const A& allocator);
+    HllArray(uint8_t lgConfigK, target_hll_type tgtHllType, const A& allocator);
     explicit HllArray(const HllArray& other, target_hll_type tgtHllType);
 
     static HllArray* newHll(const void* bytes, size_t len, const A& allocator);

hll/include/HllSketch-internal.hpp

@@ -45,7 +45,7 @@ template<typename A>
 hll_sketch_alloc<A>::hll_sketch_alloc(uint8_t lg_config_k, target_hll_type tgt_type, bool start_full_size, const A& allocator) {
   HllUtil<A>::checkLgK(lg_config_k);
   if (start_full_size) {
-    sketch_impl = HllSketchImplFactory<A>::newHll(lg_config_k, tgt_type, start_full_size, allocator);
+    sketch_impl = HllSketchImplFactory<A>::newHll(lg_config_k, tgt_type, allocator);
   } else {
     typedef typename std::allocator_traits<A>::template rebind_alloc<CouponList<A>> clAlloc;
     sketch_impl = new (clAlloc(allocator).allocate(1)) CouponList<A>(lg_config_k, tgt_type, hll_mode::LIST, allocator);
@@ -107,10 +107,8 @@ hll_sketch_alloc<A>& hll_sketch_alloc<A>::operator=(hll_sketch_alloc<A>&& other)
 }
 
 template<typename A>
-void hll_sketch_alloc<A>::reset() {
-  // TODO: need to allow starting from a full-sized sketch
-  //       (either here or in other implementation)
-  sketch_impl = sketch_impl->reset();
+void hll_sketch_alloc<A>::reset(bool full_size) {
+  sketch_impl = sketch_impl->reset(full_size);
 }
 
 template<typename A>